Problem with GMSH boundary markers

1

Hi,
I have noticed something strange with GMSH boundary markers. They work with Dirichlet BCs but not with ds. I have the following GMSH file made by GMSH 3.0.4:

// Gmsh project created on Sat Aug 20 11:37:10 2022
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {0.000025, 0, 0, 1.0};
//+
Point(3) = {-0.000025, 0, 0, 1.0};
//+
Point(4) = {0, 0.000015, 0, 1.0};
//+
Point(5) = {0, -0.000015, 0, 1.0};
//+
Point(6) = {0.00002, 0, 0, 1.0};
//+
Point(7) = {0.0000225, 0, 0, 1.0};
//+
Point(8) = {0.0000175, 0, 0, 1.0};
//+
Point(9) = {0.00002, 0.0000025, 0, 1.0};
//+
Point(10) = {0.00002, -0.0000025, 0, 1.0};
//+
Ellipse(1) = {4, 1, 2, 3};
//+
Ellipse(2) = {4, 1, 3, 2};
//+
Ellipse(3) = {5, 1, 3, 2};
//+
Ellipse(4) = {5, 1, 2, 3};
//+
Circle(5) = {7, 6, 9};
//+
Circle(6) = {9, 6, 8};
//+
Circle(7) = {8, 6, 10};
//+
Circle(8) = {10, 6, 7};
//+
Line Loop(1) = {2, -3, 4, -1};
//+
Line Loop(2) = {6, 7, 8, 5};
//+
Plane Surface(1) = {1, 2};
//+
Plane Surface(2) = {2};
//+
Physical Line(1) = {2, 1, 4, 3};
//+
Physical Line(2) = {7, 6, 5, 8};
//+
Physical Surface(3) = {1};
//+
Physical Surface(4) = {2};
//+

Then I use dolfin-convert to get the corresponding XML files. Then I use the MWE below to convert them to xdmf and check the boundaries.

from dolfin import *

###############################################################
#Reading the .xml files
###############################################################
mesh = Mesh("Mesh.xml")
Volume = MeshFunction("size_t", mesh, "Mesh_physical_region.xml")
bnd_mesh = MeshFunction("size_t", mesh, "Mesh_facet_region.xml")
###############################################################

###############################################################
#Converting to .xmdf and Saving
###############################################################
xdmf = XDMFFile(mesh.mpi_comm(),"Mesh.xdmf")
xdmf.write(mesh)
xdmf.write(Volume)
xdmf = XDMFFile(mesh.mpi_comm(),"boundaries.xdmf")
xdmf.write(bnd_mesh)
xdmf.close()
###############################################################

###############################################################
#Reading the mesh
###############################################################
mesh= Mesh()
xdmf = XDMFFile(mesh.mpi_comm(), "Mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("Mesh.xdmf") as infile:
    infile.read(mvc, "f")
Volume = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc2 = MeshValueCollection("size_t", mesh, 1)
with XDMFFile("boundaries.xdmf") as infile:
    infile.read(mvc2, "f")
bnd_mesh = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
###############################################################

###############################################################
#Defining the measures
###############################################################
dx = dx(metadata={'quadrature_degree': 2},
        subdomain_data=Volume)
ds = ds(metadata={'quadrature_degree': 2},
        subdomain_data=bnd_mesh)
###############################################################

###############################################################
#Testing Dirichlet
###############################################################
S = FunctionSpace(mesh,'P',1)
f = Function(S)
bc1 = DirichletBC(S, Constant(1), bnd_mesh,1)
bc2 = DirichletBC(S, Constant(2), bnd_mesh,2)
bcs = [bc1,bc2]
for bc in bcs:
    bc.apply(f.vector())
File('test/f.pvd')<<f
###############################################################


###############################################################
#Testing ds
###############################################################
R = FunctionSpace(mesh,'R',0)
a_1 = project(Constant(1),R)
Domain_vol = project(assemble(a_1*ds(1)),R)
print(Domain_vol.vector().get_local())

R = FunctionSpace(mesh,'R',0)
a_1 = project(Constant(1),R)
Domain_vol = project(assemble(a_1*ds(2)),R)
print(Domain_vol.vector().get_local())
###############################################################

This does the value assignments using the Dirichlet BCs for both boundaries correctly. Also the first assemble runs without any issues. But for the second one I get the following error:

*** Error:   Unable to assemble system.
*** Reason:  expected a linear form for L.
*** Where:   This error was encountered inside SystemAssembler.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  f0a90b8fffc958ed612484cd3465e59994b695f9
*** -------------------------------------------------------------------------

I am certain something is wrong with ds, but I cannot figure it out. I would really appreciate it if you could help me with this.

1 Like
2

I found the answer to my problem. So I am going to explain it here in case someone is interested.
It turned out there was nothing wrong with GMSH or conversions. The problem was with the defined measures. Outer boundaries can be integrated over using a ds measure. However, for interior boundaries, you need to use dS. So the correct measure definition for the above code looks something like this:

dx = Measure('dx', domain=mesh, subdomain_data=Volume)
ds = Measure('ds', domain=mesh, subdomain_data=bnd_mesh, subdomain_id = 1)
ds_int = Measure('dS', domain=mesh, subdomain_data=bnd_mesh, subdomain_id = 2)

And then the following codes give the right answers:

R = FunctionSpace(mesh,'R',0)
a_1 = project(Constant(1),R)
Domain_vol = project(assemble(a_1*ds(1)),R)
print(Domain_vol.vector().get_local())

R = FunctionSpace(mesh,'R',0)
a_1 = project(Constant(1),R)
Domain_vol = project(assemble(a_1*ds_int(2)),R)
print(Domain_vol.vector().get_local())