Marking subdomains in a mesh imported with meshio

General
1

I have a complex geometry. This geometry has not regular shape, so I can not access to different boundaries and domains with specifying the coordinate.
I mesh this geometry in GMSH with defining domains and boundaries. I have two PDEs to solve, PDE1 for domain 1, and PDE2 for domain2. My problem is I can not define the integrand for each subdomain. With this code, I got the following results:
The area of subdomain 1 is 0.000000e+00

which means, the subdomain is not defined.

from fenics import *
import meshio
msh = meshio.read("convertMesh/mchip_2D.msh")
def create_mesh(mesh, cell_type, prune_z=False):
     cells = numpy.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])
     cell_data = numpy.hstack([mesh.cell_data_dict["gmsh:physical"][key]
                         for key in mesh.cell_data_dict["gmsh:physical"].keys() if key==cell_type])
     # Remove z-coordinates from mesh if we have a 2D cell and all points have the same third coordinate
     points= mesh.points
     if prune_z:
         points = points[:,:2]
     mesh_new = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
     return mesh_new
 
 
 line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
 triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
 
 meshio.write("convertMesh/mchip_2D_boundary.xdmf", line_mesh)
 meshio.write("convertMesh/mchip_2D_surface.xdmf", triangle_mesh)
 
 mesh = Mesh()
 with XDMFFile("convertMesh/mchip_2D_surface.xdmf") as infile:
     infile.read(mesh)
 
 mvc_1d = MeshValueCollection("size_t", mesh, 1)
 
 with XDMFFile("convertMesh/mchip_2D_boundary.xdmf") as infile:
     infile.read(mvc_1d, "name_to_read")
 
 mf_1d = cpp.mesh.MeshFunctionSizet(mesh, mvc_1d)
 
 mvc_2d = MeshValueCollection("size_t", mesh, 2)
 
 with XDMFFile("convertMesh/mchip_2D_surface.xdmf") as infile:
     infile.read(mvc_2d, "name_to_read")
     
 mf_2d = cpp.mesh.MeshFunctionSizet(mesh, mvc_2d)
 
 class Porous(SubDomain):
     def inside(self, x, on_boundary):
         return mf_2d == 70
 class Fluid(SubDomain):
     def inside(self, x, on_boundary):
         return mf_2d == 71
 
 domains = MeshFunction("size_t", mesh, mesh.topology().dim(), mesh.domains())
 domains.set_all(10)
 
 porous = Porous()
 porous.mark(domains, 1)
 fluid = Fluid()
 fluid.mark(domains, 2)
 
 mf = MeshFunction('size_t',mesh,2, domains)
 
 dx = Measure('dx', domain=mesh, subdomain_data=domains)
 
 s1 = assemble(Constant(1)*dx(1))
 print("The area of subdomain 1 is {:2e}".format(s1))
2

First, please inspect the various xdmf files visually in Paraview or Visit.
Make sure that the surfaces are tagget corrently in

You are also doing something very weird in

as this is not accessed when you do integration. If you want to use the tags 1 or 2, the subdomain data in dx has to be domains.

I also don’t get why you generate domains, as you are overwriting it in

Please clean up your code, and use proper formatting, i.e.

```python
# Add python code here

```
3

Dear dokken, thanks for your reply.
The surfaces are tagged correctly.
The code is now cleaned up. Please give a look.
Any help can advance my work.
Regards

4

You are still:

  1. Creating domains and overwriting it.
  2. Using mf_2d instead of domains as input to dx.
5

The code is modified as you said, but I the results is 0.00 again.
I don’t know how to define the subdomains for porous and fluid domains.

6

The problem is solved! I had to put in subdomain_data mf_2d:

dx = Measure("dx",domain=mesh, subdomain_data=mf_2d)

This lines are not necessary:

class Porous(SubDomain):
     def inside(self, x, on_boundary):
         return mf_2d == 70
 class Fluid(SubDomain):
     def inside(self, x, on_boundary):
         return mf_2d == 71
 
 domains = MeshFunction("size_t", mesh, mesh.topology().dim(), mesh.domains())
 domains.set_all(10)
 
 porous = Porous()
 porous.mark(domains, 1)
 fluid = Fluid()
 fluid.mark(domains, 2)

and in dx(id), just put the id of domain defined in .msh file:

 s1 = assemble(Constant(1)*dx(70))
 print("The area of subdomain 1 is {:2e}".format(s1))