From 11db1346e05059d59438da720a049ff50b2c3826 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 2 May 2024 10:35:41 -0500 Subject: [PATCH 1/2] mesh.io.from_vertices_and_simplices: use new force_positive_orientation --- meshmode/mesh/io.py | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 84c0c863..2dbe618f 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -399,18 +399,10 @@ def from_vertices_and_simplices( grp = make_group_from_vertices(vertices, simplices, order) - if fix_orientation: - if grp.dim != vertices.shape[0]: - raise ValueError("can only fix orientation of volume meshes") - - from meshmode.mesh.processing import ( - find_volume_mesh_element_group_orientation, flip_element_group) - orient = find_volume_mesh_element_group_orientation(vertices, grp) - grp = flip_element_group(vertices, grp, orient < 0) - return make_mesh( vertices=vertices, groups=[grp], - is_conforming=True) + is_conforming=True, + force_positive_orientation=fix_orientation) # }}} From 51ccec9c3bc32ea983a3d9bba9b65d208388c6b9 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 2 May 2024 10:36:42 -0500 Subject: [PATCH 2/2] Add a meshio reader Closes gh-405 --- .test-conda-env-py3.yml | 3 + meshmode/mesh/io.py | 267 +++++++++++++++++++++++++++++++++++++++- setup.cfg | 3 + setup.py | 1 + test/coarse.bdf | 204 ++++++++++++++++++++++++++++++ test/test_mesh.py | 5 + 6 files changed, 481 insertions(+), 2 deletions(-) create mode 100644 test/coarse.bdf diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index f06be3c1..dfd861c3 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -31,6 +31,9 @@ dependencies: # for xdmf/hdf5 visualizer - h5py=*=mpi_openmpi* +# To test meshio reader +- meshio + # Only needed to make pylint succeed - matplotlib-base diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 2dbe618f..ac2e4d6f 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -1,4 +1,7 @@ -__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" +__copyright__ = """ +Copyright (C) 2014 University of Illinois Board of Trustees +Copyright (C) 2020 Miche Jansen (for the meshio cell node locations TikZ) +""" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -20,13 +23,20 @@ THE SOFTWARE. """ +from dataclasses import dataclass +from os import PathLike +from typing import BinaryIO, List, Optional, TextIO, Type, Union + import numpy as np from gmsh_interop.reader import ( # noqa: F401 FileSource, GmshMeshReceiverBase, GmshSimplexElementBase, GmshTensorProductElementBase, LiteralSource, ScriptSource, ScriptWithFilesSource) -from meshmode.mesh import Mesh +from meshmode.mesh import ( + Mesh, MeshElementGroup, SimplexElementGroup, TensorProductElementGroup, + _ModepyElementGroup) +from meshmode.mesh.tools import find_point_permutation __doc__ = """ @@ -40,6 +50,7 @@ .. autofunction:: from_meshpy .. autofunction:: from_vertices_and_simplices .. autofunction:: to_json +.. autofunction:: read_via_meshio """ @@ -455,4 +466,256 @@ def nodal_adjacency_to_json(mesh): # }}} +# {{{ read_via_meshio + +@dataclass(frozen=True) +class _MeshIOCellTypeInfo: + egrp_cls: Type[_ModepyElementGroup] + dim: int + order: int + # in bi-unit coordinates + unit_nodes: str + + +_MESHIO_CELL_TYPE_INFO = { + # https://github.com/inducer/meshio/blob/b2ee99842e119901349fdeee06b5bf61e01f450a/doc/cell_types.tex + "line": _MeshIOCellTypeInfo( + egrp_cls=SimplexElementGroup, + dim=1, + order=1, + unit_nodes=r""" + \cnode(n0,0,0,0,0,below right); + \cnode(n1,2,0,0,1,below right); + """, + ), + "line3": _MeshIOCellTypeInfo( + egrp_cls=SimplexElementGroup, + dim=1, + order=2, + unit_nodes=r""" + \cnode(n0,0,0,0,0,below right); + \cnode(n1,2,0,0,1,below right); + \cnode(n2,1,0,0,2,below right); + """, + ), + "triangle": _MeshIOCellTypeInfo( + egrp_cls=SimplexElementGroup, + dim=2, + order=1, + unit_nodes=r""" + \cnode(n0,0,0,0,0,below right); + \cnode(n1,2,0,0,1,below right); + \cnode(n2,0,2,0,2,right); + """, + ), + "triangle6": _MeshIOCellTypeInfo( + egrp_cls=SimplexElementGroup, + dim=2, + order=2, + unit_nodes=r""" + \cnode(n0,0,0,0,0,below right); + \cnode(n1,2,0,0,1,below right); + \cnode(n2,0,2,0,2,right); + \cnode(n3,1,0,0,3,below right); + \cnode(n4,1,1,0,4,right); + \cnode(n5,0,1,0,5,below right); + """, + ), + "quad": _MeshIOCellTypeInfo( + egrp_cls=TensorProductElementGroup, + dim=2, + order=1, + unit_nodes=r""" + \cnode(n0,0,0,0,0,below right); + \cnode(n1,2,0,0,1,below right); + \cnode(n2,2,2,0,2,below right); + \cnode(n3,0,2,0,3,below right); + """, + ), + "quad9": _MeshIOCellTypeInfo( + egrp_cls=TensorProductElementGroup, + dim=2, + order=2, + unit_nodes=r""" + \cnode(n0,0,0,0,0,below right); + \cnode(n1,2,0,0,1,below right); + \cnode(n2,2,2,0,2,below right); + \cnode(n3,0,2,0,3,below right); + \cnode(n4,1,0,0,4,below right); + \cnode(n5,2,1,0,5,below right); + \cnode(n6,1,2,0,6,below right); + \cnode(n7,0,1,0,7,below right); + \cnode(n8,1,1,0,8,below right); + """, + ), + "tetra": _MeshIOCellTypeInfo( + egrp_cls=SimplexElementGroup, + dim=3, + order=1, + unit_nodes=r""" + \cnode(n0,0,0,0,0,below right); + \cnode(n1,2,0,0,1,below right); + \cnode(n2,0,2,0,2,below right); + \cnode(n3,0,0,2,3,right); + """, + ), + "tetra10": _MeshIOCellTypeInfo( + egrp_cls=SimplexElementGroup, + dim=3, + order=2, + unit_nodes=r""" + \cnode(n0,0,0,0,0,below right); + \cnode(n1,2,0,0,1,below right); + \cnode(n2,0,2,0,2,below right); + \cnode(n3,0,0,2,3,right); + \cnode(n4,1,0,0,4,below right); + \cnode(n5,1,1,0,5,below right); + \cnode(n6,0,1,0,6,below right); + \cnode(n7,0,0,1,7,below right); + \cnode(n8,1,0,1,8,below right); + \cnode(n9,0,1,1,9,right); + """, + ), + "hexahedron": _MeshIOCellTypeInfo( + egrp_cls=TensorProductElementGroup, + dim=3, + order=1, + unit_nodes=r""" + \cnode(n0,0,0,0,0,below right); + \cnode(n1,2,0,0,1,below right); + \cnode(n2,2,2,0,2,below right); + \cnode(n3,0,2,0,3,below right); + \cnode(n4,0,0,2,4,below right); + \cnode(n5,2,0,2,5,below right); + \cnode(n6,2,2,2,6,below right); + \cnode(n7,0,2,2,7,below right); + """, + ), + "hexahedron27": _MeshIOCellTypeInfo( + egrp_cls=TensorProductElementGroup, + dim=3, + order=2, + unit_nodes=r"""" + \cnode(n0,0,0,0,0,below right); + \cnode(n1,2,0,0,1,below right); + \cnode(n2,2,2,0,2,below right); + \cnode(n3,0,2,0,3,below right); + \cnode(n4,0,0,2,4,below right); + \cnode(n5,2,0,2,5,below right); + \cnode(n6,2,2,2,6,below right); + \cnode(n7,0,2,2,7,below right); + \cnode(n8,1,0,0,8,below right); + \cnode(n9,2,1,0,9,below right); + \cnode(n10,1,2,0,10,below right); + \cnode(n11,0,1,0,11,below right); + \cnode(n12,1,0,2,12,below right); + \cnode(n13,2,1,2,13,below right); + \cnode(n14,1,2,2,14,below right); + \cnode(n15,0,1,2,15,below right); + \cnode(n16,0,0,1,16,below right); + \cnode(n17,2,0,1,17,below right); + \cnode(n18,2,2,1,18,below right); + \cnode(n19,0,2,1,19,below right); + \cnode(n20,0,1,1,20,below right); + \cnode(n21,2,1,1,21,below right); + \cnode(n22,1,0,1,22,below right); + \cnode(n23,1,2,1,23,below right); + \cnode(n24,1,1,0,24,below right); + \cnode(n25,1,1,2,25,below right); + \cnode(n26,1,1,1,26,below right); + """, + ) +} + + +def _mio_parse_unit_nodes(dim: int, unode_str: str) -> np.ndarray: + import re + cnode_re = re.compile(r"^\\cnode\((.+)\);$") + + result: List[List[int]] = [] + for ln in unode_str.split("\n"): + ln = ln.strip() + if not ln: + continue + match = cnode_re.match(ln) + assert match + node_id, c0, c1, c2, node_nr, _alignment = match.group(1).split(",") + assert node_id == f"n{len(result)}" + assert int(node_nr) == len(result) + coords = [int(c) for c in [c0, c1, c2]] + assert all(ci == 0 for ci in coords[dim:]) + coords = coords[:dim] + result.append([c-1 for c in coords]) + + return np.array(result).T.copy() + + +def _mio_cell_block_to_grp( + mio_mesh, cblock, force_ambient_dim: Optional[int] + ) -> MeshElementGroup: + try: + info = _MESHIO_CELL_TYPE_INFO[cblock.type] + except KeyError: + raise NotImplementedError(f"meshio cell block type '{cblock.type}'") + + unit_nodes = _mio_parse_unit_nodes(info.dim, info.unit_nodes) + + shape = info.egrp_cls._modepy_shape_cls(info.dim) + import modepy as mp + ref_unit_vertices = mp.unit_vertices_for_shape(shape) + unit_vertex_indices = find_point_permutation(ref_unit_vertices, unit_nodes) + assert unit_vertex_indices is not None + + assert cblock.data.shape[1] == unit_nodes.shape[1] + + nodes = mio_mesh.points[cblock.data].transpose(2, 0, 1).copy() + if force_ambient_dim is not None: + assert (nodes[force_ambient_dim:] == 0).all() + nodes = nodes[:force_ambient_dim].copy() + + return info.egrp_cls.make_group( + order=info.order, + vertex_indices=cblock.data[:, unit_vertex_indices].astype(np.int32), + nodes=nodes, + unit_nodes=unit_nodes + ) + + +def read_via_meshio( + file: Union[str, PathLike, BinaryIO, TextIO], + file_format: Optional[str] = None, + force_ambient_dim: Optional[int] = None, + ) -> Mesh: + from meshio import read + mio_mesh = read(file, file_format) + + max_dim = max(cblock.dim for cblock in mio_mesh.cells) + + vertices = mio_mesh.points.T.copy() + if force_ambient_dim is not None: + assert (vertices[force_ambient_dim:] == 0).all() + vertices = vertices[:force_ambient_dim].copy() + + vol_groups = [ + _mio_cell_block_to_grp( + mio_mesh, cb, force_ambient_dim=force_ambient_dim) + for cb in mio_mesh.cells + if cb.dim == max_dim] + + # FIXME: This is heuristic. + if len(vol_groups) == 1: + is_conforming = True + else: + is_conforming = max_dim < 3 + + from meshmode.mesh import make_mesh + return make_mesh( + vertices=vertices, + groups=vol_groups, + is_conforming=is_conforming, + ) + +# }}} + + # vim: foldmethod=marker diff --git a/setup.cfg b/setup.cfg index aba7af88..c58b7bf4 100644 --- a/setup.cfg +++ b/setup.cfg @@ -96,3 +96,6 @@ ignore_missing_imports = True [mypy-pytential.*] ignore_missing_imports = True + +[mypy-meshio.*] +ignore_missing_imports = True diff --git a/setup.py b/setup.py index c8a1d7e8..c736ab17 100644 --- a/setup.py +++ b/setup.py @@ -52,6 +52,7 @@ def main(): ], extras_require={ "visualization": ["h5py"], + "meshio": ["meshio"], "test": ["pytest>=2.3"], }, ) diff --git a/test/coarse.bdf b/test/coarse.bdf new file mode 100644 index 00000000..d0b36f2d --- /dev/null +++ b/test/coarse.bdf @@ -0,0 +1,204 @@ + +$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +$ +$ CUBIT version 2023.11 NX Nastran Translator +$ +$ File: coarse.bdf +$ Time Stamp: 26-Apr-24 at 09:46:50 +$ +$ +$ PLEASE CHECK YOUR MODEL FOR UNITS CONSISTENCY. +$ +$ It should be noted that load ID's from CUBIT may NOT correspond to Nastran SID's +$ The SID's for the load and restraint sets start at one and increment by one:i.e.,1,2,3,4... +$ +$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +$ +$ +$ ------------------------- +$ Executive Control Section +$ ------------------------- +$ +SOL 101 +CEND +$ +$ +$ -------------------- +$ Case Control Section +$ -------------------- +$ +ECHO = SORT +$ +$ +$ Name: Initial +$ +$ +$ Name: Default Set +$ +SUBCASE = 1 +$ +LABEL = Default Set +$ +$ ----------------- +$ Bulk Data Section +$ ----------------- +$ +BEGIN BULK +$ +$ Params +$ +$ +$ Node cards +$ +GRID 33 0 0.5 0.5 0 +GRID 17 0 0.3 0.5 0 +GRID 1 0 0.3 0.3 0 +GRID 32 0 0.5 0.3 0 +GRID 18 0 0.1 0.5 0 +GRID 2 0 0.1 0.3 0 +GRID 19 0 -0.1 0.5 0 +GRID 3 0 -0.1 0.3 0 +GRID 20 0 -0.3 0.5 0 +GRID 4 0 -0.3 0.3 0 +GRID 34 0 -0.5 0.5 0 +GRID 21 0 -0.5 0.3 0 +GRID 5 0 0.3 0.1 0 +GRID 31 0 0.5 0.1 0 +GRID 6 0 0.1 0.1 0 +GRID 7 0 -0.1 0.1 0 +GRID 8 0 -0.3 0.1 0 +GRID 22 0 -0.5 0.1 0 +GRID 9 0 0.3 -0.1 0 +GRID 30 0 0.5 -0.1 0 +GRID 10 0 0.1 -0.1 0 +GRID 11 0 -0.1 -0.1 0 +GRID 12 0 -0.3 -0.1 0 +GRID 23 0 -0.5 -0.1 0 +GRID 13 0 0.3 -0.3 0 +GRID 29 0 0.5 -0.3 0 +GRID 14 0 0.1 -0.3 0 +GRID 15 0 -0.1 -0.3 0 +GRID 16 0 -0.3 -0.3 0 +GRID 24 0 -0.5 -0.3 0 +GRID 28 0 0.3 -0.5 0 +GRID 36 0 0.5 -0.5 0 +GRID 27 0 0.1 -0.5 0 +GRID 26 0 -0.1 -0.5 0 +GRID 25 0 -0.3 -0.5 0 +GRID 35 0 -0.5 -0.5 0 +$ +$ Element cards +$ +$ +$ Name: fluid_block +$ +CQUAD4 1 1 33 17 1 32 +CQUAD4 2 1 17 18 2 1 +CQUAD4 3 1 18 19 3 2 +CQUAD4 4 1 19 20 4 3 +CQUAD4 5 1 20 34 21 4 +CQUAD4 6 1 32 1 5 31 +CQUAD4 7 1 1 2 6 5 +CQUAD4 8 1 2 3 7 6 +CQUAD4 9 1 3 4 8 7 +CQUAD4 10 1 4 21 22 8 +CQUAD4 11 1 31 5 9 30 +CQUAD4 12 1 5 6 10 9 +CQUAD4 13 1 6 7 11 10 +CQUAD4 14 1 7 8 12 11 +CQUAD4 15 1 8 22 23 12 +CQUAD4 16 1 30 9 13 29 +CQUAD4 17 1 9 10 14 13 +CQUAD4 18 1 10 11 15 14 +CQUAD4 19 1 11 12 16 15 +CQUAD4 20 1 12 23 24 16 +CQUAD4 21 1 29 13 28 36 +CQUAD4 22 1 13 14 27 28 +CQUAD4 23 1 14 15 26 27 +CQUAD4 24 1 15 16 25 26 +CQUAD4 25 1 16 24 35 25 +$ +$ Name: inflow_block +$ +CROD 26 2 35 25 +CROD 27 2 25 26 +CROD 28 2 26 27 +CROD 29 2 27 28 +CROD 30 2 28 36 +$ +$ Name: wall_block +$ +CROD 31 3 34 21 +CROD 32 3 21 22 +CROD 33 3 22 23 +CROD 34 3 23 24 +CROD 35 3 24 35 +CROD 36 3 33 17 +CROD 37 3 17 18 +CROD 38 3 18 19 +CROD 39 3 19 20 +CROD 40 3 20 34 +CROD 41 3 36 29 +CROD 42 3 29 30 +CROD 43 3 30 31 +CROD 44 3 31 32 +CROD 45 3 32 33 +$ +$ Name: all_boundary_block +$ +CROD 26 4 35 25 +CROD 27 4 25 26 +CROD 28 4 26 27 +CROD 29 4 27 28 +CROD 30 4 28 36 +CROD 36 4 33 17 +CROD 37 4 17 18 +CROD 38 4 18 19 +CROD 39 4 19 20 +CROD 40 4 20 34 +CROD 31 4 34 21 +CROD 32 4 21 22 +CROD 33 4 22 23 +CROD 34 4 23 24 +CROD 35 4 24 35 +CROD 41 4 36 29 +CROD 42 4 29 30 +CROD 43 4 30 31 +CROD 44 4 31 32 +CROD 45 4 32 33 +$ +$ Property cards +$ +$ +$ Name: fluid_block +$ +PSHELL 1 100 1 +$ +$ Name: inflow_block +$ +PROD 2 100 0 0 0 0 +$ +$ Name: wall_block +$ +PROD 3 100 0 0 0 0 +$ +$ Name: all_boundary_block +$ +PROD 4 100 0 0 0 0 +$ +$ Material cards +$ +$ +$ Name: Default-Steel +$ +MAT1* 100 206800 80155.039 0.29 +* 7e-06 1.2e-05 +$ +$ Restraint cards +$ +$ +$ Load cards +$ +$ +$ +ENDDATA \ No newline at end of file diff --git a/test/test_mesh.py b/test/test_mesh.py index 22cd615e..f99d637c 100644 --- a/test/test_mesh.py +++ b/test/test_mesh.py @@ -1467,6 +1467,11 @@ def test_urchin(): mgen.generate_urchin(3, 2, 4, 1e-4) +def test_read_nastran(): + pytest.importorskip("meshio") + mio.read_via_meshio(thisdir / "coarse.bdf", force_ambient_dim=2) + + if __name__ == "__main__": import sys if len(sys.argv) > 1: