Correction of erroneous copy of pre-split tags to new split meshes

docs/requirements.txt

@@ -1,5 +1,6 @@
+sphinx >= 7.4.7
 sphinx_rtd_theme
-sphinx-argparse
+sphinx-argparse >= 0.5.2
 sphinx-design
 # Running CLI programs and capture outputs
 sphinxcontrib-programoutput>=0.17

geos-mesh/src/geos/mesh/conversion/main.py

@@ -1,6 +1,7 @@
 import argparse
 import logging
 import sys
+from geos.mesh.conversion import abaqus_converter
 
 
 def build_abaqus_converter_input_parser() -> argparse.ArgumentParser:
@@ -25,8 +26,6 @@ def main() -> None:
         output (str): Output mesh file name
         -v/--verbose (flag): Increase verbosity level
     """
-    from geosx_mesh_tools import abaqus_converter
-
     # Parse the user arguments
     parser = build_abaqus_converter_input_parser()
     args = parser.parse_args()

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

@@ -0,0 +1 @@
+# Empty

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

@@ -23,7 +23,7 @@
     vtkXMLMultiBlockDataReader, )
 from vtkmodules.util.numpy_support import (
     vtk_to_numpy, )
-from vtk_utils import (
+from .vtk_utils import (
     vtk_iter, )
 
 

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

@@ -1,5 +1,6 @@
 from collections import defaultdict
 from dataclasses import dataclass
+from sys import exit
 import logging
 from typing import (
     Collection,
@@ -19,6 +20,7 @@
 import networkx
 import numpy
 
+from vtk import vtkDataArray
 from vtkmodules.vtkCommonCore import (
     vtkIdList,
     vtkPoints,
@@ -252,45 +254,156 @@ def __call__( self, index: int ) -> int:
     return result
 
 
-def __copy_fields( old_mesh: vtkUnstructuredGrid, new_mesh: vtkUnstructuredGrid,
-                   collocated_nodes: Sequence[ int ] ) -> None:
+def cells_points_coordinates( mesh: vtkUnstructuredGrid ) -> dict[ int, list[ tuple[ float ] ] ]:
+    """Map each cell id to a list of its points coordinates.
+
+    Args:
+        mesh (vtkUnstructuredGrid): An unstructured mesh.
+
+    Returns:
+        dict[int, list[tuple[float]]]: {cell_id1: [(pt0_x, pt0_y, pt0_z), ...], ..., cell_idN: [...]}
+        for a mesh containing N cells.
     """
-    Copies the fields from the old mesh to the new one.
-    Point data will be duplicated for collocated nodes.
-    :param old_mesh: The mesh before the split.
-    :param new_mesh: The mesh after the split. Will receive the fields in place.
-    :param collocated_nodes: New index to old index.
-    :return: None
+    cells_nodes_coordinates: dict[ int, list[ tuple[ float ] ] ] = {}
+    for i in range( mesh.GetNumberOfCells() ):
+        cell: vtkCell = mesh.GetCell( i )
+        cells_nodes_coordinates[ i ] = []
+        cell_points = cell.GetPoints()
+        for v in range( cell.GetNumberOfPoints() ):
+            node_coordinates: tuple[ float ] = cell_points.GetPoint( v )
+            truncated_coordinates: list[ float ] = [ float( '%.3f' % ( coord ) ) for coord in node_coordinates ]
+            cells_nodes_coordinates[ i ].append( tuple( truncated_coordinates ) )
+    return cells_nodes_coordinates
+
+
+def link_new_cells_id_with_old_cells_id( old_mesh: vtkUnstructuredGrid,
+                                         new_mesh: vtkUnstructuredGrid ) -> dict[ int, int ]:
+    """After mapping each cell id to a list of its points coordinates for the old and new mesh,
+    it is possible to determine the link between old and new cell id by coordinates position.
+
+    Args:
+        old_mesh (vtkUnstructuredGrid): An unstructured mesh before splitting.
+        new_mesh (vtkUnstructuredGrid): An unstructured mesh after splitting the old_mesh.
+
+    Returns:
+        dict[int, int]: { new_cell_id: old_cell_id }
+    """
+    new_cell_ids_with_old_cell_ids: dict[ int, int ] = {}
+    cpc_old_mesh: dict[ int, list[ tuple[ float ] ] ] = cells_points_coordinates( old_mesh )
+    cpc_new_mesh: dict[ int, list[ tuple[ float ] ] ] = cells_points_coordinates( new_mesh )
+    for new_cell_id, new_coords in cpc_new_mesh.items():
+        for old_cell_id, old_coords in cpc_old_mesh.items():
+            if all( elem in new_coords for elem in old_coords ):
+                new_cell_ids_with_old_cell_ids[ new_cell_id ] = old_cell_id
+                break
+    return new_cell_ids_with_old_cell_ids
+
+
+def __copy_cell_data( old_mesh: vtkUnstructuredGrid, new_mesh: vtkUnstructuredGrid, is_fracture_mesh: bool = False ):
+    """Copy the cell data arrays from the old mesh to the new mesh.
+    If the new mesh is a fracture_mesh, the fracture_mesh has less cells than the old mesh.
+    Therefore, copying the entire old cell arrays to the new one will assignate invalid
+    indexing of the values of these arrays.
+    So here, we consider the new indexing of cells to add an array with the valid size of elements
+    and correct indexes.
+
+    Args:
+        old_mesh (vtkUnstructuredGrid): An unstructured mesh before splitting.
+        new_mesh (vtkUnstructuredGrid): An unstructured mesh after splitting the old_mesh.
+        is_fracture_mesh (bool, optional): True if dealing with a new_mesh being the fracture mesh.
+        Defaults to False.
     """
     # Copying the cell data.
     # The cells are the same, just their nodes support have changed.
+    if is_fracture_mesh:
+        new_to_old_cells: dict[ int, int ] = link_new_cells_id_with_old_cells_id( old_mesh, new_mesh )
     input_cell_data = old_mesh.GetCellData()
     for i in range( input_cell_data.GetNumberOfArrays() ):
-        input_array = input_cell_data.GetArray( i )
+        input_array: vtkDataArray = input_cell_data.GetArray( i )
         logging.info( f"Copying cell field \"{input_array.GetName()}\"." )
-        new_mesh.GetCellData().AddArray( input_array )
+        if is_fracture_mesh:
+            array_name: str = input_array.GetName()
+            old_values: list = vtk_to_numpy( input_array ).tolist()
+            useful_values: list = []
+            for old_cell_id in new_to_old_cells.values():
+                useful_values.append( old_values[ old_cell_id ] )
+            useful_input_array = numpy_to_vtk( numpy.array( useful_values ) )
+            useful_input_array.SetName( array_name )
+            new_mesh.GetCellData().AddArray( useful_input_array )
+        else:
+            new_mesh.GetCellData().AddArray( input_array )
+
+
+# TODO consider the fracture_mesh case ?
+def __copy_field_data( old_mesh: vtkUnstructuredGrid, new_mesh: vtkUnstructuredGrid ):
+    """Copy the field data arrays from the old mesh to the new mesh.
+    Does not consider the fracture_mesh case.
 
