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

#1

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.

1 Like
Conversion from .msh to .xml
#2

See this gist, before converting to XDMF you’d need to “prune” all lower-dimensional elements, as currently FEniCS does not support “mixed” mesh composing from several element types (may they belong to a same dimension or not!) According to their recent released roadmap, the support for mixed elements is for long term…

#3

This can be done using meshio and writing the explicit data to xdmf file.

import meshio
msh = meshio.read("test.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"]}}))

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)
File("dolfinmesh.pvd").write(mesh)
File("dolfincellfunc.pvd").write(mf)
3 Likes
Converting simple 2D mesh from gmsh to dolfin
Should I be using XML in 2019?
Run with h5 file
Unable to convert from .msh to .xml
Create mesh with gmsh, convert it and import it in dolfin, problems
#4

Thank you very much! If that’s not too much asked, how can I proceed and do the equivalent of the lines

subdomains = MeshFunction(“size_t”, mesh, file_physical_region.xml)
boundary_parts = MeshFunction(“size_t”, mesh, file_facet_region.xml)

?
I cannot simply replace file_physical_region.xml by dolfinmesh.pvd and file_facet_region.xml by dolfincellfunc.pvd.

#5

Hmm but isn’t it strange that I do not need any manual pruning when I use dolfin-convert to xml? As far as I understand, the elements of the mesh shouldn’t depend on the file extension… I guess I’m missing something.

#6

Well dolfin_convert is specially designed for FEniCS so I suppose these lower dimensional elements are taken care of during conversion. meshio is a generic mesh convertor so we try to conserve all original mesh information.

Le ven. 22 mars 2019 à 16:58, raboynics via FEniCS Project fenicsproject1@discoursemail.com a écrit :

#7

As you can see in the code snippet above, I first read the boundary data to the file mf.xdmf. this is equivalent to facet_region.xml.
For the cellfunction:

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"]}}))
mvc = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("cf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

1 Like
#8

If I understand well, what I call subdomains corresponds to your mvc, while what I call boundary_parts corresponds to your cf.
But when I use your code and then subdomains = mvc and boundary_parts = cf, I get nonsensical results.
For example a non-zero surface area assemble(1 * ds(12, domain=mesh)) returns 0 now, instead of the correct area. I suspect it’s because in my original code I have mesh = Mesh("file.xml") while in my new code I have mesh = Mesh(). I do not know which argument I can put inside the function Mesh(). I tried the XDMF file, but this returns an error. So I still don’t know how to make the full transition to XDMF.

Furthermore I am only able to reach that far on an Ubuntu machine, where dolfin has been compiled with hdf5 support. On Arch Linux (at home), I get an error message when I use your snippet codes as described in an earlier question here, and on the fenics support mailing list. I do not know how to compile dolfin with hdf5 support, it doesn’t do it by default even when the hdf5 package(s) are installed.

Solve PDE on 2D non-flat surface
#9

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
1 Like
Defining subdomains of 3D mesh
#10

Thank you very much for the information and code snippet. I can finally get it to work (on Ubuntu). However there is a slight discrepancy with the result of solving a weak form problem I get with the file.xml from dolfin-convert, and I do not know which one is more accurate (xdmf or xml). Can I simply compare the number of mesh elements in the two cases and the ones with a higher number should in principle be the most accurate?

More precisely, for a calculation I get 0.01031139055668843 with the xdmf mesh and 0.010311376092968772 with the mesh from the xml file. I have put in bold where they match.

Edit: Hmm, that’s extremely strange!
print(mesh.num_edges(), mesh.num_cells(), mesh.num_facets()) yields exactly the same numbers for both the xml and xdmf meshes. I am not sure what’s causing the solution to a complicated weak form to differ, since I copied and pasted the whole code and modified only the mesh part.

#11

Hello!

I’m having a similar problem. I have mesh.msh file generated with

gmsh mesh.geo -2 -format msh2 -o mesh.msh

And when I read it a try to create a new meshio Mesh as

mesh = meshio.read(‘mesh.msh’)
mesh.points = mesh.points[:, :2] # remove z-values to force 2d
mesh = meshio.Mesh(points=mesh.points, cells={“triangle”: mesh.cells[“triangle”]})

I got KeyValue error ‘triangle’ since mesh.cells are only ‘lines’ not ‘triangle’.

Could anybody please help me with that?

Thanks!

#12

Could you please post your GEO file?

#13

// the square
Point(1) = {0, 0, 0, 0.1E+3};
Point(2) = {0, 11E+3, 0, 0.1E+3};
Point(3) = {11E+3, 11E+3, 0, 0.1E+3};
Point(4) = {11E+3, 0, 0, 0.1E+3};
Line(1) = {1, 4};
Line(2) = {4, 3};
Line(3) = {3, 2};
Line(4) = {2, 1};
// the circle in the centre
Point(5) = {5.5E+3, 5.5E+3, 0, 0.1E+3};
Point(6) = {1E+3, 5.5E+3, 0, 0.1E+3};
Circle(5) = {6, 5, 6};
// creating the surfaces
Line Loop(6) = {2, 3, 4, 1};
Line Loop(7) = {5};
Plane Surface(8) = {6, 7};
Plane Surface(9) = {7};
//extruding them
Extrude {0, 0, 0.55E+3} {
Surface{9, 8}; Layers{10};
}//+
Physical Volume(0) = {1, 2};
//+
Physical Surface(1) = {16};
//+
Physical Surface(2) = {9};

#14

Couple of notes, are you doing a 2D mesh ? When i Look at you mesh in gmsh, it strikes me as you would like to have a 3D mesh?

#15

I am doing 3D mesh, how can i use to write and read it for fenics?
I got errors

#16

I am aplying Dirichilet boundary on Physical surface 1 and Newman in surface 2.
I am solving a poisson system like in transient.

#17

Could you supply me with your error messages, as I your mesh file, with gmsh -3 mesh3d.geo, and

import meshio

msh = meshio.read("mesh3d.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)

works for me