Redo #143

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],
+    )

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(

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,58 +19,91 @@
     # 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(
             "The given file does not point to file. The given path is: "
             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):