Problem converting file.msh (from gmsh) to file.xml

I’ve created a file.geo with gmsh, and then generated a 3D mesh, without error nor warning. That mesh looks fine when viewed in gmsh, however dolfin-convert fails to convert it to xml.
Here’s the file.geo:

SetFactory("OpenCASCADE");

// The box
Box(1) = {0, 0, 0, 2, 1, 0.038};

// The rectangular surface over the surface of the box.
Rectangle(7) = {0, 0, 0.038, 0.1, 0.2, 0};

// Define the Physical Surfaces where the boundary conditions will be applied.
Physical Surface("rectangle") = {7};
Physical Surface("top_box") = {6};
Physical Surface("bottom_box") = {5};
Physical Surface("left_box") = {1};
Physical Surface("front_box") = {3};
Physical Surface("behind_box") = {4};
Physical Surface("right_box") = {2};
Physical Volume("box") = {1};

And here’s what “dolfin-convert -i gmsh file.msh file.xml” returns:

Converting from Gmsh format (.msh, .gmsh) to DOLFIN XML format
Expecting 182 vertices
Found all vertices
Expecting 456 cells
Found all cells
Traceback (most recent call last):
File “/usr/bin/dolfin-convert”, line 132, in
main(sys.argv[1:])
File “/usr/bin/dolfin-convert”, line 79, in main
meshconvert.convert2xml(ifilename, ofilename, iformat=iformat)
File “/usr/lib/python3/dist-packages/dolfin_utils/meshconvert/meshconvert.py”, line 1301, in convert2xml
convert(ifilename, XmlHandler(ofilename), iformat=iformat)
File “/usr/lib/python3/dist-packages/dolfin_utils/meshconvert/meshconvert.py”, line 1322, in convert
gmsh2xml(ifilename, handler)
File “/usr/lib/python3/dist-packages/dolfin_utils/meshconvert/meshconvert.py”, line 494, in gmsh2xml
index = nodes_as_facets[tuple(nodes)]
KeyError: (0, 1, 62)

Which is essentially exactly the same problem described here back in 2016.

Does someone have an idea on what’s wrong? Thank you!

I suppose your intention was to have the mesh conform to the boundary of the rectangle 7. However, with your geo code the mesh won’t have this property as can be seen by visual inspection in gmsh and this is also the reason why dolfin-convert raises an error. For conforming mesh consider e.g.

SetFactory("OpenCASCADE");

Box(1) = {0, 0, 0, 2, 1, 0.038};
Box(2) = {0, 0, 0, 0.1, 0.2, 0.038};

BooleanFragments{ Volume{1}; Delete; }{ Volume{2}; Delete; }

Physical Volume("box") = {2, 3};

The resulting msh file is converted without problems (gmsh 3.0.6, fenics 2017.2.0).

I do not really understand what you mean by “having the mesh conform to the boundary of the rectangle 7”. I would like to have a parallelepipede (a box), and I want to apply a boundary condition on a small rectangular (rectangle 7) surface (not a volume!) over that box.

I have fixed the error in my post by reordering the “Physical Surface” lines, i.e. from {1} to {7}.
However this now produces the error:

dolfin-convert fenics.msh fenics.xml
Converting from Gmsh format (.msh, .gmsh) to DOLFIN XML format
Traceback (most recent call last):
File “/usr/bin/dolfin-convert”, line 132, in
main(sys.argv[1:])
File “/usr/bin/dolfin-convert”, line 79, in main
meshconvert.convert2xml(ifilename, ofilename, iformat=iformat)
File “/usr/lib/python3.7/site-packages/dolfin_utils/meshconvert/meshconvert.py”, line 1301, in convert2xml
convert(ifilename, XmlHandler(ofilename), iformat=iformat)
File “/usr/lib/python3.7/site-packages/dolfin_utils/meshconvert/meshconvert.py”, line 1322, in convert
gmsh2xml(ifilename, handler)
File “/usr/lib/python3.7/site-packages/dolfin_utils/meshconvert/meshconvert.py”, line 271, in gmsh2xml
num_elements = int(line)
ValueError: invalid literal for int() with base 10: ‘8 695\n’

These are facets on the surface of the mesh produced by your code. Your issue
is due to yellow triangles not respecting the rectangle boundary; i.e. what I meant
by (the lack of ) conformity.

The code I posted is the shortest one to get conforming mesh for the rectangle. You can improve it to avoid the extra volume but the conformity idea is necessary, c.f.
the new error that you are seeing.

I see, I didn’t know this was the problem, thank you!
I am really not sure of a correct way to get rid of the extra volume though. In the end I wish to have a single rectangular surface on one side only.

One way to get rid of it would be to define the geometry in terms lines, plane surfaces, etc, that is, without the CAD primitives.

I see…
Here’s my current failed attempt to create a simple parallelepipede mesh:

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

// bottom motif

