Moving submesh diverges from its designated trajectory in a two-phase problem

[chenyongxin]

I have a two phase system for ALE-FSI problem. I extract a submesh from parent mesh to make immersed solid. In every step, the parent mesh is moved by appending a displacement function Delta_x as

_mesh.geometry.x[:, :tdim] += Delta_x.x.array.reshape((-1, tdim))

The submesh is also needed to be updated accordingly. Here I create a displacement function Delta_x_s on the submesh, and Delta_x_s is interpolated from Delta_x as suggested by the following.

@dokken

@chenyongxin

3. My practical application is FSI + ALE, so the mesh moves in every step. If with a submesh, any idea to associate the movement of a submesh to the parent mesh?

Then you should interpolate any displacement in the parent mesh onto the submesh and move the mesh afterwards

I have tried the aforementioned step by either matching meshes or non-matching meshes interpolation, but the submesh is wrongly updated as the attached figure shows. I sense that it is due to the DOF-vertex relation is mismatched on submesh.

MWE with DOLFINx v0.9.0:

import os
os.environ['OMP_NUM_THREADS'] = '1'
from dolfinx import io, fem, mesh
from mpi4py import MPI
import numpy as np
import ufl
from dolfinx.fem.petsc import LinearProblem

comm = MPI.COMM_WORLD
prefix = "output"
mesh_order = 1
_mesh = mesh.create_unit_square(comm, 16, 16)
tdim = 2
fdim = tdim - 1
_mesh.geometry.x[:, :tdim] -= 0.5
_mesh.geometry.x[:, :tdim] *= 2

# Create mesh and functions
V = fem.functionspace(_mesh, ("Lagrange", mesh_order, (tdim,)))
DG0 = fem.functionspace(_mesh, ("DG", 0))
mask = fem.Function(DG0)

Omega_s = mesh.locate_entities(_mesh, tdim, lambda x: np.logical_and(np.abs(x[0]) < 0.25 + 1e-5,
                                                                     np.abs(x[1]) < 0.25 + 1e-5))
_mesh.topology.create_connectivity(tdim, tdim)
_mesh.topology.create_connectivity(fdim, tdim)
dofs = fem.locate_dofs_topological(DG0, tdim, Omega_s)
mask.x.array[dofs] = 1

disp = fem.Function(V)
wall_function = fem.Function(V)
Delta_x = fem.Function(V)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Create submesh
submesh, solid_cell_map, _, _ = mesh.create_submesh(_mesh, tdim, Omega_s)
V_s = fem.functionspace(submesh, ("Lagrange", mesh_order, (tdim, )))
Delta_x_s = fem.Function(V_s)

boundary_facets = mesh.exterior_facet_indices(_mesh.topology)
bc_solid = fem.dirichletbc(disp, fem.locate_dofs_topological(V, tdim, Omega_s))
bc_wall = fem.dirichletbc(wall_function, fem.locate_dofs_topological(V, fdim, boundary_facets))
bcs = [bc_solid, bc_wall]

N = 20
delta_alpha = - np.pi / 4 / N
x = ufl.SpatialCoordinate(_mesh)
center = ufl.as_vector([0.0, 0.0])
cos = ufl.cos(delta_alpha)
sin = ufl.sin(delta_alpha)
R = ufl.as_matrix([[cos, -sin],    # rotation matrix
                   [sin,  cos]])
I = ufl.Identity(tdim)

# Weak form for ALE mapping 
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = ufl.inner(fem.Constant(_mesh, [0.] * tdim), v) * ufl.dx

# Visualize initial result
with io.XDMFFile(MPI.COMM_WORLD, f"{prefix}/mesh_0.xdmf", "w") as file:
    file.write_mesh(_mesh)
    file.write_function(mask, 0)
    
with io.XDMFFile(MPI.COMM_WORLD, f"{prefix}/submesh_0.xdmf", "w") as file:
    file.write_mesh(submesh)
    file.write_function(Delta_x, 0)

