Nanson's formula + quadrature point cloud visualization

example.py

from mpi4py import MPI
import numpy as np
import dolfinx
import ufl
import basix.ufl

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (2,)))
u = dolfinx.fem.Function(V)
u.interpolate(lambda x: (x[0] + 0.1 * np.sin(2 * np.pi * x[1]), np.zeros(x.shape[1])))


# Nanson's formula
x = ufl.SpatialCoordinate(mesh)
I = ufl.Identity(2)
F = I + ufl.grad(u)
J = ufl.det(F)
n = ufl.FacetNormal(mesh)
nh = J * ufl.inv(F.T) * n


mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)


def move_to_facet_quadrature(ufl_expr, mesh, sub_facets, scheme="default", degree=1):
    fdim = mesh.topology.dim - 1
    # Create submesh
    bndry_mesh, entity_map, _, _ = dolfinx.mesh.create_submesh(mesh, fdim, sub_facets)
    # Create quadrature space on submesh
    q_el = basix.ufl.quadrature_element(
        bndry_mesh.basix_cell(), ufl_expr.ufl_shape, scheme, degree
    )
    Q = dolfinx.fem.functionspace(bndry_mesh, q_el)

    # Compute where to evaluate expression per submesh cell
    try:
        integration_entities = dolfinx.fem.compute_integration_domains(
            dolfinx.fem.IntegralType.exterior_facet, mesh.topology, entity_map, fdim
        )
    except TypeError:
        integration_entities = dolfinx.fem.compute_integration_domains(
            dolfinx.fem.IntegralType.exterior_facet, mesh.topology, entity_map
        )

    compiled_expr = dolfinx.fem.Expression(ufl_expr, Q.element.interpolation_points())

    # Evaluate expression
    q = dolfinx.fem.Function(Q)
    q.x.array[:] = compiled_expr.eval(mesh, integration_entities).reshape(-1)
    return q


nh_facet = move_to_facet_quadrature(nh, mesh, exterior_facets)


import scifem

with scifem.xdmf.XDMFFile("facets_on_deformed_boundary.xdmf", [nh_facet]) as xdmf:
    xdmf.write(0.0)

# 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 = scifem.vertex_to_dofmap(V)
tdim = mesh.topology.dim
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)

Move submesh with respect to parent mesh: https://fenicsproject.discourse.group/t/moving-submesh-diverges-from-its-designated-trajectory-in-a-two-phase-problem/16792/2?u=dokken