Hello,
I am working on a project using a module that constrains me to “downgrade” to dolfin (instead of using dolfinx for the rest of my project). I created a mesh thanks to the python API of gmsh, and would like to read it in dolfin. So I did my research:
-
first I need to convert my .msh file into .h5. For that I use a script found here 3D problem (sphere) Navier-Stokes - #23 by dokken.
-
then I try to read the mesh but I got the error:
HDF5-DIAG: Error detected in HDF5 (1.12.2) MPI-process 0:
#000: H5L.c line 943 in H5Lexists(): unable to get link info
major: Links
minor: Can't get value
#001: H5VLcallback.c line 5173 in H5VL_link_specific(): unable to execute link specific callback
major: Virtual Object Layer
minor: Can't operate on object
#002: H5VLcallback.c line 5136 in H5VL__link_specific(): unable to execute link specific callback
major: Virtual Object Layer
minor: Can't operate on object
#003: H5VLnative_link.c line 329 in H5VL__native_link_specific(): unable to specific link info
major: Links
minor: Object not found
#004: H5L.c line 3082 in H5L__exists(): path doesn't exist
major: Links
minor: Object already exists
#005: H5Gtraverse.c line 837 in H5G_traverse(): internal path traversal failed
major: Symbol table
minor: Object not found
#006: H5Gtraverse.c line 729 in H5G__traverse_real(): component not found
major: Symbol table
minor: Object not found
HDF5-DIAG: Error detected in HDF5 (1.12.2) MPI-process 0:
#000: H5O.c line 128 in H5Oopen(): unable to open object
major: Object header
minor: Can't open object
#001: H5VLcallback.c line 5386 in H5VL_object_open(): object open failed
major: Virtual Object Layer
minor: Can't open object
#002: H5VLcallback.c line 5353 in H5VL__object_open(): object open failed
major: Virtual Object Layer
minor: Can't open object
#003: H5VLnative_object.c line 58 in H5VL__native_object_open(): unable to open object by name
major: Object header
minor: Can't open object
#004: H5Oint.c line 625 in H5O_open_name(): object not found
major: Object header
minor: Object not found
#005: H5Gloc.c line 442 in H5G_loc_find(): can't find object
major: Symbol table
minor: Object not found
#006: H5Gtraverse.c line 837 in H5G_traverse(): internal path traversal failed
major: Symbol table
minor: Object not found
#007: H5Gtraverse.c line 729 in H5G__traverse_real(): component not found
major: Symbol table
minor: Object not found
HDF5-DIAG: Error detected in HDF5 (1.12.2) MPI-process 0:
#000: H5A.c line 1668 in H5Aexists(): invalid object identifier
major: Invalid arguments to routine
minor: Inappropriate type
#001: H5VLint.c line 1749 in H5VL_vol_object(): invalid identifier type to function
major: Invalid arguments to routine
minor: Inappropriate type
HDF5-DIAG: Error detected in HDF5 (1.12.2) MPI-process 0:
#000: H5O.c line 1245 in H5Oclose(): not a valid file object ID (dataset, group, or datatype)
major: Invalid arguments to routine
minor: Unable to release object
Traceback (most recent call last):
File "/home/remi/Bureau/Python/FEniCS/Mesh/convert.py", line 30, in <module>
h5file.read(mesh, "mesh", False)
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 convert string to cell type.
*** Reason: Unknown cell type ("").
*** Where: This error was encountered inside CellType.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: 3ea2183fbfe7277de9f16cbe1a9ffaab133ba1fa
*** -------------------------------------------------------------------------
Here is the script that generates the mesh:
############################################
#------------ Mesh Creation ---------------#
############################################
import gmsh
import time
import sys
import numpy as np
import math
# Dimension of the geometry
gdim = 2
# Markers (for physical groups, used in pinball_mesh.py)
wall_marker, inlet_marker, outlet_marker = range(1, 4)
obstacle_markers = list(range(4, 7))
# radius of the 3 cylinders
r = 0.5
# coordinates of the centers of the three cylinders (front, top and bottom)
c_x, c_y = [-3*math.sqrt(3)*r/2, 0, 0], [0, 3*r/2,-3*r/2]
# Calculating the time that mesh generation will take
tic = time.perf_counter()
gmsh.initialize()
gdim = 2
gmsh.model.add("fluidic_pinball")
# target mesh size for any point
# lc = 1
lc = 1
############################################
#---------- Geometry creation -------------#
############################################
## 0 dim : Points definition
# Rectangle delimiting the flow
gmsh.model.geo.addPoint(-12*r, -12*r, 0, lc, 1)
gmsh.model.geo.addPoint(40*r, -12*r, 0, lc, 2)
gmsh.model.geo.addPoint(40*r, 12*r, 0, lc, 3)
gmsh.model.geo.addPoint(-12*r, 12*r, 0, lc, 4)
# Points to create the circles for the three cylinders (4 points for each cylinder)
for i in range(3):
gmsh.model.geo.addPoint(c_x[i], c_y[i], 0, lc, 5 + 5*i) # center
gmsh.model.geo.addPoint(c_x[i] + r, c_y[i], 0, lc, 6 + 5*i)
gmsh.model.geo.addPoint(c_x[i], c_y[i] + r, 0, lc, 7 + 5*i)
gmsh.model.geo.addPoint(c_x[i] - r, c_y[i], 0, lc, 8 + 5*i)
gmsh.model.geo.addPoint(c_x[i], c_y[i] - r, 0, lc, 9 + 5*i)
# Points to create the polygon (pentagon in this case) that will include the cylinders and the wake
# In that pentagon, the mesh will be refined
gmsh.model.geo.addPoint(c_x[0] - 2*r, c_y[0], 0, lc, 20)
gmsh.model.geo.addPoint(c_x[1] - 2*r, c_y[1] + 2*r, 0, lc, 21)
gmsh.model.geo.addPoint(40*r, c_y[1] + 2*r, 0, lc, 22)
gmsh.model.geo.addPoint(40*r, c_y[2] - 2*r, 0, lc, 23)
gmsh.model.geo.addPoint(c_x[2] - 2*r, c_y[2] - 2*r, 0, lc, 24)
# level line 1
gmsh.model.geo.addPoint(3, c_y[1]+2*r, 0, lc, 25)
gmsh.model.geo.addPoint(3, c_y[2]-2*r, 0, lc, 26)
# level line 2
gmsh.model.geo.addPoint(10, c_y[1]+2*r, 0, lc, 27)
gmsh.model.geo.addPoint(10, c_y[2]-2*r, 0, lc, 28)
## 1 dim: Lines and Curves
# Lines for the rectangle
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 23, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
gmsh.model.geo.addLine(22 , 3, 5)
# Lines for the pentagone (level 1)
gmsh.model.geo.addLine(20, 21, 6)
gmsh.model.geo.addLine(21, 25, 7)
gmsh.model.geo.addLine(25, 26, 8)
gmsh.model.geo.addLine(26, 24, 9)
gmsh.model.geo.addLine(24, 20, 10)
# Lines for the wake rectangles
gmsh.model.geo.addLine(25, 27, 11)
gmsh.model.geo.addLine(27, 28, 12)
gmsh.model.geo.addLine(28, 26, 13)
gmsh.model.geo.addLine(27, 22, 14)
gmsh.model.geo.addLine(22, 23, 15)
gmsh.model.geo.addLine(23, 28, 16)
# Creation of the circles for the three cylinders (4 CircleArcs for each circle)
for i in range(3):
gmsh.model.geo.addCircleArc(6 + 5*i, 5 + 5*i, 7 + 5*i, tag=17 + 4*i)
gmsh.model.geo.addCircleArc(7 + 5*i, 5 + 5*i, 8 + 5*i, tag=18 + 4*i)
gmsh.model.geo.addCircleArc(8 + 5*i, 5 + 5*i, 9 + 5*i, tag=19 + 4*i)
gmsh.model.geo.addCircleArc(9 + 5*i, 5 + 5*i, 6 + 5*i, tag=20 + 4*i)
############################################
#-------------- Plane surfaces ------------#
############################################
## 2 dim
# Creation of the plane surface: rectangle - wake
gmsh.model.geo.addCurveLoop([1, 2, 16, 13, 9, 10, 6, 7, 11, 14, 5, 3, 4], 1)
gmsh.model.geo.addPlaneSurface([1], 1)
# Creation of the plane surface: pentagon - cylinders
gmsh.model.geo.addCurveLoop(list(range(6, 11)), 2)
# Creation of a curve for the polygon
curved_loop_list = [2]
# Creation of curves loops for the three cylinders
for i in range(3):
gmsh.model.geo.addCurveLoop([17 + 4*i, 18 + 4*i, 19 + 4*i, 20 + 4*i], 3 + i)
curved_loop_list.append(3 + i)
gmsh.model.geo.addPlaneSurface(curved_loop_list, 2)
# second level
gmsh.model.geo.addCurveLoop([11, 12, 13, -8], 6)
gmsh.model.geo.addPlaneSurface([6], 3)
# third level
gmsh.model.geo.addCurveLoop([14, 15, 16, -12], 7)
gmsh.model.geo.addPlaneSurface([7], 4)
gmsh.model.geo.synchronize()
############################################
#---------------- Fields ------------------#
############################################
## Fields
# Distance field to calculate distance from the three cylinders
gmsh.model.mesh.field.add("Distance", 1)
gmsh.model.mesh.field.setNumbers(1, "CurvesList", list(range(17, 31)))
gmsh.model.mesh.field.setNumber(1, "Sampling", 100)
# Threshold field to refine the mesh around the three cylinders
gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "InField", 1)
gmsh.model.mesh.field.setNumber(2, "SizeMin", lc / 20)
gmsh.model.mesh.field.setNumber(2, "SizeMax", lc)
gmsh.model.mesh.field.setNumber(2, "DistMin", r)
gmsh.model.mesh.field.setNumber(2, "DistMax", 2*r)
# Constant mesh size inside polygon
gmsh.model.mesh.field.add("Constant", 3)
gmsh.model.mesh.field.setNumber(3, "VIn", lc/12)
gmsh.model.mesh.field.setNumber(3, "VOut", lc)
gmsh.model.mesh.field.setNumbers(3, "SurfacesList", [2])
# Constant mesh size inside polygon
gmsh.model.mesh.field.add("Constant", 4)
gmsh.model.mesh.field.setNumber(4, "VIn", lc/7)
gmsh.model.mesh.field.setNumber(4, "VOut", lc)
gmsh.model.mesh.field.setNumbers(4, "SurfacesList", [3])
# Constant mesh size inside polygon
gmsh.model.mesh.field.add("Constant", 5)
gmsh.model.mesh.field.setNumber(5, "VIn", lc/4)
gmsh.model.mesh.field.setNumber(5, "VOut", lc)
gmsh.model.mesh.field.setNumbers(5, "SurfacesList", [4])
# # Distance field to calculate distance from the outlet
# gmsh.model.mesh.field.add("Distance", 4)
# gmsh.model.mesh.field.setNumbers(4, "CurvesList", [15, 19, 18, 22])
# gmsh.model.mesh.field.setNumber(4, "Sampling", 100)
# # Threshold field to refine the mesh from cylinders to outlet
# gmsh.model.mesh.field.add("Threshold", 5)
# gmsh.model.mesh.field.setNumber(5, "InField", 4)
# gmsh.model.mesh.field.setNumber(5, "SizeMin", lc/20)
# gmsh.model.mesh.field.setNumber(5, "SizeMax", lc)
# gmsh.model.mesh.field.setNumber(5, "DistMin", r)
# gmsh.model.mesh.field.setNumber(5, "DistMax", 20*r)
# Final field as in t10.py, take the minimum mesh size for all Threshold fields.
gmsh.model.mesh.field.add("Min", 6)
gmsh.model.mesh.field.setNumbers(6, "FieldsList", [2, 3, 4, 5])
gmsh.model.mesh.field.setAsBackgroundMesh(6)
############################################
#------------ Physical Groups -------------#
############################################
## Create Physical Groups
gmsh.model.addPhysicalGroup(1, [1, 3], name="walls") # walls: up and down sides of the rectangle
gmsh.model.addPhysicalGroup(1, [4], name="inlet") # inlet: left side of the rectangle
gmsh.model.addPhysicalGroup(1, [2, 15, 5], name="outlet") # outlet: right side of the rectangle
gmsh.model.addPhysicalGroup(1, [17, 18, 19, 20], name="front_cylinder") # front obstacle (left most)
gmsh.model.addPhysicalGroup(1, [21, 22, 23, 24], name="top_cylinder") # top obstacle (the up most)
gmsh.model.addPhysicalGroup(1, [25, 26, 27, 28], name="bottom_cylinder") # bottom cylinder (the down most)
gmsh.model.addPhysicalGroup(2, [1, 2, 3, 4], name="control_volume", tag=7) # Reunion of the two plane surfaces: where fluid will flow
gmsh.model.geo.synchronize()
############################################
#----- Mesh Saving & Post-processing ------#
############################################
# Generate and write to file
gmsh.model.mesh.generate(gdim)
gmsh.write("mesh_pinball.msh")
# gmsh.write("mesh_pinball.xml")
# gmsh.write("mesh_pinball.xdmf")
# gmsh.write("pinball_mesh.vtk")
# Launch the GUI to see the results:
if '-nopopup' not in sys.argv:
gmsh.fltk.run()
And here is the script that converts the msh file into h5 format, and try to read it in dolfin, thanks to HDF5File (I need to read the mesh with this class). I found it here :
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_pinball.msh')
line_mesh = create_mesh(msh, "line", data_name="Boundary", prune_z=True)
meshio.write("mesh_pinball_1D.xdmf", line_mesh)
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("mesh_pinball_2D.xdmf", triangle_mesh)
from dolfin import *
comm = MPI.comm_world
mesh = Mesh()
mesh_file = "mesh_pinball_2D.h5"
with HDF5File(comm, mesh_file, "r") as h5file:
h5file.read(mesh, "mesh", False)
facet = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
h5file.read(facet, "facet")
Actually, I don’t really know what is the second argument "mesh" of HDF5File.read, but with an ancient mesh (I only have the h5 file, but have no idea how it was generated), the script works fine (I have directly the h5 file, so no conversion is needed). But when I change "mesh" to any other string, it throws the same error.