# Time iterate
for step in range(N):
    
    # Set up boundary conditions
    expr = fem.Expression((R - I) * (x - center), V.element.interpolation_points())
    disp.interpolate(expr)
    disp.x.scatter_forward()
    
    # Solve ALE mapping
    if comm.rank == 0:
        print(f"Step {step + 1} / {N}", flush=True)
    problem = LinearProblem(a, L, bcs = bcs, u = Delta_x,
                petsc_options={"ksp_type": "preonly", "pc_type": "lu", 
                               "pc_factor_mat_solver_type": "mumps"})
        
    problem.solve()
    Delta_x.x.scatter_forward()
    
    _mesh.geometry.x[:, :tdim] += Delta_x.x.array.reshape((-1, tdim))
    
    # Move submesh by interpolating displacement
    # 1. Matching mesh interpolation
    Delta_x_s.x.array[:] = 0.0
    Delta_x_s.interpolate(Delta_x,
                          cells0 = solid_cell_map, 
                          cells1 = np.arange(len(solid_cell_map)))
    Delta_x_s.x.scatter_forward()
    submesh.geometry.x[:, :tdim] += Delta_x_s.x.array.reshape((-1, tdim))
    
    # 2. Non-matching mesh interpolation
    """
    dest_mesh_cell_map = V_s.mesh.topology.index_map(tdim)
    num_cells_on_proc = dest_mesh_cell_map.size_local + dest_mesh_cell_map.num_ghosts
    cells = np.arange(num_cells_on_proc, dtype = np.int32)
    interpolation_data = fem.create_interpolation_data(V_s, V, cells, padding=1e-8)
    Delta_x_s.interpolate_nonmatching(Delta_x, cells, interpolation_data=interpolation_data)
    submesh.geometry.x[:, :tdim] += Delta_x_s.x.array.reshape((-1, tdim))
    """
        
    with io.XDMFFile(MPI.COMM_WORLD, f"{prefix}/mesh_{step + 1}.xdmf", "w") as file:
        file.write_mesh(_mesh)
        file.write_function(mask, step + 1)
    
    with io.XDMFFile(MPI.COMM_WORLD, f"{prefix}/submesh_{step + 1}.xdmf", "w") as file:
        file.write_mesh(submesh)
        file.write_function(Delta_x_s, step + 1)

[dokken]

This is not the right thing to do, as this assumes alot of things about the ordering of degrees of freedom.

Here is a minimal example that maps a displacement field from a mesh onto a submesh:

"""
Map degrees of freedom to the mesh geometry nodes
and corresponding submesh.

Author: Jørgen S. Dokken and Henrik N.T. Finsberg
SPDX-License-Identifier: MIT
"""

from mpi4py import MPI
import dolfinx
import numpy as np


def vertex_to_dof_map_vectorized(V):
    """Create a map from the vertices of the mesh to the corresponding degree of freedom."""
    mesh = V.mesh
    num_vertices_per_cell = dolfinx.cpp.mesh.cell_num_entities(
        mesh.topology.cell_type, 0
    )

    dof_layout2 = np.empty((num_vertices_per_cell,), dtype=np.int32)
    for i in range(num_vertices_per_cell):
        var = V.dofmap.dof_layout.entity_dofs(0, i)
        assert len(var) == 1
        dof_layout2[i] = var[0]

    num_vertices = (
        mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
    )

    c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
    assert (c_to_v.offsets[1:] - c_to_v.offsets[:-1] == c_to_v.offsets[1]).all(), (
        "Single cell type supported"
    )

    vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)
    vertex_to_dof_map[c_to_v.array] = V.dofmap.list[:, dof_layout2].reshape(-1)
    return vertex_to_dof_map


