I have read that the support from dolfin-convert is almost nill and will soon be deprecated. This means I have to change the way I read the mesh file created with Gmsh in my codes, but so far I am unable to do it. My goal is to make the transition as smooth as possible. As of now I create a file.geo with Gmsh, convert it to a file.msh using the command gmsh file.geo -format msh2 -3. Then I use dolfin-convert file.msh file.xml which creates 2 files I can read with dolfin, and everything is easy and fine. More precisely, I do subdomains = MeshFunction(“size_t”, mesh, file_physical_region.xml) boundary_parts = MeshFunction(“size_t”, mesh, file_facet_region.xml) This is what I would like to achieve with meshio. From what I have read, meshio supports some features of dolfin-convert regarding the xml file format, but it isn’t quite similar so it is impossible to reproduce what dolfin-convert does in the above example. Furthermore, XDMF is instead the suggested alternative. So, let’s take the following file.geo: SetFactory("OpenCASCADE"); Mesh.RandomFactor=1.0e-6; lc=0.00017; a = 0.002; b = 0.001; c = 0.00003; Point(1) = {-a/2, -b/2, -c/2, lc}; Point(2) = {-a/2, b/2, -c/2, lc}; Point(3) = {a/2, b/2, -c/2, lc}; Point(4) = {a/2, -b/2, -c/2, lc}; Line(1) = {1,2}; Line(2) = {2,3}; Line(3) = {3,4}; Line(4) = {4,1}; Line Loop(5) = {1,2,3,4}; Point(7) = {-a/10 + a/4.4,-b/10-b/3, -c/2, lc}; Point(8) = {-a/10 + a/4.4,b/10-b/3,-c/2, lc}; Point(9) = {a/10 + a/4.4, b/10-b/3,-c/2, lc}; Point(10) = {a/10 + a/4.4,-b/10-b/3,-c/2, lc}; Line(11) = {7,8}; Line(12) = {8,9}; Line(13) = {9,10}; Line(14) = {10,7}; Line Loop(15) = {11,12,13,14}; Plane Surface(16) = {15}; Point(17) = {-a/16+a/4,-b/16+b/3.2, -c/2, lc}; Point(18) = {-a/16+a/4,b/16+b/3.2,-c/2, lc}; Point(19) = {a/16+a/4, b/16+b/3.2,-c/2, lc}; Point(20) = {a/16+a/4,-b/16+b/3.2,-c/2, lc}; Line(21) = {17,18}; Line(22) = {18,19}; Line(23) = {19,20}; Line(24) = {20,17}; Line Loop(25) = {21,22,23,24}; Plane Surface(26) = {25}; Point(27) = {-a/20-a/2.5,-b/3, -c/2, lc}; Point(28) = {-a/20-a/2.5,b/3,-c/2, lc}; Point(29) = {a/20-a/2.5, b/3,-c/2, lc}; Point(30) = {a/20-a/2.5,-b/3,-c/2, lc}; Line(31) = {27,28}; Line(32) = {28,29}; Line(33) = {29,30}; Line(34) = {30,27}; Line Loop(35) = {31,32,33,34}; Plane Surface(36) = {35}; Point(37) = {-a/4,-b/3, -c/2, lc}; Point(38) = {-a/4,b/3,-c/2, lc}; Point(39) = {-a/5, b/3,-c/2, lc}; Point(40) = {-a/5,-b/3,-c/2, lc}; Line(41) = {37,38}; Line(42) = {38,39}; Line(43) = {39,40}; Line(44) = {40,37}; Line Loop(45) = {41,42,43,44}; Plane Surface(46) = {45}; Plane Surface(47)={5,25,15,35,45}; Physical Surface("bottom")={47}; Physical Surface("bot_1")={16}; Physical Surface("bot_2")={26}; Physical Surface("bot_3")={36}; Physical Surface("bot_4")={46}; surfaceVector[] = Extrude {0, 0, c} { Surface{47,26,16,36,46}; Layers{31}; }; Physical Surface("top") = {surfaceVector[0]}; Physical Volume("parallelepiped") = {surfaceVector[1]}; Physical Surface("left") = {surfaceVector[2]}; Physical Surface("back") = {surfaceVector[3]}; Physical Surface("right") = {surfaceVector[4]}; Physical Surface("front") = {surfaceVector[5]}; Physical Surface("top_2") = {surfaceVector[22]}; Physical Volume("stuff_2") = {surfaceVector[23]}; Physical Surface("top_1") = {surfaceVector[28]}; Physical Volume("stuff_1") = {surfaceVector[29]}; Physical Surface("top_3") = {surfaceVector[34]}; Physical Volume("stuff_3") = {surfaceVector[35]}; Physical Surface("top_4") = {surfaceVector[40]}; Physical Volume("stuff_4") = {surfaceVector[41]}; Then we create the file.msh using gmsh file.geo -3. Then we use meshio-convert file.msh file.xdmf. We now have a file.h5 and a file.xdmf. And this is where I am stuck, I do not know how to read the mesh and assign the boundaries as I used to do with the file.xml that was produced with dolfin-convert. I tried: from dolfin import * mesh = Mesh() with XDMFFile("file.xdmf") as infile: infile.read(mesh) But this yields RuntimeError Traceback (most recent call last) in 8 with XDMFFile(“file.xdmf”) as infile: ----> 9 infile.read(mesh) 10 RuntimeError: *** ------------------------------------------------------------------------- *** DOLFIN encountered an error. If you are not able to resolve this issue *** using the information listed below, you can ask for help at *** fenics-support@googlegroups.com *** Remember to include the error message listed below and, if possible, *** include a minimal running example to reproduce the error. *** ------------------------------------------------------------------------- *** Error: Unable to recognise cell type. *** Reason: Unknown value “mixed”. *** Where: This error was encountered inside XDMFFile.cpp. *** Process: 0 *** DOLFIN version: 2018.1.0 *** Git changeset: unknown *** ------------------------------------------------------------------------- So it seems unable to recognize the cell type. I get exactly the same error if I use the deprecated -format msh2 flag in the gmsh command to produce the file.msh. Thus, overall, I do not know how to reproduce the simple commands that work with the file.xml produced with dolfin-convert, using meshio. Note that I had posted a similiar question there but I eventually ran into errors and so this lead to no solution.

