What is wrong with my mesh?

[SvenH]

Hi,
I have created a mesh with gmsh and I am having trouble importing it into my script.
Dolfin has trouble reading the facets xdmf file and shows the following error message:

*** Error:   Unable to find entity in map.
*** Reason:  Error reading MeshValueCollection.
*** Where:   This error was encountered inside HDF5File.cpp. 

this is my .geo file:

// Gmsh project created on Tue Jan 18 10:52:11 2022
SetFactory("OpenCASCADE");
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(3) = {0.007421, 0.001296, 0, 1.0};
//+
Point(5) = {0.007421, 0.002798, 0, 1.0};
//+
Point(6) = {0.007421, 0.003812, 0, 1.0};
//+
Point(7) = {0.007421, 0.002748, 0, 1.0};
//+
Point(8) = {0.006405, 0.001273, 0, 1.0};
//+
Point(9) = {0.006405, 0.002727, 0, 1.0};
//+
Point(10) = {0.006405, 0.002783, 0, 1.0};
//+
Point(11) = {0.006405, 0.003791, 0, 1.0};
//+
Point(12) = {0.005553, 0.001266, 0, 1.0};
//+
Point(13) = {0.005553, 0.002718, 0, 1.0};
//+
Point(14) = {0.005553, 0.002769, 0, 1.0};
//+
Point(15) = {0.005553, 0.003788, 0, 1.0};
//+
Point(16) = {0.004349, 0.001228, 0, 1.0};
//+
Point(17) = {0.004349, 0.002597, 0, 1.0};
//+
Point(18) = {0.004349, 0.002660, 0, 1.0};
//+
Point(19) = {0.004349, 0.003704, 0, 1.0};
//+
Point(20) = {0.004120, 0.000980, 0, 1.0};
//+
Point(21) = {0.004120, 0.002534, 0, 1.0};
//+
Point(22) = {0.004120, 0.002596, 0, 1.0};
//+
Point(23) = {0.004120, 0.003669, 0, 1.0};
//+
Point(24) = {0.003964, 0.00-0.000688, 0, 1.0};
//+
Point(25) = {0.003964, 0.002482, 0, 1.0};
//+
Point(26) = {0.003964, 0.002542, 0, 1.0};
//+
Point(27) = {0.003964, 0.003631, 0, 1.0};
//+
Point(28) = {0.003535, -0.000730, 0, 1.0};
//+
Point(29) = {0.003535, 0.002272, 0, 1.0};
//+
Point(30) = {0.003535, 0.002336, 0, 1.0};
//+
Point(31) = {0.003535, 0.003501, 0, 1.0};
//+
Point(32) = {0.003227, -0.000693, 0, 1.0};
//+
Point(33) = {0.003227, 0.002033, 0, 1.0};
//+
Point(34) = {0.003227, 0.002101, 0, 1.0};
//+
Point(35) = {0.003227, 0.003370, 0, 1.0};
//+
Point(36) = {0.002993, 0.001490, 0, 1.0};
//+
Point(37) = {0.002993, 0.003231, 0, 1.0};
//+
Point(38) = {0.002993, -0.000661, 0, 1.0};
//+
Point(39) = {0.002790, 0.003081, 0, 1.0};
//+
Point(40) = {0.002670, -0.000627, 0, 1.0};
//+
Point(41) = {0.002670, 0.000256, 0, 1.0};
//+
Point(42) = {0.002395, 0, 0, 1.0};
//+
Point(43) = {0, -0.000660, 0, 1.0};
//+
Line(1) = {1, 42};
//+
Point(44) = {0, 0, 0, 1.0};
//+
Point(45) = {0, 0, -0, 1.0};
//+
Point(46) = {0, 0, -0, 1.0};
//+
Point(47) = {0, 0, -0, 1.0};
//+
Point(48) = {0, 0, -0, 1.0};
//+
Point(49) = {0, 0, -0, 1.0};
//+
Point(50) = {0, 0, -0, 1.0};
//+
Point(51) = {0, 0, -0, 1.0};
//+
Point(52) = {0, 0, -0, 1.0};
//+
Point(53) = {0, 0, -0, 1.0};
//+
Point(54) = {0, 0, -0, 1.0};
//+
Point(55) = {0, 0, -0, 1.0};
//+
Point(56) = {0, 0, -0, 1.0};
//+
Point(57) = {0, 0, -0, 1.0};
//+
Point(58) = {0, 0, -0, 1.0};
//+
Point(59) = {0, 0, -0, 1.0};
//+
Point(60) = {0, 0, -0, 1.0};
//+
Line(2) = {41, 39};
//+
Spline(3) = {39, 37, 35, 31};
//+
Spline(4) = {31, 27, 23, 19};
//+
Spline(5) = {19, 15, 11, 6};
//+
Line(6) = {6, 5};
//+
Spline(7) = {5, 10, 14};
//+
Spline(8) = {14, 18, 22};
//+
Spline(9) = {22, 26, 30, 34, 36};
//+
Spline(10) = {36, 33, 29, 25};
//+
Spline(11) = {25, 21, 17, 13};
//+
Spline(12) = {13, 9, 7};
//+
Line(13) = {3, 7};
//+
Spline(14) = {3, 8, 12};
//+
Line(15) = {12, 16};
//+
Line(16) = {16, 20};
//+
Line(17) = {20, 24};
//+
Line(18) = {24, 28};
//+
Line(19) = {28, 32};
//+
Line(20) = {32, 38};
//+
Line(21) = {38, 40};
//+
Line(22) = {40, 43};
//+
Point(61) = {0, 0, -0, 1.0};
//+
Point(62) = {0, 0, -0, 1.0};
//+
Point(63) = {0, 0, -0, 1.0};
//+
Point(64) = {0, 0, -0, 1.0};
//+
Point(65) = {0, 0, -0, 1.0};
//+
Point(66) = {0.002598, 0.000107, 0, 1.0};
//+
Spline(23) = {41, 66, 42};
//+
Point(67) = {0, -0.000284, 0, 1.0};
//+
Point(68) = {0.000698, -0.000284, 0, 1.0};
//+
Point(69) = {0.001298, -0.000322, 0, 1.0};
//+
Point(70) = {0.001522, -0.000314, 0, 1.0};
//+
Point(71) = {0.003118, -0.000029, 0, 1.0};
//+
Point(72) = {0.003166, 0.000059, 0, 1.0};
//+
Point(73) = {0.003118, 0.000710, 0, 1.0};
//+
Line(25) = {67, 68};
//+
Spline(26) = {68, 69, 70};
//+
Spline(30) = {36, 73, 72, 71};
//+
Line(31) = {71, 70};
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Curve{1}; Curve{23}; Curve{2}; Curve{3}; Curve{4}; Curve{5}; Curve{6}; Curve{7}; Curve{8}; Curve{9}; Curve{30}; Curve{31}; Curve{26}; Curve{25}; Curve{22}; Curve{21}; Curve{20}; Curve{19}; Curve{18}; Curve{17}; Curve{16}; Curve{15}; Curve{14}; Curve{13}; Curve{12}; Curve{11}; Curve{10}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{66}; Point{37}; Point{35}; Point{27}; Point{23}; Point{11}; Point{33}; Point{29}; Point{25}; Point{21}; Point{9}; Point{8}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{34}; Point{30}; Point{26}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{14}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{14}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{10}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{15}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{70}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{69}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{71}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{72}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{73}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{26}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{22}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{18}; }
}
//+
Symmetry {1, 0, 0, 0} {
  Duplicata { Point{17}; }
}
Curve Loop(1) = {1,23,2,3,4,5,6,7,8,9,30,31,26,25,45,44,43,42,41,40,39,38,37,36,35,34,33,32};
Plane Surface (1) = {1};
Curve Loop(2) = {25,26,31,30,10,11,12,13,14,15,16,17,18,19,20,21,22,46,47,48,49,50,51,52,53,54,55,56,57,58,42,43,44,45};
Plane Surface(2) = {2};
//+
Physical Curve(1) = {35, 37, 36, 34, 32, 33, 1, 23, 2, 3, 5, 4};
//+
Physical Curve(2) = {38};
//+
Physical Curve(3) = {6};
//+
Physical Curve(4) = {8, 10, 9, 7, 11, 12};
//+
Physical Curve(5) = {13};
//+
Physical Curve(6) = {14, 15, 16, 17, 18, 19, 20, 21, 22, 46, 47, 48, 49, 50, 51, 52, 53, 54};
//+
Physical Curve(7) = {55};
//+
Physical Curve(8) = {56, 57, 58, 41, 40, 39};
//+
Physical Curve(11) = {42, 43, 44, 45, 25, 26, 31, 30};
//+
Physical Surface(9) = {1};
//+
Physical Surface(10) = {2};

