refactoring of the reference mesh (shape mesher without pyvista)

Author: otvam

docs/scripts/generate_example.py

@@ -39,15 +39,15 @@ def run_image(folder_example):
     data_viewer["domain"]["data_options"]["plot_theme"]["background_color"] = "white"
     data_viewer["domain"]["data_options"]["plot_theme"]["axis_add"] = False
     data_viewer["domain"]["data_options"]["plot_view"]["grid_plot"] = False
-    data_viewer["domain"]["data_options"]["plot_view"]["geom_plot"] = False
+    data_viewer["domain"]["data_options"]["plot_view"]["voxel_plot"] = False
     data_viewer["domain"]["data_options"]["plot_view"]["point_plot"] = False
     data_viewer["domain"]["data_plot"]["title"] = None
 
     # run the mesher
     data_voxel = pypeec.run_mesher_data(data_geometry)
 
     # run the viewer
-    pypeec.run_viewer_data(data_voxel, data_viewer, tag_plot=["domain"])
+    pypeec.run_viewer_data(data_voxel, data_viewer, tag_plot=["domain"], plot_mode="qt")
 
 
 if __name__ == "__main__":

examples/config/cfg_data_pyvista.yaml

@@ -28,10 +28,10 @@
     "grid_thickness": 1.0              # line thickness for the complete voxel geometry
     "grid_color": "black"              # line opacity for the complete voxel geometry
     "grid_opacity": 0.1                # line color for the complete voxel geometry
-    "geom_plot": true                  # plot (or not) the non-empty voxels as wireframe
-    "geom_thickness": 1.0              # line thickness for the non-empty voxels
-    "geom_color": "black"              # line color for the non-empty voxels
-    "geom_opacity": 0.1                # line opacity for the non-empty voxels
+    "voxel_plot": true                 # plot (or not) the non-empty voxels as wireframe
+    "voxel_thickness": 1.0             # line thickness for the non-empty voxels
+    "voxel_color": "black"             # line color for the non-empty voxels
+    "voxel_opacity": 0.1               # line opacity for the non-empty voxels
     "point_plot": true                 # plot (or not) the provided point cloud as dots
     "point_color": "red"               # color of the point cloud dots
     "point_size": 8.0                  # size of the point cloud dots

pypeec/data/schema_plot_pyvista.yaml

@@ -68,10 +68,10 @@
             - "grid_thickness"
             - "grid_color"
             - "grid_opacity"
-            - "geom_plot"
-            - "geom_thickness"
-            - "geom_color"
-            - "geom_opacity"
+            - "voxel_plot"
+            - "voxel_thickness"
+            - "voxel_color"
+            - "voxel_opacity"
             - "point_plot"
             - "point_color"
             - "point_size"
@@ -92,15 +92,15 @@
                 "type": "number"
                 "minimum": 0
                 "maximum": 1
-            "geom_plot":
+            "voxel_plot":
                 "type": "boolean"
-            "geom_thickness":
+            "voxel_thickness":
                 "type": "number"
                 "minimum": 0
-            "geom_color":
+            "voxel_color":
                 "type": "string"
                 "minLength": 1
-            "geom_opacity":
+            "voxel_opacity":
                 "type": "number"
                 "minimum": 0
                 "maximum": 1

pypeec/lib_mesher/mesher_png.py

@@ -178,7 +178,7 @@ def get_mesh(param, domain_color, layer_stack):
     size = param["size"]
 
     # no reference geometry, direct voxelization
-    reference = None
+    geom_def = []
 
     # get the image size
     (nx, ny) = size
@@ -211,4 +211,4 @@ def get_mesh(param, domain_color, layer_stack):
     # assemble
     n = [nx, ny, nz]
 
-    return n, d, c, domain_def, reference
+    return n, d, c, domain_def, geom_def

pypeec/lib_mesher/mesher_shape.py

@@ -21,18 +21,13 @@
 __copyright__ = "Thomas Guillod - Dartmouth College"
 __license__ = "Mozilla Public License Version 2.0"
 
-import vtk
 import warnings
 import scilogger
 import numpy as np
-import pyvista as pv
 import shapely as sha
 import rasterio.features as raf
 import rasterio.transform as rat
 
-# prevent VTK to mess up the output
-vtk.vtkObject.GlobalWarningDisplayOff()
-
 # prevent problematic linear transform to trigger warnings
 warnings.filterwarnings("ignore", module="shapely")
 warnings.filterwarnings("ignore", module="rasterio.features")
@@ -44,7 +39,7 @@
 
 def _get_boundary_polygon(bnd, z_min):
     """
-    Convert a Shapely boundary into a PyVista polygon.
+    Convert a Shapely boundary into a line mesh.
     """
 
     # check that the boundary is closed
