Problem with nested submeshes

dolfinx
,
1

Hello everyone! I have been working recently in a magneto-mechanical problem in which I have several permanent magnets pointing to a solid shere (partially hollow, for BC imposition) which is sensible to magnetic stimuli. The whole scenario, including the magnets, the solid and the surrounding air, has been meshed and tagged. I wanted to tackle the problem in the most efficient way, so I decided to use submeshes, following the guidelines given in Coupling PDEs of multiple dimensions. Therefore, the whole scenario would be the parent_mesh, where the magnetic potential field is computed; the solid_submesh, where the displacement field is solved coupled to the magnetic potential field, and the inner_solid_surface_submesh, where a lagrange multiplier is imposed to force radial-only displacements.

The issue is that NonlinearProblem requires a list of entity_maps that map from the submeshes to the parent mesh and the nested submesh maps from the inner_solid_surface_submesh to the solid_submesh. I have tried to create en entity_map with “dolfinx.mesh.entity_map” that maps from inner_solid_surface_submesh to parent_mesh, but that function does not create the cpp object required for the forms compilation. Is there a way to fix this?

Here is a MWE (still lengthy due to multiphysics coupling):

from dolfinx import fem, default_scalar_type
import dolfinx.fem.petsc
import dolfinx
from dolfinx.common import Timer
from dolfinx.io.gmsh import model_to_mesh
import scifem

import ufl
from mpi4py import MPI

import gmsh
import numpy as np

def get_entity_map(entity_map: dolfinx.mesh.EntityMap, inverse: bool = False) -> np.typing.NDArray[np.int32]:
    """Get an entity map from the sub-topology to the topology.

    Args:
        entity_map: An `EntityMap` object or a numpy array representing the mapping.
        inverse: If `True`, return the inverse mapping.
    Returns:
        Mapped indices of entities.
    """
    sub_top = entity_map.sub_topology
    assert isinstance(sub_top, dolfinx.mesh.Topology)
    sub_map = sub_top.index_map(entity_map.dim)
    indices = np.arange(sub_map.size_local + sub_map.num_ghosts, dtype=np.int32)
    return entity_map.sub_topology_to_topology(indices, inverse=inverse)

def transfer_meshtags_to_submesh(
    mesh: dolfinx.mesh.Mesh,
    entity_tag: dolfinx.mesh.MeshTags,
    submesh: dolfinx.mesh.Mesh,
    vertex_entity_map: dolfinx.mesh.EntityMap,
    cell_entity_map: dolfinx.mesh.EntityMap,
) -> tuple[dolfinx.mesh.MeshTags, np.typing.NDArray[np.int32]]:
    """
    Transfer a meshtag from a parent mesh to a sub-mesh.

    Args:
        mesh: Mesh containing the meshtags
        entity_tag: The meshtags object to transfer
        submesh: The submesh to transfer the `entity_tag` to
        sub_to_vertex_map: Map from each vertex in `submesh` to the corresponding
            vertex in the `mesh`
        sub_cell_to_parent: Map from each cell in the `submesh` to the corresponding
            entity in the `mesh`
    Returns:
        The entity tag defined on the submesh, and a map from the entities in the
        `submesh` to the entities in the `mesh`.

    """

    sub_cell_to_parent = get_entity_map(cell_entity_map, inverse=False)
    sub_vertex_to_parent = get_entity_map(vertex_entity_map, inverse=False)

    num_cells = (
        mesh.topology.index_map(mesh.topology.dim).size_local + mesh.topology.index_map(mesh.topology.dim).num_ghosts
    )
    mesh_to_submesh = np.full(num_cells, -1, dtype=np.int32)
    mesh_to_submesh[sub_cell_to_parent] = np.arange(len(sub_cell_to_parent), dtype=np.int32)

    submesh.topology.create_connectivity(entity_tag.dim, 0)

    num_child_entities = (
        submesh.topology.index_map(entity_tag.dim).size_local + submesh.topology.index_map(entity_tag.dim).num_ghosts
    )
    submesh.topology.create_connectivity(submesh.topology.dim, entity_tag.dim)

    c_c_to_e = submesh.topology.connectivity(submesh.topology.dim, entity_tag.dim)
    c_e_to_v = submesh.topology.connectivity(entity_tag.dim, 0)

    child_markers = np.full(num_child_entities, 0, dtype=np.int32)

    mesh.topology.create_connectivity(entity_tag.dim, 0)
    mesh.topology.create_connectivity(entity_tag.dim, mesh.topology.dim)
    p_f_to_v = mesh.topology.connectivity(entity_tag.dim, 0)
    p_f_to_c = mesh.topology.connectivity(entity_tag.dim, mesh.topology.dim)
    sub_to_parent_entity_map = np.full(num_child_entities, -1, dtype=np.int32)
    for facet, value in zip(entity_tag.indices, entity_tag.values):
        facet_found = False
        for cell in p_f_to_c.links(facet):
            if facet_found:
                break
            if (child_cell := mesh_to_submesh[cell]) != -1:
                for child_facet in c_c_to_e.links(child_cell):
                    child_vertices = c_e_to_v.links(child_facet)
                    child_vertices_as_parent = sub_vertex_to_parent[child_vertices]
                    is_facet = np.isin(child_vertices_as_parent, p_f_to_v.links(facet)).all()
                    if is_facet:
                        child_markers[child_facet] = value
                        facet_found = True
                        sub_to_parent_entity_map[child_facet] = facet
    tags = dolfinx.mesh.meshtags(
        submesh,
        entity_tag.dim,
        np.arange(num_child_entities, dtype=np.int32),
        child_markers,
    )
    tags.name = entity_tag.name
    return tags, sub_to_parent_entity_map

