Read complex geometry using meshio

mesh
1

Hi everyone,

I’m trying to read in a slightly complex geometry from the following GMSH file, in which i declared only small parts of it as boundary, the rest is my interior:

$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
3
2 2 "Neumann"
2 3 "Robin"
3 1 "Volume"
$EndPhysicalNames
$Nodes
2191
1 1.45935 0.948814 -0.69911
2 1.44345 0.948814 -0.69911
3 1.29746 0.23 -0.63
4 1.34746 0.23 -0.68
...
2190 3.34704 0.23 0.712405
2191 3.36103 0.23 0.742412
$EndNodes
$Elements
7503
1 2 2 3 3 55 60 62
2 2 2 3 3 55 63 60
3 2 2 3 3 63 53 59
4 2 2 3 3 63 59 60
5 2 2 3 3 3 5 55
6 2 2 3 3 5 63 55
7 2 2 3 3 5 6 53
8 2 2 3 3 5 53 63
9 2 2 3 3 75 78 169
10 2 2 3 3 78 234 469
11 2 2 3 3 78 74 77
12 2 2 3 3 77 234 78
13 2 2 3 3 234 77 170
14 2 2 3 3 77 73 170
15 2 2 3 3 169 234 170
16 2 2 3 3 169 170 172
17 2 2 2 2 1620 1621 1622
18 2 2 2 2 1621 1876 1622
19 2 2 2 2 1725 1876 1621
20 4 2 1 1 10 23 1 11
21 4 2 1 1 23 1 11 35
22 4 2 1 1 35 23 1 2
23 4 2 1 1 35 37 23 2
...
7501 4 2 1 1 2185 2188 2190 2189
7502 4 2 1 1 2187 2185 2190 2186
7503 4 2 1 1 2188 2191 2190 2189
$EndElements

with the following python code

import dolfin as df
pathname = './CarInteriorMesh/'
filename = 'Mesh'  # file to be opened

msh = meshio.read(pathname + filename + '.msh')

meshio.write(pathname + filename + '_xdmf.xdmf',
             meshio.Mesh(points=msh.points,
                             cells={'tetra': msh.cells['tetra']},
                             cell_data={
                                 'tetra': {'name_to_read': msh.cell_data['tetra']['gmsh:physical']}}))

meshio.write(pathname + filename + '_Boundary_xdmf.xdmf',
             meshio.Mesh(points=msh.points,
                         cells={'triangle': msh.cells['triangle']},
                         cell_data={
                            'triangle': {'name_to_read': msh.cell_data['triangle']['gmsh:physical']}}))

mesh = df.Mesh()

with df.XDMFFile(pathname + filename + '_xdmf.xdmf') as infile:
    infile.read(mesh)
mvc = df.MeshValueCollection('size_t', mesh, 3)
with df.XDMFFile(pathname + filename + '_Boundary_xdmf.xdmf') as infile:
    infile.read(mvc, 'name_to_read')
mf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)
with df.XDMFFile(pathname + filename + '_xdmf.xdmf') as infile:
    infile.read(mvc, 'name_to_read')
cf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_excited = df.Measure('ds', domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_admittance = df.Measure('ds', domain=mesh, subdomain_data=mf, subdomain_id=3)
dx_volume = df.Measure('dx', domain=mesh, subdomain_data=cf, subdomain_id=1)

Edit: The error occurs in

 infile.read(mvc, 'name_to_read')

and i get the following error:

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

I just don’t understand why FEniCS/Dolfin is not able to find an entity, i have declared everything i need in my gmsh-file.
I am able to read other GMSH-generated meshes without problems but this one won’t work.

Thanks in advance, any help is highly appreciated!

2

For the boundary data, you Need to create a separate MeshValueCollection, with argument 2, instead of 3.

2 Likes
3

When i write

mesh = df.Mesh()

with df.XDMFFile(pathname + filename + '_xdmf.xdmf') as infile:
    infile.read(mesh)
mvc = df.MeshValueCollection('size_t', mesh, 3)
mvc_b = df.MeshValueCollection('size_t', mesh, 2)
with df.XDMFFile(pathname + filename + '_Boundary_xdmf.xdmf') as infile:
    infile.read(mvc_b, 'name_to_read')
mf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc_b)
with df.XDMFFile(pathname + filename + '_xdmf.xdmf') as infile:
    infile.read(mvc, 'name_to_read')
cf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)

i still get the same error message

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

Could you supply the full error message, right now this does not tell which line causes the error.
The syntax for reading your mesh and cell/facet functions should be the same as in Complex PDE on GMSH generated mesh

1 Like
5

Code:

import dolfin as df
import meshio

pathname = './Mesh/'
filename = 'Mesh'
dim = 3

msh = meshio.read(pathname + filename + '.msh')

