Fix elements orderings to use cells volume

docs/geos-mesh.rst

@@ -123,26 +123,35 @@ This can be convenient if you cannot regenerate the mesh.
 .. code-block::
 
       $ python src/geos/mesh/doctor/mesh_doctor.py fix_elements_orderings --help
-      usage: mesh_doctor.py fix_elements_orderings [-h] [--Hexahedron 1,6,5,4,7,0,2,3] [--Prism5 8,2,0,7,6,9,5,1,4,3]
-                                                   [--Prism6 11,2,8,10,5,0,9,7,6,1,4,3] [--Pyramid 3,4,0,2,1]
-                                                   [--Tetrahedron 2,0,3,1] [--Voxel 1,6,5,4,7,0,2,3]
-                                                   [--Wedge 3,5,4,0,2,1] --output OUTPUT [--data-mode binary, ascii]
+      usage: mesh_doctor.py fix_elements_orderings [-h] [--Tetrahedron 2,0,3,1] [--Pyramid 3,4,0,2,1] [--Wedge 3,5,4,0,2,1]
+                                                   [--Hexahedron 1,6,5,4,7,0,2,3] [--Prism5 8,2,0,7,6,9,5,1,4,3] [--Prism6 11,2,8,10,5,0,9,7,6,1,4,3]
+                                                   --volume_to_reorder all,positive,negative --output OUTPUT [--data-mode binary, ascii]
 
       options:
-      -h, --help              show this help message and exit
+      -h, --help                  show this help message and exit
+      --Tetrahedron 2,0,3,1       [list of integers]: node permutation for "Tetrahedron".
+      --Pyramid 3,4,0,2,1         [list of integers]: node permutation for "Pyramid".
+      --Wedge 3,5,4,0,2,1         [list of integers]: node permutation for "Wedge".
       --Hexahedron 1,6,5,4,7,0,2,3
-                              [list of integers]: node permutation for "Hexahedron".
+                                  [list of integers]: node permutation for "Hexahedron".
       --Prism5 8,2,0,7,6,9,5,1,4,3
-                              [list of integers]: node permutation for "Prism5".
+                                  [list of integers]: node permutation for "Prism5".
       --Prism6 11,2,8,10,5,0,9,7,6,1,4,3