petsc_options = {
    "snes_type": "newtonls",
    "snes_linesearch_type": "basic",
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "snes_atol": 0.0,
    "snes_rtol": 1e-12,
    "snes_max_it": 100,
    "ksp_error_if_not_converged": True,
    "snes_error_if_not_converged": True,
    "snes_monitor": None,
}

def create_mesh(
    comm: MPI.Comm,
    order: int = 1,
    rank: int = 0,
    partitioner = None,
    solid_radius: float = 0.5
) -> dolfinx.mesh.Mesh:
    
    if comm.rank == rank:
        gmsh.initialize()
        gmsh.model.add("Model")
        box = gmsh.model.occ.addBox(-5,-5,-5, 10,10,10)
        outer_sphere = gmsh.model.occ.addSphere(0,0,0, solid_radius)
        inner_sphere = gmsh.model.occ.addSphere(0,0,0, 0.6*solid_radius)
        gmsh.model.occ.synchronize()

        air = gmsh.model.occ.cut([(3, box)], [(3, outer_sphere)], removeObject=True, removeTool=False)
        hollow_solid = gmsh.model.occ.cut([(3, outer_sphere)], [(3, inner_sphere)], removeObject=True, removeTool=False)
        volume_PG = {"air": [box, inner_sphere],
                     "solid": outer_sphere}
        gmsh.model.occ.synchronize()

        box_surfs = gmsh.model.get_boundary([(3, box)], oriented=False)
        box_surfs_bb = [bbs[3]-bbs[0] for bbs in (gmsh.model.get_bounding_box(*inds) for inds in box_surfs)]
        top_surf = np.array(box_surfs)[np.array(box_surfs_bb) <= 1e-6][:,1]
        surface_PG = {"box_surface": top_surf,
                      "solid_inner": gmsh.model.get_boundary([(3, inner_sphere)], oriented=False)[0][1]}

        for i, (name, value) in zip(range(len(volume_PG)), volume_PG.items()):
            gmsh.model.addPhysicalGroup(3, np.unique(value), name=name)
        for name, value in surface_PG.items():
            gmsh.model.addPhysicalGroup(2, np.unique(value), name=name)

        gmsh.model.occ.synchronize()
        ConstantField = gmsh.model.mesh.field.add("MathEval")
        general_value = 2.0
        inside_solid = solid_radius/3.0
        gmsh.model.mesh.field.setString(ConstantField, "F", f"{general_value}")
        solidField = gmsh.model.mesh.field.add("Constant")
        gmsh.model.mesh.field.setNumber(solidField, "VIn",  inside_solid)
        gmsh.model.mesh.field.setNumber(solidField, "VOut", general_value)
        gmsh.model.mesh.field.setNumbers(solidField, "VolumesList", [volume_PG["solid"]])

        DistanceField = gmsh.model.mesh.field.add("Distance")
        distance_list = gmsh.model.get_boundary([(3, volume_PG["solid"])], oriented=False)
        distance_list = [dimtag[1] for dimtag in distance_list]
        gmsh.model.mesh.field.setNumbers(DistanceField, "SurfacesList", distance_list)
        gmsh.model.mesh.field.setNumber(DistanceField, "Sampling", 100)
        ThresholdField = gmsh.model.mesh.field.add("Threshold")
        gmsh.model.mesh.field.setNumber(ThresholdField, "InField", DistanceField)
        gmsh.model.mesh.field.setNumber(ThresholdField, "SizeMin", inside_solid)
        gmsh.model.mesh.field.setNumber(ThresholdField, "SizeMax", general_value)
        gmsh.model.mesh.field.setNumber(ThresholdField, "DistMin", 0.0)
        gmsh.model.mesh.field.setNumber(ThresholdField, "DistMax", 10*solid_radius)

        FieldsRefineMesh = [ConstantField, solidField, ThresholdField]
        BackgroundMesh = gmsh.model.mesh.field.add("Min")
        gmsh.model.mesh.field.setNumbers(BackgroundMesh, "FieldsList", FieldsRefineMesh)

        gmsh.model.mesh.field.setAsBackgroundMesh(BackgroundMesh)
        gmsh.option.setNumber("Mesh.Algorithm", 6)
        gmsh.option.setNumber("Mesh.Algorithm3D", 1)
        gmsh.option.setNumber("Mesh.ElementOrder", order)
        gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 1)
        gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 2)
        gmsh.option.setNumber("Mesh.HighOrderOptimize", 2)
        gmsh.option.setNumber("Mesh.ColorCarousel", 2)
        gmsh.option.setNumber("General.Verbosity", 0)
        gmsh.model.mesh.removeDuplicateNodes()
        gmsh.model.mesh.generate(3)
        mesh_data = model_to_mesh(gmsh.model, comm, rank, gdim=3, partitioner=partitioner)
        gmsh.finalize()
    else:
        mesh_data = model_to_mesh(gmsh.model, comm, rank, gdim=3, partitioner=partitioner)
    return mesh_data, mesh_data.mesh, mesh_data.physical_groups, mesh_data.cell_tags, mesh_data.facet_tags

