Hi guys!
Some time ago I managed to write some code that can read MED files into FEniCS. Now I am refactoring that code to work with FEniCSx, but I am running into issues. This is what I have so far:
import typing
import pathlib
import meshio
import enum
import numpy as np
import tempfile
import mpi4py
import dolfinx
def load_raw(path: typing.Union[pathlib.Path, str]) -> meshio.Mesh:
return meshio.read(path)
CellTagsT = typing.Dict[int, typing.List[str]]
def calculate_tags_offset(cell_tags: CellTagsT) -> int:
tags_offset = 0
min_tag = min(k for k in cell_tags.keys())
if min_tag < 0:
tags_offset = -min_tag
return tags_offset
def get_cell_tags(mesh: meshio.Mesh) -> typing.Tuple[CellTagsT]:
cell_tags = mesh.cell_tags
tags_offset = calculate_tags_offset(cell_tags=cell_tags)
return {
k + tags_offset: cell_tags[k]
for k in cell_tags
}
class CellType(enum.Enum):
HEXAHEDRON = 'hexahedron'
INTERVAL = 'line'
PRISM = 'prism'
PYRAMID = 'pyramid'
QUADRILATERAL = 'quad'
TETRAHEDRON = 'tetra'
TRIANGLE = 'triangle'
def convert_entity(
comm: mpi4py.MPI,
mesh: meshio.Mesh,
cell_type: CellType,
tag_key: str,
prune_mask: np.ndarray,
coordinate_scaling: float
) -> typing.Tuple[dolfinx.mesh.Mesh, dolfinx.mesh.MeshTagsMetaClass]:
cells = mesh.get_cells_type(cell_type=cell_type.value)
cell_data = mesh.get_cell_data(name=tag_key, cell_type=cell_type.value)
points = mesh.points[:, prune_mask] * coordinate_scaling
tags_offset = calculate_tags_offset(cell_tags=mesh.cell_tags)
entity_mesh = meshio.Mesh(
points=points,
cells={cell_type.value: cells},
cell_data={'marker': [cell_data + tags_offset]}
)
with tempfile.TemporaryDirectory() as xdmf_dir:
target_xdmf_dir = pathlib.Path(xdmf_dir)
tmp_path = target_xdmf_dir.joinpath(cell_type.value + '_mesh.xdmf')
# entity_mesh.write(tmp_path)
meshio.write(tmp_path, entity_mesh)
with dolfinx.io.XDMFFile(comm=comm, filename=tmp_path, file_mode='r') as f:
out_mesh = f.read_mesh(name='Grid')
out_tags = f.read_meshtags(out_mesh, name='Grid')
out_mesh.topology.create_connectivity(
*[mesh.topology.dim - n for n in range(mesh.topology.dim)]
)
return out_mesh, out_tags
convert_entity(
comm=mpi4py.MPI.COMM_WORLD,
mesh=load_raw('/home/admin/test.med'),
cell_type=CellType.TETRAHEDRON,
tag_key='cell_tags',
prune_mask=np.zeros(3, dtype=bool),
coordinate_scaling=1.0
)
A bit long, but I guess the only important bit is the convert_entity
function. Basically, it is all based off:
However, when I use dolfinx.io.XDMFFile
’s read_mesh
method I get this error:
Traceback (most recent call last):
File "/usr/lib/python3.10/runpy.py", line 196, in _run_module_as_main
return _run_code(code, main_globals, None,
File "/usr/lib/python3.10/runpy.py", line 86, in _run_code
exec(code, run_globals)
File "/home/admin/.vscode-server/extensions/ms-python.python-2023.4.1/pythonFiles/lib/python/debugpy/adapter/../../debugpy/launcher/../../debugpy/__main__.py", line 39, in <module>
cli.main()
File "/home/admin/.vscode-server/extensions/ms-python.python-2023.4.1/pythonFiles/lib/python/debugpy/adapter/../../debugpy/launcher/../../debugpy/../debugpy/server/cli.py", line 430, in main
run()
File "/home/admin/.vscode-server/extensions/ms-python.python-2023.4.1/pythonFiles/lib/python/debugpy/adapter/../../debugpy/launcher/../../debugpy/../debugpy/server/cli.py", line 284, in run_file
runpy.run_path(target, run_name="__main__")
File "/home/admin/.vscode-server/extensions/ms-python.python-2023.4.1/pythonFiles/lib/python/debugpy/_vendored/pydevd/_pydevd_bundle/pydevd_runpy.py", line 321, in run_path
return _run_module_code(code, init_globals, run_name,
File "/home/admin/.vscode-server/extensions/ms-python.python-2023.4.1/pythonFiles/lib/python/debugpy/_vendored/pydevd/_pydevd_bundle/pydevd_runpy.py", line 135, in _run_module_code
_run_code(code, mod_globals, init_globals,
File "/home/admin/.vscode-server/extensions/ms-python.python-2023.4.1/pythonFiles/lib/python/debugpy/_vendored/pydevd/_pydevd_bundle/pydevd_runpy.py", line 124, in _run_code
exec(code, run_globals)
File "/home/admin/projects/acoupy/acoupy_meshutil/acoupy_meshutil/util.py", line 81, in <module>
convert_entity(
File "/home/admin/projects/acoupy/acoupy_meshutil/acoupy_meshutil/util.py", line 72, in convert_entity
out_mesh = f.read_mesh(name='Grid')
File "/usr/lib/python3.10/site-packages/dolfinx/io/utils.py", line 163, in read_mesh
x = super().read_geometry_data(name, xpath)
RuntimeError: Cannot determine geometric dimension. GeometryType "" in XDMF file is unknown or unsupported
If I hack the XDMF file and I change it to <Geometry GeometryType="XYZ">
things get even worse?
python: /build/dolfinx/src/dolfinx-0.6.0/cpp/dolfinx/io/xdmf_mesh.cpp:264: std::pair<std::vector<double>, std::array<long unsigned int, 2> > dolfinx::io::xdmf_mesh::read_geometry_data(MPI_Comm, hid_t, const pugi::xml_node&): Assertion `gdims[1] == (int)gdim' failed.
Loguru caught a signal: SIGABRT
Stack trace:
53 0x5606f6449045 _start + 37
52 0x7fcf143ee84a __libc_start_main + 138
51 0x7fcf143ee790 /usr/lib/libc.so.6(+0x23790) [0x7fcf143ee790]
50 0x7fcf141e336b Py_BytesMain + 59
49 0x7fcf14096057 /usr/lib/libpython3.10.so.1.0(+0x96057) [0x7fcf14096057]
48 0x7fcf14216617 /usr/lib/libpython3.10.so.1.0(+0x216617) [0x7fcf14216617]
47 0x7fcf141562d9 _PyFunction_Vectorcall + 121
46 0x7fcf141454d6 _PyEval_EvalFrameDefault + 838
45 0x7fcf141562d9 _PyFunction_Vectorcall + 121
44 0x7fcf141454d6 _PyEval_EvalFrameDefault + 838
43 0x7fcf141564cf /usr/lib/libpython3.10.so.1.0(+0x1564cf) [0x7fcf141564cf]
42 0x7fcf141f7f1b /usr/lib/libpython3.10.so.1.0(+0x1f7f1b) [0x7fcf141f7f1b]
41 0x7fcf141f1b94 PyEval_EvalCode + 148
40 0x7fcf14144120 /usr/lib/libpython3.10.so.1.0(+0x144120) [0x7fcf14144120]
39 0x7fcf1414a03f _PyEval_EvalFrameDefault + 20143
38 0x7fcf141562d9 _PyFunction_Vectorcall + 121
37 0x7fcf141454d6 _PyEval_EvalFrameDefault + 838
36 0x7fcf141562d9 _PyFunction_Vectorcall + 121
35 0x7fcf1414b047 _PyEval_EvalFrameDefault + 24247
34 0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
33 0x7fcf141562d9 _PyFunction_Vectorcall + 121
32 0x7fcf1414b047 _PyEval_EvalFrameDefault + 24247
31 0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
30 0x7fcf141562d9 _PyFunction_Vectorcall + 121
29 0x7fcf1414c696 _PyEval_EvalFrameDefault + 29958
28 0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
27 0x7fcf141562d9 _PyFunction_Vectorcall + 121
26 0x7fcf1414c696 _PyEval_EvalFrameDefault + 29958
25 0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
24 0x7fcf141564cf /usr/lib/libpython3.10.so.1.0(+0x1564cf) [0x7fcf141564cf]
23 0x7fcf141f7f1b /usr/lib/libpython3.10.so.1.0(+0x1f7f1b) [0x7fcf141f7f1b]
22 0x7fcf141f1b94 PyEval_EvalCode + 148
21 0x7fcf14144120 /usr/lib/libpython3.10.so.1.0(+0x144120) [0x7fcf14144120]
20 0x7fcf1414b047 _PyEval_EvalFrameDefault + 24247
19 0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
18 0x7fcf141562d9 _PyFunction_Vectorcall + 121
17 0x7fcf1414b047 _PyEval_EvalFrameDefault + 24247
16 0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
15 0x7fcf14161446 /usr/lib/libpython3.10.so.1.0(+0x161446) [0x7fcf14161446]
14 0x7fcf1414b3e5 _PyEval_EvalFrameDefault + 25173
13 0x7fcf141f0232 /usr/lib/libpython3.10.so.1.0(+0x1f0232) [0x7fcf141f0232]
12 0x7fcf141616a6 /usr/lib/libpython3.10.so.1.0(+0x1616a6) [0x7fcf141616a6]
11 0x7fcf1414f4eb _PyObject_MakeTpCall + 683
10 0x7fcf14155e21 /usr/lib/libpython3.10.so.1.0(+0x155e21) [0x7fcf14155e21]
9 0x7fcf01674129 /usr/lib/python3.10/site-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x74129) [0x7fcf01674129]
8 0x7fcf017984fb /usr/lib/python3.10/site-packages/dolfinx/cpp.cpython-310-x86_64-linux-gnu.so(+0x1984fb) [0x7fcf017984fb]
7 0x7fcf01537210 dolfinx::io::XDMFFile::read_geometry_data(std::string, std::string) const + 624
6 0x7fcf0154d62e dolfinx::io::xdmf_mesh::read_geometry_data(ompi_communicator_t*, long, pugi::xml_node const&) + 5598
5 0x7fcf143fc9f6 /usr/lib/libc.so.6(+0x319f6) [0x7fcf143fc9f6]
4 0x7fcf143ed45c /usr/lib/libc.so.6(+0x2245c) [0x7fcf143ed45c]
3 0x7fcf143ed53d abort + 215
2 0x7fcf14403ea8 gsignal + 24
1 0x7fcf144528ec /usr/lib/libc.so.6(+0x878ec) [0x7fcf144528ec]
0 0x7fcf14403f50 /usr/lib/libc.so.6(+0x38f50) [0x7fcf14403f50]
2023-03-11 16:07:34.629 ( 44.973s) [main ] :0 FATL| Signal: SIGABRT
It makes me think the mesh is the problem but it is similar to others I successfully loaded into FEnICS. The mesh is here.
Any idea?