I am using the following code to import the mesh:

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, "Boundary")

I’m thankful for any suggestions for what might be wrong with my mesh

Edit 1
Code for mesh conversion:

import meshio
import numpy as np

# import mesh file from gmesh
msh = meshio.read('mesh.msh')

# make
line_cells = []
for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif  cell.type == "line":
        if len(line_cells) == 0:
            line_cells = cell.data
        else:
            line_cells = np.vstack([line_cells, cell.data])

line_data = []
for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "line":
        if len(line_data) == 0:
            line_data = msh.cell_data_dict["gmsh:physical"][key]
        else:
            line_data = np.vstack([line_data, msh.cell_data_dict["gmsh:physical"][key]])
    elif key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]

triangle_mesh = meshio.Mesh(points=msh.points[:, :2], cells={"triangle": triangle_cells})
line_mesh =meshio.Mesh(points=msh.points[:, :2],
                           cells=[("line", line_cells)],
                           cell_data={"Boundary":[line_data]})
meshio.write("mesh.xdmf", triangle_mesh)

meshio.xdmf.write("mf.xdmf", line_mesh)

[dokken]

You need to add the script/code you used to convert the geo file to xdmf as well for anyone to properly assist you.

[SvenH]

Thanks, I added the code at the bottom of the original post. I create the mesh and export the .msh file with the gmsh GUI.