def assemble_Jacobian_nest(form, funcs, W):
    trials = ufl.TrialFunctions(W)
    Jacobian = []
    for form_i in form:
        Jacobian_i = []
        for func_j, trial_j in zip(funcs, trials):
            Jacobian_i += [ufl.derivative(form_i, func_j, trial_j)]
        Jacobian += [Jacobian_i]
    return Jacobian

def main():

    ####### Problem definition #######
    mesh_comm = MPI.COMM_WORLD
    model_rank = MPI.COMM_WORLD.rank
    TOLERANCE = float(np.finfo(default_scalar_type).resolution)
    # Mesh generation
    with Timer("Mesh generation"):
        print('\nCREATING MESH...')
        file_data, domain, physicalGroups, cell_markers, facet_markers = create_mesh(mesh_comm, 2, model_rank, solid_radius=0.5)
        domain_s, cell_map_s, vertex_map_s, node_map_s = \
            dolfinx.mesh.create_submesh(domain, domain.topology.dim, cell_markers.find(physicalGroups["solid"].tag))
        solid_tags, solid_map = transfer_meshtags_to_submesh(domain, facet_markers, domain_s, vertex_map_s, cell_map_s)
        domain_inter, cell_map_inter, vertex_map_inter, node_map_inter = \
            dolfinx.mesh.create_submesh(domain_s, domain_s.topology.dim-1, solid_tags.find(physicalGroups["solid_inner"].tag))
        inter_to_domain_map = dolfinx.mesh.entity_map(domain.topology, domain_inter.topology, domain_inter.topology.dim, get_entity_map(cell_map_inter, inverse=False))
        quadrature_degree = 4
        dx = ufl.Measure(integral_type='dx',
                         domain=domain,
                         subdomain_data=cell_markers,
                         metadata={'quadrature_degree': quadrature_degree})
        ds = ufl.Measure("ds",
                         domain=domain,
                         subdomain_data=facet_markers,
                         metadata={'quadrature_degree': quadrature_degree})
    
    ####### Function spaces and functions #######
    # Create function spaces
    with Timer("Function spaces"):
        Vu = fem.functionspace(domain_s, ("Lagrange", 2, (domain.geometry.dim,)))
        Vv = fem.functionspace(domain, ("Lagrange", 2))
        Vlm = fem.functionspace(domain_inter, ("Lagrange", 2, (domain.geometry.dim,)))
        Vlmu = scifem.create_real_functionspace(domain)
        W = ufl.MixedFunctionSpace(Vu, Vv, Vlm, Vlmu)

    # Functions:
    (du, dv, dlm, dlmu) = ufl.TestFunctions(W)
    u = fem.Function(Vu, name='displacement')
    v = fem.Function(Vv, name='mag potential')
    lm = fem.Function(Vlm, name='lm')
    lmu = fem.Function(Vlmu, name='lm_u')
    funcs = [u, v, lm, lmu]
    x = ufl.SpatialCoordinate(domain)

    ####### Boundary conditions #######
    bcs = []

    ####### Formulation coupled problem  #######
    # Define derived fields
    I = ufl.variable(ufl.Identity(3))
    F = ufl.variable(ufl.grad(u) + I)
    Finv = ufl.inv(F)
    C = ufl.dot(F.T, F)
    J = ufl.det(F)
    I1c = (J**(-2/3))*ufl.tr(C)
    H = -ufl.grad(v)

    mu = fem.Constant(domain_s, 2000.0)
    k = fem.Constant(domain_s, 1e5)
    psi_solid = (mu/2) * (I1c - 3) + (k/2) * (ufl.ln(J))**2

    # Mechanical term
    mechanical = ufl.diff(psi_solid, F)

    # Maxwell term
    mu0 = fem.Constant(domain, 4*np.pi*1e-7)
    I7 = ufl.inner(ufl.outer(H, H), ufl.dot(Finv, Finv.T))
    dI7dF = -2 * ufl.outer(ufl.dot(Finv.T, H), ufl.dot(ufl.dot(H, Finv), Finv.T))
    maxwell = -(mu0/2) * J * (I7*(Finv.T) + dI7dF)

    # Magnetization term
    mur = fem.Constant(domain, 1.2)
    ms  = fem.Constant(domain, 2.4)
    chi = mur - 1 + TOLERANCE
    I7_sqrt = ufl.sqrt(I7 + TOLERANCE)
    term = ((3*chi)/ms) * I7_sqrt
    magnetization = -((mu0*(ms**2))/(3*chi)) * J * (Finv.T) * (ufl.ln(ufl.sinh(term)) - ufl.ln(term)) +\
                    -((mu0*ms)/2) * J * ((1/ufl.tanh(term))-(1/term)) * (1/I7_sqrt) * dI7dF

    # Total
    Psolid = mechanical + maxwell + magnetization

    form_ = ufl.inner(Psolid, ufl.grad(du)) * dx(physicalGroups["solid"].tag)\
          + ufl.inner(ufl.cross(x, u), dlm) * ds(physicalGroups["solid_inner"].tag)\
          + ufl.inner(ufl.cross(x, lm), du) * ds(physicalGroups["solid_inner"].tag)\
          + ufl.inner(v, dlmu) * dx\
          + ufl.inner(lmu, dv) * dx
    B_dict = {}
    for key, value in physicalGroups.items():
        if value.dim != 3:
            continue
        elif key.find("air") != -1:
            B_term = mu0 * H
        elif key.find("solid") != -1:
            # Maxwell term
            dI7dH = 2 * ufl.dot(ufl.dot(Finv, Finv.T), H)
            maxwell_B = (mu0/2) * J * dI7dH
            # Magnetization term (Mukherjee 2020)
            magnetization_B = ((mu0*ms)/2) * J * ((1/ufl.tanh(term))-(1/term)) * (1/I7_sqrt) * dI7dH
            B_term = maxwell_B + magnetization_B
        else:
            raise ValueError
        B_dict[key] = B_term
        form_ += - ufl.dot(B_term, ufl.grad(dv)) * dx(value.tag)
    
    Bimposed = fem.Constant(domain, np.array([2.0, 0.0, 0.0], dtype=default_scalar_type))
    form_ += - ufl.inner(Bimposed, ufl.grad(dv)) * ds(physicalGroups["box_surface"].tag)

    # EXTRACTING RESIDUAL AND JACOBIAN:
    residual = ufl.extract_blocks(form_)
    if model_rank == 0:
        print(f"- Residual extracted ({len(residual)} blocks)")
    Jacobian = assemble_Jacobian_nest(residual, funcs, W)

    # SET UP NONLINEAR PROBLEM:
    solver_fsi = dolfinx.fem.petsc.NonlinearProblem(residual,
                                                    u=funcs,
                                                    J=Jacobian,
                                                    kind='nest',
                                                    bcs=bcs,
                                                    entity_maps=[cell_map_s, inter_to_domain_map],
                                                    petsc_options_prefix="Example_",
                                                    petsc_options=petsc_options,
    )

    print('\nSTARTING SIMULATION:')
    with Timer("solving"):
        solver_fsi.solve()

