hello,
I’m having a problem reading a 3d mesh created with GMSH using msh2xdmf.py (GitHub - floiseau/msh2xdmf: Converter from MSH mesh format to XDMF mesh format) or with Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #182 by gangadhara. Here’s the error I get: .
Traceback (most recent call last):
File "solver.py", line 87, in <module>
mesh, boundaries_mf, association_table = import_mesh(
^^^^^^^^^^^^
File "msh2xdmf.py", line 228, in import_mesh
infile.read(boundaries_mvc, 'boundaries')
File ".conda/lib/python3.12/site-packages/fenics_adjoint/types/io.py", line 30, in XDMFFile_read
output = __XDMFFile_read__(self, *args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError:
*** -------------------------------------------------------------------------
*** 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: 828b3965c3ebcd128da494c596c78b17c0a303da
*** -------------------------------------------------------------------------
I’ve created a very simple cube mesh to try debugging, but I get the same problem. In 2D it works fine, so I think I’m missing something. The msh2xmdf script produces cube_domain.xdmf and cube_boundaries.xdmf as well as cube_ass_table.ini which contains the Physical Surface defined in the .geo. In paraview all the files seems to be ok. I’ve attached the .geo and msh2xdmf scripts to make things easier.
// Définition des points du cube
Point(1) = {0, 0, 0, 1.0};
Point(2) = {1, 0, 0, 1.0};
Point(3) = {1, 1, 0, 1.0};
Point(4) = {0, 1, 0, 1.0};
Point(5) = {0, 0, 1, 1.0};
Point(6) = {1, 0, 1, 1.0};
Point(7) = {1, 1, 1, 1.0};
Point(8) = {0, 1, 1, 1.0};
// Définition des arêtes
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Line(9) = {1, 5};
Line(10) = {2, 6};
Line(11) = {3, 7};
Line(12) = {4, 8};
Curve Loop(1) = {1, 2, 3, 4}; // Face inférieure
Plane Surface(1) = {1};
Curve Loop(2) = {5, 6, 7, 8}; // Face supérieure
Plane Surface(2) = {2};
Curve Loop(3) = {1, 10, -5, -9}; // Face latérale avant
Plane Surface(3) = {3};
Curve Loop(4) = {2, 11, -6, -10}; // Face latérale droite
Plane Surface(4) = {4};
Curve Loop(5) = {3, 12, -7, -11}; // Face latérale arrière
Plane Surface(5) = {5};
Curve Loop(6) = {4, 9, -8, -12}; // Face latérale gauche
Plane Surface(6) = {6};
// Définition du volume
Surface Loop(1) = {1, 2, 3, 4, 5, 6};
Volume(1) = {1};
// Maillage transfinite
Transfinite Curve {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12} = 10 Using Progression 1.0;
Transfinite Surface {1};
Transfinite Surface {2};
Transfinite Surface {3};
Transfinite Surface {4};
Transfinite Surface {5};
Transfinite Surface {6};
Transfinite Volume {1};
//+
Physical Surface("inlet", 1) = {6};
//+
Physical Surface("wall", 2) = {5, 2, 3, 1};
//+
Physical Surface("outlet", 3) = {4};
//+
Physical Volume("fluid", 4) = {1};
import argparse
import meshio
import os
import mpi4py
import numpy as np
from configparser import ConfigParser
try:
from dolfin import XDMFFile, Mesh, MeshValueCollection
from dolfin.cpp.mesh import MeshFunctionSizet
except ImportError:
print("Could not import dolfin. Continuing without Dolfin support.")
def msh2xdmf(mesh_name, dim=2, directory="."):
"""
Function converting a MSH mesh into XDMF files.
The XDMF files are:
- "domain.xdmf": the domain;
- "boundaries.xdmf": the boundaries physical groups from GMSH;
"""
# Get the mesh name has prefix
prefix = mesh_name.split('.')[0]
# Read the input mesh
msh = meshio.read("{}/{}".format(directory, mesh_name))
# Generate the domain XDMF file
export_domain(msh, dim, directory, prefix)
# Generate the boundaries XDMF file
export_boundaries(msh, dim, directory, prefix)
# Export association table
export_association_table(msh, prefix, directory)
def export_domain(msh, dim, directory, prefix):
"""
Export the domain XDMF file as well as the subdomains values.
"""
# Set cell type
if dim == 2:
cell_type = "triangle"
elif dim == 3:
cell_type = "tetra"
# Generate the cell block for the domain cells
for cell in msh.cells:
print("cell.type: ",cell.type) # Vérifie les types de cellules dans ton mesh
data_array = [cell.data for cell in msh.cells if cell.type == cell_type]
if len(data_array) == 0:
print("WARNING: No domain physical group found.")
return
else:
data = np.concatenate(data_array)
cells = [
meshio.CellBlock(
cell_type=cell_type,
data=data,
)
]
# Generate the domain cells data (for the subdomains)
try:
cell_data = {
"subdomains": [
np.concatenate(
[
msh.cell_data["gmsh:physical"][i]
for i, cellBlock in enumerate(msh.cells)
if cellBlock.type == cell_type
]
)
]
}
except KeyError:
raise ValueError(
"""
No physical group found for the domain.
Define the domain physical group.
- if dim=2, the domain is a surface
- if dim=3, the domain is a volume
"""
)
# Generate a meshio Mesh for the domain
domain = meshio.Mesh(
points=msh.points[:, :dim],
cells=cells,
cell_data=cell_data,
)
# Export the XDMF mesh of the domain
meshio.write(
"{}/{}_{}".format(directory, prefix, "domain.xdmf"),
domain,
file_format="xdmf"
)
def export_boundaries(msh, dim, directory, prefix):
"""
Export the boundaries XDMF file.
"""
# Set the cell type
if dim == 2:
cell_type = "line"
elif dim == 3:
cell_type = "triangle"
# Generate the cell block for the boundaries cells
data_array = [cell.data for cell in msh.cells if cell.type == cell_type]
if len(data_array) == 0:
print("WARNING: No boundary physical group found.")
return
else:
data = np.concatenate(data_array)
boundaries_cells = [
meshio.CellBlock(
cell_type=cell_type,
data=data,
)
]
# Generate the boundaries cells data
cell_data = {
"boundaries": [
np.concatenate(
[
msh.cell_data["gmsh:physical"][i]
for i, cellBlock in enumerate(msh.cells)
if cellBlock.type == cell_type
]
)
]
}
# Generate the meshio Mesh for the boundaries physical groups
boundaries = meshio.Mesh(
points=msh.points[:, :dim],
cells=boundaries_cells,
cell_data=cell_data,
)
# Export the XDMF mesh of the lines boundaries
meshio.write(
"{}/{}_{}".format(directory, prefix, "boundaries.xdmf"),
boundaries,
file_format="xdmf"
)
def export_association_table(msh, prefix='mesh', directory='.', verbose=True):
"""
Display the association between the physical group label and the mesh
value.
"""
# Create association table
association_table = {}
# Display the correspondance
formatter = "|{:^20}|{:^20}|"
topbot = "+{:-^41}+".format("")
separator = "+{:-^20}+{:-^20}+".format("", "")
# Display
if verbose:
print('\n' + topbot)
print(formatter.format("GMSH label", "MeshFunction value"))
print(separator)
for label, arrays in msh.cell_sets.items():
# Get the index of the array in arrays
for i, array in enumerate(arrays):
if array.size != 0:
index = i
# Added check to make sure that the association table
# doesn't try to import irrelevant information.
if label != "gmsh:bounding_entities":
value = msh.cell_data["gmsh:physical"][index][0]
# Store the association table in a dictionnary
association_table[label] = value
# Display the association
if verbose:
print(formatter.format(label, value))
if verbose:
print(topbot)
# Export the association table
file_content = ConfigParser()
file_content["ASSOCIATION TABLE"] = association_table
file_name = "{}/{}_{}".format(directory, prefix, "association_table.ini")
with open(file_name, 'w') as f:
file_content.write(f)
def import_mesh(
prefix="mesh",
subdomains=False,
dim=2 ,
directory=".",
):
"""Function importing a dolfin mesh.
Arguments:
prefix (str, optional): mesh files prefix (eg. my_mesh.msh,
my_mesh_domain.xdmf, my_mesh_bondaries.xdmf). Defaults to "mesh".
subdomains (bool, optional): True if there are subdomains. Defaults to
False.
dim (int, optional): dimension of the domain. Defaults to 2.
directory (str, optional): directory of the mesh files. Defaults to ".".
Output:
- dolfin Mesh object containing the domain;
- dolfin MeshFunction object containing the physical lines (dim=2) or
surfaces (dim=3) defined in the msh file and the sub-domains;
- association table
"""
# Set the file name
domain = "{}_domain.xdmf".format(prefix)
boundaries = "{}_boundaries.xdmf".format(prefix)
# create 2 xdmf files if not converted before
if not os.path.exists("{}/{}".format(directory, domain)) or \
not os.path.exists("{}/{}".format(directory, boundaries)):
msh2xdmf("{}.msh".format(prefix), dim=dim, directory=directory)
# Import the converted domain
mesh = Mesh()
with XDMFFile("{}/{}".format(directory, domain)) as infile:
infile.read(mesh)
# Import the boundaries
boundaries_mvc = MeshValueCollection("size_t", mesh, dim=dim)
with XDMFFile("{}/{}".format(directory, boundaries)) as infile:
infile.read(boundaries_mvc, 'boundaries')
boundaries_mf = MeshFunctionSizet(mesh, boundaries_mvc)
# Import the subdomains
if subdomains:
subdomains_mvc = MeshValueCollection("size_t", mesh, dim=dim)
with XDMFFile("{}/{}".format(directory, domain)) as infile:
infile.read(subdomains_mvc, 'subdomains')
subdomains_mf = MeshFunctionSizet(mesh, subdomains_mvc)
# Import the association table
association_table_name = "{}/{}_{}".format(
directory, prefix, "association_table.ini")
file_content = ConfigParser()
file_content.read(association_table_name)
association_table = dict(file_content["ASSOCIATION TABLE"])
# Convert the value from string to int
for key, value in association_table.items():
association_table[key] = int(value)
# Return the Mesh and the MeshFunction objects
if not subdomains:
return mesh, boundaries_mf, association_table
else:
return mesh, boundaries_mf, subdomains_mf, association_table
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument(
"msh_file",
help="input .msh file",
type=str,
)
parser.add_argument(
"-d",
"--dimension",
help="dimension of the domain",
type=int,
default=2,
)
args = parser.parse_args()
# Get current directory
current_directory = os.getcwd()
# Conert the mesh
msh2xdmf(args.msh_file, args.dimension, directory=current_directory)
Thank you for yout help and your time !
