Gmsh mesh in freefem

I’ve been trying to export a gmsh mesh to freefem including the physical groups (for labeling purpose). Does anyone have any pointers on how to do that?
I have 3 physical groups: inlet, outlet, wall (which includes all the walls in the geometry). I’m exporting the file as .mesh file, and I can import it into freefem using mesh Th(nameoffile). However, I can’t seem to find the labels.
I’ve read in another thread about using labels(Th), but that apparently only finds 2 labels.

An update: there’s a difference between using gmsh 2.7.1 and 4.7.1. It looks like it’s almost working in 4.7.1 (i.e. uploading the mesh into freefem and asking for its labels returns 3 labels).
On the other hand, what should be a hole in the mesh is now meshed, too. The outside mesh is red, while the inside mesh is orange, and I couldn’t figure out what it means.

If you want to reproduce this, here is the .geo file.

In gmsh 4.7.1: save as .msh, version 2 ASCII. If you save all elements you get a correct mesh, but incorrect labels. If you don’t save all elements you get an incorrect mesh but (presumably) correct labels.
In gmsh 2.7.1: save as .msh, version 2 ASCII. No matter whether you select to ignore physical groups or not, the mesh is correct but the labels are not.

I meshed with an element size factor of 50 in both cases.
Import it into freefem using:

load “gmsh”
mesh Th = gmshload(“tesla.msh”);
plot(Th);
int[int] labs = labels(Th);
cout << “LABELS \n” << labs << endl;

As said earlier, the results of plot and labels vary based on which version of gmsh you use and which option you select. As a cheatsheet:
2.7.1 => correct mesh, incorrect labels
4.7.1, all elements => correct mesh, incorrect labels
4.7.1, unselect all elements => incorrect mesh, presumably correct labels.

Does anyone know how to get both a correct mesh and the correct labels?

I’d suggest you use the .mesh (Inria MEDIT) format. This is partially covered in this tutorial.

1 Like

I have a couple of problems with MEDIT files. First of all, the mesh is 3D, and I don’t really know how to handle it (in the definition of the problem). Secondly, if I don’t export all elements with a MEDIT file, I get an error. Specifically, “WARNING!!! The mesh file just contains a set of vertices”.
If I do export all elements, though, I get the incorrect labels once again.

Even if you follow the different steps from the tutorial?

I just tried adding the readmesh3 (which seemed to be the only difference from my code to his) and I still get the same problem.

load “msh3”
mesh3 Th(“tesla.mesh”);
Th = readmesh3(“tesla.mesh”);

I think I can reproduce your issue. This maybe be due to the “Reverse Surface” in your .geo. Not sure though, have you tried a simplified geometry?

I initially made the geometry using gmsh GUI and added the reverse surface manually to the .geo because I was getting an “assertion error: Area > 0” in freefem.
Even then, if I do remove the reverse surface and make the .msh in 4.7.1, I either get the assertion error (if I export all elements) or keep getting the same wrong geometry from earlier.
I remember trying out a holed square to check out if gmsh was having some problems with the not simply connected surface, but the meshing worked alright. I don’t think I had labels on, however.
I tried it out again just to be sure and, yeah, I get the same problem.

On the bright side, by trying it on a simple square, I noticed that the mesh as read by freefem is different from the one in gmsh even for simply connected surfaces, so it’s most likely a problem with gmshload.

It seems to be necessary to create an element set as well.

Physical Surface(“channel”) = {31};

Then when I export the mesh in gmsh-format version 2 ascii without the include all elements option and import it in Freefem++ then I get the three labels 2, 3, 4. Label 1 is assigned to channel in the msh-file.

By the way in gmsh the elements inherit their numbering scheme (clockwise or counterclockwise) from the underlying geometry. In Freefem++ node numbers have to be given in counterclockwise order (as with most FEM-programs I know) otherwise you get the “assertion error: Area > 0”.
This means for Freefem++ the elements have a negative area. Such elements are also called inside out in some commercial finite element programs.

2 Likes

The physical surface seems to solve the problem. Thanks to both of you!

I figured the assertion error was due to the number scheme, but it felt off to me that gmsh didn’t automatically reverse the surface to make it compliant with what most programs need.

Hello !
I’m also trying to import 2d meshes from gmsh in FreeFem and found this thread.
I tried to follow the solution proposed here but did not succeed.

I first started from the “tesla.geo” and added the line
“Physical Surface(“channel”) = {31};” as suggested.
But I am not able to export under format “gmsh-format version 2”.
This seems to be a legacy format which is not handled any more by gmsh… With the version I have (gmsh 4.8.4), this format is not proposed in the “export” dialog box.

With default output format I get the error
“Element of type 17 is not considered in Freefem++”

I also tried using Inria MEDIT format but this does not work as well:
mshptg bug in number of points 0 > 0 == max nb points
Sorry Error build delaunay triangle error = 1
current line = 2
Exec error : Error mshptg8_

Also tried with a much simpler mesh as found here but did not succeed.

Can someone help ?
It would be very useful to provide a simple example of how to import a 2D mesh with gmshload.
I have seen in the examples/ folder a case with gmshload3 but none with gmshload

Thanks for help !!
David