meshio.write(pathname + filename + '_xdmf.xdmf',
             meshio.Mesh(points=msh.points,
                         cells={'tetra': msh.cells['tetra']},
                         cell_data={
                             'tetra': {'name_to_read': msh.cell_data['tetra']['gmsh:physical']}}))

meshio.write(pathname + filename + '_Boundary_xdmf.xdmf',
             meshio.Mesh(points=msh.points,
                         cells={'triangle': msh.cells['triangle']},
                         cell_data={
                             'triangle': {'name_to_read': msh.cell_data['triangle']['gmsh:physical']}}))

mesh = df.Mesh()
# print(type(mesh))

with df.XDMFFile(pathname + filename + '_xdmf.xdmf') as infile:
    infile.read(mesh)
mvc = df.MeshValueCollection('size_t', mesh, dim)
mvc_boundary = df.MeshValueCollection('size_t', mesh, dim - 1)
# mvc.dim()
with df.XDMFFile(pathname + filename + '_Boundary_xdmf.xdmf') as infile:
    infile.read(mvc_boundary, 'name_to_read')  # this is where the error occurs
mf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc_boundary)
with df.XDMFFile(pathname + filename + '_xdmf.xdmf') as infile:
    infile.read(mvc, 'name_to_read')
cf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_excited = df.Measure('ds', domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_admittance = df.Measure('ds', domain=mesh, subdomain_data=mf, subdomain_id=3)
dx_volume = df.Measure('dx', domain=mesh, subdomain_data=cf, subdomain_id=1)

Error message:

    infile.read(mvc_boundary, 'name_to_read')
RuntimeError: 
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to find entity in map.
*** Reason:  Error reading MeshValueCollection.
*** Where:   This error was encountered inside HDF5File.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  473fc681d80577c29aa757c3b36ffec0b3f067e9
*** -------------------------------------------------------------------------
6

Hi again,
even when i’m writing my code like this, i still get the same error message as above (i commented the line where the error occurs).

import dolfin as df
import meshio
import matplotlib.pyplot as plt

pathname = './Mesh/'
filename = 'Mesh'
dim = 3

msh = meshio.read(pathname + filename + '.msh')

meshio.write(pathname + filename + '_xdmf.xdmf',
             meshio.Mesh(points=msh.points,
                         cells={'tetra': msh.cells['tetra']},
                         cell_data={
                             'tetra': {'name_to_read': msh.cell_data['tetra']['gmsh:physical']}}))

meshio.write(pathname + filename + '_Boundary_xdmf.xdmf',
             meshio.Mesh(points=msh.points,
                         cells={'triangle': msh.cells['triangle']},
                         cell_data={
                             'triangle': {'name_to_read': msh.cell_data['triangle']['gmsh:physical']}}))

mesh = df.Mesh()

with df.XDMFFile(pathname + filename + '_xdmf.xdmf') as infile:
    infile.read(mesh)
mvc = df.MeshValueCollection('size_t', mesh)
mvc_boundary = df.MeshValueCollection('size_t', mesh, dim - 1)
with df.XDMFFile(pathname + filename + '_Boundary_xdmf.xdmf') as infile:
     infile.read(mvc_boundary, 'name_to_read')  # this is where the error occurs
mf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc_2 = df.MeshValueCollection("size_t", mesh, dim)
with df.XDMFFile(pathname + filename + '_xdmf.xdmf') as infile:
    infile.read(mvc_2, 'name_to_read')
cf = df.cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_excited = df.Measure('ds', domain=mesh, subdomain_data=mf, subdomain_id=2)
ds_admittance = df.Measure('ds', domain=mesh, subdomain_data=mf, subdomain_id=3)
dx_volume = df.Measure('dx', domain=mesh, subdomain_data=cf, subdomain_id=1)

df.plot(mesh)
plt.show()

I really don’t understand what i’m doing wrong.

7

Since you have not supplied the full mesh-file, i really can’t Do anything other than pointing you to your previous thread, Where you were able to load facet and cell functions. Have you tried running that code with both the mesh that that i worked for and the one that you are working on now?

1 Like
8

The thing is, that not the whole surface of the mesh is written in the .msh file as you can see here:

$Elements
7503
1 2 2 3 3 55 60 62
2 2 2 3 3 55 63 60
3 2 2 3 3 63 53 59
4 2 2 3 3 63 59 60
5 2 2 3 3 3 5 55
6 2 2 3 3 5 63 55
7 2 2 3 3 5 6 53
8 2 2 3 3 5 53 63
9 2 2 3 3 75 78 169
10 2 2 3 3 78 234 469
11 2 2 3 3 78 74 77
12 2 2 3 3 77 234 78
13 2 2 3 3 234 77 170
14 2 2 3 3 77 73 170
15 2 2 3 3 169 234 170
16 2 2 3 3 169 170 172
17 2 2 2 2 1620 1621 1622
18 2 2 2 2 1621 1876 1622
19 2 2 2 2 1725 1876 1621
20 4 2 1 1 10 23 1 11
21 4 2 1 1 23 1 11 35
22 4 2 1 1 35 23 1 2
23 4 2 1 1 35 37 23 2
...
7501 4 2 1 1 2185 2188 2190 2189
7502 4 2 1 1 2187 2185 2190 2186
7503 4 2 1 1 2188 2191 2190 2189
$EndElements

As you can see, there are only 19 facet elements out of 7503 total elements, where 7484 elements are volume elements, since i’m not able to generate facet data for the entire geometry.
So basically, only a small part of the boundary is available in the msh file.
May this cause the problem in FEniCS?

9

Why is that? Your mesh file should contain the full mesh.

1 Like
10

I only have the original file as .nas file format with only volume information, but no information of the boundary, when opening the .nas file in the GMSH GUI and saving the file as .msh file format, consequently, there is only volume information available.
I was now trying to find only a few nodes manually on the boundary and mark them with different boundary types in the .msh file.

So it’s not possible to read the mesh in FEniCS correctly if i don’t have the whole boundary as information?

Thank you for your help anyway, i really appreciate it!

11

I have not worked with nas files myself, so i dont know how They work. However, if you do not have the whole boundary information, i think reading a facet-function will be tricky. I suggest asking the meshio-developers if They know how to handle those files. Or the Gmsh-developers.

1 Like
12

OK thank you very much anyway for your help, i will try it like as you said!

Thanks again!

13

I’m experiencing the same error:

*** Error: Unable to find entity in map.
*** Reason: Error reading MeshValueCollection.
*** Where: This error was encountered inside HDF5File.cpp.
*** Process: 0
*** DOLFIN version: 2019.1.0
*** Git changeset: c5b9b269f4a6455a739109e3a66e036b5b8412f5

Here’s my gmsh geo file:

//+
SetFactory(“OpenCASCADE”);
Circle(1) = {0, 0, 0, 1, 0, 2Pi};
Circle(2) = {0, 0, 0, 0.9, 0, 2Pi};
//+
Extrude {0, 0, 5} {
Curve{1};
}
Extrude {0, 0, 5} {
Curve{2};
}
//+
Curve Loop(3) = {4};
//+
Curve Loop(4) = {6};
//+
Plane Surface(3) = {3, 4};
//+
Curve Loop(5) = {1};
//+
Curve Loop(6) = {2};
//+
Plane Surface(4) = {5, 6};
//+
Surface Loop(1) = {2, 3, 1, 4};
//+
Volume(1) = {1};
//+
Physical Surface(“Bottom”) = {3};
//+
Physical Surface(“Top”) = {4};
//+
Physical Surface(“Inner”) = {2};
//+
Physical Surface(“Outer”) = {1};
//+
Physical Volume(“Main”) = {1};
//+
Physical Surface(“Bottom”) += {3};
//+
Physical Surface(“Inner”) += {4, 2};
//+
Physical Surface(“Outer”) += {1};

I’m using this meshio to read in the msh file as follows:

import meshio
from dolfin import *
import numpy as np

msh=meshio.read("Pipe.msh")

triangle_cells=[]
tetra_cells=[]

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

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]
    elif key == "tetra":
        tetra_data = msh.cell_data_dict["gmsh:physical"][key]
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells})

