Gmsh import error

mesh
1

Hi,

with Gmesh I build a 2D-mesh using:

SetFactory("OpenCASCADE");

ln1 = 0.2; 
ln2 = 3.0;

lc1 = 0.04;
lc3 = 0.005;
lc4 = 0.01;

Point(1) = {0,0,0,lc3}; 
Point(2) = {ln1,0.0,0,lc1};
Point(3) = {ln1,0.6,0,lc1};
Point(4) = {ln1,ln2,0,lc1};     
Point(5) = {0,ln2,0,lc1};

Line(1) = {1,2}; 
Line(2) = {2,3};  
Line(3) = {3,4};  
Line(4) = {4,5}; 
Line(5) = {5,1}; 

Curve Loop(5) = {1,2,3,4,5}; 
Plane Surface(6) = {5};

Field[1] = Distance;
Field[1].CurvesList = {1,2};
Field[1].NumPointsPerCurve = ln2/lc3; 

Field[2] = Distance;
Field[2].CurvesList = {3};
Field[2].NumPointsPerCurve = ln2/lc4; 

Field[3] = Threshold;
Field[3].InField = 1;
Field[3].SizeMin = lc3;
Field[3].SizeMax = lc1;
Field[3].DistMin = 0.01;
Field[3].DistMax = 0.1;

Field[4] = Threshold;
Field[4].InField = 2;
Field[4].SizeMin = lc4;
Field[4].SizeMax = lc1;
Field[4].DistMin = 0.01;
Field[4].DistMax = 0.2;

Field[5] = Min;
Field[5].FieldsList = {3,4};
Background Field = 5;

Mesh.MeshSizeExtendFromBoundary = 0;
Mesh.MeshSizeFromPoints = 0;
Mesh.MeshSizeFromCurvature = 0;

and further

meshio-convert --input-format gmsh --output-format xdmf mymesh.gmsh mymesh.xdmf

which should result in a mesh like the attached file.

mesh
mesh1442×161 26.9 KB

Doing the following:

mesh = Mesh()
filename = "mymesh.xdmf"
f = XDMFFile(MPI.comm_world, filename)
f.read(mesh)
plot(mesh)

gives the error:

*** Error:   Unable to recognise cell type.
*** Reason:  Unknown value "mixed".
*** Where:   This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown

What did I do wrong and how can I just replace

mesh = RectangleMesh(p0,p1,npX,npY,'left')

by the Gmsh mesh?

Thank’s ins advance.

2

Please search the on the forum for your error message before posting a new question.
If you search for

Unknown value “mixed”.

The follow thread comes up.

3

Sorry, I saw this thread but was confused by “mesh.xml to mesh.xdmf” which was the opposite direction for me. I check this out. Thank you.

4

I added now physical line and surface commands (as described in the other post) even if I’m not totally sure what I’m doing (sorry):
…
Physical Line(10) = {1};
Physical Line(11) = {2};
Physical Line(12) = {3};
Physical Line(13) = {4};
Physical Line(14) = {5};
Physical Surface(15) = {6};
…

I still get problems with xdmf. I’m doing:

meshio-convert --input-format gmsh --output-format xdmf --prune-z-0 mymesh.gmsh mymesh.xdmf
meshio-convert --input-format gmsh --output-format dolfin-xml --prune-z-0 mymesh.gmsh mymesh.xml

Even if there is a warning on dolfin output it works now to do:
mesh = Mesh("mymesh.xml")

While reading xml finally is ok, this

mesh = Mesh()
filename = "mymesh.xdmf"
f = XDMFFile(MPI.comm_world, filename)
f.read(mesh)

results in the described error message Unknown value "mixed". In the xdmf file there is a setting

<Topology TopologyType="Mixed" ...

If I change Mixed to Triangle, than the read command crashes.

So finally I can live with the dolfin-xml export but don’t feel comfortable because the dolfin-xml type seems to be deprecated.

5

It seems like you have not read the post I linked you to. You cannot use the meshio-convert script out of the box for XDMF, and you have to use the meshio python interface. See the second section of: Mesh generation and conversion with GMSH and PYGMSH | Jørgen S. Dokken
where I read in a mesh:

import meshio
mesh_from_file = meshio.read("mesh.msh")

extract data for a given cell type

def create_mesh(mesh, cell_type, prune_z=False):
    cells = numpy.vstack([cell.data for cell in mesh.cells if cell.type==cell_type])
    cell_data = numpy.hstack([mesh.cell_data_dict["gmsh:physical"][key]
                         for key in mesh.cell_data_dict["gmsh:physical"].keys() if key==cell_type])
    # Remove z-coordinates from mesh if we have a 2D cell and all points have the same third coordinate
    points= mesh.points
    if prune_z:
        points = points[:,:2]
    mesh_new = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return mesh_new

and write a separate mesh file for the cell and facet data

line_mesh = create_mesh(mesh_from_file, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)
6

Ok, thank you. As not native english speaker it is sometimes difficult to get the main point of a post.