a = 2;
b = 1;
c = 0.038;

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};

// 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(40) = {2,3};
Line(41) = {3,4};
Line(42) = {4,5};
Line(43) = {5,2};

Line(44) = {6,7};
Line(45) = {7,8};
Line(46) = {8,9};
Line(47) = {9,6};

//+
Line(48) = {2, 6};
//+
Line(49) = {7, 3};
//+
Line(50) = {4, 8};
//+
Line(51) = {5, 9};


Line Loop(60) = {40,41,42,43};
Line Loop(61) = {44,45,46,47};

// left
Line Loop(62) = {48,44,49,-40};

// front
Line Loop(63) = {51,47,-48,-43};

// back
Line Loop(64) = {50,-45,49,41};

// right
Line Loop(65) = {50,46,-51,-42};

Plane Surface(70) = {60};
Plane Surface(71) = {61};
Plane Surface(72) = {62};
Plane Surface(73) = {63};
Plane Surface(74) = {64};
Plane Surface(75) = {65};

Physical Surface("bot_box") = {70};
Physical Surface("top_box") = {71};
Physical Surface("left_box") = {72};
Physical Surface("front_box") = {73};
Physical Surface("back_box") = {74};
Physical Surface("right_box") = {75};//+
Extrude {0, 0, c} {
  Surface{70}; Layers{1}; Recombine;
}
//+
Physical Volume("the_box") = {1};

Do you have an idea why it doesn’t work? I’m getting errors like Error : Could not find extruded vertex (-0.3540811831507598, -0.1075709531782967, 0.019) and Error : Vertex 0x6060a30 not MEdgeVertex when creating a 3d mesh. Note that I haven’t implemented the rectangular shape yet…

What gmsh version are you using? Anyways, consider adapting the following code

size = 1;

// Box spec
x0 = 0;
y0 = 0;
z0 = 0;
dx = 1;
dy = 2;
dz = 3;

// Rectangle
dxx = 0.2;
dyy = 0.3;
//
//   8----------7
//11/--10      /| 
// 5---/------6 |
// |   9      | |
// | 4--------|-3           /y
// |/         |/           /
// 1----------2            ------
//                              x 


// Base
Point(1) = {x0, y0, z0, size};
Point(2) = {x0+dx, y0, z0, size};
Point(3) = {x0+dx, y0+dy, z0, size};
Point(4) = {x0, y0+dy, z0, size};

Line(12) = {1, 2};
Line(23) = {2, 3};
Line(34) = {3, 4};
Line(41) = {4, 1};

// Shifted Base
Point(5) = {x0, y0, z0+dz, size};
Point(6) = {x0+dx, y0, z0+dz, size};
Point(7) = {x0+dx, y0+dy, z0+dz, size};
Point(8) = {x0, y0+dy, z0+dz, size};

// Aux rectangle points
Point(9) = {x0+dxx, y0, z0+dz, size};
Point(10) = {x0+dxx, y0+dyy, z0+dz, size};
Point(11) = {x0, y0+dyy, z0+dz, size};

// Verticals
Line(15) = {1, 5};
Line(26) = {2, 6};
Line(37) = {3, 7};
Line(48) = {4, 8};

// Top lines
Line(59) = {5, 9};
Line(96) = {9, 6};
Line(67) = {6, 7};
Line(78) = {7, 8};
Line(811) = {8, 11};
Line(115) = {11, 5};
Line(910) = {9, 10};
Line(1011) = {10, 11};

Line Loop(1) = {12, 23, 34, 41};  // Might need curve loop here instead

Line Loop(2) = {26, 67, -37, -23};
Plane Surface(1) = {2};

Line Loop(3) = {34, 48, -78, -37};
Plane Surface(2) = {3};

Plane Surface(3) = {1};

Line Loop(4) = {48, 811, 115, -15, -41};
Plane Surface(4) = {4};

Line Loop(5) = {15, 59, 96, -26, -12};
Plane Surface(5) = {5};

Line Loop(6) = {910, 1011, 115, 59};
Plane Surface(6) = {6};

Line Loop(7) = {78, 811, -1011, -910, 96, 67};
Plane Surface(7) = {7};

Surface Loop(1) = {7, 2, 3, 5, 4, 6, 1};
Volume(1) = {1};
// TODO Physical Surface(x) = {y};
//             Physical Volume(xx) = {yy};

Oh, I didn’t know I could define volumes with Surface Loop().
Back to my parallelepipede example (without the rectangle), I can now generate a mesh that looks totally fine with gmsh, but dolfin still fails to convert it properly.
My code is:

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

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

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};

// 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(40) = {2,3};
Line(41) = {3,4};
Line(42) = {4,5};
Line(43) = {5,2};

Line(44) = {6,7};
Line(45) = {7,8};
Line(46) = {8,9};
Line(47) = {9,6};

