In this answer of you, you made an example where you differentiated the sides. It is all based on a first construction of the entities of the interior_facet integration domain. However, this array is empty for the given gmsh example. Maybe its due to the construction with two boxes, so it is somehow regarded as exterior facet, even though this would still contradict the logic since the side is not neccessarily picked correctly.
Here is the full code, where the above function is copied in the code for debugging:
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 = 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')
mesh, cell_tags, facet_tags = dolfinx.io.gmshio.read_from_msh("doublebox.msh", MPI.COMM_WORLD, 0, gdim=3)
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)
expr = ufl.dot(ufl.grad(u), n)
# for renaming due to copying the function
facets = facet_for_evaluation
meshdom = mesh
# function
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 of function
coordinates = coord_expr.eval(mesh, boundary_entities)
values =bndry_expr.eval(mesh, boundary_entities)
# handling of interior facets
interior_indices = dolfinx.fem.compute_integration_domains(
dolfinx.fem.IntegralType.interior_facet, mesh.topology, facets,
mesh.topology.dim - 1)
#print(coordinates)
#print(values)
The interior_indices variable gives an empty array. Hence, the picking of the side is impossible. The handling as exterior facets even though the facet with tag 901 is the interior connection of the two created cubes.