I think you can use the ‘guess from extension’ option, and save the file like ‘tesla.msh2’ and choose then the version 2 ASCII option, or you can also run gmsh from command line.

Thanks for the tip Aszaboa ! indeed it works this way.

I spent some time doing a number of elementary tests, starting from the examples in the gmsh tutorial directory. My diagnostics is as follows:

  • Works for simple geometries (such as example t1.geo), using gmsh-format version 2 ascii (version 4 and INRIA format don’t work).
  • Assigning labels effectively works using the solution given in this thread.
  • On the other hand, I was not able to succeed with more complicated geometries such as the “t4.geo” from the gmsh tutorial. Importing the mesh seems to work but plotting crashes.

Overall, I see a solution for what I expected to do, but the facts that (i) the solution relies on a legacy format, (ii) gmshload is not documented in the freefem doc and (iii) there is no example using this in the freefem repository make me fear that the solution may not be perennal…

I’m thinking that probably the most satisfying solution would be to suggest to gmsh developers to allow directly exportation under freefem .msh format.

What do you think ?

If you want to use version 4 and/or binary format, you could use

load "PETSc"
DM dm("your_mesh.msh");
mesh Th = dm;

It relies on PETSc DMPlex reader, which is much more up-to-date than FreeFEM gmshload (and support other formats, e.g., Salome, Exodus, HDF5…).

1 Like

Thanks Pierre ! I’ll try this !

Dear forum,
I am trying to import a 2D mesh from GMSH 4.9.3 (Windows version) to FF++ v. 4.10 (Windows version). I attach the .geo code:

//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {0, 0, 0.5, 1.0};
//+
Point(3) = {0.5, 0, 0.5, 1.0};
//+
Point(4) = {0.5, 0, 0, 1.0};
//+
Line(1) = {1, 4};
//+
Line(2) = {4, 3};
//+
Line(3) = {3, 2};
//+
Line(4) = {2, 1};
//+
Curve Loop(1) = {1, 2, 3, 4};
//+
Plane Surface(1) = {1};
//+
Physical Curve("Inlet", 1) = {1};
//+
Physical Curve("Outlet", 2) = {3};
//+
Physical Curve("Upper", 3) = {2};
//+
Physical Curve("Lower", 4) = {4};
//+
Physical Surface("Fluid", 5) = {1};

To export the .msh file I use: File->Export-> Test.msh file. For the MSH Options I selected Version 2 ASCII and I unchecked “Save all elements” and “Save parametric coordinates”. When I try to plot the mesh in FF++ with

load "gmsh"

mesh Th=gmshload("Test.msh");
plot(Th);

I get the “Assertion fail : (area > 0)” error message.

What else I have tried:

  1. Checked “Save parametric coordinates” when saving the mesh but I get the error:
    “ParametricNodes is not considered yet in freefem++”
  2. Changed in the .geo file the line
    “Curve Loop(1) = {1, 2, 3, 4};” to “Curve Loop(1) = {-1, -2, -3, -4};”, reloaded the script and exported the mesh again (with both “Save all elements” and “Save parametric coordinates” checked, then both unchecked, with only 1 of them checked etc. but still couldn’t import the .msh file).

Is this method not supported anymore in the newer version or is there something which I miss? In the exported .msh file I didn’t change/removed anything. If I should change to older version which combination of GMSH and FF++ versions would work?

Best regards,
Robert

Without checking myself:

It seems your mesh is in the x-z plane.
Could you try using the x-y plane instead.

I have a suspicion that freefem assumes 2D-meshes to lie in the x-y plane.

Good luck.

Dear gero,
Thank you for the interest. I tried with a geometry/mesh in the x-y plane and still couldn’t make it run. I attach the new .geo and .msh files and scripts.

// Gmsh project created on Wed Feb 02 17:46:08 2022
//+
Point(1) = {-0, 0, 0, 1.0};
//+
Point(1) = {-0, 0, 0, 1.0};
//+
Point(2) = {0, 0.5, 0, 1.0};
//+
Point(3) = {0.5, 0.5, 0, 1.0};
//+
Point(4) = {0.5, 0, 0, 1.0};
//+
Line(1) = {2, 3};
//+
Line(2) = {3, 4};
//+
Line(3) = {4, 1};
//+
Line(4) = {1, 2};
//+
Curve Loop(1) = {4, 1, 2, 3};
//+
Plane Surface(1) = {1};
//+
Physical Curve("Inlet", 1) = {4};
//+
Physical Curve("Outlet", 2) = {2};
//+
Physical Curve("Upper", 3) = {1};
//+
Physical Curve("Lower", 4) = {3};
$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
4
1 1 "Inlet"
1 2 "Outlet"
1 3 "Upper"
1 4 "Lower"
$EndPhysicalNames
$Nodes
4
1 -0 0 0
2 0 0.5 0
3 0.5 0.5 0
4 0.5 0 0
$EndNodes
$Elements
4
1 1 2 3 1 2 3
2 1 2 2 2 3 4
3 1 2 4 3 4 1
4 1 2 1 4 1 2
$EndElements

Layer edit: It works, I changed the Physical Curves tag numbers, probably there was a conflict between Physical Curve(“Upper”, 3) = {1}; and Plane Surface(1) = {1};
Best regards,
Robert