3d .msh file to Fenics

1

Hi,
I have created 3d domain in gmsh and imported it to fenics with help of @dokken. Now, the problem is when I verified mesh file paraview, only one subdomain is imported . But domain contains 4 subdomains.
I am attaching .geo file and fenics code.

// Gmsh project created by Jørgen S. Dokken 2020
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 1, 1, .285};
Box(2) = {0.4, 0.4, 0.285, .2, .2, .03};
Box(3) = {0., 0, 0.315, 1, 1, .01};
Box(4) = {0.0, 0.0, 0.325, 1, 1, 1};
v() = BooleanFragments{ Volume{1}; Delete;}{ Volume{2,3,4}; Delete; };
Physical Volume("air", 1) ={#v()};
Physical Volume("hBN", 2) ={#v()-1};
Physical Volume("cavity", 3) ={#v()-2};
Physical Volume("substrate", 4) ={#v()-3};
Physical Surface("backgate", 5) = {29};
Physical Surface("patterned_gate", 6) = {18, 23};

import meshio

#msh = meshio.read("revised_3d_model_22.msh")
msh = meshio.read("bokada.msh")
for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]
    elif key == "tetra":
        tetra_data = msh.cell_data_dict["gmsh:physical"][key]
for cell in msh.cells:
    if cell.type == "tetra":
        tetra_cells = cell.data
    elif cell.type == "triangle":
        triangle_cells = cell.data
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},
                         cell_data={"name_to_read":[tetra_data]})
triangle_mesh =meshio.Mesh(points=msh.points,
                           cells=[("triangle", triangle_cells)],
                           cell_data={"name_to_read":[triangle_data]})
meshio.write("plate.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh)

from dolfin import *
set_log_level(LogLevel.ERROR)

mesh = Mesh()
with XDMFFile("plate.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)


mvc2 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("plate.xdmf") as infile:
    infile.read(mvc2, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
2

When I run your code, I encounter an error:

Traceback (most recent call last):

  File "<filename>", line 23, in <module>
    tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},

  File "/home/connor/anaconda3/envs/fenics/lib/python3.8/site-packages/meshio/_mesh.py", line 77, in __init__
    raise ValueError(

ValueError: Incompatible cell data. Cell block 0 has length 1102, but corresponding cell data name_to_read item has length 2029.

This is because your elements are divided among multiple blocks:

In[102]: msh.cells
Out[102]: 
[<meshio CellBlock, type: triangle, num cells: 128>,
 <meshio CellBlock, type: triangle, num cells: 128>,
 <meshio CellBlock, type: triangle, num cells: 90>,
 <meshio CellBlock, type: tetra, num cells: 527>,
 <meshio CellBlock, type: tetra, num cells: 49>,
 <meshio CellBlock, type: tetra, num cells: 351>,
 <meshio CellBlock, type: tetra, num cells: 1102>]

By using a for loop, i.e.

for cell in msh.cells:
    if cell.type == "tetra":
        tetra_cells = cell.data
    elif cell.type == "triangle":
        triangle_cells = cell.data

tetra_cells will only hold the contents of the last CellBlock of type tetra. By using the following, I find that all subdomains are exported:

import meshio

#msh = meshio.read("revised_3d_model_22.msh")
msh = meshio.read("bokada.msh")

triangle_data = msh.cell_data_dict["gmsh:physical"]["triangle"]
tetra_data = msh.cell_data_dict["gmsh:physical"]["tetra"]

triangle_cells = msh.cells_dict["triangle"]
tetra_cells = msh.cells_dict["tetra"]

tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells},
                         cell_data={"name_to_read":[tetra_data]})
triangle_mesh =meshio.Mesh(points=msh.points,
                           cells=[("triangle", triangle_cells)],
                           cell_data={"name_to_read":[triangle_data]})
meshio.write("plate.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh)

plate

1 Like
3

Thank you for your help

1 Like