//+
Line(48) = {2, 6};
//+
Line(49) = {7, 3};
//+
Line(50) = {4, 8};
//+
Line(51) = {5, 9};
Line Loop(60) = {40,41,42,43};
Line Loop(61) = {44,45,46,47};

// left
Line Loop(62) = {48,44,49,-40};

// front
Line Loop(63) = {51,47,-48,-43};

// back
Line Loop(64) = {50,-45,49,41};

// right
Line Loop(65) = {50,46,-51,-42};

Plane Surface(70) = {60};
Plane Surface(71) = {61};
Plane Surface(72) = {62};
Plane Surface(73) = {63};
Plane Surface(74) = {64};
Plane Surface(75) = {65};

Surface Loop(80) = {70, 71,72,73,74,75};
Physical Volume("the_box") = {80};

This time there is no apparent problem with the mesh, yet the error I get is still:

dolfin-convert parallelepipede_fenics.msh parallelepipede_fenics.xml
Converting from Gmsh format (.msh, .gmsh) to DOLFIN XML format
Traceback (most recent call last):
File “/usr/bin/dolfin-convert”, line 132, in
main(sys.argv[1:])
File “/usr/bin/dolfin-convert”, line 79, in main
meshconvert.convert2xml(ifilename, ofilename, iformat=iformat)
File “/usr/lib/python3.7/site-packages/dolfin_utils/meshconvert/meshconvert.py”, line 1301, in convert2xml
convert(ifilename, XmlHandler(ofilename), iformat=iformat)
File “/usr/lib/python3.7/site-packages/dolfin_utils/meshconvert/meshconvert.py”, line 1322, in convert
gmsh2xml(ifilename, handler)
File “/usr/lib/python3.7/site-packages/dolfin_utils/meshconvert/meshconvert.py”, line 271, in gmsh2xml
num_elements = int(line)
ValueError: invalid literal for int() with base 10: ‘26 1310\n’

By the way, I use gmsh 4.1.0.

Now, do you have an idea why my code fails for a simple parallelepipede?

Can you test with an older gmsh version? As dolfin-convert is no longer maintained, (meshio is recommended replacement) there might be missing support for the new numbering in Gmsh 4. There might be python3 problems as well.

1 Like

Hi,

are you exporting your mesh using msh version 2 (ASCII)? I think that dolfin-convert only works with that version of msh format.

Cheers,

Rodolfo

2 Likes

Thank you very much to both of you! I didn’t know about meshio and the obsolete format. I just tried with the msh version 2, and dolfin-convert has now no problem to convert the .msh to .xml for the parallelepipede. I still have to try to add the rectangular surface.

I also installed meshio and gave it a try, it worked fine although returned a warning about the xml file format being deprecated.

Unfortunately I haven’t been able to adapt your code when the rectangular surface is not at a corner.
I like the idea of using BooleanFragments, but for some reason I do not get the same nice result than you get.

I do get a nice mesh that is adapted to the rectangular surface, but I also get the default one that doesn’t conform with the rectangular surface. I could not, for the life of me, get rid of that extra annoying meshing. Below is my attached code (note that I’d like to use extrusion):

SetFactory(“OpenCASCADE”);
Mesh.RandomFactor=1.0e-6;
lc=0.8;

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

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};
Plane Surface(6) = {5};

// rectangular surface
Point(7) = {-a/16+a/4,-b/16, -c/2};
Point(8) = {-a/16+a/4,b/16,-c/2};
Point(9) = {a/16+a/4, b/16,-c/2};
Point(10) = {a/16+a/4,-b/16,-c/2};

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};

surfaceVector = Extrude {0, 0, c} {
Surface{6};
Layers{1};
Recombine;
};
/* surfaceVector contains in the following order:
[0] - front surface (opposed to source surface)
[1] - extruded volume
[2] - bottom surface (belonging to 1st line in “Line Loop (6)”)
[3] - right surface (belonging to 2nd line in “Line Loop (6)”)
[4] - top surface (belonging to 3rd line in “Line Loop (6)”)
[5] - left surface (belonging to 4th line in “Line Loop (6)”) */

Physical Surface(17)=BooleanFragments{ Surface{16}; Delete; }{ Surface{6}; Delete; };

//Physical Surface("front") = surfaceVector[0];
Physical Volume("internal") = surfaceVector[1];
//Physical Surface("bottom") = surfaceVector[2];
Physical Surface("right") = surfaceVector[3];
//Physical Surface("top") = surfaceVector[4];
Physical Surface("left") = surfaceVector[5];
Physical Surface("back") = {6}; // from Plane Surface (6) ...

I’ve tried several variants, by not defining any Physical Surface, yet the meshing would be the same. By removing surfaces, but gmsh would complain it cannot modify the surfaces, etc.

To me it seems like the BooleanFragments line produces a correct mesh, but for some reason it doesn’t remove or delete what it should, despite that I specify “delete”.

I’m almost there, it would be marvelous if this finally worked…