diff --git a/examples/show_gmsh.py b/examples/show_gmsh.py index 0a30d6aa8..ce906de05 100644 --- a/examples/show_gmsh.py +++ b/examples/show_gmsh.py @@ -7,16 +7,17 @@ import load_sample_file -from gustaf import io +from gustaf import Vertices, io, show if __name__ == "__main__": mesh_file_tri = pathlib.Path("faces/tri/2DChannelTria.msh") mesh_file_quad = pathlib.Path("faces/quad/2DChannelQuad.msh") + mesh_file_tetra = pathlib.Path("volumes/tet/3DBrickTet.msh") base_samples_path = pathlib.Path("samples") - if not base_samples_path.exists(): - load_sample_file.load_sample_file(str(mesh_file_tri)) - load_sample_file.load_sample_file(str(mesh_file_quad)) + load_sample_file.load_sample_file(str(mesh_file_tri)) + load_sample_file.load_sample_file(str(mesh_file_quad)) + load_sample_file.load_sample_file(str(mesh_file_tetra)) # load the .msh file directly with the correct io module (meshio) loaded_mesh_tri = io.meshio.load(base_samples_path / mesh_file_tri) @@ -33,3 +34,44 @@ loaded_mesh_default = io.load(base_samples_path / mesh_file_tri) loaded_mesh_default.show() + + loaded_meshes = io.load( + base_samples_path / mesh_file_tri, + return_only_one_mesh=False, + set_boundary=True, + ) + boundary_node_indices = loaded_meshes[0].BC["edges-nodes"] + boundary_vertices = Vertices( + loaded_meshes[0].vertices[boundary_node_indices] + ) + + show.show_vedo( + ["faces", loaded_meshes[0]], + ["boundaries: edges", loaded_meshes[1]], + ["vertices", boundary_vertices], + ) + + rect_meshes = io.load( + base_samples_path / mesh_file_tetra, + return_only_one_mesh=False, + set_boundary=True, + ) + + show.show_vedo( + *[ + [name, mesh] + for name, mesh in zip(["volume", "face", "line"], rect_meshes) + ] + ) + + volume_nodes = Vertices( + rect_meshes[0].vertices[rect_meshes[0].BC["tri-nodes"]] + ) + faces_nodes = Vertices( + rect_meshes[0].vertices[rect_meshes[1].BC["edges-nodes"]] + ) + + show.show_vedo( + ["volumes: boundary-nodes", rect_meshes[0], volume_nodes], + ["faces: boundary-nodes", rect_meshes[1], faces_nodes], + ) diff --git a/gustaf/io/default.py b/gustaf/io/default.py index b15b16952..18dce089d 100644 --- a/gustaf/io/default.py +++ b/gustaf/io/default.py @@ -3,7 +3,7 @@ from gustaf.io import meshio, mfem, mixd -def load(fname): +def load(fname, **kwargs): """Load function for all supported file formats. This function tries to guess the correct io module for the given file. @@ -27,7 +27,7 @@ def load(fname): fname = pathlib.Path(fname).resolve() if fname.suffix in extensions_to_load_functions: - return extensions_to_load_functions[fname.suffix](fname) + return extensions_to_load_functions[fname.suffix](fname, **kwargs) else: raise ValueError( diff --git a/gustaf/io/meshio.py b/gustaf/io/meshio.py index 232373838..245b4a792 100644 --- a/gustaf/io/meshio.py +++ b/gustaf/io/meshio.py @@ -7,6 +7,7 @@ import numpy as np +from gustaf.edges import Edges from gustaf.faces import Faces from gustaf.helpers.raise_if import ModuleImportRaiser from gustaf.volumes import Volumes @@ -18,26 +19,34 @@ # from meshio import Mesh as MeshioMesh -def load(fname): +def load(fname, set_boundary=True): """Load mesh in meshio format. Loads vertices and their connectivity. - Currently cannot process boundary. - Note - ----- - This is more or less a direct copy from the original gustav implementation. - A lot of the meshio information are lost. When boundaries and multi-patch - definitions are added this needs to be revisited and extended. + The function reads all gustaf-processable element types from meshio, which + are sorted by the dimensionality of their element types. Different + elements can occur at the same dimensionality, e.g., in mixed meshes, + or as sub-topologies, e.g., as boundaries or interfaces. + + If set_boundary is set true, all (d-1)-dimensional submesh node identifiers + are assigned to MESH_TYPES.BC without check. + + Return_only_one_mesh ensures compatibility with previous gustaf versions. + + Loading physical groups is not supported yet. Parameters ------------ fname: str | pathlib.Path + Input filename + set_boundary : bool, optional, default is False + Assign nodes of lower-dimensional meshes mesh.BC + return_only_one_mesh : bool, optional, default is True + Return only the highest-dimensional mesh, ensures compartibility. Returns -------- - MESH_TYPES + MESH_TYPES or list(MESH_TYPES) """ - mesh_type = Faces - fname = pathlib.Path(fname) if not (fname.exists() and fname.is_file()): raise ValueError( @@ -45,31 +54,56 @@ def load(fname): f"{fname.resolve()}" ) + # Read mesh meshio_mesh: meshio.Mesh = meshio.read(fname) vertices = meshio_mesh.points - # check for 2D mesh - # Try for triangle grid - cells = meshio_mesh.get_cells_type("triangle") - # If no triangle elements, try for square - if len(cells) == 0: - cells = meshio_mesh.get_cells_type("quad") - if not len(cells) > 0: - # 3D mesh - mesh_type = Volumes - cells = meshio_mesh.get_cells_type("tetra") - if len(cells) == 0: - cells = meshio_mesh.get_cells_type("hexahedron") - - for i, cell in enumerate(cells): - if i == 0: - elements = cell.data - else: - elements = np.vstack((elements, cell.data)) - - mesh = mesh_type(vertices=vertices, elements=elements) + # Maps meshio element types to gustaf types and dimensionality + # Meshio vertex elements are ommited, as they do not follow gustaf logics + meshio_mapping = { + "hexahedron": [Volumes, 3], + "tetra": [Volumes, 3], + "quad": [Faces, 2], + "triangle": [Faces, 2], + "line": [Edges, 1], + } + meshes = [ + [ + meshio_mapping[key][0]( + vertices=vertices, elements=meshio_mesh.cells_dict[key] + ), + meshio_mapping[key][1], + ] + for key in meshio_mapping.keys() + if key in meshio_mesh.cells_dict.keys() + ] + # Sort by decreasing dimensionality + meshes.sort(key=lambda x: -x[1]) + + # Raises exception if the mesh file contains gustaf-incompartible types + for elem_type in meshio_mesh.cells_dict.keys(): + if elem_type not in list(meshio_mapping.keys()) + ["vertex"]: + # Be aware that gustaf currently only supports linear element + # types. + raise NotImplementedError( + f"Sorry, {elem_type} elements currently are not supported." + ) - return mesh + if set_boundary is True: + # mesh.BC is only defined for volumes and faces + for i in range(len(meshes)): + # Indicator for boundary + if meshes[i][1] in [3, 2]: + for j in range(i, len(meshes)): + # Considers only (d-1)-dimensional subspaces + if meshes[j][1] == meshes[i][1] - 1: + meshes[i][0].BC[ + f"{meshes[j][0].whatami}-nodes" + ] = np.unique(meshes[j][0].elements) + + return_value = [mesh[0] for mesh in meshes] + + return return_value def export(mesh, fname, submeshes=None, **kwargs):