+    Args:
+        old_mesh (vtkUnstructuredGrid): An unstructured mesh before splitting.
+        new_mesh (vtkUnstructuredGrid): An unstructured mesh after splitting the old_mesh.
+    """
     # Copying field data. This data is a priori not related to geometry.
     input_field_data = old_mesh.GetFieldData()
     for i in range( input_field_data.GetNumberOfArrays() ):
         input_array = input_field_data.GetArray( i )
         logging.info( f"Copying field data \"{input_array.GetName()}\"." )
         new_mesh.GetFieldData().AddArray( input_array )
 
+
+# TODO consider the fracture_mesh case ?
+def __copy_point_data( old_mesh: vtkUnstructuredGrid, new_mesh: vtkUnstructuredGrid,
+                       collocated_nodes: Sequence[ int ] ):
+    """Copy the point data arrays from the old mesh to the new mesh.
+    Does not consider the fracture_mesh case.
+
+    Args:
+        old_mesh (vtkUnstructuredGrid): An unstructured mesh before splitting.
+        new_mesh (vtkUnstructuredGrid): An unstructured mesh after splitting the old_mesh.
+    """
     # Copying the point data.
     input_point_data = old_mesh.GetPointData()
     for i in range( input_point_data.GetNumberOfArrays() ):
         input_array = input_point_data.GetArray( i )
-        logging.info( f"Copying point field \"{input_array.GetName()}\"" )
+        logging.info( f"Copying point field \"{input_array.GetName()}\"." )
         tmp = input_array.NewInstance()
         tmp.SetName( input_array.GetName() )
         tmp.SetNumberOfComponents( input_array.GetNumberOfComponents() )
         tmp.SetNumberOfTuples( new_mesh.GetNumberOfPoints() )
         for p in range( tmp.GetNumberOfTuples() ):
-            tmp.SetTuple( p, input_array.GetTuple( collocated_nodes[ p ] ) )
+            collocated_nodes_p = collocated_nodes[ p ]
+            if isinstance( collocated_nodes_p, numpy.integer ):
+                tmp.SetTuple( p, input_array.GetTuple( collocated_nodes_p ) )
+            else:  # when using collocated_nodes from fracture mesh, it ia an array
+                tmp.SetTuple( p, input_array.GetTuple( collocated_nodes_p[ 0 ] ) )
         new_mesh.GetPointData().AddArray( tmp )
 
 
+def __copy_fields( old_mesh: vtkUnstructuredGrid,
+                   new_mesh: vtkUnstructuredGrid,
+                   collocated_nodes: Sequence[ int ],
+                   is_fracture_mesh: bool = False ) -> None:
+    """
+    Copies the fields from the old mesh to the new one.
+    Point data will be duplicated for collocated nodes.
+    :param old_mesh: The mesh before the split.
+    :param new_mesh: The mesh after the split. Will receive the fields in place.
+    :param collocated_nodes: New index to old index.
+    :param is_fracture_mesh: True when copying fields for a fracture mesh, False instead.
+    :return: None
+    """
+    # Mesh cannot contains global ids before splitting.
+    if vtk_utils.has_invalid_field( old_mesh, [ "GLOBAL_IDS_POINTS", "GLOBAL_IDS_CELLS" ] ):
+        logging.critical( f"The mesh cannot contain global ids for neither cells nor points." )
+        logging.critical( f"The correct procedure is to split the mesh and then generate global" +
+                          "ids for new split meshes. Dying ..." )
+        exit( 1 )
+
+    __copy_cell_data( old_mesh, new_mesh, is_fracture_mesh )
+    __copy_field_data( old_mesh, new_mesh )
+    __copy_point_data( old_mesh, new_mesh, collocated_nodes )
+
+
 def __perform_split( old_mesh: vtkUnstructuredGrid,
                      cell_to_node_mapping: Mapping[ int, Mapping[ int, int ] ] ) -> vtkUnstructuredGrid:
     """