I dont really now how to help you further. If you use the docker to install fenics, and then run the following: docker run -ti --rm -v $(pwd):/home/fenics/shared quay.io/fenicsproject/stable:latest sudo pip install meshio lxml h5py I used gmsh 3.0.6 to generate the msh file. Then running the following script: import meshio msh = meshio.read("file.msh") meshio.write("mesh.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells["tetra"]})) meshio.write("mf.xdmf", meshio.Mesh(points=msh.points, cells={"triangle": msh.cells["triangle"]}, cell_data={"triangle": {"name_to_read": msh.cell_data["triangle"]["gmsh:physical"]}})) meshio.write("cf.xdmf", meshio.Mesh( points=msh.points, cells={"tetra": msh.cells["tetra"]}, cell_data={"tetra": {"name_to_read": msh.cell_data["tetra"]["gmsh:physical"]}})) from dolfin import * mesh = Mesh() with XDMFFile("mesh.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) mvc = MeshValueCollection("size_t", mesh, 3) with XDMFFile("cf.xdmf") as infile: infile.read(mvc, "name_to_read") cf = cpp.mesh.MeshFunctionSizet(mesh, mvc) ds_custom = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=12) print(assemble(1*ds_custom)) I get a non-zero output: 3.125e-08

Yes, the tetra_data is a simple list of shape (numcells,). To illustrate this, you could use this example: import numpy as np import pygmsh import meshio def create_mesh_gmsh(): geom = pygmsh.built_in.Geometry() rect = geom.add_rectangle(0.0, 2.0, 0.0, 1.0, 0.0, lcar=0.2) geom.add_physical([rect.line_loop.lines[0], rect.line_loop.lines[2]], 1) geom.add_physical([rect.line_loop.lines[1]], 2) geom.add_physical([rect.line_loop.lines[3]], 3) geom.add_physical([rect.surface], 4) # Generate mesh mesh = pygmsh.generate_mesh(geom, dim=2, prune_z_0=True) cells = np.vstack(np.array([cells.data for cells in mesh.cells if cells.type == "triangle"])) triangle_mesh = meshio.Mesh(points=mesh.points, cells=[("triangle", cells)]) facet_cells = np.vstack(np.array([cells.data for cells in mesh.cells if cells.type == "line"])) facet_data = mesh.cell_data_dict["gmsh:physical"]["line"] facet_mesh = meshio.Mesh(points=mesh.points, cells=[("line", facet_cells)], cell_data={"name_to_read": [facet_data]}) # Write mesh meshio.xdmf.write("mesh.xdmf", triangle_mesh) meshio.xdmf.write("facet_mesh.xdmf", facet_mesh) Here mesh.cells contains 5 CellBlock, four of them with facet data.