if __name__ == "__main__":
    main()

The error is due to the “artificial” entity_map and it is the following:

AttributeError: 'dolfinx.cpp.mesh.EntityMap' object has no attribute '_cpp_object'

Note: Instead of using lagrange multipliers, one could use dolfinx_mpc for imposing no-slip condition on two tangential directions to the surface, needing only the solid_submesh without the nested surface submesh, but I’m getting a “core dumped” error (need to test it deeper).

2

Is there a specific reason of making the domain_inter submesh on the domain_s rather than on the domain? I think this would simplify your setup alot.

To fix the map, you can call:

        inter_to_domain_map = dolfinx.mesh.EntityMap(dolfinx.mesh.entity_map(domain.topology, domain_inter.topology, domain_inter.topology.dim, get_entity_map(cell_map_inter, inverse=False)))

rather than

This minor API change you have to here should really be fixed upstream.

3

Thank you very much @dokken . The main reason behind my setup is the use of functionspaces created on parent and sub-parent meshes. In my case, the Lagrange Multiplier is created at the inter_solid_surface_submesh and needs to operate together with the dispalcement functionspace, which is created at the solid_submesh. For clarification:

  • Parent mesh \rightarrow \Omega_{parent}
  • Submeshes \rightarrow \Gamma_{inter} \subset \Omega_{solid} \subset \Omega_{parent}
  • Magnetic potential field \rightarrow v: \Omega_{parent} \to \mathbb{R}
  • Displacement field \rightarrow u: \Omega_{solid} \to \mathbb{R}^3
  • Lagrange multiplier field \rightarrow lm: \Gamma_{inter} \to \mathbb{R}^3