[dokken]

The issue with your mesh is that several of the boundaries you have marked is not part of a volume mesh, see attached images:

Note that you can easily inspect your mesh and boundary xdmf file using Paraview.

[SvenH]

The small gaps without mesh on the left and right side are intended. I didn’t think this would be a problem. Thank you, I will look into this.

[dokken]

The point is that there is no mesh (cells) attached to those facets, which will cause the reading of boundary markers to fail, as it cannot find any cell with those facets.

[SvenH]

Sorry, I’m not quite sure whether I understand correctly. This is what the mesh looks like on my end in paraview:

which at least looks like what I want it to.

[dokken]

Using a slightly different conversion syntax, I can reproduce your domain with XDMF.

import meshio
import numpy as np


def create_mesh(mesh: meshio.Mesh, cell_type: str, data_name: str = "name_to_read",
                prune_z: bool = False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
                           data_name: [cell_data]})
    return out_mesh


# import and convert msh file to xdmf
msh = meshio.read('mesh.msh')
line_mesh = create_mesh(msh, "line", data_name="Boundary", prune_z=True)
meshio.write("mf.xdmf", line_mesh)
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

The issue is that you have tagged facets that are not actually in the mesh.
This can be shown by calling:

points_in_facets = np.unique(np.sort(msh.get_cells_type("line").reshape(-1)))
points_in_cells = np.sort(
    np.unique(msh.get_cells_type("triangle").reshape(-1)))
vertices_in_facets_not_cells = np.isin(
    points_in_facets, points_in_cells, invert=True)
print(msh.points[np.flatnonzero(vertices_in_facets_not_cells)])

All the vertices printed here are vertices that are found in your facet mesh, but not in your actual mesh.

The problem is probably how you duplicate points of your mesh. However, as I do not really use the GMSH geo format (I rather use the gmsh python API), I don’t know a workaround for your problem.

[SvenH]

Thank you for your time, Dokken. The problem arose from the transform option called “symmetry”. Apparently, you are not supposed to use curves as entities, just points and then connect them manually. It works perfectly fine now.

[Rafael_Ferro]

