Surface evaluation with correct normal of gmsh mesh

General
, , ,
1

Hello,

I’m working on the post-processing of a 3D electrial field (Nedelec elements) problem. I use perfectly matched layers on the surrounding, hence I want to evaluate the electrical field on an interior surface. The mesh is generated with gmsh and physical are defined within. Therefore, the selection of the corresponding surface for evaluation is possible with the defined facet tag.

With the surface evaluation of a simple problem provided by Dokken in a similar issue and the interior surface handling I generated the following code:

from mpi4py import MPI
import dolfinx
import ufl
import basix
import numpy as np
import numpy.typing as npt
import gmsh
import sys

def evaluate_expression_at_facet_midpoints(
    expr: ufl.core.expr.Expr, mesh: dolfinx.mesh.Mesh, facets: npt.NDArray[np.int32]
) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.inexact]]:
    # Get coordinate element for facets of cell
    facet_types = basix.cell.subentity_types(mesh.basix_cell())[fdim]
    unique_facet_types = np.unique(facet_types)
    assert len(unique_facet_types) == 1, (
        "All facets must have the same type for interpolation."
    )
    facet_type = facet_types[0]
    facet_geom = basix.cell.geometry(facet_type)
    facet_midpoint = np.sum(facet_geom, axis=0) / facet_geom.shape[0]

    # Compile code for evaluation
    bndry_expr = dolfinx.fem.Expression(expr, facet_midpoint)

    # Map facet indices to (cell, local_facet) pairs
    try:
        boundary_entities = dolfinx.fem.compute_integration_domains(
            dolfinx.fem.IntegralType.exterior_facet, mesh.topology, facets
        )
    except TypeError:
        boundary_entities = dolfinx.fem.compute_integration_domains(
            dolfinx.fem.IntegralType.exterior_facet,
            mesh.topology,
            facets,
            mesh.topology.dim - 1,
        )

    x = ufl.SpatialCoordinate(mesh)
    coord_expr = dolfinx.fem.Expression(x, facet_midpoint)

    return coord_expr.eval(mesh, boundary_entities), bndry_expr.eval(
        mesh, boundary_entities
    )

gmsh.initialize(sys.argv)
dim=3

box = gmsh.model.occ.addBox(0.0, 0.0, 0.0, 1.0, 1.0 , 1.0)
box2 = gmsh.model.occ.addBox(1.0, 0.0, 0.0, 1.0, 1.0 , 1.0)
gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(dim, [box], tag=801, name='box')
gmsh.model.addPhysicalGroup(dim, [box2], tag=802, name='box2')

surfaces = gmsh.model.getEntities(dim=dim-1)
comlist = []
othersurf = []
for surf in surfaces:
    com = gmsh.model.occ.getCenterOfMass(surf[0], surf[1])
    if com in comlist:# to avoid double tags of overlapping bodies
        continue
    comlist.append(com)
    if np.allclose(com, [1.0, 0.5, 0.5]):
        gmsh.model.addPhysicalGroup(surf[0], [surf[1]], tag=901, name='evalsurf')
    else:
        othersurf.append(surf[1])
dummy = gmsh.model.addPhysicalGroup(dim-1, othersurf, tag=902, name='other')