Regarding your suggestion, I have tried to create both submeshes from the parent mesh, but a problem arises (core dumped) with the form when trying to mix both displacement and Lagrange multipliers in the same expression. Is there a way of mixing fields coming from different submeshes as in my case (i.e., mixing lm and u fields if both come from non-related submeshes)? That would simplify my setup as you mentioned.

4

As illustrated in: Primal mixed-domain formulation — FEniCS-in-the-wild
one can have sub-meshes of a subset of facets on the parent mesh that interacts with a space created on a submesh of the parent grid.

I would suggest having a look at that tutorial, and see if it helps to clarify your issue.

5

Thank you very much for your suggestion @dokken ! The guides on FEniCS-in-the-wild have been very useful. However, I’m stumbling upon an obstacle over and over again related to the solution you told me. Maybe you can easily show me the way to go.

I have discovered a mistake on my formulation. Since I’m definig the displacement field only on a submesh, I should consistently impose maxwell tractions on the boundary for the formulation to be correct. The issue comes here. There are multiple examples on FEniCS-in-the-wild about defining the interface between two submeshes, but in my case I have the parent mesh, \Omega, and the submesh, \Omega_{solid}, which is a subset of \Omega. The example in Mixed mixed-domain formulation is somewhat similar, except for the fact that I have to satisfy the following equality on the submesh boundary, \partial\Omega_{solid}, with v (from both sides) and u and \delta u (from the solid side, i.e. the submesh):

