Trouble with Multilayer boundaries in 2D mesh

Hello everyone,

I will appreciate it if anyone can help me to understand what is the problem.

I am reading a mesh like this.

I can read the mesh and select the boundaries and I can get physical groups.
I used this function to create subdomains and boundaries:

def make2Dmesh(msh):
    triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
    meshio.write(os.path.join(main_path, "Mesh/mesh.xdmf"), triangle_mesh)

    triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
    meshio.write(os.path.join(main_path, "Mesh/cf.xdmf"), triangle_mesh)

    line_mesh = create_mesh(msh, "line", prune_z=True)
    meshio.write(os.path.join(main_path, "Mesh/mf.xdmf"), line_mesh)

    inmesh = Mesh(MPI.comm_self)#Mesh()
    # Mesh information
    with XDMFFile(MPI.comm_self, os.path.join(main_path, 'Mesh/mesh.xdmf')) as f:
        f.read(inmesh)

    # Boundary line information
    mvc = MeshValueCollection("size_t", inmesh, 1)
    with XDMFFile(MPI.comm_self, os.path.join(main_path, 'Mesh/mf.xdmf')) as f_ph:
        f_ph.read(mvc, 'name_to_read')
    boundaries = MeshFunctionSizet(inmesh, mvc)

    # Material information
    mvc2 = MeshValueCollection("size_t", inmesh, inmesh.geometric_dimension())
    with XDMFFile(MPI.comm_self, os.path.join(main_path, 'Mesh/cf.xdmf')) as f_ph:
        f_ph.read(mvc2, 'name_to_read')
    materials = MeshFunctionSizet(inmesh, mvc2)

    return inmesh, materials, boundaries

The problem that I encountered is that when I use:

dx = Measure('dx', domain=msh, subdomain_data=materials)
ds = Measure('ds', domain=msh, subdomain_data=boundaries)

I have all information regarding “dx” but “ds” only has value on the exterior boundary curve.
I used this to check the values:

for i in range(20):
    if boundaries.where_equal(i):
        print(f"boundary: {i}")
        print("Circumference is: ", assemble(Constant(1) * ds(i)))

for i in range(20):
    if materials.where_equal(i):
        print(f"subdomain: {i}")
        print("Area is: ", assemble(10*dx(I)))

and the output is:

boundary: 8
Circumference is:  0.0
boundary: 9
Circumference is:  0.0
boundary: 10
Circumference is:  0.0
boundary: 11
Circumference is:  0.4083935358272394
****************************************
subdomain: 12
Area is:  0.0031358538980296065
subdomain: 13
Area is:  0.047112009142460345
subdomain: 14
Area is:  0.013351739776037257
subdomain: 15
Area is:  0.0691151254391912

Thank you in advance.

Hi @Ali_Rajabi,

The problem is that you are using ds (with a lowercase ‘s’) instead of dS (with a capital ‘S’). ds integrates only over exterior facets, while dS integrates over all facets, as described in the UFL manual: Form language — Unified Form Language (UFL) 2021.1.0 documentation

Thank you very much. It solved my problem.
I have another question.

I want to obtain the FacetNormal on one of the inner subdomains. Of course, my geometry will not remain this simple. I did something like this, but I am not sure if I did it correctly,

in my code Mesh=msh, Domains = materials,
material - 13 is the dark orange circle and I want “n” on the boundary of the dark orange and cyan circle.

innerMesh = SubMesh(msh, materials, 13)
n = FacetNormal(innerMesh)

Now I should have normal to the boundary of that special subdomain.

Is my procedure correct?
Thank you in advance.

Hi @Ali_Rajabi,

I believe this would work but it is difficult to tell without a minimum working example. Note that FacetNormal can be used without a SubMesh, as detailed in Integrating over an interior surface - #4 by MiroK. You can control which way the FacetNormal points by using the domain markers and a restriction (either ('+') or ('-')).

Thanks for your kind help. I will check the link.