Skip to content

Medit higher order elements #1526

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
matnoel wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Medit higher order elements #1526

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5bdb753
#1511 Higher-order Medit elements
matnoel Mar 16, 2025
a48ed24
#1511 Higher-order Medit elements
matnoel Mar 26, 2025
c68ce77
#1511 Higher-order Medit elements
matnoel Apr 9, 2025
cef267f
#1517 Update gmsh/common.py
matnoel Apr 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,4 +23,4 @@ doc/_build/
.coverage
.tox/
foo.vtk
.vscode/
.vscode/
1 change: 1 addition & 0 deletions src/meshio/_mesh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
"quad9": 2,
"tetra10": 3,
"hexahedron27": 3,
"wedge15": 3,
"wedge18": 3,
"pyramid14": 3,
"vertex": 0,
Expand Down
71 changes: 29 additions & 42 deletions src/meshio/gmsh/common.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,57 +160,44 @@ def _read_data(f, tag, data_dict, data_size, is_ascii):
}
_meshio_to_gmsh_type = {v: k for k, v in _gmsh_to_meshio_type.items()}

_gmsh_to_meshio = {
# fmt: off
"tetra10": [0, 1, 2, 3, 4, 5, 6, 7, 9, 8],
"hexahedron20": [
0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 13,
9, 16, 18, 19, 17, 10, 12, 14, 15,
], # https://vtk.org/doc/release/4.2/html/classvtkQuadraticHexahedron.html and https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering
"hexahedron27": [
0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 13,
9, 16, 18, 19, 17, 10, 12, 14, 15,
22, 23, 21, 24, 20, 25, 26,
],
"wedge15": [
0, 1, 2, 3, 4, 5, 6, 9, 7, 12, 14, 13, 8, 10, 11
], # http://davis.lbl.gov/Manuals/VTK-4.5/classvtkQuadraticWedge.html and https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering
"wedge18": [
0, 1, 2, 3, 4, 5, 6, 9, 7, 12, 14, 13, 8, 10, 11, 15, 17, 16
],
"pyramid13": [0, 1, 2, 3, 4, 5, 8, 10, 6, 7, 9, 11, 12],
# fmt: on
}

_meshio_to_gmsh = {cell_type: [order.index(i) for i in range(len(order))]
for cell_type, order in _gmsh_to_meshio.items()}

def _gmsh_to_meshio_order(cell_type: str, idx: ArrayLike) -> np.ndarray:
# Gmsh cells are mostly ordered like VTK, with a few exceptions:
meshio_ordering = {
# fmt: off
"tetra10": [0, 1, 2, 3, 4, 5, 6, 7, 9, 8],
"hexahedron20": [
0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 13,
9, 16, 18, 19, 17, 10, 12, 14, 15,
], # https://vtk.org/doc/release/4.2/html/classvtkQuadraticHexahedron.html and https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering
"hexahedron27": [
0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 13,
9, 16, 18, 19, 17, 10, 12, 14, 15,
22, 23, 21, 24, 20, 25, 26,
],
"wedge15": [
0, 1, 2, 3, 4, 5, 6, 9, 7, 12, 14, 13, 8, 10, 11
], # http://davis.lbl.gov/Manuals/VTK-4.5/classvtkQuadraticWedge.html and https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering
"pyramid13": [0, 1, 2, 3, 4, 5, 8, 10, 6, 7, 9, 11, 12],
# fmt: on
}
idx = np.asarray(idx)
if cell_type not in meshio_ordering:
if cell_type not in _gmsh_to_meshio:
return idx
return idx[:, meshio_ordering[cell_type]]

return idx[:, _gmsh_to_meshio[cell_type]]

