Exact Procedure to work with a .geo file on FEniCS (with meshio)?

mesh
1

Hi,
I would like to work on “real image” with FEniCS. So I try gmsh with simple geometry to understand the process.
If I understood correctly the other threads of the forum, there are three steps : convert the geo file into a msh file ; convert the msh into a xdmf file with meshio ; read on FEniCS.
For the first step, should I use the line * gmsh -2 mesh.geo* with pygmsh?
For the second step (meshio), I’m using a msh file find on the internet with the lines import meshio and msh = meshio.read(“circle_2d.msh”) but this give an assertion error.
If it is not too much asked, will it be possible to send me the complete code (python) to achieve these three steps?

Thank you in advance,
AM

2

Hi,
Does this help:

1 Like
3

Hi Dokken,

Thank you for your answer.

When I use pygmsh (on my FEniCS environment), your geometry (that I called testDokken.geo) and the line gmsh -2 testDokken.geo, I get an error message (invalid synthax).

What do I need to import for this line to work?

AM

4

Could you please supply the full error message? Also, Which Gmsh version are you using?

5

First of all, I am really sorry to bother you with a subject that is not related to fenics.

I have the 4.3.6 version of gmsh (pygmsh) and the 2.3.3 version of meshio.

Also, the message error is an invalid syntax on the line gmsh -2 testDokken.geo

AM

6

Dear AMM,
Which version of gmsh is installed on your computer. this can be checked by typing gmsh --version in your terminal.
With the following geo file:

// the square
Point(1) = {0, 0, 0,0.1};
Point(2) = {0, 1, 0,0.1};
Point(3) = {1, 1, 0,0.1};
Point(4) = {1, 0, 0,0.1};
Line(1) = {1, 4};
Line(2) = {4, 3};
Line(3) = {3, 2};
Line(4) = {2, 1};
// creating the surfaces
Line Loop(6) = {2, 3, 4, 1};
Plane Surface(8) = {6};

I am able to call gmsh -2 filename.geo for both gmsh 3.0.6 and 4.2.3

7

Dokken,

The Lab did not reinstall gmsh on all computers when they change the ubuntu version : the problem comes from there.

Thanks for your answers and again sorry for making you waste your time on such a trivial question.

8

Hi dokken,

So, I ran a code (taken on a thread of the forum) but I get an error message :

 meshio.write("mf.xdmf", meshio.Mesh(points=msh.points, cells={"line": msh.cells["line"]},
 KeyError: 'line'

I tried with “triangle” instead of “line” (although my geometry is 2d) and it recognizes the syntaxe.

Do you know where this error may come from?

The complete code :
msh = meshio.read(“EssaiDokken.msh”)
meshio.write(“mesh.xdmf”, meshio.Mesh(points=msh.points, cells={“triangle”: msh.cells[“triangle”]}))
meshio.write(“mf.xdmf”, meshio.Mesh(points=msh.points, cells={“line”: msh.cells[“line”]},
cell_data={“line”: {“name_to_read”: msh.cell_data[“line”][“gmsh:physical”]}}))

mesh = Mesh()
with XDMFFile(“msh.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 = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)
File(“dolfinmesh.pvd”).write(mesh)
File(“dolfincellfunc.pvd”).write(mf)

9

Hi,
Which mesh are you using?
You can only write things as a facet-function if you have defined physical lines in gmsh.

10

Dokken,

I use as geometry one of your example (I can not find the exact thread anymore) and I “mesh” with the command gmsh -2 mygeometry.geo -format msh2

The geometry :

// Gmsh project created on Tue Mar 26 15:00:10 2019
//+
// THe square
Point(1) = {0, 0, 0,0.1};
Point(2) = {0, 1, 0,0.1};
Point(3) = {1, 1, 0,0.1};
Point(4) = {1, 0, 0,0.1};
Line(1) = {1, 4};
Line(2) = {4, 3};
Line(3) = {3, 2};
Line(4) = {2, 1};
// creating the surfaces
Line Loop(6) = {2, 3, 4, 1};
Plane Surface(8) = {6};

Physical Surface(2) = {8};

11

To be able to read line data, you Need to add ‘Physical Line(5) = {1};’ . This adds a physical facet marker with value 5 on Line 1.

12

Hi dokken,

It now recognizes the “line” synthax.

From here, I try my luck but do not feel obliged to answer this question since it no longer concerns FEniCS.
In fact I can not display my mesh with “plot (mesh)” and he gives me an error message (see screenshot in attachment).
Do you know the cause of this?

Thank you

Edit : Paraview does the job well.
Thanks again for all your answers!