Create mesh with gmsh, convert it and import it in dolfin, problems

I’ve created a mesh from a file.geo using gmsh. I used dolfin-convert to convert it to xml, with hopes that I could use it with dolfin. However dolfin-convert does not produce the file_facet_region.xml, essential for the boundary conditions. If only I could make it produce that file, my problem would be solved.

So I tried meshio-convert, which produced two files: file_gmsh:geometrical.xml and file_gmsh:physical.xml, but it seems that they are not equivalent to the dolfin-convert files, and using them in dolfin as I use the dolfin-converted files return errors like

*** Error: Unable to read mesh function from XML file.
*** Reason: Type mismatch reading XML MeshFunction. MeshFunction type is “uint”, but file type is “int”.

I tried to fix the error by modifying the .xml file and changing “uint” to “int”, to no avail, i.e. I get the exact same error message despite the file containing only “int”, and no “uint”.

I then tried to use meshio-convert to produce xdmf files, but it produced an xdmf and an h5 file, which I do not know how to read/use with dolfin properly.

I’m used to do:

subdomains = MeshFunction("size_t", mesh, physical_name_file)
boundary_parts = MeshFunction("size_t", mesh, facet_name_file)

So here I am, lost on how to use the mesh I produced with gmsh in dolfin.

Here’s my file.geo:

Mesh.RandomFactor=1.0e-6;
lc=0.008;
Mesh.CharacteristicLengthFromPoints = 0;
Mesh.CharacteristicLengthExtendFromBoundary = 1;
Mesh.CharacteristicLengthFromCurvature = 0;
Mesh.CharacteristicLengthMin = lc/20;
Mesh.CharacteristicLengthMax = lc;

// bottom motif
a = 2;
b = 1;
c = 0.03;

Point(2) = {-a/2, -b/2, -c/2};
Point(3) = {-a/2, b/2, -c/2};
Point(4) = {a/2, b/2, -c/2};
Point(5) = {a/2, -b/2, -c/2};

Line(100) = {2,3};
Line(101) = {3,4};
Line(102) = {4,5};
Line(103) = {5,2};

Line Loop(200) = {100,101,102,103};
Plane Surface(300) = {200};

// top motif
Point(6) = {-a/2, -b/2, c/2};
Point(7) = {-a/2, b/2, c/2};
Point(8) = {a/2, b/2, c/2};
Point(9) = {a/2, -b/2, c/2};

Line(104) = {6,7};
Line(105) = {7,8};
Line(106) = {8,9};
Line(107) = {9,6};

Line Loop(201) = {104,105,106,107};
Plane Surface(301) = {201};

// vertical lines between the 2 surfaces
Line(108) = {2, 6};
Line(109) = {7, 3};
Line(110) = {4, 8};
Line(111) = {5, 9};

// left
Line Loop(202) = {108,104,109,-100};
Plane Surface(302) = {202};

// front
Line Loop(203) = {103,108,-107,-111};
Plane Surface(303) = {203};

// right
Line Loop(204) = {106,-111,-102,110};
Plane Surface(304) = {204};

// back
Line Loop(205) = {-105,109,101,110};
Plane Surface(305) = {205};

// auxiliary points over top and bottom surface
//rectangle 1
Point(112) = {-a/8,-b/8, c/2};
Point(113) = {-a/8,b/8,c/2};
Point(114) = {a/8, b/8,c/2};
Point(115) = {a/8,-b/8,c/2};

Line(112) = {112,113};
Line(113) = {113,114};
Line(114) = {114,115};
Line(115) = {115,112};

Line Loop(206) = {112,113,114,115};
Plane Surface(306) = {206};

// rectangle 2
Point(116) = {-a/16+a/4,-b/16, -c/2};
Point(117) = {-a/16+a/4,b/16,-c/2};
Point(118) = {a/16+a/4, b/16,-c/2};
Point(119) = {a/16+a/4,-b/16,-c/2};

Line(116) = {116,117};
Line(117) = {117,118};
Line(118) = {118,119};
Line(119) = {119,116};

Line Loop(207) = {116,117,118,119};
Plane Surface(307) = {207};

// rectangle 3
Point(120) = {-a/16+a/4,-b/16, c/2};
Point(121) = {-a/16+a/4,b/16,c/2};
Point(122) = {a/16+a/4, b/16,c/2};
Point(123) = {a/16+a/4,-b/16,c/2};

Line(120) = {120,121};
Line(121) = {121,122};
Line(122) = {122,123};
Line(123) = {123,120};