@@ -366,7 +479,7 @@ def __perform_split( old_mesh: vtkUnstructuredGrid,
     return new_mesh
 
 
-def __generate_fracture_mesh( mesh_points: vtkPoints, fracture_info: FractureInfo,
+def __generate_fracture_mesh( old_mesh: vtkUnstructuredGrid, fracture_info: FractureInfo,
                               cell_to_node_mapping: Mapping[ int, Mapping[ int, int ] ] ) -> vtkUnstructuredGrid:
     """
     Generates the mesh of the fracture.
@@ -377,6 +490,7 @@ def __generate_fracture_mesh( mesh_points: vtkPoints, fracture_info: FractureInf
     """
     logging.info( "Generating the meshes" )
 
+    mesh_points: vtkPoints = old_mesh.GetPoints()
     is_node_duplicated = numpy.zeros( mesh_points.GetNumberOfPoints(), dtype=bool )  # defaults to False
     for node_mapping in cell_to_node_mapping.values():
         for i, o in node_mapping.items():
@@ -399,11 +513,12 @@ def __generate_fracture_mesh( mesh_points: vtkPoints, fracture_info: FractureInf
         # for dfns in discarded_face_nodes:
         #     tmp.append(", ".join(map(str, dfns)))
         msg: str = "(" + '), ('.join( map( lambda dfns: ", ".join( map( str, dfns ) ), discarded_face_nodes ) ) + ")"
-        # logging.info(f"The {len(tmp)} faces made of nodes ({'), ('.join(tmp)}) were/was discarded from the fracture mesh because none of their/its nodes were duplicated.")
-        # print(f"The {len(tmp)} faces made of nodes ({'), ('.join(tmp)}) were/was discarded from the fracture mesh because none of their/its nodes were duplicated.")
-        print(
-            f"The faces made of nodes [{msg}] were/was discarded from the fracture mesh because none of their/its nodes were duplicated."
-        )
+        # logging.info(f"The {len(tmp)} faces made of nodes ({'), ('.join(tmp)}) were/was discarded"
+        #              + "from the fracture mesh because none of their/its nodes were duplicated.")
+        # print(f"The {len(tmp)} faces made of nodes ({'), ('.join(tmp)}) were/was discarded"
+        #              + "from the fracture mesh because none of their/its nodes were duplicated.")
+        print( f"The faces made of nodes [{msg}] were/was discarded" +
+               "from the fracture mesh because none of their/its nodes were duplicated." )
 
     fracture_nodes_tmp = numpy.ones( mesh_points.GetNumberOfPoints(), dtype=int ) * -1
     for ns in face_nodes:
@@ -448,6 +563,9 @@ def __generate_fracture_mesh( mesh_points: vtkPoints, fracture_info: FractureInf
     if polygons.GetNumberOfCells() > 0:
         fracture_mesh.SetCells( [ VTK_POLYGON ] * polygons.GetNumberOfCells(), polygons )
     fracture_mesh.GetPointData().AddArray( array )
+
+    __copy_fields( old_mesh, fracture_mesh, collocated_nodes, True )
+
     return fracture_mesh
 
 
@@ -458,7 +576,7 @@ def __split_mesh_on_fracture( mesh: vtkUnstructuredGrid,
     cell_to_node_mapping: Mapping[ int, Mapping[ int, int ] ] = __identify_split( mesh.GetNumberOfPoints(),
                                                                                   cell_to_cell, fracture.node_to_cells )
     output_mesh: vtkUnstructuredGrid = __perform_split( mesh, cell_to_node_mapping )
-    fractured_mesh: vtkUnstructuredGrid = __generate_fracture_mesh( mesh.GetPoints(), fracture, cell_to_node_mapping )
+    fractured_mesh: vtkUnstructuredGrid = __generate_fracture_mesh( mesh, fracture, cell_to_node_mapping )
     return output_mesh, fractured_mesh
 
 

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

@@ -50,6 +50,38 @@ def vtk_iter( l ) -> Iterator[ Any ]:
             yield l.GetCellType( i )
 
 
+def has_invalid_field( mesh: vtkUnstructuredGrid, invalid_fields: list[ str ] ) -> bool:
+    """Checks if a mesh contains at least a data arrays within its cell, field or point data
+    having a certain name. If so, returns True, else False.
+
+    Args:
+        mesh (vtkUnstructuredGrid): An unstructured mesh.
+        invalid_fields (list[str]): Field name of an array in any data from the data.
+
+    Returns:
+        bool: True if one field found, else False.
+    """
+    # Check the cell data fields
+    cell_data = mesh.GetCellData()
+    for i in range( cell_data.GetNumberOfArrays() ):
+        if cell_data.GetArrayName( i ) in invalid_fields:
+            logging.error( f"The mesh contains an invalid cell field name '{cell_data.GetArrayName( i )}'." )
+            return True
+    # Check the field data fields
+    field_data = mesh.GetFieldData()
+    for i in range( field_data.GetNumberOfArrays() ):
+        if field_data.GetArrayName( i ) in invalid_fields:
+            logging.error( f"The mesh contains an invalid field name '{field_data.GetArrayName( i )}'." )
+            return True
+    # Check the point data fields
+    point_data = mesh.GetPointData()
+    for i in range( point_data.GetNumberOfArrays() ):
+        if point_data.GetArrayName( i ) in invalid_fields:
+            logging.error( f"The mesh contains an invalid point field name '{point_data.GetArrayName( i )}'." )
+            return True
+    return False
+
+
 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,6 +115,9 @@ 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.
     """
+    if not os.path.exists( vtk_input_file ):
+        logging.critical( f"Invalid file path. COuld not read \"{vtk_input_file}\". Dying..." )
+        sys.exit( 1 )
     file_extension = os.path.splitext( vtk_input_file )[ -1 ]
     extension_to_reader = { ".vtk": __read_vtk, ".vtu": __read_vtu }
     # Testing first the reader that should match