if __name__ == "__main__":
    # Create parent mesh and mock displacement
    mesh = dolfinx.mesh.create_rectangle(MPI.COMM_WORLD, [[-1, -1], [1, 1]], [30, 30])
    V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim,)))
    u = dolfinx.fem.Function(V)
    u.interpolate(lambda x: (0.1 * np.sin(np.pi * x[0]), 0.1 * np.sin(np.pi * x[1])))

    # Create a submesh
    tdim = mesh.topology.dim
    fdim = tdim - 1
    tol = 1e3 * np.finfo(dolfinx.default_scalar_type).eps
    Omega_s = dolfinx.mesh.locate_entities(
        mesh,
        tdim,
        lambda x: np.logical_and(
            np.abs(x[0]) < 0.25 + 1e-5, np.abs(x[1]) < 0.25 + 1e-5
        ),
    )

    # Create submesh
    submesh, solid_cell_map, vertex_map, node_map = dolfinx.mesh.create_submesh(
        mesh, tdim, Omega_s
    )

    # Displace parent mesh
    num_parent_vertices = (
        mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
    )
    # Create the dof -> vertex map and vertex->node map
    vertex_to_dof_map = vertex_to_dof_map_vectorized(V)
    mesh.topology.create_connectivity(0, tdim)
    vertex_to_node_map = dolfinx.mesh.entities_to_geometry(
        mesh, 0, np.arange(num_parent_vertices, dtype=np.int32)
    ).reshape(-1)
    u_values = u.x.array.reshape(-1, mesh.geometry.dim)

    mesh.geometry.x[vertex_to_node_map, :tdim] += u_values[vertex_to_dof_map, :]

    with dolfinx.io.XDMFFile(mesh.comm, "moved_mesh.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)

    # Displace submesh
    with dolfinx.io.XDMFFile(mesh.comm, "submesh.xdmf", "w") as xdmf:
        xdmf.write_mesh(submesh)
    num_parent_nodes = (
        mesh.geometry.index_map().size_local + mesh.geometry.index_map().num_ghosts
    )
    parent_to_sub_node = np.full(num_parent_nodes, -1, dtype=np.int32)
    parent_to_sub_node[node_map] = np.arange(
        submesh.geometry.index_map().size_local
        + submesh.geometry.index_map().num_ghosts,
        dtype=np.int32,
    )
    # Map u_values to the vertices of the mesh
    u_values_at_parent_nodes = u_values[vertex_to_dof_map, :][vertex_to_node_map]
    # Map u_values to the submesh
    u_values_at_sub_nodes = u_values_at_parent_nodes[node_map]
    # Deform mesh
    submesh.geometry.x[:, :tdim] += u_values_at_sub_nodes
    with dolfinx.io.XDMFFile(mesh.comm, "moved_submesh.xdmf", "w") as xdmf:
        xdmf.write_mesh(submesh)

yielding

where the blue cells are the deformed parent mesh, the white the unperturbed submesh and green the perturbed submesh.

[Stein]

I trust your question is answered in @dokken’s above response, so more on a meta level: are you sure this is actually what you want to do?
The whole point of ALE is that you pull-back the moving domain to a static reference domain, and describe the mesh-displacement by its own solution field. This combines the Lagrangian and Eulerian perspectives, and hence ‘ALE’.

I imagine that repeatedly moving mesh-nodes is detrimental to FEniCSx’ performance, as many chached quantities would have to be recomputed every timestep. (Please correct me if I’m wrong.)

See mae-207-fea-for-coupled-problems/fsi/fitted-fsi-example.py at master · david-kamensky/mae-207-fea-for-coupled-problems · GitHub for a nice example (legacy fenics). You’ll find that no mesh displacement operations are used.

[chenyongxin]

@Stein, thanks a lot for your kind reminder. To answer your concern about my implementation of ALE, I think the following two points are equivalent:

1. Solve ALE mapping equation, e.g. (\nabla u, \nabla v)_{\Omega_f^0} on the reference mesh, as what you described. The solution u is the total displacement. One can implement this by either
   1.1 involving Jacobian matrix J to the weak forms to transfer terms between reference static mesh and current moving mesh. An example is Simula lab’s code TurtleFSI. See for instance: GitHub - KVSlab/turtleFSI: Monolithic Fluid-Structure Interaction (FSI) solver, or
   1.2 move mesh backward and forward in one step, and assemble system on different domains, i.e. for ()_{\Omega_f^n} + ()_{\Omega_s^0}, but without J.
2. Solve ALE mapping equation, e.g. (\nabla u, \nabla v)_{\Omega_f^n}, on the current mesh. So the solution u is displacement increment \Delta x in the curret step.

I used 2 in my minimal working example. That may not follow what ALE really describes, hence I would like to call it updated Lagrangian description, see for instance:

A unified-field monolithic fictitious domain-finite element method for...

[chenyongxin]

@Stein

I imagine that repeatedly moving mesh-nodes is detrimental to FEniCSx’ performance, as many chached quantities would have to be recomputed every timestep. (Please correct me if I’m wrong.)

From my perspective, this cost should be marginal as it only executes simple addition on the mesh coordinate. No matter the ALE is implmented by ‘Jacobian formulation’ or move mesh method, the linear system should be the same, hence yields similar cost in time.
@dokken Could you please also comment on this?

[chenyongxin]

Thanks, @dokken. I see that you use vertex (topology) to bridge node (geometry) and DOF (fem). My mistake is because the vertex-dof map of submesh is not always ordered as that of parent mesh.

It also works for my case with the following:

submesh.geometry.x[vertex_to_node_map_s, :tdim] += Delta_x_s.x.array.reshape((-1, tdim))[vertex_to_dof_map_s, :]

A step further: how to deal with 2nd order mesh and its submesh?

@dokken

num_vertices = (
        mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
    )

The above one always returns 3 for either 1st or 2nd order mesh. The middle points on the facets cannot be correctly indexed by those functions for topology.