Importing 3d .msh file to fenics

1

Dear sir,
I have imported 3d .msh file to fenics, but I am not sure whether it is correct or wrong. Please verify my code and mention mistakes if any. Here I am attaching the .geo file and a bit of fenics code.
Please help me

// Gmsh project created on Mon Nov  1 20:43:38 2021
SetFactory("OpenCASCADE");
//+
Cylinder(1) = {0, 0, 0, 0, 0, 0.08, 0.03, 2*Pi};
//+
Cylinder(2) = {0, 0, 0.01, 0, 0, 0.07, 0.000895, 2*Pi};
//+
BooleanDifference{ Volume{1}; Delete; }{ Volume{2}; Delete; }
//+
Cylinder(2) = {0, 0, 0.01, 0, 0, 0.07, 0.000895, 2*Pi};
//+
Cylinder(3) = {0, 0, 0.011, 0, 0, 0.069, 0.000595, 2*Pi};
//+
BooleanDifference{ Volume{2}; Delete; }{ Volume{3}; Delete; }
//+
Cone(3) = {0, 0, 0.01, 0, 0, -0.002, 0.000895, 0.0, 2*Pi};
//+
Cylinder(4) = {0, 0, 0.016, 0, 0, 0.001, 0.000595, 2*Pi};
//+
Cylinder(5) = {0, 0, 0.016, 0, 0, 0.001, 0.000470, 2*Pi};
//+
BooleanDifference{ Volume{4}; Delete; }{ Volume{5}; Delete; }
//+
Cylinder(5) = {0, 0, 0.0115, 0, 0, 0.0685, 0.000470, 2*Pi};
//+
Cylinder(6) = {0, 0, 0.011, 0, 0, 0.069, 0.000135, 2*Pi};
//+
BooleanDifference{ Volume{5}; Delete; }{ Volume{6}; Delete; }
//+
Physical Volume(1) = {1};
//+
Physical Volume(2) = {3, 2};
//+
Physical Volume(3) = {4};
//+
Physical Volume(4) = {5};
//+
Physical Surface(5) = {7, 9, 8, 17};
//+
Physical Surface(6) = {29};
//+
Physical Surface(7) = {13, 15, 30, 28, 31, 27, 25, 24, 26};
//+
Physical Curve(8) = {34, 31, 13};
//+
Characteristic Length {6, 5} = 0.3;
//+
Characteristic Length {13, 14, 4, 10, 22, 24, 20, 18, 19, 17, 3, 9, 21, 23} = 0.003;
//+
Characteristic Length {6, 5} = 0.03;
//+
Characteristic Length {13, 4, 14, 10, 22, 24, 18, 20, 17, 19, 23, 21, 9, 3} = 0.003;
//+
Characteristic Length {6, 5} = 0.3;
//+
Characteristic Length {13, 4, 14, 10, 3, 9} = 0.003;
//+
Characteristic Length {23, 21, 17, 19, 18, 20, 24, 22} = 0.0003;
//+
Characteristic Length {13, 14, 4, 10, 22, 24, 20, 18, 19, 17, 23, 21, 9, 3} = 0.0003;
//+
Characteristic Length {6, 5} = 0.03;

and fenics code

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]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh


#msh = meshio.read("3d_mwamodel.msh")
msh = meshio.read("3d_mwa_geometry_new_4.msh")
#print(msh.cell_data_dict)

triangle_mesh = create_mesh(msh, "triangle", True)
tetra_mesh = create_mesh(msh, "tetra", True)
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("mf.xdmf", triangle_mesh) 

from dolfin import *
parameters["allow_extrapolation"] = True 
parameters["form_compiler"]["optimize"] = True 

mesh=Mesh()
mvc = MeshValueCollection("size_t", mesh,3)
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, 2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")   
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc) 
ds = Measure("ds", domain=mesh, subdomain_data=mf)
2

As this post is almost a duplicate of: How tom import 3d domain to fenics
I am going to close this post. Please continue the discussion in the previously created post rather than creating a new one.

1 Like