diff --git a/.gitignore b/.gitignore index 513e3e81..99341663 100644 --- a/.gitignore +++ b/.gitignore @@ -23,4 +23,4 @@ doc/_build/ .coverage .tox/ foo.vtk -.vscode/ +.vscode/ \ No newline at end of file diff --git a/src/meshio/_mesh.py b/src/meshio/_mesh.py index 741831d9..2bc21495 100644 --- a/src/meshio/_mesh.py +++ b/src/meshio/_mesh.py @@ -21,6 +21,7 @@ "quad9": 2, "tetra10": 3, "hexahedron27": 3, + "wedge15": 3, "wedge18": 3, "pyramid14": 3, "vertex": 0, diff --git a/src/meshio/gmsh/common.py b/src/meshio/gmsh/common.py index 717e82ee..0456eba9 100644 --- a/src/meshio/gmsh/common.py +++ b/src/meshio/gmsh/common.py @@ -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): diff --git a/src/meshio/medit/_medit.py b/src/meshio/medit/_medit.py index 272bd16e..5ba754f8 100644 --- a/src/meshio/medit/_medit.py +++ b/src/meshio/medit/_medit.py @@ -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: @@ -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 @@ -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 @@ -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) @@ -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 @@ -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]) @@ -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()) @@ -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()) @@ -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) @@ -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] @@ -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: @@ -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) @@ -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. "