\textbf{P}_{solid}(v_-, u) \cdot \textbf{N} = \textbf{P}_{air}(v_+) \cdot \textbf{N} \;\;\; \textrm{on} \;\;\; \partial\Omega_{solid}
\forall v \in V(\Omega), \forall u \in Q(\Omega_{solid}),

where \textbf{P} is a rank-2 stress tensor and \textbf{N} is the normal to \partial\Omega_{solid}\subset\Omega. Should I be using “ds” or “dS”?

PD: I hope you don’t mind that I’ve added this last question to the same forum thread.

6

You would have to use the dS measure to capture the jump of u or sigma on both sides.
Here is a small example that kind of illustrates what you would like to do:

from mpi4py import MPI
import dolfinx
import numpy as np
import scifem
import ufl

mesh  = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, ghost_mode=dolfinx.mesh.GhostMode.shared_facet)
tol = 10 *np.finfo(mesh.geometry.x.dtype).eps

def sub_cells(x):
    return x[0]<=0.5 + tol

tdim = mesh.topology.dim
num_cells = mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
left_marker = 1
right_marker = 2
marker = np.full(num_cells, left_marker, dtype=np.int32)
marker[dolfinx.mesh.locate_entities(mesh, tdim, sub_cells)] = right_marker
ct = dolfinx.mesh.meshtags(mesh, tdim, np.arange(num_cells, dtype=np.int32), marker)

# Tag exterior boundary and interface
interface = scifem.find_interface(ct, left_marker, right_marker)
mesh.topology.create_connectivity(tdim - 1, tdim)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
facet_map = mesh.topology.index_map(tdim - 1)
num_facets_local = facet_map.size_local + facet_map.num_ghosts
facets = np.arange(num_facets_local, dtype=np.int32)
interface_marker = 4
boundary_marker = 5
marker = np.full_like(facets, -1, dtype=np.int32)
marker[interface] = interface_marker
marker[exterior_facets] = boundary_marker
marker_filter = np.flatnonzero(marker != -1).astype(np.int32)
ft = dolfinx.mesh.meshtags(mesh, tdim - 1, marker_filter, marker[marker_filter])

# Create interface submesh
gamma = scifem.find_interface(ct, left_marker, 2)
Gamma, interface_to_parent, _, _, _ = scifem.extract_submesh(
    mesh, ft, interface_marker
)

# Extract submesh
submesh, interface_map, _, _, _ = scifem.extract_submesh(mesh, ct, left_marker)

# Compute integration data where the lowest marker comes first.
ordered_integration_data = scifem.compute_interface_data(ct, ft.find(interface_marker))
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, )))
Q = dolfinx.fem.functionspace(submesh, ("Lagrange", 1, (mesh.geometry.dim, )))
L = dolfinx.fem.functionspace(Gamma, ("Lagrange", 2, (mesh.geometry.dim, )))
# Restrict anything on "Q" with the correct marker
q_restriction = "+" if left_marker < right_marker else "-"

lambda_ = dolfinx.fem.Constant(mesh, 1.2e4)
mu = dolfinx.fem.Constant(mesh, 0.4)

u = ufl.TrialFunction(V)
def epsilon(u):
    return ufl.sym(
        ufl.grad(u)
    )


def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)



q = ufl.TestFunction(Q)
q_r = q(q_restriction)
lmbda = dolfinx.fem.Function(L)
lmbda.interpolate(lambda x: (1+x[0],1+x[1]))
lmbda.x.scatter_forward()
dS = ufl.Measure("dS", domain=mesh, subdomain_data=[(interface_marker, ordered_integration_data.flatten())], subdomain_id=interface_marker)
n = ufl.FacetNormal(mesh)
term = ufl.dot(lmbda, lmbda)*ufl.dot(ufl.jump(sigma(u), n), q_r)*dS
compiled = dolfinx.fem.form(term, entity_maps=[interface_map, interface_to_parent])

A = dolfinx.fem.petsc.assemble_matrix(compiled)
A.assemble()

print(A.norm(0), A.norm(2))
7

I have followed that example and tried to extrapolate it to my exact problem without succes. Here is a slight modification of the MWE you provided so that it focuses on the exact issue:

from mpi4py import MPI
import dolfinx
import numpy as np
import scifem
import ufl

