Mesh moving / ALE in DOLFINx — example or official API

[jialu_gao]

Hi DOLFINx community,

I’m working on a simulation requiring moving mesh / mesh deformation (ALE or Lagrangian description) using DOLFINx, but I found that many older example codes use dolfin.ALE.move(...) from legacy DOLFIN and are incompatible with recent DOLFINx versions.

I would kindly ask:

1. Is there an official API or recommended workflow in DOLFINx (version ≥ 0.x) for mesh movement in the sense of ALE (arbitrary Lagrangian-Eulerian) or purely Lagrangian mesh deformation?

2. Are there official example scripts or demos showcasing mesh movement (especially with non-trivial geometry, boundary displacement, or deformable domains)?

3. If not yet fully integrated, what is the community recommended workaround (for example, manually updating mesh.geometry.x)?

4. Any guidance regarding mesh quality / smoothing after large deformations (for example, harmonic smoothing, mesh relaxation) in DOLFINx would also be appreciated.

My current scenario:

• Mesh imported from Gmsh with a cut geometry (e.g., rectangle minus circle)
• I want to apply prescribed boundary displacement on a subset of facets, update interior nodes accordingly
• I’ve tried directly updating mesh.geometry.x[:] = old + u_eval and it seems to work for small deformation, but I’m worried about cell quality and complexity for larger deformations

Any pointers, example code, or official plans would be great. Thank you for your time and help!

Best regards,
jialu gao

[Stein]

@jialu_gao

Are there official example scripts or demos showcasing mesh movement (especially with non-trivial geometry, boundary displacement, or deformable domains)?

I don’t believe there is an explicit ALE example available, I’m afraid.

@jialu_gao

Is there an official API or recommended workflow in DOLFINx (version ≥ 0.x) for mesh movement in the sense of ALE (arbitrary Lagrangian-Eulerian) or purely Lagrangian mesh deformation?

These are two separate questions. For Lagrangian deformation, you indeed simply can update mesh.geometry.x. But in my experience, you do not want to move the actual mesh nodes. Instead, you solve for an additional mesh displacement field. You’d re-define your gradient and divergence operators to take into account this additional mesh displacement.

@jialu_gao

Any guidance regarding mesh quality / smoothing after large deformations (for example, harmonic smoothing, mesh relaxation) in DOLFINx would also be appreciated.

The equation that you solve for the mesh displacement can involve terms that counteract severe mesh distortion (e.g. ‘‘jacobian-based stiffness’’).

I really like this old demo by David: mae-207-fea-for-coupled-problems/fsi/fitted-fsi-example.py at master · david-kamensky/mae-207-fea-for-coupled-problems · GitHub

But, it’s FEniCS, not FEniCSx… I started translating it a while back but must have a bug somewhere as I wasn’t getting the Turek-Hron benchmark results correctly. If that’s what you’re interested in I could post that here. (I’d have to clean it up though, so it’d take some effort).

[jialu_gao]

Thank you very much for your detailed explanation, Stein — this is extremely helpful!

[federicobetti]

Hi @Stein , this post caught my attention as I am also implementing a monolithic FEniCSx solver for the Turek Hron benchmark and I am struggling to get the correct results in terms of drag, lift coefficients and tip displacement. To be fair, the latter two are correct, the drag is approximately correct but the lift coefficient is completely off for some reason. I took as a starting point the turtleFSI repository (in legacy dolfin). May I ask if you are encountering similar issues?

[Stein]

To be honest, my ‘flag’ doesn’t even start to oscillate It just bounces a couple of times before leveling off. I’m sure it’s some trivial error in the material properties or boundary conditions or so. I could aford myself only a couple of hours to work on it, and I spend most my time trying out FEniCSx’s then-new multi domain capabilities. Maybe I’ll find some time later this week to clean up the code into something stand-alone.

[federicobetti]

Hi @Stein and thanks for your reply. I have a different problem at the moment in terms of results, but may I ask if you are also using a theta method to solve the equations monolithically?

[Stein]

Yes, simple theta method.

[federicobetti]

Thank you very much for your reply @stein. I launched this morning a pure CFD case (without flag domain) and my results are now matching those of Turek Hron for all vcases, which makes me think that there is something wrong in how I am computing the drag and lift coefficient when the boundary with the flag is a dS measure instead of a ds one.