Hi everyone, I have followed all the suggestion you are giving using the last version of meshio. However, the code doesn’t recognize 3 of the 4 boundary conditions. (*** Warning: Found no facets matching domain for boundary condition.) I am trying to compute a mesh in 2D, following the code: msh = meshio.read("presa.msh") for cell in msh.cells: if cell.type == "triangle": triangle_cells = cell.data elif cell.type == "line": line_cells = cell.data for key in msh.cell_data_dict["gmsh:physical"].keys(): if key == "line": line_data = msh.cell_data_dict["gmsh:physical"][key] elif key == "triangle": triangle_data = msh.cell_data_dict["gmsh:physical"][key] triangle_mesh = meshio.Mesh(points=msh.points, cells={"triangle": triangle_cells}) line_mesh =meshio.Mesh(points=msh.points, cells=[("line", line_cells)], cell_data={"name_to_read":[line_data]}) meshio.write("mesh.xdmf", triangle_mesh) meshio.write("mf.xdmf", line_mesh) mesh = Mesh() with XDMFFile("mesh.xdmf") as infile: infile.read(mesh) File("Dolfin_circle_mesh.pvd").write(mesh) mvc = MeshValueCollection("size_t", mesh, 1) with XDMFFile("mf.xdmf") as infile: infile.read(mvc, "name_to_read") mf = cpp.mesh.MeshFunctionSizet(mesh, mvc) File("Dolfin_circle_facets.pvd").write(mf). V = FunctionSpace(mesh, 'P', 1) bc1 = DirichletBC(V, Constant(0.0), mf, 1) bc2 = DirichletBC(V, Constant(10.0), mf, 2) bc3 = DirichletBC(V, Constant(6.0), mf, 3) bc4 = DirichletBC(V, Constant(3.0), mf, 4) Additionally, I am using the next .geo file to generate the mesh file: //#gmsh file Point(6) = {0, 0, 0, 1}; Point(7) = {2, 0, 0, 1}; Point(8) = {0, 2, 0, 1}; Point(9) = {2, 2, 0, 1}; Line(10) = {6, 7}; Line(20) = {7, 9}; Line(30) = {9, 8}; Line(40) = {8, 6}; Line Loop(1) = {10, 20, 30, 40}; Plane Surface(1) = {1}; Physical Surface(0) = {1}; Physical Line(1) = {40}; Physical Line(2) = {20}; Physical Line(4) = {30}; Physical Line(3) = {10}; Could you please help me with the problem? I am not sure if it is because of the version of meshio, or I am doing something wrong here. Thank you so much.

Try naming all of your Lines as Physical Line, and not Physical Curve.

Hi Sebastian. Yes, I already did that modification, but I have the same result. I think that my problem is not there. Thank you