msize = 0.1
points = gmsh.model.getEntities(dim=dim-3)
for p in points:
    coords = gmsh.model.getValue(p[0], p[1], [])
    if np.allclose(coords, [0.0,0.0,0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    elif np.allclose(coords, [1.0,0.0,0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    else:
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)

gmsh.model.getPhysicalGroups()


gmsh.model.mesh.generate(dim-1)
gmsh.model.mesh.generate(dim)
gmsh.model.mesh.refine()

gmsh.write('doublebox.msh')

# read mesh
mesh, cell_tags, facet_tags = dolfinx.io.gmshio.read_from_msh("doublebox.msh", MPI.COMM_WORLD, 0, gdim=3)

# fe problem
V = dolfinx.fem.functionspace(mesh, ("N1curl", 3))
u = dolfinx.fem.Function(V)
u.interpolate(
    lambda x: (x[1] * x[0] ** 2 + x[1] ** 2, np.zeros(x[0].shape), np.zeros(x[0].shape))
)

n = ufl.FacetNormal(mesh)
fdim = mesh.topology.dim - 1

facet_for_evaluation = facet_tags.find(901)
# rename for later use in function
facets = facet_for_evaluation
meshdom = mesh

facet_types = basix.cell.subentity_types(mesh.basix_cell())[fdim]
unique_facet_types = np.unique(facet_types)
assert len(unique_facet_types) == 1, (
    "All facets must have the same type for interpolation."
)
facet_type = facet_types[0]
facet_geom = basix.cell.geometry(facet_type)
facet_midpoint = np.sum(facet_geom, axis=0) / facet_geom.shape[0]


mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
interior_entities = dolfinx.fem.compute_integration_domains(
    dolfinx.fem.IntegralType.interior_facet, mesh.topology, facets,
    mesh.topology.dim - 1)

connected_cells = interior_entities.reshape(-1, 4)[:, [0, 2]]
switch = cell_tags.values[connected_cells[:,0]] < cell_tags.values[connected_cells[:,1]]
if len(switch)> 0 and  any(switch):
    tmp_entities = integration_entities.reshape(-1, 4)
    tmp_entities[switch] = tmp_entities[switch][:, [2, 3, 0, 1]]
    interior_entities = tmp_entities.flatten()

dS_custom = ufl.Measure("dS", domain=mesh, subdomain_data=[(1111, np.array(interior_entities, dtype=np.int32))])
vect = ufl.dot(ufl.grad(u("+")),n("+"))

# dolfinx.fem.form(vect*dS_custom(1111)) not working due to required scalar in integration
v0 = dolfinx.fem.form(vect[0]*dS_custom(1111))
v1 = dolfinx.fem.form(vect[1]*dS_custom(1111))
v2 = dolfinx.fem.form(vect[2]*dS_custom(1111))

# Compile code for evaluation
bndry_expr0 = dolfinx.fem.Expression(v0, facet_midpoint, comm=MPI.COMM_WORLD)
bndry_expr1 = dolfinx.fem.Expression(v1, facet_midpoint, comm=MPI.COMM_WORLD)
bndry_expr2 = dolfinx.fem.Expression(v2, facet_midpoint, comm=MPI.COMM_WORLD)

# Map facet indices to (cell, local_facet) pairs
try:
    boundary_entities = dolfinx.fem.compute_integration_domains(
        dolfinx.fem.IntegralType.interior_facet, mesh.topology, facet_for_evaluation
    )
except TypeError:
    boundary_entities = dolfinx.fem.compute_integration_domains(
        dolfinx.fem.IntegralType.interior_facet,
        mesh.topology,
        facet_for_evaluation,
        mesh.topology.dim - 1,
    )

x = ufl.SpatialCoordinate(mesh)
coord_expr = dolfinx.fem.Expression(x, facet_midpoint)


coordinates2 = coord_expr.eval(mesh, boundary_entities)
values20 = bndry_expr0.eval(mesh, boundary_entities)
values21 = bndry_expr1.eval(mesh, boundary_entities)
values22 = bndry_expr2.eval(mesh, boundary_entities)

print(coordinates2)
print(values20)
print(values21)
print(values22)

There are a couple of issues with that. First of all, the interior_entities in line 120 remain an empty array. So further processing is not useful. However, the defined measure dS_custom is then defined on that region and planned to be multiplied to the gradient projection onto the normal vector, resulting in a moment. The integration requires scalar expressions, hence the moment is split in its components. But when the dolfinx.fem.form is then used in the dolfinx.fem.Expression in line 140 I get a TypeError: <class ‘tuple’>. If the mpi-communication option is not given, it yields the error AttributeError: ‘Form’ object has no attribute '_ufl_is_terminal_' .

Why is the interior_entities variable empty in the first place and how do I handle the custom measure properly?

Without the side-handling of the surface:

vect = ufl.dot(ufl.grad(u),n)

bndry_expr = dolfinx.fem.Expression(vect, facet_midpoint, comm=MPI.COMM_WORLD)

the program is running through, though the output is empty due to the empty array interior_entities.

2

The issue quite clear.
You would like to treat an “interior” facet as an exterior facet, to get a “one-sided” evaluation at that boundary.

There are two ways of handling this.

by finding the local (cell, local_facet_index)-pairs that correspond to one side of your integral.
This is for instance illustrated in: Add one-sided integration example as example of custom integration · Issue #158 · jorgensd/dolfinx-tutorial · GitHub
with the code:

# Compute integration entities for a single facet of a cell.
# Each facet is represented as a tuple (cell_index, local_facet_index), where cell_index is local to process
# local_facet_index is the local indexing of a facet for a given cell
integration_entities = []
for i, facet in enumerate(facets_to_integrate):
    # Only loop over facets owned by the process to avoid duplicate integration
    if facet >= facet_map.size_local:
        continue
    # Find cells connected to facet
    cells = f_to_c.links(facet)
    # Get value of cells
    marked_cells = ct.values[cells]
    # Get the cell marked with 2
    correct_cell = np.flatnonzero(marked_cells == 2)

    assert len(correct_cell) == 1
    # Get local index of facet
    local_facets = c_to_f.links(cells[correct_cell[0]])
    local_index = np.flatnonzero(local_facets == facet)
    assert len(local_index) == 1

    # Append integration entities
    integration_entities.append(cells[correct_cell[0]])
    integration_entities.append(local_index[0])
3

The function evaluate_expression_at_facet_midpoints is not used in the code and was by mistake included, sorry. In reality, the dolfinx.fem.IntegralType.interior_facet is used.

I adapted the procedure of the git issue, resulting in

from mpi4py import MPI
import dolfinx
import ufl
import basix
import numpy as np
import numpy.typing as npt
import gmsh
import sys


gmsh.initialize(sys.argv)
dim=3

box = gmsh.model.occ.addBox(0.0, 0.0, 0.0, 1.0, 1.0 , 1.0)
box2 = gmsh.model.occ.addBox(1.0, 0.0, 0.0, 1.0, 1.0 , 1.0)
gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(dim, [box], tag=801, name='box')
gmsh.model.addPhysicalGroup(dim, [box2], tag=802, name='box2')

surfaces = gmsh.model.getEntities(dim=dim-1)
comlist = []
othersurf = []
for surf in surfaces:
    com = gmsh.model.occ.getCenterOfMass(surf[0], surf[1])
    if com in comlist:# to avoid double tags of overlapping bodies
        continue
    comlist.append(com)
    if np.allclose(com, [1.0, 0.5, 0.5]):
        gmsh.model.addPhysicalGroup(surf[0], [surf[1]], tag=901, name='evalsurf')
    else:
        othersurf.append(surf[1])
dummy = gmsh.model.addPhysicalGroup(dim-1, othersurf, tag=902, name='other')

msize = 1.0
points = gmsh.model.getEntities(dim=dim-3)
for p in points:
    coords = gmsh.model.getValue(p[0], p[1], [])
    if np.allclose(coords, [0.0,0.0,0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    elif np.allclose(coords, [1.0,0.0,0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    else:
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)

gmsh.model.getPhysicalGroups()


gmsh.model.mesh.generate(dim-1)
gmsh.model.mesh.generate(dim)
gmsh.model.mesh.refine()

gmsh.write('doublebox.msh')

# read mesh
mesh, cell_tags, facet_tags = dolfinx.io.gmshio.read_from_msh("doublebox.msh", MPI.COMM_WORLD, 0, gdim=3)

# fe problem
V = dolfinx.fem.functionspace(mesh, ("N1curl", 3))
u = dolfinx.fem.Function(V)
u.interpolate(
    lambda x: (x[1] * x[0] ** 2 + x[1] ** 2, np.zeros(x[0].shape), np.zeros(x[0].shape))
)

n = ufl.FacetNormal(mesh)
fdim = mesh.topology.dim - 1

facet_for_evaluation = facet_tags.find(901)


#
# the adapted code of the git issue
#
f_to_c = mesh.topology.connectivity(fdim, mesh.topology.dim)
c_to_f = mesh.topology.connectivity(mesh.topology.dim, fdim)
# Compute integration entities for a single facet of a cell.
# Each facet is represented as a tuple (cell_index, local_facet_index), where cell_index is local to process
# local_facet_index is the local indexing of a facet for a given cell
integration_entities = []
for i, facet in enumerate(facet_for_evaluation):
    # Only loop over facets owned by the process to avoid duplicate integration
    if facet >= mesh.topology.index_map(fdim).size_local:
        continue
    # Find cells connected to facet
    cells = f_to_c.links(facet)
    # Get value of cells
    marked_cells = cell_tags.values[cells]
    # Get the cell marked with 2
    correct_cell = np.flatnonzero(marked_cells == 801)

    assert len(correct_cell) == 1
    # Get local index of facet
    local_facets = c_to_f.links(cells[correct_cell[0]])
    local_index = np.flatnonzero(local_facets == facet)
    assert len(local_index) == 1

    # Append integration entities
    integration_entities.append(cells[correct_cell[0]])
    integration_entities.append(local_index[0])


#
#
#


facet_types = basix.cell.subentity_types(mesh.basix_cell())[fdim]
unique_facet_types = np.unique(facet_types)
assert len(unique_facet_types) == 1, (
    "All facets must have the same type for interpolation."
)
facet_type = facet_types[0]
facet_geom = basix.cell.geometry(facet_type)
facet_midpoint = np.sum(facet_geom, axis=0) / facet_geom.shape[0]




dS_custom = ufl.Measure("dS", domain=mesh, subdomain_data=[(1111, np.array(integration_entities, dtype=np.int32))])
vect = ufl.dot(ufl.grad(u("+")),n("+"))

# dolfinx.fem.form(vect*dS_custom(1111)) not working due to required scalar in integration
# ve = dolfinx.fem.form(vect*dS_custom(1111))
v0 = dolfinx.fem.form(vect[0]*dS_custom(1111))
v1 = dolfinx.fem.form(vect[1]*dS_custom(1111))
v2 = dolfinx.fem.form(vect[2]*dS_custom(1111))

# Compile code for evaluation
bndry_expr0 = dolfinx.fem.Expression(v0, facet_midpoint, comm=MPI.COMM_WORLD)
bndry_expr1 = dolfinx.fem.Expression(v1, facet_midpoint, comm=MPI.COMM_WORLD)
bndry_expr2 = dolfinx.fem.Expression(v2, facet_midpoint, comm=MPI.COMM_WORLD)

# Map facet indices to (cell, local_facet) pairs
try:
    boundary_entities = dolfinx.fem.compute_integration_domains(
        dolfinx.fem.IntegralType.interior_facet, mesh.topology, facet_for_evaluation
    )
except TypeError:
    boundary_entities = dolfinx.fem.compute_integration_domains(
        dolfinx.fem.IntegralType.interior_facet,
        mesh.topology,
        facet_for_evaluation,
        mesh.topology.dim - 1,
    )

x = ufl.SpatialCoordinate(mesh)
coord_expr = dolfinx.fem.Expression(x, facet_midpoint)


coordinates2 = coord_expr.eval(mesh, boundary_entities)
values20 = bndry_expr0.eval(mesh, boundary_entities)
values21 = bndry_expr1.eval(mesh, boundary_entities)
values22 = bndry_expr2.eval(mesh, boundary_entities)

print(coordinates2)
print(values20)
print(values21)
print(values22)

and still facing the problem, that vect is a vector, hence the multiplication with the measure dS_custom has to be split, realized as variables v0, v1 and v2. The TypeError: <class ‘tuple’>, as described more thorougly above, remains.

Additionally, I made the mesh coarser, resulting in 16 facets on the surface to evaluate. The variable facet_for_evaluation gives exactly 16 results, unfortunately the integration_entities are of length 32, regarding just unique numbers end up in 18 facets. The numbers of both arrays are not equal, not even a single. How comes this difference? Is the whole procedure even required, since the facets can be identified by the tags in the first place and the positive normal is chosen via ufl.dot(ufl.grad(u(“+”)),n(“+”)), which is independent how the facets are found?

Edit: coarser mesh for better interpretation and comparison of tag-find and cell-connected facets.

4

This is exactly what I was trying to get at in the last comment.
A facet, in itself is just described by an index, say j. However, a facet with relation to some cell c_i, is described as a tuple (c_i, l_i) where l_i is the local number of the facet from the viewpoint of the cell.

The code in the issue

Illustrates how to extract one of these tuples, with a consistent orientation/view of the facet from a cell, such that the normal points «outward».

In general, the “+” and “-” restrictions are arbitrary in dolfin (and dolfinx) if not packed in a specific way. This aligns with the general mathematical papers on these methods, as they work on abstract grids, where a specific orientation doesn’t make sense.
What you will note, is that in most variational formulations, one always have tuples of n_F*avg(a)jump(b) or jump(a)jump(b) which means that the orientation doesn’t matter.

Thus, to evaluate something at a facet, in a DOLFINx Expression, one supply a (c_i,l_i) tuple, to ensure that the normals are pointing out of c_i.

5

Okay, I get that. Then the entries in the list integration_entities are understandable and correct. They are not stored as tuples, they are subsequent as [c_1, l_1, c_2, l_2, c_3, l_3, …] for a list of tuples one would need the additional line integration_entities=list(zip(integration_entities[::2], integration_entities[1::2])), which results in [(c_1, l_1), (c_2, l_2), (c_3, l_3), …]. But I think the list of tuples is not desired from the later functions. I just tried it due to the TypeError: <class ‘tuple’> I’m getting from bndry_expr0 = dolfinx.fem.Expression(v0, facet_midpoint, comm=MPI.COMM_WORLD), which I still can’t wrap my head around from where its originating. It still remains and exactly looks like this:

bndry_expr0 = dolfinx.fem.Expression(v0, facet_midpoint, comm=MPI.COMM_WORLD)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File “/usr/lib/petscdir/petsc-complex/lib/python3/dist-packages/dolfinx/fem/function.py”, line 151, in __init__
self._ufcx_expression, module, self._code = jit.ffcx_jit(
^^^^^^^^^^^^^
File “/usr/lib/petscdir/petsc-complex/lib/python3/dist-packages/dolfinx/jit.py”, line 62, in mpi_jit
return local_jit(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^
File “/usr/lib/petscdir/petsc-complex/lib/python3/dist-packages/dolfinx/jit.py”, line 218, in ffcx_jit
raise TypeError(type(ufl_object))
TypeError: <class ‘tuple’>

Do I understand it correctly, that I still need

dS_custom = ufl.Measure(“dS”, domain=mesh, subdomain_data=[(1111, np.array(integration_entities, dtype=np.int32))])

vect = ufl.dot(ufl.grad(u(“+”)),n(“+”))

v0 = dolfinx.fem.form(vect[0]*dS_custom(1111))
v1 = dolfinx.fem.form(vect[1]*dS_custom(1111))
v2 = dolfinx.fem.form(vect[2]*dS_custom(1111))

the positive signs (“+”) in the vect array and to split the product with the measure in scalar quantitites for three expressions to evaluate?

And what’s the TypeError about, what would I need to fix?

6

Let us break this down into two different sections.
One question is asking:

  • If I have a dS (interior facet measure), how do I ensure that the n("+") normal is pointing out of a given cell. For this, you would order your integration entities as for instance done in: presentations/icosahom/code/emi_primal_mixed.py at baa2f0deb68c406cc24d9e016cc7d68ffeff3eac · jorgensd/presentations · GitHub
    This code also shows you how to attach it to the mesh
    ordered_integration_data = compute_interface_data(ct, ft.find(interface_marker))
dGamma = Measure(
  "dS",
  domain=mesh,
  subdomain_data=[(2, ordered_integration_data.flatten())],
  subdomain_id=2,
)
  • If you want to use ds measure (usually known exterior facet, i.e. the integral kernel only knows about one (cell, facet_index) pair), or evaluate an expression on an “interior facet” but in a consistent direction, see how to pack the data at: GitHub · Where software is built

Depending on what version of DOLFINx you are using one should reshape the entities when sending them into expression.eval,
i.e.

    if Version(dolfinx.__version__) > Version("0.9.0"):
        data = expr.eval(mesh, integration_entities)
    else:
        data = expr.eval(mesh, integration_entities.flatten())
        shaped_data = data.reshape(expr.value_size, -1)
7

I’m not even at the expression evaluation process, so the second question is not yet to tackle.

I implemented the referred function. Interestingly, it returns an empty array. The compute_integration_domains(dolfinx.fem.IntegralType.interior_facet, …)-function is non-existant for the gmsh-mesh. There are 384 cells in my mesh with 16 facets, so both objects are present and not empty. This is before the ordering, hence it should work as it is. Why is the result of compute_integration_domains empty, whats missing?

8

You need to start providing reproducible examples of what you are trying to do.
You say that you have implemented compute_integration_domains, but you have not provided me the call signature you are using (i.e. what goes into this function in your code)?
The compute integration domains function I referred to earlier, i.e.

Are you using the one I provided in the link:

or are you using the one that is part of the DOLFINx API.

So let us assume that this is the code you are using:

from mpi4py import MPI
import dolfinx
import ufl
import basix
import numpy as np
import numpy.typing as npt
import gmsh
import sys


gmsh.initialize(sys.argv)
dim = 3

box = gmsh.model.occ.addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0)
box2 = gmsh.model.occ.addBox(1.0, 0.0, 0.0, 1.0, 1.0, 1.0)
gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(dim, [box], tag=801, name="box")
gmsh.model.addPhysicalGroup(dim, [box2], tag=802, name="box2")

surfaces = gmsh.model.getEntities(dim=dim - 1)
comlist = []
othersurf = []
for surf in surfaces:
    com = gmsh.model.occ.getCenterOfMass(surf[0], surf[1])
    if com in comlist:  # to avoid double tags of overlapping bodies
        continue
    comlist.append(com)
    if np.allclose(com, [1.0, 0.5, 0.5]):
        gmsh.model.addPhysicalGroup(surf[0], [surf[1]], tag=901, name="evalsurf")
    else:
        othersurf.append(surf[1])
dummy = gmsh.model.addPhysicalGroup(dim - 1, othersurf, tag=902, name="other")

msize = 1.0
points = gmsh.model.getEntities(dim=dim - 3)
for p in points:
    coords = gmsh.model.getValue(p[0], p[1], [])
    if np.allclose(coords, [0.0, 0.0, 0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    elif np.allclose(coords, [1.0, 0.0, 0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    else:
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)

gmsh.model.getPhysicalGroups()


gmsh.model.mesh.generate(dim - 1)
gmsh.model.mesh.generate(dim)
gmsh.model.mesh.refine()

gmsh.write("doublebox.msh")

# read mesh
try:
    mesh, cell_tags, facet_tags = dolfinx.io.gmshio.read_from_msh(
        "doublebox.msh", MPI.COMM_WORLD, 0, gdim=3
    )
except ValueError:
    mesh_data = dolfinx.io.gmshio.read_from_msh(
        "doublebox.msh", MPI.COMM_WORLD, 0, gdim=3
    )
    mesh = mesh_data.mesh
    cell_tags = mesh_data.cell_tags
    facet_tags = mesh_data.facet_tags
# fe problem
V = dolfinx.fem.functionspace(mesh, ("N1curl", 3))
u = dolfinx.fem.Function(V)
u.interpolate(
    lambda x: (x[1] * x[0] ** 2 + x[1] ** 2, np.zeros(x[0].shape), np.zeros(x[0].shape))
)

n = ufl.FacetNormal(mesh)
fdim = mesh.topology.dim - 1

facet_for_evaluation = facet_tags.find(901)

from packaging.version import Version

# Future compatibilty check
integration_args: tuple[int] | tuple
if Version("0.10.0") <= Version(dolfinx.__version__):
    integration_args = ()
else:
    fdim = cell_tags.dim - 1
    integration_args = (fdim,)
idata = dolfinx.cpp.fem.compute_integration_domains(
    dolfinx.fem.IntegralType.interior_facet,
    cell_tags.topology,
    facet_for_evaluation,
    *integration_args,
)
print(facet_for_evaluation, idata)

This returns

[193 194 196 202 205 206 209 257 328 332 333 342 344 349 420 425] []

which indicates that DOLFINx don’t think these are interior facets.
Then, if we look at the GMSH script, it doesn’t seem like you ever enforce that this should be a single interface between the two cubes.

This will make two distinct surfaces in GMSH.
To resolve this, we can use fragment in GMSH, as explained in: Meshes from external sources — FEniCS Workshop
This means that if you change the GMSH mesh generation logic to:

from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
import gmsh
import sys


gmsh.initialize(sys.argv)
dim = 3

box = gmsh.model.occ.addBox(0.0, 0.0, 0.0, 1.0, 1.0, 1.0)
box2 = gmsh.model.occ.addBox(1.0, 0.0, 0.0, 1.0, 1.0, 1.0)
gmsh.model.occ.synchronize()
whole_domain, map_to_input = gmsh.model.occ.fragment([(3, box)], [(3, box2)])
gmsh.model.occ.synchronize()

box = [idx for (dim, idx) in map_to_input[1] if dim == 3]
box2 = [idx for (dim, idx) in map_to_input[0] if dim == 3 and idx not in box]
gmsh.model.addPhysicalGroup(dim, box, tag=801, name="box")
gmsh.model.addPhysicalGroup(dim, box2, tag=802, name="box2")

surfaces = gmsh.model.getEntities(dim=dim - 1)
comlist = []
othersurf = []
for surf in surfaces:
    com = gmsh.model.occ.getCenterOfMass(surf[0], surf[1])
    if com in comlist:  # to avoid double tags of overlapping bodies
        continue
    comlist.append(com)
    if np.allclose(com, [1.0, 0.5, 0.5]):
        gmsh.model.addPhysicalGroup(surf[0], [surf[1]], tag=901, name="evalsurf")
    else:
        othersurf.append(surf[1])
dummy = gmsh.model.addPhysicalGroup(dim - 1, othersurf, tag=902, name="other")

msize = 1.0
points = gmsh.model.getEntities(dim=dim - 3)
for p in points:
    coords = gmsh.model.getValue(p[0], p[1], [])
    if np.allclose(coords, [0.0, 0.0, 0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    elif np.allclose(coords, [1.0, 0.0, 0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    else:
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)

gmsh.model.getPhysicalGroups()


gmsh.model.mesh.generate(dim - 1)
gmsh.model.mesh.generate(dim)
gmsh.model.mesh.refine()

gmsh.write("doublebox.msh")

and run the rest of the code for loading into dolfinx and getting interior facets, you get the output

[362 363 374 375 379 382 445 446 457 459 471 475 479 481 500 544] [136   3 159   3 137   1 160   3 152   1 174   1 153   3 175   3 161   2
 184   1 176   1 198   3 171   3 192   3 172   1 193   3 177   1 199   3
 179   2 201   1 185   1 208   3 200   1 220   3 202   1 222   3 212   2
 231   1 214   1 234   3 230   1 252   3]

which I have described the meaning of in the previous post.

9

Thank you very much, I was not aware of the gmsh requiring the fragment. You were right, I am using the function you provided in the presentation repository and as you wrote in the beginning, I got the empty return of the idata. Since the return of the corrected gmsh is of length 64, I interpret it as both side information. Hence, when I want to get the information of a single side, then my previous function

f_to_c = mesh.topology.connectivity(fdim, mesh.topology.dim)
c_to_f = mesh.topology.connectivity(mesh.topology.dim, fdim)
# Compute integration entities for a single facet of a cell.
# Each facet is represented as a tuple (cell_index, local_facet_index), where cell_index is local to process
# local_facet_index is the local indexing of a facet for a given cell
integration_entities = []
for i, facet in enumerate(facet_for_evaluation):
    # Only loop over facets owned by the process to avoid duplicate integration
    if facet >= mesh.topology.index_map(fdim).size_local:
        continue
    # Find cells connected to facet
    cells = f_to_c.links(facet)
    # Get value of cells
    marked_cells = cell_tags.values[cells]
    # Get the cell marked with 2
    correct_cell = np.flatnonzero(marked_cells == 802)

    assert len(correct_cell) == 1
    # Get local index of facet
    local_facets = c_to_f.links(cells[correct_cell[0]])
    local_index = np.flatnonzero(local_facets == facet)
    assert len(local_index) == 1

    # Append integration entities
    integration_entities.append(cells[correct_cell[0]])
    integration_entities.append(local_index[0])

using the physical volume tag 801 (box) for left side and 802 (box2) for right side evaluation, should be sufficient. Is this correct? Or is this too much, since the surface normals n(“-”) orn(“+”), respectively, are handling this automatically, which is used later in the expression?

Independent on which of the both function I’m using, the described TypeError: <class ‘tuple’> remains. The corresponding measures

dS_idata = ufl.Measure("dS", domain=mesh, subdomain_data=[(2222, np.array(idata, dtype=np.int32))])
dS_integration_entities = ufl.Measure("dS", domain=mesh, subdomain_data=[(1111, np.array(integration_entities, dtype=np.int32))])

are defined, The vector vect = ufl.dot(ufl.grad(u(“+”)),n(“+”)) is of shape (3,), hence requires the split in scalar components
v0 = dolfinx.fem.form(vect[0]*dS_integration_entities(1111))
v1 = dolfinx.fem.form(vect[1]*dS_integration_entities(1111))
v2 = dolfinx.fem.form(vect[2]*dS_integration_entities(1111)),
or
vi0 = dolfinx.fem.form(vect[0]*dS_idata(2222))
vi1 = dolfinx.fem.form(vect[1]*dS_idata(2222))
vi2 = dolfinx.fem.form(vect[2]*dS_idata(2222)).

The TypeError occurs at bndry_expr0 = dolfinx.fem.Expression(v0, facet_midpoint, comm=MPI.COMM_WORLD) and correspondingly bndry_expri0 = dolfinx.fem.Expression(vi0, facet_midpoint, comm=MPI.COMM_WORLD).

Is there something wrong with the scalar split or the definition of the measure(s)?

Here the full code:

from mpi4py import MPI
import dolfinx
import ufl
import basix
import numpy as np
import numpy.typing as npt
import gmsh
import sys

def evaluate_expression_at_facet_midpoints(
    expr: ufl.core.expr.Expr, mesh: dolfinx.mesh.Mesh, facets: npt.NDArray[np.int32]
) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.inexact]]:
    # Get coordinate element for facets of cell
    facet_types = basix.cell.subentity_types(mesh.basix_cell())[fdim]
    unique_facet_types = np.unique(facet_types)
    assert len(unique_facet_types) == 1, (
        "All facets must have the same type for interpolation."
    )
    facet_type = facet_types[0]
    facet_geom = basix.cell.geometry(facet_type)
    facet_midpoint = np.sum(facet_geom, axis=0) / facet_geom.shape[0]

    # Compile code for evaluation
    bndry_expr = dolfinx.fem.Expression(expr, facet_midpoint)

    # Map facet indices to (cell, local_facet) pairs
    try:
        boundary_entities = dolfinx.fem.compute_integration_domains(
            dolfinx.fem.IntegralType.exterior_facet, mesh.topology, facets
        )
    except TypeError:
        boundary_entities = dolfinx.fem.compute_integration_domains(
            dolfinx.fem.IntegralType.exterior_facet,
            mesh.topology,
            facets,
            mesh.topology.dim - 1,
        )

    x = ufl.SpatialCoordinate(mesh)
    coord_expr = dolfinx.fem.Expression(x, facet_midpoint)

    return coord_expr.eval(mesh, boundary_entities), bndry_expr.eval(
        mesh, boundary_entities
    )

gmsh.initialize(sys.argv)
dim=3

box = gmsh.model.occ.addBox(0.0, 0.0, 0.0, 1.0, 1.0 , 1.0)
box2 = gmsh.model.occ.addBox(1.0, 0.0, 0.0, 1.0, 1.0 , 1.0)
gmsh.model.occ.synchronize()
whole_domain, map_to_input = gmsh.model.occ.fragment([(3, box)], [(3, box2)])
gmsh.model.occ.synchronize()

box = [idx for (dim, idx) in map_to_input[1] if dim == 3]
box2 = [idx for (dim, idx) in map_to_input[0] if dim == 3 and idx not in box]

gmsh.model.addPhysicalGroup(dim, box, tag=801, name='box')
gmsh.model.addPhysicalGroup(dim, box2, tag=802, name='box2')

surfaces = gmsh.model.getEntities(dim=dim-1)
comlist = []
othersurf = []
for surf in surfaces:
    com = gmsh.model.occ.getCenterOfMass(surf[0], surf[1])
    if com in comlist:# to avoid double tags of overlapping bodies
        continue
    comlist.append(com)
    if np.allclose(com, [1.0, 0.5, 0.5]):
        gmsh.model.addPhysicalGroup(surf[0], [surf[1]], tag=901, name='evalsurf')
    else:
        othersurf.append(surf[1])
dummy = gmsh.model.addPhysicalGroup(dim-1, othersurf, tag=902, name='other')

msize = 1.0#0.1
points = gmsh.model.getEntities(dim=dim-3)
for p in points:
    coords = gmsh.model.getValue(p[0], p[1], [])
    if np.allclose(coords, [0.0,0.0,0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    elif np.allclose(coords, [1.0,0.0,0.0]):
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)
    else:
        gmsh.model.mesh.setSize([(p[0], p[1])], size=msize)

gmsh.model.getPhysicalGroups()


gmsh.model.mesh.generate(dim-1)
gmsh.model.mesh.generate(dim)
gmsh.model.mesh.refine()

gmsh.write('doublebox.msh')

# read mesh
mesh, cell_tags, facet_tags = dolfinx.io.gmshio.read_from_msh("doublebox.msh", MPI.COMM_WORLD, 0, gdim=3)

# fe problem
V = dolfinx.fem.functionspace(mesh, ("N1curl", 3))
u = dolfinx.fem.Function(V)
u.interpolate(
    lambda x: (x[1] * x[0] ** 2 + x[1] ** 2, np.zeros(x[0].shape), np.zeros(x[0].shape))
)

n = ufl.FacetNormal(mesh)
fdim = mesh.topology.dim - 1

facet_for_evaluation = facet_tags.find(901)


f_to_c = mesh.topology.connectivity(fdim, mesh.topology.dim)
c_to_f = mesh.topology.connectivity(mesh.topology.dim, fdim)
# Compute integration entities for a single facet of a cell.
# Each facet is represented as a tuple (cell_index, local_facet_index), where cell_index is local to process
# local_facet_index is the local indexing of a facet for a given cell
integration_entities = []
for i, facet in enumerate(facet_for_evaluation):
    # Only loop over facets owned by the process to avoid duplicate integration
    if facet >= mesh.topology.index_map(fdim).size_local:
        continue
    # Find cells connected to facet
    cells = f_to_c.links(facet)
    # Get value of cells
    marked_cells = cell_tags.values[cells]
    # Get the cell marked with 2
    correct_cell = np.flatnonzero(marked_cells == 802)

    assert len(correct_cell) == 1
    # Get local index of facet
    local_facets = c_to_f.links(cells[correct_cell[0]])
    local_index = np.flatnonzero(local_facets == facet)
    assert len(local_index) == 1

    # Append integration entities
    integration_entities.append(cells[correct_cell[0]])
    integration_entities.append(local_index[0])

#integration_entities=list(zip(integration_entities[::2], integration_entities[1::2]))
dS_integration_entities = ufl.Measure("dS", domain=mesh, subdomain_data=[(1111, np.array(integration_entities, dtype=np.int32))])

def compute_interface_data(
    cell_tags: dolfinx.mesh.MeshTags, facet_indices: npt.NDArray[np.int32]
) -> npt.NDArray[np.int32]:
    """
    Compute interior facet integrals that are consistently ordered according to the `cell_tags`,
    such that the data `(cell0, facet_idx0, cell1, facet_idx1)` is ordered such that
    `cell_tags[cell0]`<`cell_tags[cell1]`, i.e the cell with the lowest cell marker is considered the
    "+" restriction".

    Args:
        cell_tags: MeshTags that must contain an integer marker for all cells adjacent to the `facet_indices`
        facet_indices: List of facets (local index) that are on the interface.
    Returns:
        The integration data.
    """
    # Future compatibilty check
    integration_args: tuple[int] | tuple
    if 10 <= int(dolfinx.__version__.split(".")[1]):
        integration_args = ()
    else:
        fdim = cell_tags.dim - 1
        integration_args = (fdim,)
    idata = dolfinx.cpp.fem.compute_integration_domains(
        dolfinx.fem.IntegralType.interior_facet,
        cell_tags.topology,
        facet_indices,
        *integration_args,
    )
    # ordered_idata = idata.reshape(-1, 4).copy()
    # switch = (
    #     cell_tags.values[ordered_idata[:, 0]] > cell_tags.values[ordered_idata[:, 2]]
    # )
    # if True in switch:
    #     ordered_idata[switch, :] = ordered_idata[switch][:, [2, 3, 0, 1]]
    return idata#ordered_idata

ordered_integration_data = compute_interface_data(cell_tags, facet_for_evaluation)
dS_idata = ufl.Measure("dS", domain=mesh, subdomain_data=[(2222, np.array(ordered_integration_data, dtype=np.int32))])


facet_types = basix.cell.subentity_types(mesh.basix_cell())[fdim]
unique_facet_types = np.unique(facet_types)
assert len(unique_facet_types) == 1, (
    "All facets must have the same type for interpolation."
)
facet_type = facet_types[0]
facet_geom = basix.cell.geometry(facet_type)
facet_midpoint = np.sum(facet_geom, axis=0) / facet_geom.shape[0]

vect = ufl.dot(ufl.grad(u("+")),n("+"))

# dolfinx.fem.form(vect*dS_custom(1111)) not working due to required scalar in integration
# ve = dolfinx.fem.form(vect*dS_custom(1111))
v0 = dolfinx.fem.form(vect[0]*dS_integration_entities(1111))
v1 = dolfinx.fem.form(vect[1]*dS_integration_entities(1111))
v2 = dolfinx.fem.form(vect[2]*dS_integration_entities(1111))
vi0 = dolfinx.fem.form(vect[0]*dS_idata(2222))
vi1 = dolfinx.fem.form(vect[1]*dS_idata(2222))
vi2 = dolfinx.fem.form(vect[2]*dS_idata(2222))

# Compile code for evaluation
bndry_expr0 = dolfinx.fem.Expression(v0, facet_midpoint, comm=MPI.COMM_WORLD)
bndry_expr1 = dolfinx.fem.Expression(v1, facet_midpoint, comm=MPI.COMM_WORLD)
bndry_expr2 = dolfinx.fem.Expression(v2, facet_midpoint, comm=MPI.COMM_WORLD)
bndry_expr0 = dolfinx.fem.Expression(vi0, facet_midpoint, comm=MPI.COMM_WORLD)
bndry_expr1 = dolfinx.fem.Expression(vi1, facet_midpoint, comm=MPI.COMM_WORLD)
bndry_expr2 = dolfinx.fem.Expression(vi2, facet_midpoint, comm=MPI.COMM_WORLD)

# Map facet indices to (cell, local_facet) pairs
try:
    boundary_entities = dolfinx.fem.compute_integration_domains(
        dolfinx.fem.IntegralType.interior_facet, mesh.topology, facet_for_evaluation
    )
except TypeError:
    boundary_entities = dolfinx.fem.compute_integration_domains(
        dolfinx.fem.IntegralType.interior_facet,
        mesh.topology,
        facet_for_evaluation,
        mesh.topology.dim - 1,
    )

x = ufl.SpatialCoordinate(mesh)
coord_expr = dolfinx.fem.Expression(x, facet_midpoint)


coordinates2 = coord_expr.eval(mesh, boundary_entities)
values20 = bndry_expr0.eval(mesh, boundary_entities)
values21 = bndry_expr1.eval(mesh, boundary_entities)
values22 = bndry_expr2.eval(mesh, boundary_entities)

print(coordinates2)
print(values20)
print(values21)
print(values22)
10

The point of this code, which originates from

clearly states that you should use a “ds” measure, not the “dS” measure, which alleviates the need for a restriction. i.e.

ds_integration_entities = ufl.Measure(
    "ds",
    domain=mesh,
    subdomain_data=[(1111, np.array(integration_entities, dtype=np.int32))],
)

we also modify the code for evaluating the integral to

vect = ufl.dot(ufl.grad(u), n)

# dolfinx.fem.form(vect*dS_custom(1111)) not working due to required scalar in integration
# ve = dolfinx.fem.form(vect*dS_custom(1111))
v0 = dolfinx.fem.form(vect[0] * ds_integration_entities(1111))
v1 = dolfinx.fem.form(vect[1] * ds_integration_entities(1111))
v2 = dolfinx.fem.form(vect[2] * ds_integration_entities(1111))

Similarly, I would have modified the dS integral to

ordered_integration_data = compute_interface_data(cell_tags, facet_for_evaluation)
dS_idata = ufl.Measure(
    "dS",
    domain=mesh,
    subdomain_data=[(2222, np.array(ordered_integration_data, dtype=np.int32))],
)
# other code
#....
#
vi0 = dolfinx.fem.form(vect[0]("+") * dS_idata(2222))
vi1 = dolfinx.fem.form(vect[1]("+") * dS_idata(2222))
vi2 = dolfinx.fem.form(vect[2]("+") * dS_idata(2222))

However, I think you are missing the point completely.
The initial question was how to evaluate an expression at the boundary. Now you have introduced integrals. If we want to evaluate v0, v1, ..., vi0, ... you would call dolfinx.fem.assemble_scalar, i.e.

for v in [v0, v1, v2, vi0, vi1, vi2]:
    print(dolfinx.fem.assemble_scalar(v))

What you actually want to do, is to call

bndry_expr0 = dolfinx.fem.Expression(vect[0], facet_midpoint, comm=MPI.COMM_WORLD)
bndry_expr1 = dolfinx.fem.Expression(vect[1], facet_midpoint, comm=MPI.COMM_WORLD)
bndry_expr2 = dolfinx.fem.Expression(vect[2], facet_midpoint, comm=MPI.COMM_WORLD)

coordinates2 = coord_expr.eval(mesh, np.array(integration_entities, dtype=np.int32))
values20 = bndry_expr0.eval(mesh, np.array(integration_entities, dtype=np.int32))
values21 = bndry_expr1.eval(mesh, np.array(integration_entities, dtype=np.int32))
values22 = bndry_expr2.eval(mesh, np.array(integration_entities, dtype=np.int32))
11

Yes, you are completely right. I missed the point and got confused in between. The ds or dS measure is not required at all. The expression evaluation is what I needed and then even

expr = dolfinx.fem.Expression(vect, facet_midpoint, comm=MPI.COMM_WORLD)
expr.eval(mesh, boundary_entities)

is possible, allowing the vectorial evaluation instead of splitting in scalars.
Thank you very much for the explanations.