Error in extracting surface and boundary markers for a compiled subdomain from its parent mesh

1

Here is the problem. I have a geometry imported from gmsh. There are 4 subdomains. Let’s call them mesh1, mesh2, mesh3 and mesh4. After importing the mesh, I define my boundary and surface markers. Then I use the surface markers to create a compiled subdomain of mesh2 + mesh3 + mesh4. The problem is creating surface and boundary measures for the compiled subdomain using surface and boundary markers of the parent mesh.

The geometry looks like

Screenshot (116)
Screenshot (116)1352×761 67 KB 

# Importing mesh from gmsh and defining surface and boundary markers
msh = meshio.read("mesh.msh")

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)

# Compiling subdomains

combined_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
combined_subdomains.array()[cf.array()==5] = 1 
combined_subdomains.array()[cf.array()==6] = 1
combined_subdomains.array()[cf.array()==7] = 1
mesh_ima = SubMesh(mesh, combined_subdomains, 1)

# Defining measures for the compiled subdomain

mvc_ima = MeshValueCollection("size_t", mesh_ima, mesh_ima.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mesh_ima)
   infile.read(mvc_ima, "name_to_read")
cf_ima = cpp.mesh.MeshFunctionSizet(mesh_ima, mvc_ima)

mvc_ima = MeshValueCollection("size_t", mesh_ima, mesh_ima.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc_ima, "name_to_read")
mf_ima = cpp.mesh.MeshFunctionSizet(mesh_ima, mvc_ima)

ds_custom = Measure("dS", domain=mesh_ima, subdomain_data=mf_ima) 
dx_custom = Measure("dx", domain=mesh_ima, subdomain_data=cf_ima)

I believe that in ds_custom and dx_custom, the domain and subdomain data are for the compiled subdomain, but that seems not to be the case. So is there anything wrong in creating dx_custom and ds_custom?

Let me know if you need any other information.

2

As you are not supplying the geo file (or gmsh python scrupt) that generates the mesh, it is hard to give you any further assistance.

3

Hi @dokken , here is the gmsh script:

SetFactory("OpenCASCADE");


x1 = 0;
x2 = 2*Pi;
n =200;
lc1 = 0.08;
lc2 = 0.08;
lc3 = 0.08;
lc4 = 0.08;

//First interface Points and Lines for the endothelium
For i In {0:n-1}

  x = x1 + (x2 - x1) * i/n;

  r = 1.7217+0.2858*Cos (x)-0.2452*Sin (x)+0.0204*Cos (2*x)+0.0345*Sin (2*x)-0.1876*Cos (3*x)+0.0348*Sin (3*x)-0.0358*Cos (4*x)+0.0329*Sin (4*x)-0.0101*Cos (5*x)+0.1119*Sin (5*x);

  pList1[i] = newp;

  Point (pList1[i]) = {r*Cos (x),r*Sin (x), 0}; //if you want the inner boundary mesh to be more refined change lc1 to a smaller number

EndFor

//Spline creates a continuous boundary using the points above
s1 = newreg;
Spline(s1) = {pList1[{0:n-1}],1};

//Second interface Points and Lines for the IEL
For i In {n:2*n-1}

  x = x1 + (x2 - x1) * i/n;

  r = 2.4656+0.2869*Cos(x)-0.2907*Sin(x)-0.0445*Cos(2*x)-0.0129*Sin(2*x)-0.1528*Cos(3*x)+0.0131*Sin(3*x)-0.0457*Cos(4*x)-0.0665*Sin(4*x)+0.0186*Cos(5*x)-0.0081*Sin(5*x);

  pList2[i] = newp;

  Point (pList2[i]) = {r*Cos (x),r*Sin (x), 0};

EndFor

//Spline creates a continuous boundary using the points above
s2 = newreg;
Spline(s2) = {pList2[{n:2*n-1}],n+1};

//Second interface Points and Lines for the MEDIA-ADVENTITIA interface
For i In {2*n:3*n-1}

  x = x1 + (x2 - x1) * i/n;

  r = 4.9314+0.1441*Cos(x)-0.2625*Sin(x)-0.1752*Cos(2*x)-0.2901*Sin(2*x)+0.1908*Cos(3*x)+0.0297*Sin(3*x)+0.2856*Cos(4*x)-0.2104*Sin(4*x)-0.0537*Cos(5*x)-0.1570*Sin(5*x);

  pList3[i] = newp;

  Point (pList3[i]) = {r*Cos (x),r*Sin (x), 0};

EndFor

//Spline creates a continuous boundary using the points above
s3 = newreg;
Spline(s3) = {pList3[{2*n:3*n-1}],2*n+1};

//Second interface Points and Lines for the ADVENTITIA outermost boundary
For i In {3*n:4*n-1}

  x = x1 + (x2 - x1) * i/n;

  r = 7.9536+0.0443*Cos(x)-0.5983*Sin(x)-0.0620*Cos(2*x)+0.1528*Sin(2*x)+0.4119*Cos(3*x)-0.1249*Sin(3*x)-0.1130*Cos(4*x)-0.2179*Sin(4*x)+0.0910*Cos(5*x)+0.1176*Sin(5*x);

  pList4[i] = newp;

  Point (pList4[i]) = {r*Cos (x),r*Sin (x), 0};

EndFor

//Spline creates a continuous boundary using the points above
s4 = newreg;
Spline(s4) = {pList4[{3*n:4*n-1}],3*n+1};
//+
Physical Curve("C1", 1) = {1};
//+
Physical Curve("C2", 2) = {2};
//+
Physical Curve("C3", 3) = {3};
//+
Physical Curve("C4", 4) = {4};
//+
Curve Loop(1) = {2};
//+
Curve Loop(2) = {1};
//+
Plane Surface(1) = {1, 2};
//+
Curve Loop(3) = {3};
//+
Curve Loop(4) = {2};
//+
Plane Surface(2) = {3, 4};
//+
Curve Loop(5) = {4};
//+
Curve Loop(6) = {3};
//+
Plane Surface(3) = {5, 6};
//+
Physical Surface("S5", 5) = {1};
//+
Physical Surface("S6", 6) = {2};
//+
Physical Surface("S7", 7) = {3};
//+
Curve Loop(7) = {1};
//+
Plane Surface(4) = {7};
//+
Physical Surface("S8", 8) = {4};

I need some guidance on this.