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

General
2

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

image
image1392×1355 289 KB

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

3 Likes