def _meshio_to_gmsh_order(cell_type: str, idx: ArrayLike) -> np.ndarray:
# Gmsh cells are mostly ordered like VTK, with a few exceptions:
gmsh_ordering = {
# fmt: off
"tetra10": [0, 1, 2, 3, 4, 5, 6, 7, 9, 8],
"hexahedron20": [
0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 16,
9, 17, 10, 18, 19, 12, 15, 13, 14,
],
"hexahedron27": [
0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 16,
9, 17, 10, 18, 19, 12, 15, 13, 14,
24, 22, 20, 21, 23, 25, 26,
],
"wedge15": [
0, 1, 2, 3, 4, 5, 6, 8, 12, 7, 13, 14, 9, 11, 10,
],
"pyramid13": [0, 1, 2, 3, 4, 5, 8, 9, 6, 10, 7, 11, 12],
# fmt: on
}
# Gmsh cells are mostly ordered like VTK, with a few exceptions:
idx = np.asarray(idx)
if cell_type not in gmsh_ordering:
if cell_type not in _meshio_to_gmsh:
return idx
return idx[:, gmsh_ordering[cell_type]]
return idx[:, _meshio_to_gmsh[cell_type]]


def _write_physical_names(fh, field_data):
Expand Down
153 changes: 98 additions & 55 deletions src/meshio/medit/_medit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,48 @@
from .._mesh import Mesh
from ._medit_internal import medit_codes

DICT_MESHIO: dict[str, tuple[str, int, int]] = {
# "vertex": ("Corners", 0, 13), # DOT NOT UNCOMMENT
# wrong data for GmfCorners in medit_codes?
"point": ("Vertices", 0, 4),
"line": ("Edges", 2, 5),
"line3": ("EdgesP2", 3, 25),
"line4": ("EdgesP3", 4, 92),
"line5": ("EdgesP4", 5, 93),
"triangle": ("Triangles", 3, 6),
"triangle6": ("TrianglesP2", 6, 24),
"triangle10": ("TrianglesP3", 10, 90),
"triangle15": ("TrianglesP4", 15, 91),
"quad": ("Quadrilaterals", 4, 7),
# "quad8": (_, 8, _), # not on medit!
"quad9": ("QuadrilateralsQ2", 9, 27),
"tetra": ("Tetrahedra", 4, 8),
"tetra10": ("TetrahedraP2", 10, 30),
"wedge": ("Prisms", 6, 9),
# "wedge15": (_, 15, _), # not on medit!
"wedge18": ("PrismsP2", 18, 86),
"pyramid": ("Pyramids", 5, 49),
"hexahedron": ("Hexahedra", 8, 10), # Frey
# "hexahedron20": (_, 20, _), # not on medit!
"hexahedron27": ("HexahedraQ2", 27, 33)
}
"""meshio: (medit, nodes_per_elements, tag)"""
# see _medit_internal.py for tags

DICT_MEDIT: dict[str, tuple[str, int, int]] = {}
"""medit: (meshio, nodes_per_elements, tag)"""

DICT_GMFMEDIT: dict[str, tuple[str, int, int]] = {}
"""GmfMedit: (meshio, nodes_per_elements, tag)"""

for key, value in DICT_MESHIO.items():
medit, nPe, tag = value
DICT_MEDIT[medit] = (key, nPe, tag)
DICT_GMFMEDIT["Gmf"+medit] = (key, nPe, tag)

# add Dobrzynski case
DICT_MEDIT["Hexaedra"] = ("hexahedron", 8, 10)
DICT_GMFMEDIT["GmfHexaedra"] = ("hexahedron", 8, 10)

def read(filename):
with open_file(filename) as f:
Expand All @@ -25,7 +67,6 @@ def read(filename):
mesh = read_ascii_buffer(f)
return mesh