@@ -58,24 +53,29 @@ def _get_boundary_polygon(bnd, z_min):
     # get the 3D boundary
     z = np.full(len(xy), z_min, dtype=np.float64)
     z = np.expand_dims(z, axis=1)
-    xyz = np.hstack((xy, z))
+    points = np.hstack((xy, z))
 
     # remove the duplicated points
-    xyz = xyz[:-1]
+    points = points[:-1]
 
-    # get the face indices
-    lines = np.arange(len(xyz))
-    lines = np.concatenate(([len(xyz)+1], lines, [0]))
+    # get the indices
+    faces = np.empty(0, dtype=np.int64)
+    lines = np.arange(len(points), dtype=np.int64)
+    lines = np.concatenate(([len(points)+1], lines, [0]))
 
     # create the polygon
-    mesh = pv.PolyData(xyz, lines=lines)
+    mesh = {
+        "faces": np.array(faces, dtype=np.int64),
+        "lines": np.array(lines, dtype=np.int64),
+        "points": np.array(points, dtype=np.float64),
+    }
 
     return mesh
 
 
-def _get_shape_mesh(z_min, z_max, obj):
+def _get_shape_mesh(obj, z_min, z_max):
     """
-    Extrude a Shapely polygon into a PyVista mesh.
+    Extrude a Shapely polygon into line meshes.
     """
 
     # ensure that the polygon has a positive orientation
@@ -85,19 +85,19 @@ def _get_shape_mesh(z_min, z_max, obj):
     bnd = obj.exterior
     holes = obj.interiors
 
-    # create an empty mesh
-    mesh = pv.PolyData()
+    # init mesh list
+    geom_def = []
 
     # polygon for the external boundaries
-    mesh += _get_boundary_polygon(bnd, z_min)
-    mesh += _get_boundary_polygon(bnd, z_max)
+    geom_def.append(_get_boundary_polygon(bnd, z_min))
+    geom_def.append(_get_boundary_polygon(bnd, z_max))
 
     # polygon for the holes
     for bnd in holes:
-        mesh += _get_boundary_polygon(bnd, z_min)
-        mesh += _get_boundary_polygon(bnd, z_max)
+        geom_def.append(_get_boundary_polygon(bnd, z_min))
+        geom_def.append(_get_boundary_polygon(bnd, z_max))
 
-    return mesh
+    return geom_def
 
 
 def _get_shape_single(tag, shape_type, shape_data):
