First of all thank you @dokken for your prompt and detailed response and also I have to apologize for not preparing my MWE with the necessary care.
I have prepared another example below. The problem occurs when integrating over subdomains. Integrals over the whole domain remain unaffected.
As in your “subdomains” tutorial I prepare a domain with two physical regions. Before the refinement of the original mesh the integral yields 0.5 in each case as expected. After calling refine the same integral yields wildly different results (0.122 and 0.124).
Info : Meshing 1D...
Info : [ 0%] Meshing curve 1 (Line)
Info : [ 20%] Meshing curve 2 (Line)
Info : [ 30%] Meshing curve 3 (Line)
Info : [ 50%] Meshing curve 4 (Line)
Info : [ 60%] Meshing curve 5 (Line)
Info : [ 80%] Meshing curve 6 (Line)
Info : [ 90%] Meshing curve 7 (Line)
Info : Done meshing 1D (Wall 0.000915329s, CPU 0.001274s)
Info : Meshing 2D...
Info : [ 0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Info : [ 50%] Meshing surface 2 (Plane, Frontal-Delaunay)
Info : Done meshing 2D (Wall 0.00665657s, CPU 0.006914s)
Info : 103 nodes 218 elements
before refine
dx 0.9999999999999996
dx(1) 0.5
dx(2) 0.5
after refine
dx 1.0000000000000009
dx(1) 0.1224846896364883
dx(2) 0.12420233804951691
from mpi4py import MPI
import gmsh
import numpy as np
from dolfinx import __version__, git_commit_hash
from dolfinx.mesh import refine,create_rectangle
from dolfinx.fem import assemble_scalar, form
from dolfinx.io.gmshio import model_to_mesh
from ufl import Measure
# https://jsdokken.com/dolfinx-tutorial/chapter3/subdomains.html
gmsh.initialize()
proc = MPI.COMM_WORLD.rank
top_marker = 2
bottom_marker = 1
left_marker = 1
if proc == 0:
# We create one rectangle for each subdomain
gmsh.model.occ.addRectangle(0, 0, 0, 1, 0.5, tag=1)
gmsh.model.occ.addRectangle(0, 0.5, 0, 1, 0.5, tag=2)
# We fuse the two rectangles and keep the interface between them
gmsh.model.occ.fragment([(2, 1)], [(2, 2)])
gmsh.model.occ.synchronize()
# Mark the top (2) and bottom (1) rectangle
top, bottom = None, None
for surface in gmsh.model.getEntities(dim=2):
com = gmsh.model.occ.getCenterOfMass(surface[0], surface[1])
if np.allclose(com, [0.5, 0.25, 0]):
bottom = surface[1]
else:
top = surface[1]
gmsh.model.addPhysicalGroup(2, [bottom], bottom_marker)
gmsh.model.addPhysicalGroup(2, [top], top_marker)
# Tag the left boundary
left = []
for line in gmsh.model.getEntities(dim=1):
com = gmsh.model.occ.getCenterOfMass(line[0], line[1])
if np.isclose(com[0], 0):
left.append(line[1])
gmsh.model.addPhysicalGroup(1, left, left_marker)
gmsh.model.mesh.generate(2)
msh, ct, ft = model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
gmsh.finalize()
dx = Measure("dx", domain=msh, subdomain_data=ct)
print("before refine")
print("dx", assemble_scalar(form(1*dx)))
print("dx(1)", assemble_scalar(form(1*dx(1))))
print("dx(2)", assemble_scalar(form(1*dx(2))))
msh2 = refine(msh)
print("after refine")
print("dx", assemble_scalar(form(1*dx)))
print("dx(1)", assemble_scalar(form(1*dx(1))))
print("dx(2)", assemble_scalar(form(1*dx(2))))