Integrating on a submesh without mesh tags

General
1

Hi,

I have a 3D geometry where I create a submesh on a 2D facet. I need to integrate on the edges of the submesh but I do not know how it is done without facet tags. All I have on the submesh are cell tags which were facet tags on the original mesh. I do have geometrical markers for the entities I wish to integrate though. So any alternate to using tags with ds/dS will help.

Below is an MWE where I am able to integrate on all edges but not on a select one.

import dolfinx, ufl, gmsh, petsc4py, mpi4py, scifem, basix
from dolfinx.io import gmshio
import numpy as np, matplotlib.pyplot as plt
from dolfinx.fem.petsc import LinearProblem

gdim = 3
gmsh.finalize()
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
box = gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1)
gmsh.model.occ.synchronize()
all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

gmsh.model.mesh.generate(gdim)
model_rank = 0
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)

submesh_idx = (1,)
fdim = mesh.geometry.dim - 1

# create submesh
submesh, sm_to_mesh = dolfinx.mesh.create_submesh(mesh, mesh.geometry.dim-1, np.hstack([ft.find(j) for j in submesh_idx]))[0:2]
facet_imap = mesh.topology.index_map(ft.dim)
num_facets = facet_imap.size_local + facet_imap.num_ghosts
mesh_to_sm = np.full(num_facets, -1)
mesh_to_sm[sm_to_mesh] = np.arange(len(sm_to_mesh))
entity_maps = {submesh: mesh_to_sm}

# transfer facet tags to submesh
# https://fenicsproject.discourse.group/t/transfer-meshtags-to-submesh-in-dolfinx/8952/6
f_map = mesh.topology.index_map(fdim)
all_facets = f_map.size_local + f_map.num_ghosts
# Create array with zeros for all facets that are not marked
all_values = np.zeros(all_facets, dtype=np.int32)
all_values[ft.indices] = ft.values
sub_values = np.empty(len(sm_to_mesh), dtype=np.int32)
for i, entity in enumerate(sm_to_mesh):
    sub_values[i] = all_values[entity]
submesh_ct = dolfinx.mesh.meshtags(submesh, submesh.topology.dim, np.arange(
    len(sm_to_mesh), dtype=np.int32), sub_values)

dx = ufl.Measure("dx", submesh, subdomain_data=ct)
ds = ufl.Measure("ds", submesh, subdomain_data=ft)

# How to integrate only on one (say left) edge?
mesh.comm.allreduce(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(dolfinx.fem.form(1*ds)), op=mpi4py.MPI.SUM))
2

If you have markers on the original mesh, you can use:
scifem/src/scifem/mesh.py at 64b354d319506fce74cb9db706b86a56502a13cb · scientificcomputing/scifem · GitHub to move them from the parent mesh to the submesh.

1 Like