mesh  = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, ghost_mode=dolfinx.mesh.GhostMode.shared_facet)
tol = 10 *np.finfo(mesh.geometry.x.dtype).eps

def sub_cells(x):
    return x[0]<=0.5 + tol

tdim = mesh.topology.dim
num_cells = mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
left_marker = 1
right_marker = 2
marker = np.full(num_cells, left_marker, dtype=np.int32)
marker[dolfinx.mesh.locate_entities(mesh, tdim, sub_cells)] = right_marker
ct = dolfinx.mesh.meshtags(mesh, tdim, np.arange(num_cells, dtype=np.int32), marker)

# Tag exterior boundary and interface
interface = scifem.find_interface(ct, left_marker, right_marker)
mesh.topology.create_connectivity(tdim - 1, tdim)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
facet_map = mesh.topology.index_map(tdim - 1)
num_facets_local = facet_map.size_local + facet_map.num_ghosts
facets = np.arange(num_facets_local, dtype=np.int32)
interface_marker = 4
boundary_marker = 5
marker = np.full_like(facets, -1, dtype=np.int32)
marker[interface] = interface_marker
marker[exterior_facets] = boundary_marker
marker_filter = np.flatnonzero(marker != -1).astype(np.int32)
ft = dolfinx.mesh.meshtags(mesh, tdim - 1, marker_filter, marker[marker_filter])

# Create interface submesh
gamma = scifem.find_interface(ct, left_marker, 2)
Gamma, interface_to_parent, _, _, _ = scifem.extract_submesh(
    mesh, ft, interface_marker
)

# Extract submesh
submesh, interface_map, _, _, _ = scifem.extract_submesh(mesh, ct, left_marker)

# Compute integration data where the lowest marker comes first.
ordered_integration_data = scifem.compute_interface_data(ct, ft.find(interface_marker))
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, )))
Q = dolfinx.fem.functionspace(submesh, ("Lagrange", 1, (mesh.geometry.dim, )))
# Restrict anything on "Q" with the correct marker
inside = "+" if left_marker < right_marker else "-"
outside = "-" if left_marker < right_marker else "+"

G = dolfinx.fem.Constant(mesh, 0.4)
K = dolfinx.fem.Constant(mesh, 0.4e+3)

dv = ufl.TestFunction(V)
du = ufl.TestFunction(Q)
v = dolfinx.fem.Function(V)
u = dolfinx.fem.Function(Q)

def epsilon(u):
    return ufl.sym(ufl.grad(u))

def sigma(field):
    epsilon_ = epsilon(field)
    vol_epsilon = (1/3) * ufl.tr(epsilon_) * ufl.Identity(len(u))
    dev_epsilon = epsilon_ - vol_epsilon
    return 2 * G * dev_epsilon + 3 * K * vol_epsilon

def sigma_not_submesh(v):
    outside_function = 2.0 * sigma(v) # Lets suppose if similar to the one in submesh, although it is arbitrary
    return outside_function

def sigma_submesh(u, v):
    inside_function = sigma(v) + sigma(u)
    return inside_function

dS = ufl.Measure("dS", domain=mesh, subdomain_data=[(interface_marker, ordered_integration_data.flatten())], subdomain_id=interface_marker)
n = ufl.FacetNormal(mesh)

custom_jump = ufl.dot(sigma_not_submesh(v(outside)), n(outside)) + ufl.dot(sigma_submesh(u(inside), v(inside)), n(inside))
term = ufl.dot(custom_jump, du(inside)) * dS
compiled = dolfinx.fem.form(term, entity_maps=[interface_map])

The thing is that \sigma is not equally defined inside and outside of the submesh and I need to integrate it using the variation of u (\delta u) which is only defined on the submesh. The weak formulation would be:

\int_{\Gamma} ( \sigma_{inside} \cdot \textbf{n} - \sigma_{outside} \cdot \textbf{n}) \; \delta u \; \textrm{d} \Gamma = 0 \;\;\; \forall \delta u

or following ufl convention (what I think it would be):

\int_{\Gamma} ( \sigma_{inside} \cdot \textbf{n}_{inside} + \sigma_{outside} \cdot \textbf{n}_{outside}) \; \delta u \; \textrm{d} \Gamma = 0 \;\;\; \forall \delta u

