Problem of defining global measure on a mesh imported from gmsh

1

I created the following geometry in gmsh:

Screenshot (111)
Screenshot (111)1519×855 54.8 KB

Then I defined the surface and boundary measures on the whole mesh. But upon calculating the length of curve 1, it’s giving me 0. Whereas, when I extracted the submesh 3, defined a boundary and surface measure on the mesh tagged 3, and calculated the length of curve 1, it’s giving me a non-zero number.

Here is the gmsh script:

//+
SetFactory("OpenCASCADE");
//+
Point(3) = {-1.6, 0.3, 0, 1.0};
//+
Point(4) = {-1.4, -0.7, 0, 1.0};
//+
Point(5) = {-1.1, -1.3, 0, 1.0};
//+
Point(6) = {-0.3, -1.7, 0, 1.0};
//+
Point(7) = {0.8, -2.1, 0, 1.0};
//+
Point(8) = {1.7, -1.3, 0, 1.0};
//+
Point(9) = {1.7, -0.4, 0, 1.0};
//+
Point(10) = {1.9, 0.5, 0, 1.0};
//+
Point(11) = {1.6, 1.2, 0, 1.0};
//+
Point(12) = {0.9, 1.5, 0, 1.0};
//+
Point(13) = {0, 1.5, 0, 1.0};
//+
Point(14) = {-0.6, 0.9, 0, 1.0};
//+
Spline(3) = {3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 3};
//+
Circle(4) = {0, 0, 0, 3, 0, 2*Pi};
//+
Curve Loop(1) = {3};

//+
Curve Loop(2) = {4};
//+
Curve Loop(3) = {3};


//+
Hide "*";
//+
Show {
  Point{3}; Point{15}; Curve{3}; Curve{4}; Surface{1}; Surface{2}; 
}
//+
Curve Loop(4) = {3};
//+
Plane Surface(1) = {4};
//+
Curve Loop(5) = {4};
//+
Curve Loop(6) = {3};
//+
Plane Surface(2) = {5, 6};
//+
Physical Curve("C1", 1) = {3};
//+
Physical Curve("C2", 2) = {4};
//+
Physical Surface("S1", 3) = {1};
//+
Physical Surface("S2", 4) = {2};

The FEniCS script is:

from dolfin import *     
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


import numpy
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

import meshio
msh = meshio.read("MWE.msh") # Importing mesh from gmsh

triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh)
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)
dx_custom = Measure("dx", domain=mesh, subdomain_data=cf)

print(assemble(Constant(1.0)*ds_custom(1)))

#-------------------------------------------------------------------Extracting mesh with tag 3-------------------------------------------

mesh_1=SubMesh(mesh,cf,3)
boundary_marker_m1 = MeshFunction("size_t", mesh_1, mesh_1.topology().dim()-1, 0)
ncells1 = MeshFunction("size_t", mesh_1, mesh_1.topology().dim())
surface_marker_m1 = MeshFunction("size_t", mesh_1, mesh_1.topology().dim(), 0)
vmap1 = mesh_1.data().array("parent_vertex_indices", 0)
cmap1 = mesh_1.data().array("parent_cell_indices", mesh_1.topology().dim())

        
n = 0
for c in cells(mesh_1):
  parent_cell = Cell(mesh, cmap1[c.index()])
  surface_marker_m1.array()[c.index()] = cf.array()[parent_cell.index()]
  for f in facets(parent_cell):
    for g in facets(c):
      g_vertices = vmap1[g.entities(0)]
      if set(f.entities(0)) == set(g_vertices):
        boundary_marker_m1.array()[g.index()] = mf.array()[f.index()]
    n=n+1
    
ds_m1 = Measure('ds', domain=mesh_1, subdomain_data=boundary_marker_m1)
dx_m1 = Measure('dx', domain=mesh_1, subdomain_data=surface_marker_m1)

print(assemble(Constant(1)*ds_m1(1)))

which gives me the outputs:

0.0
11.425079060142222

But these numbers should be the same since they are the length of the same boundary curve.

I am confused as to why it is not registering the global definition of measure and accepting the local submesh measures. My gmsh script is not producing any error when executed.

Please let me know if you have any questions.

2

'ds' integrates only over external boundaries. 'dS' integrates over all facets. When considering your submission, 1 is an external boundary, and thus the integration succeeds. When considering your entire mesh, there are no external boundaries marked 1, and thus you are integrating over the null set; hence the integration gives 0.

Changing 'ds' to 'dS' in your Measures should give the correct result.

2 Likes
3

@conpierce8 , Ah! That’s interesting. I didn’t know that.
Thanks a lot. It works.

1 Like