@@ -426,13 +426,11 @@ def _get_voxel_size(dx, dy, dz, stack_pos, xy_max_obj, xy_min_obj, xy_max, xy_mi
 def _get_merge_shape(stack_pos, shape_obj):
     """
     Transform all the 2D vector shapes into 3D meshes.
-    Merge all the meshes into a single mesh.
-    The resulting mesh represent the original geometry.
-    This mesh can be used to assess the quality of the voxelization.
+    The resulting meshes represent the original geometry.
     """
 
-    # create an empty mesh
-    reference = pv.PolyData()
+    # init mesh list
+    geom_def = []
 
     # merge all the shapes
     for shape_obj_tmp in shape_obj:
@@ -446,14 +444,14 @@ def _get_merge_shape(stack_pos, shape_obj):
 
         # transform the shapes into meshes
         if isinstance(obj, sha.Polygon):
-            reference += _get_shape_mesh(z_min, z_max, obj)
+            geom_def += _get_shape_mesh(obj, z_min, z_max)
         elif isinstance(obj, sha.MultiPolygon):
             for obj_tmp in obj.geoms:
-                reference += _get_shape_mesh(z_min, z_max, obj_tmp)
+                geom_def += _get_shape_mesh(obj_tmp, z_min, z_max)
         else:
             raise ValueError("invalid shape type")
 
-    return reference
+    return geom_def
 
 
 def get_mesh(param, layer_stack, geometry_shape):
@@ -496,7 +494,7 @@ def get_mesh(param, layer_stack, geometry_shape):
 
         # merge the shapes representing the original geometry
         LOGGER.debug("merge the shapes")
-        reference = _get_merge_shape(stack_pos, shape_obj)
+        geom_def = _get_merge_shape(stack_pos, shape_obj)
 
     # voxelize the shapes
     LOGGER.debug("voxelize the shapes")
@@ -508,11 +506,4 @@ def get_mesh(param, layer_stack, geometry_shape):
     d = d.tolist()
     c = c.tolist()
 
-    # cast reference mesh
-    reference = {
-        "faces": np.array(reference.faces, dtype=np.int64),
-        "lines": np.array(reference.lines, dtype=np.int64),
-        "points": np.array(reference.points, dtype=np.float64),
-    }
-
-    return n, d, c, domain_def, reference
+    return n, d, c, domain_def, geom_def

pypeec/lib_mesher/mesher_stl.py

@@ -84,9 +84,11 @@ def _get_mesh_stl(domain_stl, check):
     Find the bounding box for all the meshes (minimum and maximum coordinates).
     """
 
-    # init STL mesh list
+    # init mesh list (for the voxelization)
     mesh_stl = []
-    mesh_all = []
+
+    # init mesh list (for the reference meshes)
+    geom_def = []
 
     # init the coordinate (minimum and maximum coordinates)
     xyz_min = np.full(3, +np.inf, dtype=np.float64)
@@ -113,14 +115,16 @@ def _get_mesh_stl(domain_stl, check):
             xyz_max = np.maximum(xyz_max, tmp_max)
 
             # assign the mesh
-            mesh_all.append(mesh)
             mesh_stl.append({"tag": tag, "mesh": mesh})
 
-    # merge all the meshes
-    reference = pv.MultiBlock(mesh_all)
-    reference = reference.combine().extract_surface()
+            # cast the surface mesh to a dict
+            geom_def.append({
+                "faces": np.array(mesh.faces, dtype=np.int64),
+                "lines": np.array(mesh.lines, dtype=np.int64),
+                "points": np.array(mesh.points, dtype=np.float64),
+            })
 
-    return mesh_stl, reference, xyz_min, xyz_max
+    return mesh_stl, geom_def, xyz_min, xyz_max
 
 
 def _get_voxel_size(d, xyz_max_stl, xyz_min_stl, xyz_max, xyz_min):
@@ -303,7 +307,7 @@ def get_mesh(param, domain_stl):
     with LOGGER.BlockIndent():
         # load the mesh and get the STL bounds
         LOGGER.debug("load STL files")
-        (mesh_stl, reference, xyz_min_stl, xyz_max_stl) = _get_mesh_stl(domain_stl, check)
+        (mesh_stl, geom_def, xyz_min_stl, xyz_max_stl) = _get_mesh_stl(domain_stl, check)
 
         # geometry size
         LOGGER.debug("get the voxel size")
@@ -331,11 +335,4 @@ def get_mesh(param, domain_stl):
     d = d.tolist()
     c = c.tolist()
 
-    # cast reference mesh
-    reference = {
-        "faces": np.array(reference.faces, dtype=np.int64),
-        "lines": np.array(reference.lines, dtype=np.int64),
-        "points": np.array(reference.points, dtype=np.float64),
-    }
-
-    return n, d, c, domain_def, reference
+    return n, d, c, domain_def, geom_def

pypeec/lib_mesher/mesher_voxel.py

@@ -29,7 +29,7 @@ def get_mesh(param, domain_index):
     nv = np.prod(n)
 
     # no reference geometry, direct voxelization
-    reference = None
+    geom_def = []
 
     # init new domain indices
     domain_def = {}
@@ -51,4 +51,4 @@ def get_mesh(param, domain_index):
         # add the new item
         domain_def[tag] = idx_tmp
 
-    return n, d, c, domain_def, reference
+    return n, d, c, domain_def, geom_def

pypeec/lib_plot/manage_pyvista.py

@@ -48,13 +48,13 @@ def _get_plot_view_theme(pl, grid, voxel, point, plot_view, plot_theme):
         )
 
     # plot the non-empty voxels
-    if plot_view["geom_plot"] and (voxel.n_cells > 0):
+    if plot_view["voxel_plot"] and (voxel.n_cells > 0):
         pl.add_mesh(
             voxel,
             style="wireframe",
-            color=plot_view["geom_color"],
-            opacity=plot_view["geom_opacity"],
-            line_width=plot_view["geom_thickness"],
+            color=plot_view["voxel_color"],
+            opacity=plot_view["voxel_opacity"],
+            line_width=plot_view["voxel_thickness"],
         )
 
     # plot the point cloud

pypeec/lib_plot/parse_voxel.py

@@ -64,17 +64,19 @@ def get_point(pts_cloud):
     return point
 
 
-def get_reference(reference):
+def get_reference(geom_def):
     """
     Construct a PyVista object from the reference mesh.
     """
 
-    if reference is None:
-        reference = pv.PolyData()
-    else:
-        points = reference["points"]
-        lines = reference["lines"]
-        faces = reference["faces"]
-        reference = pv.PolyData(points, lines=lines, faces=faces)
+    # create an empty mesh
+    reference = pv.PolyData()
+
+    # add the different reference mesh objects
+    for geom_tmp in geom_def:
+        points = geom_tmp["points"]
+        lines = geom_tmp["lines"]
+        faces = geom_tmp["faces"]
+        reference += pv.PolyData(points, lines=lines, faces=faces)
 
     return reference