I tried to implement this in the MWE above, but errors keep arising no matter what. Any suggestion anyone? I must be missing something.

8

I get no errors when running:

from mpi4py import MPI
import dolfinx
import numpy as np
import scifem
import ufl

mesh  = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, ghost_mode=dolfinx.mesh.GhostMode.shared_facet)
tol = 10 *np.finfo(mesh.geometry.x.dtype).eps

def sub_cells(x):
    return x[0]<=0.5 + tol

tdim = mesh.topology.dim
num_cells = mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
left_marker = 1
right_marker = 2
marker = np.full(num_cells, left_marker, dtype=np.int32)
marker[dolfinx.mesh.locate_entities(mesh, tdim, sub_cells)] = right_marker
ct = dolfinx.mesh.meshtags(mesh, tdim, np.arange(num_cells, dtype=np.int32), marker)

# Tag exterior boundary and interface
interface = scifem.find_interface(ct, left_marker, right_marker)
mesh.topology.create_connectivity(tdim - 1, tdim)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
facet_map = mesh.topology.index_map(tdim - 1)
num_facets_local = facet_map.size_local + facet_map.num_ghosts
facets = np.arange(num_facets_local, dtype=np.int32)
interface_marker = 4
boundary_marker = 5
marker = np.full_like(facets, -1, dtype=np.int32)
marker[interface] = interface_marker
marker[exterior_facets] = boundary_marker
marker_filter = np.flatnonzero(marker != -1).astype(np.int32)
ft = dolfinx.mesh.meshtags(mesh, tdim - 1, marker_filter, marker[marker_filter])

# Create interface submesh
gamma = scifem.find_interface(ct, left_marker, 2)
Gamma, interface_to_parent, _, _, _ = scifem.extract_submesh(
    mesh, ft, interface_marker
)

# Extract submesh
submesh, interface_map, _, _, _ = scifem.extract_submesh(mesh, ct, left_marker)

# Compute integration data where the lowest marker comes first.
ordered_integration_data = scifem.compute_interface_data(ct, ft.find(interface_marker))
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (mesh.geometry.dim, )))
Q = dolfinx.fem.functionspace(submesh, ("Lagrange", 1, (mesh.geometry.dim, )))
# Restrict anything on "Q" with the correct marker
inside = "+" if left_marker < right_marker else "-"
outside = "-" if left_marker < right_marker else "+"

G = dolfinx.fem.Constant(mesh, 0.4)
K = dolfinx.fem.Constant(mesh, 0.4e+3)

dv = ufl.TestFunction(V)
du = ufl.TestFunction(Q)
v = dolfinx.fem.Function(V)
u = dolfinx.fem.Function(Q)
u.interpolate(lambda x: (x[0]+x[1], x[0]**2))
v.interpolate(lambda x: (x[1], x[0]**2))

def epsilon(u):
    return ufl.sym(ufl.grad(u))

def sigma(field):
    epsilon_ = epsilon(field)
    vol_epsilon = (1/3) * ufl.tr(epsilon_) * ufl.Identity(len(u))
    dev_epsilon = epsilon_ - vol_epsilon
    return 2 * G * dev_epsilon + 3 * K * vol_epsilon

def sigma_not_submesh(v):
    outside_function = 2.0 * sigma(v) # Lets suppose if similar to the one in submesh, although it is arbitrary
    return outside_function

def sigma_submesh(u, v):
    inside_function = sigma(v) + sigma(u)
    return inside_function

dS = ufl.Measure("dS", domain=mesh, subdomain_data=[(interface_marker, ordered_integration_data.flatten())], subdomain_id=interface_marker)
n = ufl.FacetNormal(mesh)
custom_jump = ufl.dot(sigma_not_submesh(v(outside)), n(outside)) + ufl.dot(sigma_submesh(u(inside), v(inside)), n(inside))
term = ufl.dot(custom_jump, du(inside)) * dS
compiled = dolfinx.fem.form(term, entity_maps=[interface_map])

from petsc4py import PETSc
b = dolfinx.fem.petsc.assemble_vector(compiled)
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
print(b.norm(0), b.norm(1),b.norm(2))

What version of DOLFINx are you using?