The problem is that meshio defines multiple cell block for facets. The following works: import meshio import numpy as np msh = meshio.read("juanfe360.msh") line_cells = [] for cell in msh.cells: if cell.type == "triangle": triangle_cells = cell.data elif cell.type == "line": if len(line_cells) == 0: line_cells = cell.data else: line_cells = np.vstack([line_cells, cell.data]) line_data = [] for key in msh.cell_data_dict["gmsh:physical"].keys(): if key == "line": if len(line_data) == 0: line_data = msh.cell_data_dict["gmsh:physical"][key] else: line_data = np.vstack([line_data, msh.cell_data_dict["gmsh:physical"][key]]) elif key == "triangle": triangle_data = msh.cell_data_dict["gmsh:physical"][key] triangle_mesh = meshio.Mesh(points=msh.points, cells={"triangle": triangle_cells}) line_mesh =meshio.Mesh(points=msh.points, cells=[("line", line_cells)], cell_data={"name_to_read":[line_data]}) meshio.write("mesh.xdmf", triangle_mesh) meshio.xdmf.write("mf.xdmf", line_mesh) I’ve explained the same case for 3D-2D in the post above: [image] Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio mesh Yes, the tetra_data is a simple list of shape (numcells,). To illustrate this, you could use this example: import numpy as np import pygmsh import meshio def create_mesh_gmsh(): geom = pygmsh.built_in.Geometry() rect = geom.add_rectangle(0.0, 2.0, 0.0, 1.0, 0.0, lcar=0.2) geom.add_physical([rect.line_loop.lines[0], rect.line_loop.lines[2]], 1) geom.add_physical([rect.line_loop.lines[1]], 2) geom.add_physical([rect.line_loop.lines[3]], 3) geom.add_physical([rect.surface],…

I am currently trying to covert a simple a .msh to XDMF mesh with the following code from the tutorials above: msh = meshio.read("test_box.msh") for cell in msh.cells: if cell.type == "triangle": triangle_cells = cell.data elif cell.type == "tetra": tetra_cells = cell.data 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] 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("mesh.xdmf", tetra_mesh) meshio.write("mf.xdmf", triangle_mesh) mesh = Mesh() with XDMFFile("mesh.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) cf = cpp.mesh.MeshFunctionSizet(mesh, mvc) but when it tries to read in the file in infield.read(mesh) I get a ValueError: vector. I cannot figure out why, as the same code works on my colleagues computer. I have attached the beginning of my mesh file! Any help would be appreciated. mesh file: $MeshFormat 2.2 0 8 $EndMeshFormat $Nodes 358 1 0 0 1 2 0 0 0 3 0 1 1 4 0 1 0 5 1 0 1 6 1 0 0 7 1 1 1 8 1 1 0 9 0 0 0.1666666666666662 10 0 0 0.3333333333333324 11 0 0 0.4999999999999986 12 0 0 0.6666666666666646 13 0 0 0.8333333333333321 14 0 0.1666666666666662 1 15 0 0.3333333333333324 1 16 0 0.4999999999999986 1 17 0 0.6666666666666646 1 18 0 0.8333333333333321 1 19 0 1 0.1666666666666662 20 0 1 0.3333333333333324 21 0 1 0.4999999999999986 22 0 1 0.6666666666666646 23 0 1 0.8333333333333321 24 0 0.1666666666666662 0 25 0 0.3333333333333324 0 26 0 0.4999999999999986 0 27 0 0.6666666666666646 0 28 0 0.8333333333333321 0 29 1 0 0.1666666666666662 30 1 0 0.3333333333333324 31 1 0 0.4999999999999986 32 1 0 0.6666666666666646 33 1 0 0.8333333333333321 34 1 0.1666666666666662 1 35 1 0.3333333333333324 1 36 1 0.4999999999999986 1 37 1 0.6666666666666646 1 38 1 0.8333333333333321 1 39 1 1 0.1666666666666662 40 1 1 0.3333333333333324 41 1 1 0.4999999999999986 42 1 1 0.6666666666666646 43 1 1 0.8333333333333321 44 1 0.1666666666666662 0 45 1 0.3333333333333324 0 46 1 0.4999999999999986 0 47 1 0.6666666666666646 0

If it works on your colleagues computers, it is probably a version issue. So which versions of Meshio, FEniCS and gmsh are you using? And what versions are they using? Additionally, you need to supply the whole mesh file (make the coarsest possible mesh) for anyone to be able to debug it.

Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio

mesh

plugged March 10, 2020, 3:43pm 109

What are the boundary conditions you are using? Are they all Dirichlet type?

Code from Poisson subdomain problem tutorial not working

.msh to .xdmf errors, paraview crash

Topic		Replies	Views
Extracting mesh regions for gmsh mesh mesh	24	5563	December 5, 2023
Transitioning from Triangle -> GMSH -> XDMF I/O	8	1080	April 16, 2023
Need help converting GMSH to FEniCS installation	61	7333	September 28, 2022
Create mesh with gmsh, convert it and import it in dolfin, problems mesh	5	11064	April 1, 2019
How to import mixed triangle and quad meshes from gmsh to fenics mesh	22	5158	June 20, 2020

Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio

Related topics