What about a 3D mesh? I made the following code but the faces do not appear, the file is blank.

import meshio
import numpy as np


def create_mesh(mesh: meshio.Mesh, cell_type: str, data_name: str = "name_to_read",
                prune_z: bool = False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
                           data_name: [cell_data]})
    return out_mesh


# import and convert msh file to xdmf
msh = meshio.read('test.msh')
#line_mesh = create_mesh(msh, "line", data_name="Boundary", prune_z=True)
#meshio.write("test_mf.xdmf", line_mesh)

facet_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("test_facet.xdmf", facet_mesh)

triangle_mesh = create_mesh(msh, "tetra", prune_z=True)
meshio.write("test_mesh.xdmf", triangle_mesh)
```

[dokken]

@Rafael_Ferro

facet_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("test_facet.xdmf", facet_mesh)

triangle_mesh = create_mesh(msh, "tetra", prune_z=True)
meshio.write("test_mesh.xdmf", triangle_mesh)

Why are you pruning the z coordinate if your mesh is 3D?

As you haven’t supplied the msh file or the geo file, there isn’t much more people can do to help you.

[Rafael_Ferro]

I think I got it now, I changed 2 for 3 on the line:

points = mesh.points[:, :3] if prune_z else mesh.points

is this correct?

gmsh’s code is simple:

// Gmsh project created on Fri Feb 17 16:59:21 2023
SetFactory("OpenCASCADE");
//+
Box(1) = {0, 0, 0, 1, 1, 0.02};
//+
Physical Volume(13) = {1};
//+
Physical Surface(14) = {3, 4};

[dokken]

I cannot reproduce any errors with:

import meshio
import numpy as np


def create_mesh(mesh: meshio.Mesh, cell_type: str, data_name: str = "name_to_read",
                prune_z: bool = False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
                           data_name: [cell_data]})
    return out_mesh


# import and convert msh file to xdmf
msh = meshio.read('mesh.msh')

print(msh)
facet_mesh = create_mesh(msh, "triangle")
meshio.write("test_facet.xdmf", facet_mesh)

triangle_mesh = create_mesh(msh, "tetra")
meshio.write("test_mesh.xdmf", triangle_mesh)

which prints

<meshio mesh object>
  Number of points: 199
  Number of cells:
    triangle: 16
    triangle: 16
    tetra: 512
  Cell sets: gmsh:bounding_entities
  Point data: gmsh:dim_tags
  Cell data: gmsh:physical, gmsh:geometrical

Please look at all the other posts regarding gmsh/meshio

Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio

As far as I can tell you would like to have a 3D mesh (tetrahedrons). Then you also need to mark the volumes: // Gmsh project created on Tue May 18 16:34:55 2021 SetFactory("OpenCASCADE"); //+ Circle(1) = {0, 00, 00, 0.6, 0, 2*Pi}; //+ Circle(2) = {0, 00, 00, 0.4, 0, 2*Pi}; //+ Curve Loop(1) = {1}; //+ Curve Loop(2) = {2}; //+ Plane Surface(1) = {1, 2}; //+ Extrude {0, 0, 10} { Curve{2}; Curve{1}; Layers {1}; Recombine; } //+ Curve Loop(5) = {6}; //+ Curve Loop(6) = {4}; //+ Plane Surface(4) …

Need help converting GMSH to FEniCS

Please note that you have some errors in your geo-file, i.e. Correcting these, and adding Physical Surface("surf", 1) = {2}; at the end, the following file (mesh.geo) // Gmsh project created on Thu Nov 18 16:23:27 2021 SetFactory("OpenCASCADE"); Rectangle(1) = {0, 0, 0, 18, 6, 0}; Point(5) = {3, 2.8, 0, 0.05}; Point(6) = {3, 3.2, 0, 0.05}; Point(7) = {3.3464, 3.4, 0, 0.05}; Point(8) = {3.3464, 2.6, 0, 0.05}; Point(9) = {3.6928, 2.8, 0, 0.05}; Point(10) = {3.6928, 3.2, 0, 0.05}; Point(11) = …

as they contain a lot of tips and tricks.