Cant read .msh file generated from .geo file in Fenics

Errors
1

Hello, everyone,

I generated a .geo file by myself, transformed it in a .msh file and tried to plot it in Fenics. But got an error when I plotted it. Is anyone know what is the reason inside and how to fix it? These are my code and the error I got.

from dolfin import *     
import meshio
import matplotlib.pyplot as plt

# os.system('gmsh -2 -format msh2 test.geo -o test.msh')
msh = meshio.gmsh.read("test.msh")
plot(msh)
plt.savefig("msh.png")
File("mesh.pvd") << msh

image
image1143×509 18.9 KB

This is part of my .geo file

Point(1) = {0,0.12949856476352284,0,1};
Point(2) = {0.0015080736462709616,0.12917801381452632,0,1};
Point(3) = {0.00816718010823942,0.12917801381452632,0,1};
Point(4) = {0.014680769117044795,0.13056251989836964,0,1};
...
Point(1535) = {1.25,1.75,0,1};
Point(1536) = {1.75,1.75,0,1};
Point(1537) = {0.25,1.75,0,1};
Point(1538) = {0.75,1.75,0,1};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
...
Line(1572) = {181,162};
Line(1573) = {191,182};
Line(1574) = {207,192};
Line Loop(1) = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208};
...
Line Loop(75) = {182,183,184,185,186,187,188,189,190,1573};
Line Loop(76) = {192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,1574};
Plane Surface(1) = {1,3,5,7,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35};
Plane Surface(2) = {2,36};
Plane Surface(3) = {4,37,38,39};
Plane Surface(4) = {6};
...
Plane Surface(69) = {40};
Plane Surface(70) = {59};
Plane Surface(71) = {61};
Plane Surface(72) = {63};
Physical Surface(1) = {1};
Physical Surface(2) = {2};
Physical Surface(3) = {3};
Physical Surface(4) = {4};
Physical Surface(5) = {5};
Physical Surface(6) = {6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55};
Physical Surface(7) = {56,57,58,59};
Physical Surface(8) = {60,61,62,63,64,65};
Physical Surface(9) = {66,67,68};
Physical Surface(10) = {69,70,71,72};

and how the .geo file looks like this in Gmsh interface.

image
image880×884 15 KB

Thank you very much for any help!

2

You can give this a try:

msh = meshio.read("mymesh.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()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
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)
4

Thanks for your reply. I added this code with you code and got the following error

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("test.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) 

mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
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)

image

And if I change prune_z to False will get this

image
image996×224 7.6 KB

Do you know how to fix it?

Thanks

5

Hi,

Can you post the full Gmsh code please?

1 Like