def _produce_dtype(string_type, dim, itype, ftype):
"""
convert a medit_code to a dtype appropriate for building a numpy array
Expand All @@ -49,19 +90,33 @@ def _produce_dtype(string_type, dim, itype, ftype):
res += ","
return res

_medit_to_meshio = {
"hexahedron27": [
*list(range(20)), 25, 23, 22, 24, 20, 21, 26
]
}

def read_binary_buffer(f):
_meshio_to_medit = {cell_type: [order.index(i) for i in range(len(order))]
for cell_type, order in _medit_to_meshio.items()}

def _convert_cell_data(cell_type: str, data, dict_convert: dict[str, list[int]]) -> np.ndarray:

if cell_type in dict_convert.keys():
idx = dict_convert[cell_type]
data = data[:, idx]

return data

meshio_from_medit = {
"GmfVertices": ("point", None),
"GmfEdges": ("line", 2),
"GmfTriangles": ("triangle", 3),
"GmfQuadrilaterals": ("quad", 4),
"GmfTetrahedra": ("tetra", 4),
"GmfPrisms": ("wedge", 6),
"GmfPyramids": ("pyramid", 5),
"GmfHexahedra": ("hexahedron", 8),
}
def _convert_cells(cells: list[tuple[str, np.ndarray]], dict_convert: dict[str, list[int]]):

for i, cell in enumerate(cells):
cell_type, data = cell
data = _convert_cell_data(cell_type, data, dict_convert)
cells[i] = (cell_type, data)

return cells

def read_binary_buffer(f):

dim = 0
points = None
Expand Down Expand Up @@ -150,20 +205,22 @@ def read_binary_buffer(f):
field_template = field_code[2]
dtype = np.dtype(_produce_dtype(field_template, dim, itype, ftype))
out = np.asarray(np.fromfile(f, count=nitems, dtype=dtype))
if field_code[0] not in meshio_from_medit.keys():
if field_code[0] not in DICT_GMFMEDIT.keys():
warn(f"meshio doesn't know {field_code[0]} type. Skipping.")
continue
continue

elif field_code[0] == "GmfVertices":
points = out["f0"]
point_data["medit:ref"] = out["f1"]
else:
meshio_type, ncols = meshio_from_medit[field_code[0]]
meshio_type, ncols, _ = DICT_GMFMEDIT[field_code[0]]
# transform the structured array to integer array which suffices
# for the cell connectivity
out_view = out.view(itype).reshape(nitems, ncols + 1)
cells.append((meshio_type, out_view[:, :ncols] - 1))
cell_data["medit:ref"].append(out_view[:, -1])

cells = _convert_cells(cells, _medit_to_meshio)

return Mesh(points, cells, point_data=point_data, cell_data=cell_data)

Expand All @@ -173,17 +230,7 @@ def read_ascii_buffer(f):
cells = []
point_data = {}
cell_data = {"medit:ref": []}

meshio_from_medit = {
"Edges": ("line", 2),
"Triangles": ("triangle", 3),
"Quadrilaterals": ("quad", 4),
"Tetrahedra": ("tetra", 4),
"Prisms": ("wedge", 6),
"Pyramids": ("pyramid", 5),
"Hexahedra": ("hexahedron", 8), # Frey
"Hexaedra": ("hexahedron", 8), # Dobrzynski
}

points = None
dtype = None

Expand All @@ -199,12 +246,12 @@ def read_ascii_buffer(f):

items = line.split()

if not items[0].isalpha():
if (not items[0].isalpha()) and (items[0] not in DICT_MEDIT.keys()):
raise ReadError()

if items[0] == "MeshVersionFormatted":
version = items[1]
dtype = {"0": c_float, "1": c_float, "2": c_double}[version]
dtype = {"0": c_float, "1": c_float, "2": c_double, "3": c_double}[version]
elif items[0] == "Dimension":
if len(items) >= 2:
dim = int(items[1])
Expand All @@ -223,8 +270,8 @@ def read_ascii_buffer(f):
).reshape(num_verts, dim + 1)
points = out[:, :dim]
point_data["medit:ref"] = out[:, dim].astype(int)
elif items[0] in meshio_from_medit:
meshio_type, points_per_cell = meshio_from_medit[items[0]]
elif items[0] in DICT_MEDIT:
meshio_type, points_per_cell, _ = DICT_MEDIT[items[0]]
# The first value is the number of elements
num_cells = int(f.readline())

Expand Down Expand Up @@ -257,6 +304,12 @@ def read_ascii_buffer(f):
np.fromfile(
f, count=num_sub_domain_from_mesh * 4, dtype=int, sep=" "
).reshape(num_sub_domain_from_mesh, 4)
elif items[0] == "SubDomainFromGeom":
# those are just discarded
num_sub_domain_from_geom = int(f.readline())
np.fromfile(
f, count=num_sub_domain_from_geom * 4, dtype=int, sep=" "
).reshape(num_sub_domain_from_geom, 4)
elif items[0] == "VertexOnGeometricVertex":
# those are just discarded
num_vertex_on_geometric_vertex = int(f.readline())
Expand Down Expand Up @@ -289,11 +342,14 @@ def read_ascii_buffer(f):
for _ in range(num_to_pass):
f.readline()
else:
if items[0] != "End":
if items[0] not in ("End", "END"):
raise ReadError(f"Unknown keyword '{items[0]}'.")

if points is None:
raise ReadError("Expected `Vertices`")

cells = _convert_cells(cells, _medit_to_meshio)

return Mesh(points, cells, point_data=point_data, cell_data=cell_data)


Expand All @@ -303,7 +359,6 @@ def write(filename, mesh, float_fmt=".16e"):
else:
write_ascii_file(filename, mesh, float_fmt)


def write_ascii_file(filename, mesh, float_fmt=".16e"):
with open_file(filename, "wb") as fh:
version = {np.dtype(c_float): 1, np.dtype(c_double): 2}[mesh.points.dtype]
Expand Down Expand Up @@ -331,17 +386,7 @@ def write_ascii_file(filename, mesh, float_fmt=".16e"):
fmt = " ".join(["{:" + float_fmt + "}"] * d) + " {:d}\n"
for x, label in zip(mesh.points, labels):
fh.write(fmt.format(*x, label).encode())

medit_from_meshio = {
"line": ("Edges", 2),
"triangle": ("Triangles", 3),
"quad": ("Quadrilaterals", 4),
"tetra": ("Tetrahedra", 4),
"wedge": ("Prisms", 6),
"pyramid": ("Pyramids", 5),
"hexahedron": ("Hexahedra", 8),
}


# pick out cell_data
labels_key, other = _pick_first_int_data(mesh.cell_data)
if labels_key and other:
Expand All @@ -352,10 +397,13 @@ def write_ascii_file(filename, mesh, float_fmt=".16e"):
)

for k, cell_block in enumerate(mesh.cells):

cell_type = cell_block.type
data = cell_block.data
data = _convert_cell_data(cell_type, data, _meshio_to_medit)

try:
medit_name, num = medit_from_meshio[cell_type]
medit_name, num, _ = DICT_MESHIO[cell_type]
except KeyError:
msg = f"MEDIT's mesh format doesn't know {cell_type} cells. Skipping."
warn(msg)
Expand Down Expand Up @@ -464,19 +512,14 @@ def write_binary_file(f, mesh):
f"Picking {labels_key}, skipping {string}."
)

# first component is medit keyword id see _medit_internal.py
medit_from_meshio = {
"line": 5,
"triangle": 6,
"quad": 7,
"tetra": 8,
"wedge": 9,
"pyramid": 49,
"hexahedron": 10,
}
for k, cell_block in enumerate(mesh.cells):

# reorder data
cell_block.data = _convert_cell_data(
cell_block.type, cell_block.data, _meshio_to_medit)

try:
medit_key = medit_from_meshio[cell_block.type]
_, _, medit_key = DICT_MESHIO[cell_block.type]
except KeyError:
warn(
f"MEDIT's mesh format doesn't know {cell_block.type} cells. "
Expand Down