Unexpected behavior of surface integral at internal boundary

General
1

Hi,

I am working with dolfinx 0.6.0 installed through conda and presently struggling with unexpected behavior of the surface integral over internal boundaries for my use case. I know that ‘ds’ is for external boundaries while ‘dS’ is for internal ones. However, puzzlingly I am getting a non-zero integral with ‘ds’ on an internal boundary (‘dS’ being zero then).

I am not simplifying the geometry for the MWE because then I do not see this behavior. I would appreciate if someone willing to help executes the MWE below with the mesh file from this link placed in the same folder.

import dolfinx, gmsh, petsc4py, mpi4py, ufl, time
from dolfinx.io import gmshio

gmsh_file = 'perm_magnet.msh'

gmsh.initialize()
gmsh.model.add("modelo_1")
gmsh.option.setNumber("General.Terminal", 1)
in_geom = gmsh.merge("perm_magnet.msh")
mpi_rank = 0 
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, mpi_rank)

internal_surf_tag = 17 # internal surface, seen visually in gmsh

facet_tag = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1, ft.indices, ft.values)
ds = ufl.Measure('ds', domain=mesh, subdomain_data=facet_tag)
dS = ufl.Measure('dS', domain=mesh, subdomain_data=facet_tag)
n = ufl.FacetNormal(mesh)
flux_ds = dolfinx.fem.form(1 * ds(internal_surf_tag))
flux_dS = dolfinx.fem.form(1 * dS(internal_surf_tag))
print("Integral using ds:", dolfinx.fem.assemble.assemble_scalar(flux_ds))
print("Integral using dS:", dolfinx.fem.assemble.assemble_scalar(flux_dS))

I will sincerely appreciate help finding the mistake.

Thanks!

1 Like
2

Upon further look, turns out that I did not fragment the two adjacent objects in gmsh which caused this behavior.