triangle_mesh =meshio.Mesh(points=msh.points,
                           cells={"triangle": triangle_cells},
                           cell_data={"name_to_read":[triangle_data]})

meshio.write("mesh.xdmf", tetra_mesh)

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



mesh=Mesh()
file_in=XDMFFile("mesh.xdmf")
file_in.read(mesh)

mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")

Any suggestions?

14

For my specific case, the error is caused by the msh file. The geo file that I used contains physical surface groups detached from the volume’s actual surfaces. The error is disappeared with msh generated from the following geo code:

Mesh.Algorithm3D = 4; // 3D mesh algorithm (1=Delaunay, 4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D, 9=R-tree)
SetFactory("OpenCASCADE");

Mesh.CharacteristicLengthMax = 0.5;

//+
Rectangle(1) = {-1, 0, 0, 0.1, 5, 0};
//+
Extrude {{0, 1, 0}, {0, 0, 0}, Pi*2} {
  Curve{3}; Curve{2}; Curve{1}; Curve{4};
}
//+
Surface Loop(1) = {5, 2, 3, 4};
//+
Volume(1) = {1};
//+
Physical Surface("Bottom") = {2};
//+
Physical Surface("Top") = {4};
//+
Physical Surface("Inner") = {7};
//+
Physical Surface("Outer") = {6};
//+
Physical Volume("Main") = {1};