Line Loop(208) = {120,121,122,123};
Plane Surface(308) = {208};

//rectangle 4
Point(124) = {-a/22,-b/22, -c/2};
Point(125) = {-a/22,b/22,-c/2};
Point(126) = {a/22, b/22,-c/2};
Point(127) = {a/22,-b/22,-c/2};

Line(124) = {124,125};
Line(125) = {125,126};
Line(126) = {126,127};
Line(127) = {127,124};

Line Loop(209) = {124,125,126,127};
Plane Surface(309) = {209};

// Physical Surfaces to build the mesh.
Physical Surface(“bot_box”) = {300};
Physical Surface(“top_box”) = {301};
Physical Surface(“left_box”) = {302};
Physical Surface(“front_box”) = {303};
Physical Surface(“right_box”) = {304};
Physical Surface(“back_box”) = {305};

// rectangular surfaces
Physical Surface(“rect_1”) = {306};
Physical Surface(“rect_2”) = {307};
Physical Surface(“rect_3”) = {308};
Physical Surface(“rect_4”) = {309};

Surface Loop(400) ={300,301,302,304,303,305};
Physical Volume(1) = {400};

In the end I wish to apply boundary conditions on the Physical Surfaces(“rect_1”) through “rect_4”.
I would appreciate some guidance. Thank you!

1 Like

I would suggest not using xml, as this is a legacy-format.
How about using meshio to save your file as XDMF:

import meshio

geometry = meshio.read(filename)
meshio.write("mesh.xdmf", meshio.Mesh(points=geometry.points, cells={"triangle": geometry.cells["triangle"]}))
meshio.write("mf.xdmf", meshio.Mesh(points=geometry.points, cells={"line": geometry.cells["line"]},
                                                               cell_data={"line": {"name_to_read": geometry.cell_data["line"]["gmsh:physical"]}}))

This can be read into dolfin in the following way:

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(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))

1 Like

Thanks for the suggestion. Unfortunately the first code snippet returns

—> 11 meshio.write(“mf.xdmf”, meshio.Mesh(points=geometry.points, cells={“line”: geometry.cells[“line”]},
12 cell_data={“line”: {“name_to_read”: geometry.cell_data[“line”][“gmsh:physical”]}}))
13
KeyError: ‘line’

Nevertheless, the file mesh.xdmf is created, so I could proceed with the 1st part of your 2nd snippet. Unfortunately I also got an error:

RuntimeError Traceback (most recent call last)
in
1 mesh = Mesh()
2 with XDMFFile(“mesh.xdmf”) as infile:
----> 3 infile.read(mesh)
4 mvc = MeshValueCollection(“size_t”, mesh, 1)

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 get dataset.
*** Reason: DOLFIN has not been configured with HDF5.
*** Where: This error was encountered inside XDMFFile.cpp.
*** Process: 0


*** DOLFIN version: 2018.2.0.dev0
*** Git changeset: d4731dd5de3c6dbdf6c18a6e2764beb5b624ff5b
*** -------------------------------------------------------------------------

But as far as I understand, I have just recompiled dolfin-git and python-dolfin-git (I use Arch Linux, both packages are in the AUR) which should have compiled with the HDF5 configuration. I am unsure how to fix any of the 2 errors…

Hi again,

First of all, i did not realize you were using a 3D mesh, see
Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio for a detailed description on how to load cell and facet functions.

Another point is that I am not able to generate a 3D mesh with gmsh using your code. Could you please verify that your mesh generates a 3D mesh, (using i.e. gmsh -3 meshfile.geo and open the msh file in gmsh.

1 Like

Thank you. I am the author of the thread you link to, so I fixed the transition from dolfin-convert to meshio thanks to you, even though there is a very slight discrepancy in the results I get when solving a problem (with a given weak form) between the meshes from dolfin-convert and from meshio.

However in case you wonder, I am able to produce a correct file.msh using gmsh file.geo -3. I can examine it with gmsh without any problem nor error. What problem do you face?

If i run gmsh -3 on your geo file, it only meshes the 2D surfaces:

Info    : Done meshing 2D (3.864 s)
Info    : Meshing 3D...
Info    : Done meshing 3D (0 s)
Info    : 90432 vertices 182392 elements
Info    : Writing 'test_geo_3d.msh'...
Info    : Done writing 'test_geo_3d.msh'
Info    : Stopped on Mon Apr  1 09:52:14 201

I think its because the surface loop is not correct, but havent had time to investigate further.