-                              [list of integers]: node permutation for "Prism6".
-      --Pyramid 3,4,0,2,1     [list of integers]: node permutation for "Pyramid".
-      --Tetrahedron 2,0,3,1   [list of integers]: node permutation for "Tetrahedron".
-      --Voxel 1,6,5,4,7,0,2,3 [list of integers]: node permutation for "Voxel".
-      --Wedge 3,5,4,0,2,1     [list of integers]: node permutation for "Wedge".
-      --output OUTPUT         [string]: The vtk output file destination.
-      --data-mode binary, ascii
-                              [string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary.
+                                  [list of integers]: node permutation for "Prism6".
+      --volume_to_reorder all,positive,negative
+                                  [str]: Select which element volume is invalid and needs reordering. 'all' will allow reordering of nodes for every element, regarding of their volume.
+                                  'positive' or 'negative' will only reorder the element with the corresponding volume.
+      --output OUTPUT             [string]: The vtk output file destination.
+      --data-mode binary, ascii   [string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary.
+
+For example, assume that you have a mesh that contains 3 different cell types like Hexahedrons, Pyramids and Tetrahedrons.
+After checking ``element_volumes``, you found that all your Hexahedrons and half of your Tetrahedrons have a negative volume.
+To correct that, assuming you know the correct node ordering to correct these volumes, you can use this command:
+
+.. code-block::
+
+      $ python src/geos/mesh/doctor/mesh_doctor.py -i /path/to/your/mesh/file fix_elements_orderings --Hexahedron 1,6,5,4,7,0,2,3
+      --Tetrahedron 2,0,3,1 --volume_to_reorder negative --output /new/path/for/your/new/mesh/reordered/file --data-mode binary
 
 ``generate_cube``
 """""""""""""""""
@@ -340,8 +349,4 @@ API
 ^^^
 
 .. automodule:: geos.mesh.conversion.abaqus_converter
-    :members:
-
-
-
-
+    :members:

geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py

@@ -1,65 +1,173 @@
 from dataclasses import dataclass
+import numpy as np
 import logging
-from typing import (
-    List,
-    Dict,
-    Set,
-    FrozenSet,
-)
-
-from vtkmodules.vtkCommonCore import (
-    vtkIdList, )
-
-from . import vtk_utils
-from .vtk_utils import (
-    to_vtk_id_list,
-    VtkOutput,
-)
+from vtk import vtkCellSizeFilter
+from vtkmodules.vtkCommonCore import vtkIdList
+from vtkmodules.util.numpy_support import vtk_to_numpy
+from vtkmodules.vtkCommonDataModel import ( vtkDataSet, VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE,
+                                            VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM )
+from .vtk_utils import VtkOutput, to_vtk_id_list, write_mesh, read_mesh
 
 
 @dataclass( frozen=True )
 class Options:
     vtk_output: VtkOutput
-    cell_type_to_ordering: Dict[ int, List[ int ] ]
+    cell_name_to_ordering: dict[ str, list[ int ] ]
+    volume_to_reorder: str
 
 
 @dataclass( frozen=True )
 class Result:
     output: str
-    unchanged_cell_types: FrozenSet[ int ]
+    reordering_stats: dict[ str, list[ int ] ]
+
+
+GEOS_ACCEPTED_TYPES = [ VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM ]
+# the number of different nodes that needs to be entered in parsing when dealing with a specific vtk element
+NAME_TO_VTK_TYPE = {
+    "Hexahedron": VTK_HEXAHEDRON,
+    "Tetrahedron": VTK_TETRA,
+    "Pyramid": VTK_PYRAMID,
+    "Wedge": VTK_WEDGE,
+    "Prism5": VTK_PENTAGONAL_PRISM,
+    "Prism6": VTK_HEXAGONAL_PRISM
+}
+VTK_TYPE_TO_NAME = { vtk_type: name for name, vtk_type in NAME_TO_VTK_TYPE.items() }
+
+
+# Knowing the calculation of cell volumes in vtk was discussed there: https://github.com/GEOS-DEV/GEOS/issues/2253
+# Here, we do not need to have the exact volumes matching the simulation softwares results
+# because we are mostly interested in knowing the sign of the volume for the rest of the workflow.
+# Therefore, there is no need to use vtkMeshQuality for specific cell types when vtkCellSizeFilter works with all types.
+def compute_mesh_cells_volume( mesh: vtkDataSet ) -> np.array:
+    """Generates a volume array that was calculated on all cells of a mesh.
+
+    Args:
+        mesh (vtkDataSet): A vtk grid.
+
+    Returns:
+        np.array: Volume for every cell of a mesh.
+    """
+    cell_size_filter = vtkCellSizeFilter()
+    cell_size_filter.SetInputData( mesh )
+    cell_size_filter.SetComputeVolume( True )
+    cell_size_filter.Update()
+    return vtk_to_numpy( cell_size_filter.GetOutput().GetCellData().GetArray( "Volume" ) )
+
+
+def get_cell_types_and_number( mesh: vtkDataSet ) -> tuple[ list[ int ] ]:
+    """Gets the cell type for every cell of a mesh and the amount for each cell type.
+
+    Args:
+        mesh (vtkDataSet): A vtk grid.
+
+    Raises:
+        ValueError: "Invalid type '{cell_type}' for GEOS is in the mesh. Dying ..."
+
+    Returns:
+        tuple[ list[ int ] ]: ( unique_cell_types, total_per_cell_types )
+    """
+    number_cells: int = mesh.GetNumberOfCells()
+    all_cells_type: np.array = np.ones( number_cells, dtype=int )
+    for cell_id in range( number_cells ):
+        all_cells_type[ cell_id ] = mesh.GetCellType( cell_id )
+    unique_cell_types, total_per_cell_types = np.unique( all_cells_type, return_counts=True )
+    unique_cell_types, total_per_cell_types = unique_cell_types.tolist(), total_per_cell_types.tolist()
+    for cell_type in unique_cell_types:
+        if cell_type not in GEOS_ACCEPTED_TYPES:
+            err_msg: str = f"Invalid type '{cell_type}' for GEOS is in the mesh. Dying ..."
+            logging.error( err_msg )
+            raise ValueError( err_msg )
+    return ( unique_cell_types, total_per_cell_types )
+
+
+def is_cell_to_reorder( cell_volume: str, options: Options ) -> bool:
+    """Check if the volume of vtkCell qualifies the cell to be reordered.
+
+    Args:
+        cell_volume (float): The volume of a vtkCell.
+        options (Options): Options defined by the user.
+
+    Returns:
+        bool: True if cell needs to be reordered
+    """
+    if options.volume_to_reorder == "all":
+        return True
+    if cell_volume == 0.0:
+        return True
+    sign_of_cell_volume: int = int( cell_volume / abs( cell_volume ) )
+    if options.volume_to_reorder == "positive" and sign_of_cell_volume == 1:
+        return True
+    elif options.volume_to_reorder == "negative" and sign_of_cell_volume == -1:
+        return True
+    return False
+
+
+def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple:
+    """From an old mesh, creates a new one where certain cell nodes are reordered to obtain a correct volume.
+
+    Args:
+        old_mesh (vtkDataSet): A vtk grid needing nodes to be reordered.
+        options (Options): Options defined by the user.
+
+    Returns:
+        tuple: ( vtkDataSet, reordering_stats )
+    """
+    unique_cell_types, total_per_cell_types = get_cell_types_and_number( old_mesh )
+    unique_cell_names: list[ str ] = [ VTK_TYPE_TO_NAME[ vtk_type ] for vtk_type in unique_cell_types ]
+    names_with_totals: dict[ str, int ] = { n: v for n, v in zip( unique_cell_names, total_per_cell_types ) }
+    # sorted dict allow for sorted output of reordering_stats
+    names_with_totals_sorted: dict[ str, int ] = dict( sorted( names_with_totals.items() ) )
+    useful_VTK_TYPEs: list[ int ] = [ NAME_TO_VTK_TYPE[ name ] for name in options.cell_name_to_ordering.keys() ]
+    all_cells_volume: np.array = compute_mesh_cells_volume( old_mesh )
+    # a new mesh with the same data is created from the old mesh
+    new_mesh: vtkDataSet = old_mesh.NewInstance()
+    new_mesh.CopyStructure( old_mesh )
+    new_mesh.CopyAttributes( old_mesh )
+    cells = new_mesh.GetCells()
+    # Statistics on how many cells have or have not been reordered
+    reordering_stats: dict[ str, list[ any ] ] = {
+        "Types reordered": list(),
+        "Number of cells reordered": list(),
+        "Types non reordered": list( names_with_totals_sorted.keys() ),
+        "Number of cells non reordered": list( names_with_totals_sorted.values() )
+    }
+    counter_cells_reordered: dict[ str, int ] = { name: 0 for name in options.cell_name_to_ordering.keys() }
+    # sorted dict allow for sorted output of reordering_stats
+    ounter_cells_reordered_sorted: dict[ str, int ] = dict( sorted( counter_cells_reordered.items() ) )
+    # Reordering operations
+    for cell_id in range( new_mesh.GetNumberOfCells() ):
+        vtk_type: int = new_mesh.GetCellType( cell_id )
+        if vtk_type in useful_VTK_TYPEs:
+            if is_cell_to_reorder( float( all_cells_volume[ cell_id ] ), options ):
+                cell_name: str = VTK_TYPE_TO_NAME[ vtk_type ]
+                support_point_ids = vtkIdList()
+                cells.GetCellAtId( cell_id, support_point_ids )
+                new_support_point_ids: list[ int ] = list()
+                node_ordering: list[ int ] = options.cell_name_to_ordering[ cell_name ]
+                for i in range( len( node_ordering ) ):
+                    new_support_point_ids.append( support_point_ids.GetId( node_ordering[ i ] ) )
+                cells.ReplaceCellAtId( cell_id, to_vtk_id_list( new_support_point_ids ) )
+                ounter_cells_reordered_sorted[ cell_name ] += 1
+    # Calculation of stats
+    for cell_name_reordered, amount in ounter_cells_reordered_sorted.items():
+        if amount > 0:
+            reordering_stats[ "Types reordered" ].append( cell_name_reordered )
+            reordering_stats[ "Number of cells reordered" ].append( amount )
+            index_non_reordered: int = reordering_stats[ "Types non reordered" ].index( cell_name_reordered )
+            reordering_stats[ "Number of cells non reordered" ][ index_non_reordered ] -= amount
+            if reordering_stats[ "Number of cells non reordered" ][ index_non_reordered ] == 0:
+                reordering_stats[ "Types non reordered" ].pop( index_non_reordered )
+                reordering_stats[ "Number of cells non reordered" ].pop( index_non_reordered )
+    return ( new_mesh, reordering_stats )
 
 
 def __check( mesh, options: Options ) -> Result:
-    # The vtk cell type is an int and will be the key of the following mapping,
-    # that will point to the relevant permutation.
-    cell_type_to_ordering: Dict[ int, List[ int ] ] = options.cell_type_to_ordering
-    unchanged_cell_types: Set[ int ] = set()  # For logging purpose
-
-    # Preparing the output mesh by first keeping the same instance type.
-    output_mesh = mesh.NewInstance()
-    output_mesh.CopyStructure( mesh )
-    output_mesh.CopyAttributes( mesh )
-
-    # `output_mesh` now contains a full copy of the input mesh.
-    # We'll now modify the support nodes orderings in place if needed.
-    cells = output_mesh.GetCells()
-    for cell_idx in range( output_mesh.GetNumberOfCells() ):
-        cell_type: int = output_mesh.GetCell( cell_idx ).GetCellType()
-        new_ordering = cell_type_to_ordering.get( cell_type )
-        if new_ordering:
-            support_point_ids = vtkIdList()
-            cells.GetCellAtId( cell_idx, support_point_ids )
-            new_support_point_ids = []
-            for i, v in enumerate( new_ordering ):
-                new_support_point_ids.append( support_point_ids.GetId( new_ordering[ i ] ) )
-            cells.ReplaceCellAtId( cell_idx, to_vtk_id_list( new_support_point_ids ) )
-        else:
-            unchanged_cell_types.add( cell_type )
-    is_written_error = vtk_utils.write_mesh( output_mesh, options.vtk_output )
-    return Result( output=options.vtk_output.output if not is_written_error else "",
-                   unchanged_cell_types=frozenset( unchanged_cell_types ) )
+    output_mesh, reordering_stats = reorder_nodes_to_new_mesh( mesh, options )
+    write_mesh( output_mesh, options.vtk_output )
+    return Result( output=options.vtk_output.output, reordering_stats=reordering_stats )
 
 
 def check( vtk_input_file: str, options: Options ) -> Result:
-    mesh = vtk_utils.read_mesh( vtk_input_file )
+    mesh = read_mesh( vtk_input_file )
     return __check( mesh, options )

geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py

@@ -1,7 +1,7 @@
 from dataclasses import dataclass
-import os.path
+from os import path, access, R_OK, W_OK
 import logging
-import sys
+from sys import exit
 from typing import (
     Any,
     Iterator,
@@ -50,6 +50,39 @@ def vtk_iter( l ) -> Iterator[ Any ]:
             yield l.GetCellType( i )
 
 
+#  Inspired from https://stackoverflow.com/a/78870363
+def is_filepath_valid( filepath: str, is_input=True ) -> bool:
+    """Checks if a filepath can be used to read or to create a new file.
+
+    Args:
+        filepath (str): A filepath.
+        is_input (bool, optional): If the filepath is used to read a file, use True. 
+        If the filepath is used to create a new file, use False. Defaults to True.
+
+    Returns:
+        bool: False if invalid, True instead.
+    """
+    dirname = path.dirname( filepath )
+    if not path.isdir( dirname ):
+        logging.error( f"The directory '{dirname}' specified does not exist." )
+        return False
+    if is_input:
+        if not access( dirname, R_OK ):
+            logging.error( f"You do not have rights to read in directory '{dirname}' specified for the file." )
+            return False
+        if not path.exists( filepath ):
+            logging.error( f"The file specified does not exist." )
+            return False
+    else:
+        if not access( dirname, W_OK ):
+            logging.error( f"You do not have rights to write in directory '{dirname}' specified for the file." )
+            return False
+        if path.exists( filepath ):
+            logging.error( f"A file with the same name already exists. No over-writing possible." )
+            return False
+    return True
+
+
 def __read_vtk( vtk_input_file: str ) -> Optional[ vtkUnstructuredGrid ]:
     reader = vtkUnstructuredGridReader()
     logging.info( f"Testing file format \"{vtk_input_file}\" using legacy format reader..." )
@@ -83,7 +116,10 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid:
         If first guess does not work, eventually all the others reader available will be tested.
     :return: A unstructured grid.
     """
-    file_extension = os.path.splitext( vtk_input_file )[ -1 ]
+    if not is_filepath_valid( vtk_input_file, True ):
+        logging.error( f"Invalid filepath to read. Dying ..." )
+        exit( 1 )
+    file_extension = path.splitext( vtk_input_file )[ -1 ]
     extension_to_reader = { ".vtk": __read_vtk, ".vtu": __read_vtu }
     # Testing first the reader that should match
     if file_extension in extension_to_reader:
@@ -97,7 +133,7 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid:
             return output_mesh
     # No reader did work. Dying.
     logging.critical( f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\". Dying..." )
-    sys.exit( 1 )
+    exit( 1 )
 
 
 def __write_vtk( mesh: vtkUnstructuredGrid, output: str ) -> int:
@@ -125,16 +161,16 @@ def write_mesh( mesh: vtkUnstructuredGrid, vtk_output: VtkOutput ) -> int:
     :param vtk_output: Where to write. The file extension will be used to select the VTK file format.
     :return: 0 in case of success.
     """
-    if os.path.exists( vtk_output.output ):
-        logging.error( f"File \"{vtk_output.output}\" already exists, nothing done." )
-        return 1
-    file_extension = os.path.splitext( vtk_output.output )[ -1 ]
+    if not is_filepath_valid( vtk_output.output, False ):
+        logging.error( f"Invalid filepath to write. Dying ..." )
+        exit( 1 )
+    file_extension = path.splitext( vtk_output.output )[ -1 ]
     if file_extension == ".vtk":
         success_code = __write_vtk( mesh, vtk_output.output )
     elif file_extension == ".vtu":
         success_code = __write_vtu( mesh, vtk_output.output, vtk_output.is_data_mode_binary )
     else:
         # No writer found did work. Dying.
         logging.critical( f"Could not find the appropriate VTK writer for extension \"{file_extension}\". Dying..." )
-        sys.exit( 1 )
+        exit( 1 )
     return 0 if success_code else 2  # the Write member function return 1 in case of success, 0 otherwise.

geos-mesh/src/geos/mesh/doctor/mesh_doctor.py

@@ -1,17 +1,17 @@
 import sys
+import logging
+from geos.mesh.doctor.parsing import CheckHelper
+from geos.mesh.doctor.parsing.cli_parsing import parse_and_set_verbosity
+import geos.mesh.doctor.register as register
 
+min_python_version = ( 3, 7 )
 try:
-    min_python_version = ( 3, 7 )
     assert sys.version_info >= min_python_version
 except AssertionError as e:
     print( f"Please update python to at least version {'.'.join(map(str, min_python_version))}." )
     sys.exit( 1 )
 
-import logging
-
-from geos.mesh.doctor.parsing import CheckHelper
-from geos.mesh.doctor.parsing.cli_parsing import parse_and_set_verbosity
-import geos.mesh.doctor.register as register
+MESH_DOCTOR_FILEPATH = __file__
 
 
 def main():