So far I am doing it like this:

d = dvp_["n"].sub(0)
v = dvp_["n"].sub(1)
p = dvp_["n"].sub(2)


mesh = d.function_space.mesh
Dr = -dolfinx.fem.assemble_scalar(
   dolfinx.fem.form((sigma(v, p, d, mu_f) * n)[0] * ds(6))
)
Li = -dolfinx.fem.assemble_scalar(
   dolfinx.fem.form((sigma(v, p, d, mu_f) * n)[1] * ds(6))
)

Dr += -dolfinx.fem.assemble_scalar(
   dolfinx.fem.form((sigma(v("+"), p("+"), d("+"), mu_f) * n("+"))[0] * dS(5)
)
Li += -dolfinx.fem.assemble_scalar(
   dolfinx.fem.form((sigma(v("+"), p("+"), d("+"), mu_f) * n("+"))[1] * dS(5)
)

so perhaps the dS(5) computation is somehow wrong.

[dokken]

Using a “+” restriction is not correct on an interior facet (as it is arbitrary).
To do a proper, one-sided, oriented integral, see for instance:

Add one-sided integration example as example of custom integration

and

Correct way to compute internal interface fluxes

The key take-away is that: In most finite element methods, the jump of a quantity across a surface does not require a specific orientation. For instance, in DG methods, you get jump integrals across all internal facets of the grids, but mathematically the formulation doesn’t care which cell is positive and which one is negative. DOLFINx follows the same notion, that for an interior facet integral (dS) which cell comes first (+) and which one is last (-) doesn’t matter if put in the correct va…

[federicobetti]

Thanks for the input @dokken . I am now doing this

def create_internal_ds_measure(mesh, domain_tag, facet_tag):

    """Create internal ds measure for the interior facets."""

    fdim = mesh.topology.dim - 1
    tdim = mesh.topology.dim
    f_to_c = mesh.topology.connectivity(fdim, tdim)
    c_to_f = mesh.topology.connectivity(tdim, fdim)


    cell_map = mesh.topology.index_map(tdim)
    num_cells = cell_map.size_local + cell_map.num_ghosts

    facets_to_integrate = facet_tag.find(5)
    facet_map = mesh.topology.index_map(fdim)
    num_facets = facet_map.size_local + facet_map.num_ghosts

    integration_entities = []
    for i, facet in enumerate(facets_to_integrate):
        if facet >= facet_map.size_local:
            continue
        cells = f_to_c.links(facet)
        marked_cells = domain_tag.values[cells]
        correct_cell = np.flatnonzero(marked_cells == 2)

        assert len(correct_cell) == 1
        local_facets = c_to_f.links(cells[correct_cell[0]])
        local_index = np.flatnonzero(local_facets == facet)
        assert len(local_index) == 1

        integration_entities.append(cells[correct_cell[0]])
        integration_entities.append(local_index[0])

    ds_internal = ufl.Measure("ds", domain=mesh, subdomain_data=[
                    (5, np.asarray(integration_entities, dtype=np.int32).flatten())])




    return ds_internal

and I then change the call for the drag and lift calculation to

Dr = -dolfinx.fem.assemble_scalar(
dolfinx.fem.form((sigma(v, p, d, mu_f) * n)[0] * ds(6))
)
Li = -dolfinx.fem.assemble_scalar(
dolfinx.fem.form((sigma(v, p, d, mu_f) * n)[1] * ds(6))
)

Dr += -dolfinx.fem.assemble_scalar(
dolfinx.fem.form((sigma(v, p, d, mu_f) * n)[0] * ds_internal(5))
)
Li += -dolfinx.fem.assemble_scalar(
dolfinx.fem.form((sigma(v, p, d, mu_f) * n)[1] * ds_internal(5))
)

but I am getting the error

Traceback (most recent call last):
File “/root/main_fsi.py”, line 194, in 
main(args.case)
File “/root/main_fsi.py”, line 150, in main
Dr, Li = compute_drag_lift(
File “/root/src/fsi/modules/solver_utils.py”, line 124, in compute_drag_lift
dolfinx.fem.form((sigma(v, p, d, mu_f) * n)[0] * ds_internal(5))
File “/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py”, line 176, in form
return _create_form(form)
File “/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py”, line 171, in _create_form
return _form(form)
File “/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py”, line 165, in _form
return formcls(ufcx_form, V, coeffs, constants, subdomains, mesh, module.ffi, code)
File “/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py”, line 51, in init
super().init(ffi.cast(“uintptr_t”, ffi.addressof(self._ufcx_form)),
TypeError: init(): incompatible constructor arguments. The following argument types are supported:
1. dolfinx.cpp.fem.Form_float64(spaces: List[dolfinx::fem::FunctionSpace], integrals: Dict[dolfinx::fem::IntegralType, Tuple[List[Tuple[int, object]], dolfinx.cpp.mesh.MeshTags_int32]], coefficients: List[dolfinx.cpp.fem.Function_float64], constants: List[dolfinx.cpp.fem.Constant_float64], need_permutation_data: bool, mesh: dolfinx.cpp.mesh.Mesh = None)
2. dolfinx.cpp.fem.Form_float64(form: int, spaces: List[dolfinx::fem::FunctionSpace], coefficients: List[dolfinx.cpp.fem.Function_float64], constants: List[dolfinx.cpp.fem.Constant_float64], subdomains: Dict[dolfinx::fem::IntegralType, dolfinx.cpp.mesh.MeshTags_int32], mesh: dolfinx.cpp.mesh.Mesh)

Invoked with: <cdata ‘uintptr_t’ 126747154867488>, 
, [<dolfinx.cpp.fem.Function_float64 object at 0x73469ec770f0>, <dolfinx.cpp.fem.Function_float64 object at 0x73469e9a71f0>, <dolfinx.cpp.fem.Function_float64 object at 0x73469e9a76b0>], 
, {<IntegralType.cell: 0>: None, <IntegralType.exterior_facet: 1>: [(5, array([ 375,    2,  397,    2,  398,    2,  421,    2,  444,    2,  464,
2,  488,    2,  489,    2,  514,    2,  532,    2,  536,    2,
559,    2,  581,    2,  604,    2,  609,    2,  633,    2,  634,
2,  659,    2,  660,    2,  684,    2,  706,    2,  728,    2,
748,    2,  772,    2,  779,    2,  806,    2,  824,    2,  835,
2,  856,    2,  861,    2,  886,    2,  891,    2,  917,    2,
927,    2,  963,    2,  964,    2,  933,    2,  935,    2,  973,
2,  960,    2, 1011,    2, 1027,    2, 1070,    2, 1106,    2,
1110,    2, 1107,    2, 1121,    2, 1086,    2, 1125,    1, 1162,
2, 1246,    2, 1254,    1], dtype=int32))], <IntegralType.interior_facet: 2>: None, <IntegralType.vertex: 3>: None}, <dolfinx.mesh.Mesh object at 0x73469e94cfe0> 

[Stein]

Is there a way for you to create an MWE? It’s not immediately obvious to me what is wrong with your code, but for some reason, (sigma(v, p, d, mu_f) * n)[0] * ds_internal(5) is not producing something dolfinx.fem.form can work with. That’s probably somehow due to ds_internal(5), but I’d have to play around with an MWE to pinpoint what is wrong.

Your MWE doesn’t need to actually perform a solve, you can simply initialize arbitray fields.

[federicobetti]

Hi @Stein , I fixed the issue by updating to dolfinx 0.8.0 instead of the 0.6.0 version that I was previously using. And indeed, using ds_internal instead of dS are giving me results that are closer to the Turek-Hron benchmark for the lift as well.

I am now running the examples on larger meshes to understand if they will converge to the reference values, but so far it looks like this is the case.

[Stein]

@Stein

Maybe I’ll find some time later this week to clean up the code into something stand-alone.

For those landing here looking for a FEniCSx FSI script:

FEniCSx4Flow / Fx4F FSI · GitLab

I didn’t end up spending the time to make it run stand-alone. So this has to be cloned recursively:

git clone --recurse-submodules git@gitlab.com:fenicsx4flow/fx4f_fsi.git

and ran with a custom docker image:

docker run --volume="$(pwd):/app" --user $(id -u):$(id -g) steinkfstoter/fx4f:latest python3 simulation.py <arguments>

where <arguments> can be specified as (e.g.) --nx 16 --ny 16 --T 5, etc.
Log files and VTX output are then generated in a time-stamped ./output/... directory.