diff --git a/geos-mesh/src/geos/mesh/doctor/checks/check_fractures.py b/geos-mesh/src/geos/mesh/doctor/checks/check_fractures.py index 9f14aed..21695ad 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/check_fractures.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/check_fractures.py @@ -1,30 +1,14 @@ -from dataclasses import dataclass import logging - -from typing import ( - Collection, - FrozenSet, - Iterable, - Sequence, - Set, - Tuple, -) - -from tqdm import tqdm import numpy - -from vtkmodules.vtkCommonDataModel import ( - vtkUnstructuredGrid, - vtkCell, -) -from vtkmodules.vtkCommonCore import ( - vtkPoints, ) -from vtkmodules.vtkIOXML import ( - vtkXMLMultiBlockDataReader, ) -from vtkmodules.util.numpy_support import ( - vtk_to_numpy, ) -from vtk_utils import ( - vtk_iter, ) +from dataclasses import dataclass +from typing import Collection, Iterable, Sequence +from tqdm import tqdm +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkCell +from vtkmodules.vtkCommonCore import vtkPoints +from vtkmodules.vtkIOXML import vtkXMLMultiBlockDataReader +from vtkmodules.util.numpy_support import vtk_to_numpy +from geos.mesh.doctor.checks.vtk_utils import vtk_iter +from geos.mesh.doctor.checks.generate_fractures import Coordinates3D @dataclass( frozen=True ) @@ -44,7 +28,7 @@ class Result: def __read_multiblock( vtk_input_file: str, matrix_name: str, - fracture_name: str ) -> Tuple[ vtkUnstructuredGrid, vtkUnstructuredGrid ]: + fracture_name: str ) -> tuple[ vtkUnstructuredGrid, vtkUnstructuredGrid ]: reader = vtkXMLMultiBlockDataReader() reader.SetFileName( vtk_input_file ) reader.Update() @@ -73,9 +57,9 @@ def format_collocated_nodes( fracture_mesh: vtkUnstructuredGrid ) -> Sequence[ I def __check_collocated_nodes_positions( - matrix_points: Sequence[ Tuple[ float, float, float ] ], fracture_points: Sequence[ Tuple[ float, float, float ] ], - g2l: Sequence[ int ], collocated_nodes: Iterable[ Iterable[ int ] ] -) -> Collection[ Tuple[ int, Iterable[ int ], Iterable[ Tuple[ float, float, float ] ] ] ]: + matrix_points: Sequence[ Coordinates3D ], fracture_points: Sequence[ Coordinates3D ], g2l: Sequence[ int ], + collocated_nodes: Iterable[ Iterable[ int ] ] +) -> Collection[ tuple[ int, Iterable[ int ], Iterable[ Coordinates3D ] ] ]: issues = [] for li, bucket in enumerate( collocated_nodes ): matrix_nodes = ( fracture_points[ li ], ) + tuple( map( lambda gi: matrix_points[ g2l[ gi ] ], bucket ) ) @@ -98,14 +82,14 @@ def my_iter( ccc ): def __check_neighbors( matrix: vtkUnstructuredGrid, fracture: vtkUnstructuredGrid, g2l: Sequence[ int ], collocated_nodes: Sequence[ Iterable[ int ] ] ): - fracture_nodes: Set[ int ] = set() + fracture_nodes: set[ int ] = set() for bucket in collocated_nodes: for gi in bucket: fracture_nodes.add( g2l[ gi ] ) # For each face of each cell, # if all the points of the face are "made" of collocated nodes, # then this is a fracture face. - fracture_faces: Set[ FrozenSet[ int ] ] = set() + fracture_faces: set[ frozenset[ int ] ] = set() for c in range( matrix.GetNumberOfCells() ): cell: vtkCell = matrix.GetCell( c ) for f in range( cell.GetNumberOfFaces() ): @@ -116,7 +100,7 @@ def __check_neighbors( matrix: vtkUnstructuredGrid, fracture: vtkUnstructuredGri # Finding the cells for c in tqdm( range( fracture.GetNumberOfCells() ), desc="Finding neighbor cell pairs" ): cell: vtkCell = fracture.GetCell( c ) - cns: Set[ FrozenSet[ int ] ] = set() # subset of collocated_nodes + cns: set[ frozenset[ int ] ] = set() # subset of collocated_nodes point_ids = frozenset( vtk_iter( cell.GetPointIds() ) ) for point_id in point_ids: bucket = collocated_nodes[ point_id ] @@ -129,9 +113,8 @@ def __check_neighbors( matrix: vtkUnstructuredGrid, fracture: vtkUnstructuredGri if f in fracture_faces: found += 1 if found != 2: - logging.warning( - f"Something went wrong since we should have found 2 fractures faces (we found {found}) for collocated nodes {cns}." - ) + logging.warning( f"Something went wrong since we should have found 2 fractures faces (we found {found})" + + f" for collocated nodes {cns}." ) def __check( vtk_input_file: str, options: Options ) -> Result: diff --git a/geos-mesh/src/geos/mesh/doctor/checks/generate_fractures.py b/geos-mesh/src/geos/mesh/doctor/checks/generate_fractures.py index c21325d..76eaf56 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/generate_fractures.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/generate_fractures.py @@ -1,50 +1,26 @@ +import logging +import networkx +import numpy from collections import defaultdict from dataclasses import dataclass -import logging -from typing import ( - Collection, - Dict, - FrozenSet, - Iterable, - List, - Mapping, - Optional, - Set, - Sequence, - Tuple, -) from enum import Enum - from tqdm import tqdm -import networkx -import numpy - -from vtkmodules.vtkCommonCore import ( - vtkIdList, - vtkPoints, -) -from vtkmodules.vtkCommonDataModel import ( - vtkCell, - vtkCellArray, - vtkPolygon, - vtkUnstructuredGrid, - VTK_POLYGON, - VTK_POLYHEDRON, -) -from vtkmodules.util.numpy_support import ( - vtk_to_numpy, - numpy_to_vtk, -) +from typing import Collection, Iterable, Mapping, Optional, Sequence, TypeAlias +from vtk import vtkDataArray +from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints +from vtkmodules.vtkCommonDataModel import ( vtkCell, vtkCellArray, vtkPolygon, vtkUnstructuredGrid, VTK_POLYGON, + VTK_POLYHEDRON ) +from vtkmodules.util.numpy_support import vtk_to_numpy, numpy_to_vtk from vtkmodules.util.vtkConstants import VTK_ID_TYPE - -from . import vtk_utils -from .vtk_utils import ( - vtk_iter, - VtkOutput, - to_vtk_id_list, -) -from .vtk_polyhedron import ( - FaceStream, ) +from geos.mesh.doctor.checks.vtk_utils import ( VtkOutput, vtk_iter, to_vtk_id_list, read_mesh, write_mesh, + has_invalid_field ) +from geos.mesh.doctor.checks.vtk_polyhedron import FaceStream +""" +TypeAliases used in this file +""" +IDMapping: TypeAlias = Mapping[ int, int ] +CellsPointsCoords: TypeAlias = dict[ int, list[ tuple[ float ] ] ] +Coordinates3D: TypeAlias = tuple[ float ] class FracturePolicy( Enum ): @@ -56,7 +32,7 @@ class FracturePolicy( Enum ): class Options: policy: FracturePolicy field: str - field_values: FrozenSet[ int ] + field_values: frozenset[ int ] vtk_output: VtkOutput vtk_fracture_output: VtkOutput @@ -73,17 +49,17 @@ class FractureInfo: def build_node_to_cells( mesh: vtkUnstructuredGrid, - face_nodes: Iterable[ Iterable[ int ] ] ) -> Mapping[ int, Iterable[ int ] ]: - node_to_cells: Dict[ int, - Set[ int ] ] = defaultdict( set ) # TODO normally, just a list and not a set should be enough. + face_nodes: Iterable[ Iterable[ int ] ] ) -> dict[ int, Iterable[ int ] ]: + # TODO normally, just a list and not a set should be enough. + node_to_cells: dict[ int, set[ int ] ] = defaultdict( set ) - fracture_nodes: Set[ int ] = set() + fracture_nodes: set[ int ] = set() for fns in face_nodes: for n in fns: fracture_nodes.add( n ) for cell_id in tqdm( range( mesh.GetNumberOfCells() ), desc="Computing the node to cells mapping" ): - cell_points: FrozenSet[ int ] = frozenset( vtk_iter( mesh.GetCell( cell_id ).GetPointIds() ) ) + cell_points: frozenset[ int ] = frozenset( vtk_iter( mesh.GetCell( cell_id ).GetPointIds() ) ) intersection: Iterable[ int ] = cell_points & fracture_nodes for node in intersection: node_to_cells[ node ].add( cell_id ) @@ -92,8 +68,8 @@ def build_node_to_cells( mesh: vtkUnstructuredGrid, def __build_fracture_info_from_fields( mesh: vtkUnstructuredGrid, f: Sequence[ int ], - field_values: FrozenSet[ int ] ) -> FractureInfo: - cells_to_faces: Dict[ int, List[ int ] ] = defaultdict( list ) + field_values: frozenset[ int ] ) -> FractureInfo: + cells_to_faces: dict[ int, list[ int ] ] = defaultdict( list ) # For each face of each cell, we search for the unique neighbor cell (if it exists). # Then, if the 2 values of the two cells match the field requirements, # we store the cell and its local face index: this is indeed part of the surface that we'll need to be split. @@ -109,10 +85,10 @@ def __build_fracture_info_from_fields( mesh: vtkUnstructuredGrid, f: Sequence[ i for j in range( neighbor_cell_ids.GetNumberOfIds() ): # It's 0 or 1... neighbor_cell_id = neighbor_cell_ids.GetId( j ) if f[ neighbor_cell_id ] != f[ cell_id ] and f[ neighbor_cell_id ] in field_values: - cells_to_faces[ cell_id ].append( - i ) # TODO add this (cell_is, face_id) information to the fracture_info? - face_nodes: List[ Collection[ int ] ] = list() - face_nodes_hashes: Set[ FrozenSet[ int ] ] = set() # A temporary not to add multiple times the same face. + # TODO add this (cell_is, face_id) information to the fracture_info? + cells_to_faces[ cell_id ].append( i ) + face_nodes: list[ Collection[ int ] ] = list() + face_nodes_hashes: set[ frozenset[ int ] ] = set() # A temporary not to add multiple times the same face. for cell_id, faces_ids in tqdm( cells_to_faces.items(), desc="Extracting the faces of the fractures" ): cell = mesh.GetCell( cell_id ) for face_id in faces_ids: @@ -121,23 +97,23 @@ def __build_fracture_info_from_fields( mesh: vtkUnstructuredGrid, f: Sequence[ i if fnh not in face_nodes_hashes: face_nodes_hashes.add( fnh ) face_nodes.append( fn ) - node_to_cells: Mapping[ int, Iterable[ int ] ] = build_node_to_cells( mesh, face_nodes ) + node_to_cells: dict[ int, Iterable[ int ] ] = build_node_to_cells( mesh, face_nodes ) return FractureInfo( node_to_cells=node_to_cells, face_nodes=face_nodes ) def __build_fracture_info_from_internal_surfaces( mesh: vtkUnstructuredGrid, f: Sequence[ int ], - field_values: FrozenSet[ int ] ) -> FractureInfo: - node_to_cells: Dict[ int, List[ int ] ] = {} - face_nodes: List[ Collection[ int ] ] = [] + field_values: frozenset[ int ] ) -> FractureInfo: + node_to_cells: dict[ int, list[ int ] ] = defaultdict( list ) + face_nodes: list[ Collection[ int ] ] = list() for cell_id in tqdm( range( mesh.GetNumberOfCells() ), desc="Computing the face to nodes mapping" ): cell = mesh.GetCell( cell_id ) if cell.GetCellDimension() == 2: if f[ cell_id ] in field_values: - nodes = [] + nodes = list() for v in range( cell.GetNumberOfPoints() ): point_id: int = cell.GetPointId( v ) - node_to_cells[ point_id ] = [] + node_to_cells[ point_id ] = list() nodes.append( point_id ) face_nodes.append( tuple( nodes ) ) @@ -177,24 +153,24 @@ def build_cell_to_cell_graph( mesh: vtkUnstructuredGrid, fracture: FractureInfo """ # Faces are identified by their nodes. But the order of those nodes may vary while referring to the same face. # Therefore we compute some kinds of hashes of those face to easily detect if a face is part of the fracture. - tmp: List[ FrozenSet[ int ] ] = [] + tmp: list[ frozenset[ int ] ] = list() for fn in fracture.face_nodes: tmp.append( frozenset( fn ) ) - face_hashes: FrozenSet[ FrozenSet[ int ] ] = frozenset( tmp ) + face_hashes: frozenset[ frozenset[ int ] ] = frozenset( tmp ) # We extract the list of the cells that touch the fracture by at least one node. - cells: Set[ int ] = set() + cells: set[ int ] = set() for cell_ids in fracture.node_to_cells.values(): for cell_id in cell_ids: cells.add( cell_id ) # Using the last precomputed containers, we're now building the dict which connects # every face (hash) of the fracture to the cells that touch the face... - face_to_cells: Dict[ FrozenSet[ int ], List[ int ] ] = defaultdict( list ) + face_to_cells: dict[ frozenset[ int ], list[ int ] ] = defaultdict( list ) for cell_id in tqdm( cells, desc="Computing the cell to cell graph" ): cell: vtkCell = mesh.GetCell( cell_id ) for face_id in range( cell.GetNumberOfFaces() ): - face_hash: FrozenSet[ int ] = frozenset( vtk_iter( cell.GetFace( face_id ).GetPointIds() ) ) + face_hash: frozenset[ int ] = frozenset( vtk_iter( cell.GetFace( face_id ).GetPointIds() ) ) if face_hash not in face_hashes: face_to_cells[ face_hash ].append( cell_id ) @@ -208,7 +184,7 @@ def build_cell_to_cell_graph( mesh: vtkUnstructuredGrid, fracture: FractureInfo def __identify_split( num_points: int, cell_to_cell: networkx.Graph, - node_to_cells: Mapping[ int, Iterable[ int ] ] ) -> Mapping[ int, Mapping[ int, int ] ]: + node_to_cells: dict[ int, Iterable[ int ] ] ) -> dict[ int, IDMapping ]: """ For each cell, compute the node indices replacements. :param num_points: Number of points in the whole mesh (not the fracture). @@ -228,7 +204,7 @@ class NewIndex: def __init__( self, num_nodes: int ): self.__current_last_index = num_nodes - 1 - self.__seen: Set[ int ] = set() + self.__seen: set[ int ] = set() def __call__( self, index: int ) -> int: if index in self.__seen: @@ -239,10 +215,9 @@ def __call__( self, index: int ) -> int: return index build_new_index = NewIndex( num_points ) - result: Dict[ int, Dict[ int, int ] ] = defaultdict( dict ) - for node, cells in tqdm( - sorted( node_to_cells.items() ), # Iteration over `sorted` nodes to have a predictable result for tests. - desc="Identifying the node splits" ): + result: dict[ int, IDMapping ] = defaultdict( dict ) + # Iteration over `sorted` nodes to have a predictable result for tests. + for node, cells in tqdm( sorted( node_to_cells.items() ), desc="Identifying the node splits" ): for connected_cells in networkx.connected_components( cell_to_cell.subgraph( cells ) ): # Each group of connect cells need around `node` must consider the same `node`. # Separate groups must have different (duplicated) nodes. @@ -252,24 +227,102 @@ def __call__( self, index: int ) -> int: return result -def __copy_fields( old_mesh: vtkUnstructuredGrid, new_mesh: vtkUnstructuredGrid, - collocated_nodes: Sequence[ int ] ) -> None: +def truncated_coordinates_with_id( mesh: vtkUnstructuredGrid, decimals: int = 3 ) -> dict[ Coordinates3D, int ]: + """Creates a mapping of truncated points coordinates with the their point id for every point of a mesh. + + Args: + mesh (vtkUnstructuredGrid): An unstructured mesh. + decimals (int, optional): Number of decimals to keep for truncation. Defaults to 3. + + Returns: + dict[ Coordinates3D, int ]: { coords0: pt_id0, ..., coordsN: pt_idN } """ - 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 + points: vtkPoints = mesh.GetPoints() + coords_with_id: dict[ Coordinates3D, int ] = dict() + for point_id in range( points.GetNumberOfPoints() ): + pt_coords = points.GetPoint( point_id ) + truncated_pt_coords = tuple( [ round( coord, decimals ) for coord in pt_coords ] ) + coords_with_id[ truncated_pt_coords ] = point_id + return coords_with_id + + +def link_new_cells_with_old_cells_id( old_mesh: vtkUnstructuredGrid, new_mesh: vtkUnstructuredGrid ) -> IDMapping: + """After mapping each truncated point coordinates to a list of its points ids 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: + IDMapping: { new_cell_id: old_cell_id, ... } + """ + truncated_coords_with_id_old: dict[ Coordinates3D, int ] = truncated_coordinates_with_id( old_mesh ) + truncated_coords_with_id_new: dict[ Coordinates3D, int ] = truncated_coordinates_with_id( new_mesh ) + # Every new_mesh by convention in this workflow is a mesh extracted from the old_mesh. + # So the number of elements is lower or equal in new mesh than in old mesh so we'd rather iterate over it + new_pts_to_old_pts_id: IDMapping = dict() + for coords, new_pt_id in truncated_coords_with_id_new.items(): + old_pt_id: int = truncated_coords_with_id_old[ coords ] # We can do that because new_mesh is from old_mesh + new_pts_to_old_pts_id[ new_pt_id ] = old_pt_id + # Now we have a link between point ids from the new_mesh to the old_mesh + # So we can now link the cells of new_mesh to the ones of old_mesh + new_cells_with_old_cells_id: IDMapping = dict() + for new_cell_id in range( new_mesh.GetNumberOfCells() ): + new_cell_pt_ids: tuple[ int ] = tuple( vtk_iter( new_mesh.GetCell( new_cell_id ).GetPointIds() ) ) + old_cell_pt_ids: tuple[ int ] = tuple( [ new_pts_to_old_pts_id[ new_pt_id ] for new_pt_id in new_cell_pt_ids ] ) + for old_cell_id in range( old_mesh.GetNumberOfCells() ): + pt_ids: tuple[ int ] = tuple( vtk_iter( old_mesh.GetCell( old_cell_id ).GetPointIds() ) ) + if pt_ids == old_cell_pt_ids: # the old cell was identified with the new cell + new_cells_with_old_cells_id[ new_cell_id ] = old_cell_id + break + return new_cells_with_old_cells_id + + +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: IDMapping = link_new_cells_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: numpy.array = vtk_to_numpy( input_array ) + useful_values: list = 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() ): @@ -277,29 +330,58 @@ def __copy_fields( old_mesh: vtkUnstructuredGrid, new_mesh: vtkUnstructuredGrid, 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 ): + """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 ] ) ) + tmp.DeepCopy( input_array ) new_mesh.GetPointData().AddArray( tmp ) -def __perform_split( old_mesh: vtkUnstructuredGrid, - cell_to_node_mapping: Mapping[ int, Mapping[ int, int ] ] ) -> vtkUnstructuredGrid: +def __copy_fields( old_mesh: vtkUnstructuredGrid, new_mesh: vtkUnstructuredGrid, is_fracture_mesh: bool = False ): + """ + 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 has_invalid_field( old_mesh, [ "GLOBAL_IDS_POINTS", "GLOBAL_IDS_CELLS" ] ): + err_msg: str = ( + "The mesh cannot contain global ids for neither cells nor points. The correct procedure is to split" + + " the mesh and then generate global ids for new split meshes. Dying ..." ) + logging.error( err_msg ) + raise ValueError( err_msg ) + + __copy_cell_data( old_mesh, new_mesh, is_fracture_mesh ) + __copy_field_data( old_mesh, new_mesh ) + __copy_point_data( old_mesh, new_mesh ) + + +def __perform_split( old_mesh: vtkUnstructuredGrid, cell_to_node_mapping: Mapping[ int, + IDMapping ] ) -> vtkUnstructuredGrid: """ Split the main 3d mesh based on the node duplication information contained in @p cell_to_node_mapping :param old_mesh: The main 3d mesh. :param cell_to_node_mapping: For each cell, gives the nodes that must be duplicated and their new index. :return: The main 3d mesh split at the fracture location. """ - added_points: Set[ int ] = set() + added_points: set[ int ] = set() for node_mapping in cell_to_node_mapping.values(): for i, o in node_mapping.items(): if i != o: @@ -336,16 +418,16 @@ def __perform_split( old_mesh: vtkUnstructuredGrid, new_mesh.Allocate( old_mesh.GetNumberOfCells() ) for c in tqdm( range( old_mesh.GetNumberOfCells() ), desc="Performing the mesh split" ): - node_mapping: Mapping[ int, int ] = cell_to_node_mapping.get( c, {} ) + node_mapping: IDMapping = cell_to_node_mapping.get( c, {} ) cell: vtkCell = old_mesh.GetCell( c ) cell_type: int = cell.GetCellType() # For polyhedron, we'll manipulate the face stream directly. if cell_type == VTK_POLYHEDRON: face_stream = vtkIdList() old_mesh.GetFaceStream( c, face_stream ) - new_face_nodes: List[ List[ int ] ] = [] + new_face_nodes: list[ list[ int ] ] = list() for face_nodes in FaceStream.build_from_vtk_id_list( face_stream ).face_nodes: - new_point_ids = [] + new_point_ids = list() for current_point_id in face_nodes: new_point_id: int = node_mapping.get( current_point_id, current_point_id ) new_point_ids.append( new_point_id ) @@ -361,13 +443,13 @@ def __perform_split( old_mesh: vtkUnstructuredGrid, cell_point_ids.SetId( i, new_point_id ) new_mesh.InsertNextCell( cell_type, cell_point_ids ) - __copy_fields( old_mesh, new_mesh, collocated_nodes ) + __copy_fields( old_mesh, new_mesh ) return new_mesh -def __generate_fracture_mesh( mesh_points: vtkPoints, fracture_info: FractureInfo, - cell_to_node_mapping: Mapping[ int, Mapping[ int, int ] ] ) -> vtkUnstructuredGrid: +def __generate_fracture_mesh( old_mesh: vtkUnstructuredGrid, fracture_info: FractureInfo, + cell_to_node_mapping: Mapping[ int, IDMapping ] ) -> vtkUnstructuredGrid: """ Generates the mesh of the fracture. :param mesh_points: The points of the main 3d mesh. @@ -377,6 +459,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(): @@ -386,8 +469,8 @@ def __generate_fracture_mesh( mesh_points: vtkPoints, fracture_info: FractureInf # Some elements can have all their nodes not duplicated. # In this case, it's mandatory not get rid of this element # because the neighboring 3d elements won't follow. - face_nodes: List[ Collection[ int ] ] = [] - discarded_face_nodes: Set[ Iterable[ int ] ] = set() + face_nodes: list[ Collection[ int ] ] = list() + discarded_face_nodes: set[ Iterable[ int ] ] = set() for ns in fracture_info.face_nodes: if any( map( is_node_duplicated.__getitem__, ns ) ): face_nodes.append( ns ) @@ -395,15 +478,16 @@ def __generate_fracture_mesh( mesh_points: vtkPoints, fracture_info: FractureInf discarded_face_nodes.add( ns ) if discarded_face_nodes: - # tmp = [] + # tmp = list() # 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.") + logging.info( 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: @@ -413,9 +497,9 @@ def __generate_fracture_mesh( mesh_points: vtkPoints, fracture_info: FractureInf num_points: int = len( fracture_nodes ) points = vtkPoints() points.SetNumberOfPoints( num_points ) - node_3d_to_node_2d: Dict[ int, int ] = {} # Building the node mapping, from 3d mesh nodes to 2d fracture nodes. + node_3d_to_node_2d: IDMapping = {} # Building the node mapping, from 3d mesh nodes to 2d fracture nodes. for i, n in enumerate( fracture_nodes ): - coords: Tuple[ float, float, float ] = mesh_points.GetPoint( n ) + coords: Coordinates3D = mesh_points.GetPoint( n ) points.SetPoint( i, coords ) node_3d_to_node_2d[ n ] = i @@ -427,7 +511,7 @@ def __generate_fracture_mesh( mesh_points: vtkPoints, fracture_info: FractureInf polygon.GetPointIds().SetId( i, node_3d_to_node_2d[ n ] ) polygons.InsertNextCell( polygon ) - buckets: Dict[ int, Set[ int ] ] = defaultdict( set ) + buckets: dict[ int, set[ int ] ] = defaultdict( set ) for node_mapping in cell_to_node_mapping.values(): for i, o in node_mapping.items(): k: Optional[ int ] = node_3d_to_node_2d.get( min( i, o ) ) @@ -448,31 +532,34 @@ 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, True ) + return fracture_mesh def __split_mesh_on_fracture( mesh: vtkUnstructuredGrid, - options: Options ) -> Tuple[ vtkUnstructuredGrid, vtkUnstructuredGrid ]: + options: Options ) -> tuple[ vtkUnstructuredGrid, vtkUnstructuredGrid ]: fracture: FractureInfo = build_fracture_info( mesh, options ) cell_to_cell: networkx.Graph = build_cell_to_cell_graph( mesh, fracture ) - cell_to_node_mapping: Mapping[ int, Mapping[ int, int ] ] = __identify_split( mesh.GetNumberOfPoints(), - cell_to_cell, fracture.node_to_cells ) + cell_to_node_mapping: Mapping[ int, IDMapping ] = __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 def __check( mesh, options: Options ) -> Result: output_mesh, fracture_mesh = __split_mesh_on_fracture( mesh, options ) - vtk_utils.write_mesh( output_mesh, options.vtk_output ) - vtk_utils.write_mesh( fracture_mesh, options.vtk_fracture_output ) + write_mesh( output_mesh, options.vtk_output ) + write_mesh( fracture_mesh, options.vtk_fracture_output ) # TODO provide statistics about what was actually performed (size of the fracture, number of split nodes...). return Result( info="OK" ) def check( vtk_input_file: str, options: Options ) -> Result: try: - mesh = vtk_utils.read_mesh( vtk_input_file ) + mesh = read_mesh( vtk_input_file ) return __check( mesh, options ) except BaseException as e: logging.error( e ) diff --git a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py index 9beb375..3b581b7 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py @@ -1,25 +1,11 @@ -from dataclasses import dataclass import os.path import logging -import sys -from typing import ( - Any, - Iterator, - Optional, -) - -from vtkmodules.vtkCommonCore import ( - vtkIdList, ) -from vtkmodules.vtkCommonDataModel import ( - vtkUnstructuredGrid, ) -from vtkmodules.vtkIOLegacy import ( - vtkUnstructuredGridWriter, - vtkUnstructuredGridReader, -) -from vtkmodules.vtkIOXML import ( - vtkXMLUnstructuredGridReader, - vtkXMLUnstructuredGridWriter, -) +from dataclasses import dataclass +from typing import Iterator, Optional +from vtkmodules.vtkCommonCore import vtkIdList +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from vtkmodules.vtkIOLegacy import vtkUnstructuredGridWriter, vtkUnstructuredGridReader +from vtkmodules.vtkIOXML import vtkXMLUnstructuredGridReader, vtkXMLUnstructuredGridWriter @dataclass( frozen=True ) @@ -36,7 +22,7 @@ def to_vtk_id_list( data ) -> vtkIdList: return result -def vtk_iter( l ) -> Iterator[ Any ]: +def vtk_iter( l ) -> Iterator[ any ]: """ Utility function transforming a vtk "container" (e.g. vtkIdList) into an iterable to be used for building built-ins python containers. :param l: A vtk container. @@ -50,6 +36,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 +101,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. """ + if not os.path.exists( vtk_input_file ): + err_msg: str = f"Invalid file path. Could not read \"{vtk_input_file}\". Dying..." + logging.error( err_msg ) + raise ValueError( err_msg ) 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 @@ -96,8 +118,9 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid: if output_mesh: 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 ) + err_msg = f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\". Dying..." + logging.error( err_msg ) + raise ValueError( err_msg ) def __write_vtk( mesh: vtkUnstructuredGrid, output: str ) -> int: @@ -135,6 +158,7 @@ def write_mesh( mesh: vtkUnstructuredGrid, vtk_output: VtkOutput ) -> int: 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 ) + err_msg = f"Could not find the appropriate VTK writer for extension \"{file_extension}\". Dying..." + logging.error( err_msg ) + raise ValueError( err_msg ) return 0 if success_code else 2 # the Write member function return 1 in case of success, 0 otherwise. diff --git a/geos-mesh/tests/test_generate_fractures.py b/geos-mesh/tests/test_generate_fractures.py index 077d9d7..30e0e7d 100644 --- a/geos-mesh/tests/test_generate_fractures.py +++ b/geos-mesh/tests/test_generate_fractures.py @@ -1,31 +1,20 @@ -from dataclasses import dataclass - -from typing import ( - Tuple, - Iterable, - Iterator, - Sequence, -) - +import logging import numpy - import pytest - -from vtkmodules.vtkCommonDataModel import ( - vtkUnstructuredGrid, - VTK_HEXAHEDRON, - VTK_POLYHEDRON, - VTK_QUAD, -) -from vtkmodules.util.numpy_support import ( - numpy_to_vtk, ) - -from checks.vtk_utils import ( - to_vtk_id_list, ) - -from checks.check_fractures import format_collocated_nodes -from checks.generate_cube import build_rectilinear_blocks_mesh, XYZ -from checks.generate_fractures import __split_mesh_on_fracture, Options, FracturePolicy +from dataclasses import dataclass +from typing import Iterable, Iterator, Sequence, TypeAlias +from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkQuad, VTK_HEXAHEDRON, VTK_POLYHEDRON, VTK_QUAD ) +from vtkmodules.util.numpy_support import numpy_to_vtk, vtk_to_numpy +from geos.mesh.doctor.checks.vtk_utils import to_vtk_id_list +from geos.mesh.doctor.checks.check_fractures import format_collocated_nodes +from geos.mesh.doctor.checks.generate_cube import build_rectilinear_blocks_mesh, XYZ +from geos.mesh.doctor.checks.generate_fractures import ( __split_mesh_on_fracture, Options, FracturePolicy, + Coordinates3D, IDMapping ) +""" +TypeAliases used in this file +""" +FaceNodesCoords: TypeAlias = tuple[ tuple[ float ] ] +IDMatrix: TypeAlias = Sequence[ Sequence[ int ] ] @dataclass( frozen=True ) @@ -42,11 +31,11 @@ class TestCase: __test__ = False input_mesh: vtkUnstructuredGrid options: Options - collocated_nodes: Sequence[ Sequence[ int ] ] + collocated_nodes: IDMatrix result: TestResult -def __build_test_case( xs: Tuple[ numpy.ndarray, numpy.ndarray, numpy.ndarray ], +def __build_test_case( xs: tuple[ numpy.ndarray, numpy.ndarray, numpy.ndarray ], attribute: Iterable[ int ], field_values: Iterable[ int ] = None, policy: FracturePolicy = FracturePolicy.FIELD ): @@ -95,20 +84,10 @@ def __generate_test_data() -> Iterator[ TestCase ]: # Split in 3 inc = Incrementor( 27 ) - collocated_nodes: Sequence[ Sequence[ int ] ] = ( - ( 1, *inc.next( 1 ) ), - ( 3, *inc.next( 1 ) ), - ( 4, *inc.next( 2 ) ), - ( 7, *inc.next( 1 ) ), - ( 1 + 9, *inc.next( 1 ) ), - ( 3 + 9, *inc.next( 1 ) ), - ( 4 + 9, *inc.next( 2 ) ), - ( 7 + 9, *inc.next( 1 ) ), - ( 1 + 18, *inc.next( 1 ) ), - ( 3 + 18, *inc.next( 1 ) ), - ( 4 + 18, *inc.next( 2 ) ), - ( 7 + 18, *inc.next( 1 ) ), - ) + collocated_nodes: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 3, *inc.next( 1 ) ), ( 4, *inc.next( 2 ) ), + ( 7, *inc.next( 1 ) ), ( 1 + 9, *inc.next( 1 ) ), ( 3 + 9, *inc.next( 1 ) ), + ( 4 + 9, *inc.next( 2 ) ), ( 7 + 9, *inc.next( 1 ) ), ( 1 + 18, *inc.next( 1 ) ), + ( 3 + 18, *inc.next( 1 ) ), ( 4 + 18, *inc.next( 2 ) ), ( 7 + 18, *inc.next( 1 ) ) ) mesh, options = __build_test_case( ( three_nodes, three_nodes, three_nodes ), ( 0, 1, 2, 1, 0, 1, 2, 1 ) ) yield TestCase( input_mesh=mesh, options=options, @@ -117,27 +96,13 @@ def __generate_test_data() -> Iterator[ TestCase ]: # Split in 8 inc = Incrementor( 27 ) - collocated_nodes: Sequence[ Sequence[ int ] ] = ( - ( 1, *inc.next( 1 ) ), - ( 3, *inc.next( 1 ) ), - ( 4, *inc.next( 3 ) ), - ( 5, *inc.next( 1 ) ), - ( 7, *inc.next( 1 ) ), - ( 0 + 9, *inc.next( 1 ) ), - ( 1 + 9, *inc.next( 3 ) ), - ( 2 + 9, *inc.next( 1 ) ), - ( 3 + 9, *inc.next( 3 ) ), - ( 4 + 9, *inc.next( 7 ) ), - ( 5 + 9, *inc.next( 3 ) ), - ( 6 + 9, *inc.next( 1 ) ), - ( 7 + 9, *inc.next( 3 ) ), - ( 8 + 9, *inc.next( 1 ) ), - ( 1 + 18, *inc.next( 1 ) ), - ( 3 + 18, *inc.next( 1 ) ), - ( 4 + 18, *inc.next( 3 ) ), - ( 5 + 18, *inc.next( 1 ) ), - ( 7 + 18, *inc.next( 1 ) ), - ) + collocated_nodes: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 3, *inc.next( 1 ) ), ( 4, *inc.next( 3 ) ), + ( 5, *inc.next( 1 ) ), ( 7, *inc.next( 1 ) ), ( 0 + 9, *inc.next( 1 ) ), + ( 1 + 9, *inc.next( 3 ) ), ( 2 + 9, *inc.next( 1 ) ), ( 3 + 9, *inc.next( 3 ) ), + ( 4 + 9, *inc.next( 7 ) ), ( 5 + 9, *inc.next( 3 ) ), ( 6 + 9, *inc.next( 1 ) ), + ( 7 + 9, *inc.next( 3 ) ), ( 8 + 9, *inc.next( 1 ) ), ( 1 + 18, *inc.next( 1 ) ), + ( 3 + 18, *inc.next( 1 ) ), ( 4 + 18, *inc.next( 3 ) ), ( 5 + 18, *inc.next( 1 ) ), + ( 7 + 18, *inc.next( 1 ) ) ) mesh, options = __build_test_case( ( three_nodes, three_nodes, three_nodes ), range( 8 ) ) yield TestCase( input_mesh=mesh, options=options, @@ -146,14 +111,8 @@ def __generate_test_data() -> Iterator[ TestCase ]: # Straight notch inc = Incrementor( 27 ) - collocated_nodes: Sequence[ Sequence[ int ] ] = ( - ( 1, *inc.next( 1 ) ), - ( 4, ), - ( 1 + 9, *inc.next( 1 ) ), - ( 4 + 9, ), - ( 1 + 18, *inc.next( 1 ) ), - ( 4 + 18, ), - ) + collocated_nodes: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 4, ), ( 1 + 9, *inc.next( 1 ) ), ( 4 + 9, ), + ( 1 + 18, *inc.next( 1 ) ), ( 4 + 18, ) ) mesh, options = __build_test_case( ( three_nodes, three_nodes, three_nodes ), ( 0, 1, 2, 2, 0, 1, 2, 2 ), field_values=( 0, 1 ) ) yield TestCase( input_mesh=mesh, @@ -163,16 +122,8 @@ def __generate_test_data() -> Iterator[ TestCase ]: # L-shaped notch inc = Incrementor( 27 ) - collocated_nodes: Sequence[ Sequence[ int ] ] = ( - ( 1, *inc.next( 1 ) ), - ( 4, *inc.next( 1 ) ), - ( 7, *inc.next( 1 ) ), - ( 1 + 9, *inc.next( 1 ) ), - ( 4 + 9, ), - ( 7 + 9, ), - ( 1 + 18, *inc.next( 1 ) ), - ( 4 + 18, ), - ) + collocated_nodes: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 4, *inc.next( 1 ) ), ( 7, *inc.next( 1 ) ), + ( 1 + 9, *inc.next( 1 ) ), ( 4 + 9, ), ( 7 + 9, ), ( 19, *inc.next( 1 ) ), ( 22, ) ) mesh, options = __build_test_case( ( three_nodes, three_nodes, three_nodes ), ( 0, 1, 0, 1, 0, 1, 2, 2 ), field_values=( 0, 1 ) ) yield TestCase( input_mesh=mesh, @@ -182,16 +133,9 @@ def __generate_test_data() -> Iterator[ TestCase ]: # 3x1x1 split inc = Incrementor( 2 * 2 * 4 ) - collocated_nodes: Sequence[ Sequence[ int ] ] = ( - ( 1, *inc.next( 1 ) ), - ( 2, *inc.next( 1 ) ), - ( 5, *inc.next( 1 ) ), - ( 6, *inc.next( 1 ) ), - ( 1 + 8, *inc.next( 1 ) ), - ( 2 + 8, *inc.next( 1 ) ), - ( 5 + 8, *inc.next( 1 ) ), - ( 6 + 8, *inc.next( 1 ) ), - ) + collocated_nodes: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 2, *inc.next( 1 ) ), ( 5, *inc.next( 1 ) ), + ( 6, *inc.next( 1 ) ), ( 1 + 8, *inc.next( 1 ) ), ( 2 + 8, *inc.next( 1 ) ), + ( 5 + 8, *inc.next( 1 ) ), ( 6 + 8, *inc.next( 1 ) ) ) mesh, options = __build_test_case( ( four_nodes, two_nodes, two_nodes ), ( 0, 1, 2 ) ) yield TestCase( input_mesh=mesh, options=options, @@ -199,12 +143,8 @@ def __generate_test_data() -> Iterator[ TestCase ]: result=TestResult( 6 * 4, 3, 2 * 4, 2 ) ) # Discarded fracture element if no node duplication. - collocated_nodes: Sequence[ Sequence[ int ] ] = () - mesh, options = __build_test_case( ( three_nodes, four_nodes, four_nodes ), [ - 0, - ] * 8 + [ 1, 2 ] + [ - 0, - ] * 8, + collocated_nodes: IDMatrix = tuple() + mesh, options = __build_test_case( ( three_nodes, four_nodes, four_nodes ), ( 0, ) * 8 + ( 1, 2 ) + ( 0, ) * 8, field_values=( 1, 2 ) ) yield TestCase( input_mesh=mesh, options=options, @@ -213,20 +153,11 @@ def __generate_test_data() -> Iterator[ TestCase ]: # Fracture on a corner inc = Incrementor( 3 * 4 * 4 ) - collocated_nodes: Sequence[ Sequence[ int ] ] = ( - ( 1 + 12, ), - ( 4 + 12, ), - ( 7 + 12, ), - ( 1 + 12 * 2, *inc.next( 1 ) ), - ( 4 + 12 * 2, *inc.next( 1 ) ), - ( 7 + 12 * 2, ), - ( 1 + 12 * 3, *inc.next( 1 ) ), - ( 4 + 12 * 3, *inc.next( 1 ) ), - ( 7 + 12 * 3, ), - ) - mesh, options = __build_test_case( ( three_nodes, four_nodes, four_nodes ), [ - 0, - ] * 6 + [ 1, 2, 1, 2, 0, 0, 1, 2, 1, 2, 0, 0 ], + collocated_nodes: IDMatrix = ( ( 1 + 12, ), ( 4 + 12, ), ( 7 + 12, ), ( 1 + 12 * 2, *inc.next( 1 ) ), + ( 4 + 12 * 2, *inc.next( 1 ) ), ( 7 + 12 * 2, ), ( 1 + 12 * 3, *inc.next( 1 ) ), + ( 4 + 12 * 3, *inc.next( 1 ) ), ( 7 + 12 * 3, ) ) + mesh, options = __build_test_case( ( three_nodes, four_nodes, four_nodes ), + ( 0, ) * 6 + ( 1, 2, 1, 2, 0, 0, 1, 2, 1, 2, 0, 0 ), field_values=( 1, 2 ) ) yield TestCase( input_mesh=mesh, options=options, @@ -235,12 +166,8 @@ def __generate_test_data() -> Iterator[ TestCase ]: # Generate mesh with 2 hexs, one being a standard hex, the other a 42 hex. inc = Incrementor( 3 * 2 * 2 ) - collocated_nodes: Sequence[ Sequence[ int ] ] = ( - ( 1, *inc.next( 1 ) ), - ( 1 + 3, *inc.next( 1 ) ), - ( 1 + 6, *inc.next( 1 ) ), - ( 1 + 9, *inc.next( 1 ) ), - ) + collocated_nodes: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 1 + 3, *inc.next( 1 ) ), ( 1 + 6, *inc.next( 1 ) ), + ( 1 + 9, *inc.next( 1 ) ) ) mesh, options = __build_test_case( ( three_nodes, two_nodes, two_nodes ), ( 0, 1 ) ) polyhedron_mesh = vtkUnstructuredGrid() polyhedron_mesh.SetPoints( mesh.GetPoints() ) @@ -258,12 +185,8 @@ def __generate_test_data() -> Iterator[ TestCase ]: # Split in 2 using the internal fracture description inc = Incrementor( 3 * 2 * 2 ) - collocated_nodes: Sequence[ Sequence[ int ] ] = ( - ( 1, *inc.next( 1 ) ), - ( 1 + 3, *inc.next( 1 ) ), - ( 1 + 6, *inc.next( 1 ) ), - ( 1 + 9, *inc.next( 1 ) ), - ) + collocated_nodes: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 1 + 3, *inc.next( 1 ) ), ( 1 + 6, *inc.next( 1 ) ), + ( 1 + 9, *inc.next( 1 ) ) ) mesh, options = __build_test_case( ( three_nodes, two_nodes, two_nodes ), attribute=( 0, 0, 0 ), field_values=( 0, ), @@ -286,3 +209,146 @@ def test_generate_fracture( test_case: TestCase ): res = format_collocated_nodes( fracture_mesh ) assert res == test_case.collocated_nodes assert len( res ) == test_case.result.fracture_mesh_num_points + + +def add_simplified_field_for_cells( mesh: vtkUnstructuredGrid, field_name: str, field_dimension: int ): + """Reduce functionality obtained from src.geos.mesh.doctor.checks.generate_fracture.__add_fields + where the goal is to add a cell data array with incrementing values. + + Args: + mesh (vtkUnstructuredGrid): Unstructured mesh. + field_name (str): Name of the field to add to CellData + field_dimension (int): Number of components for the field. + """ + data = mesh.GetCellData() + n = mesh.GetNumberOfCells() + array = numpy.ones( ( n, field_dimension ), dtype=float ) + array = numpy.arange( 1, n * field_dimension + 1 ).reshape( n, field_dimension ) + vtk_array = numpy_to_vtk( array ) + vtk_array.SetName( field_name ) + data.AddArray( vtk_array ) + + +def find_borders_faces_rectilinear_grid( mesh: vtkUnstructuredGrid ) -> tuple[ FaceNodesCoords ]: + """ + 6+--------+7 + / /| + / / | + 4+--------+5 | + | | | + | 2+ | +3 + | | / + | |/ + 0+--------+1 + + For a vtk rectilinear grid, gives the coordinates of each of its borders face nodes. + + Args: + mesh (vtkUnstructuredGrid): Unstructured mesh. + + Returns: + tuple[QuadCoords]: For a rectilinear grid, returns a tuple of 6 elements. + """ + mesh_bounds: tuple[ float ] = mesh.GetBounds() + min_bound: Coordinates3D = [ mesh_bounds[ i ] for i in range( len( mesh_bounds ) ) if i % 2 == 0 ] + max_bound: Coordinates3D = [ mesh_bounds[ i ] for i in range( len( mesh_bounds ) ) if i % 2 == 1 ] + center: Coordinates3D = mesh.GetCenter() + face_diag: tuple[ float ] = ( ( max_bound[ 0 ] - min_bound[ 0 ] ) / 2, ( max_bound[ 1 ] - min_bound[ 1 ] ) / 2, + ( max_bound[ 2 ] - min_bound[ 2 ] ) / 2 ) + node0: Coordinates3D = ( center[ 0 ] - face_diag[ 0 ], center[ 1 ] - face_diag[ 1 ], center[ 2 ] - face_diag[ 2 ] ) + node1: Coordinates3D = ( center[ 0 ] + face_diag[ 0 ], center[ 1 ] - face_diag[ 1 ], center[ 2 ] - face_diag[ 2 ] ) + node2: Coordinates3D = ( center[ 0 ] - face_diag[ 0 ], center[ 1 ] + face_diag[ 1 ], center[ 2 ] - face_diag[ 2 ] ) + node3: Coordinates3D = ( center[ 0 ] + face_diag[ 0 ], center[ 1 ] + face_diag[ 1 ], center[ 2 ] - face_diag[ 2 ] ) + node4: Coordinates3D = ( center[ 0 ] - face_diag[ 0 ], center[ 1 ] - face_diag[ 1 ], center[ 2 ] + face_diag[ 2 ] ) + node5: Coordinates3D = ( center[ 0 ] + face_diag[ 0 ], center[ 1 ] - face_diag[ 1 ], center[ 2 ] + face_diag[ 2 ] ) + node6: Coordinates3D = ( center[ 0 ] - face_diag[ 0 ], center[ 1 ] + face_diag[ 1 ], center[ 2 ] + face_diag[ 2 ] ) + node7: Coordinates3D = ( center[ 0 ] + face_diag[ 0 ], center[ 1 ] + face_diag[ 1 ], center[ 2 ] + face_diag[ 2 ] ) + faces: tuple[ FaceNodesCoords ] = ( ( node0, node1, node3, node2 ), ( node4, node5, node7, node6 ), + ( node0, node2, node6, node4 ), ( node1, node3, node7, node5 ), + ( node0, node1, node5, node4 ), ( node2, node3, node7, node6 ) ) + return faces + + +def add_quad( mesh: vtkUnstructuredGrid, face: FaceNodesCoords ): + """Adds a quad cell to each border of an unstructured mesh. + + Args: + mesh (vtkUnstructuredGrid): Unstructured mesh. + """ + points_coords = mesh.GetPoints().GetData() + quad: vtkQuad = vtkQuad() + ids_association: IDMapping = {} + for i in range( mesh.GetNumberOfPoints() ): + for j in range( len( face ) ): + if points_coords.GetTuple( i ) == face[ j ]: + ids_association[ i ] = j + break + if len( ids_association ) == 4: + break + + for vertice_id, quad_coord_index in ids_association.items(): + quad.GetPoints().InsertNextPoint( face[ quad_coord_index ] ) + quad.GetPointIds().SetId( quad_coord_index, vertice_id ) + + mesh.InsertNextCell( quad.GetCellType(), quad.GetPointIds() ) + + +def test_copy_fields_when_splitting_mesh(): + """This test is designed to check the __copy_fields method from generate_fractures, + that will be called when using __split_mesh_on_fracture method from generate_fractures. + """ + # Generating the rectilinear grid and its quads on all borders + x: numpy.array = numpy.array( [ 0, 1, 2 ] ) + y: numpy.array = numpy.array( [ 0, 1 ] ) + z: numpy.array = numpy.array( [ 0, 1 ] ) + xyzs: XYZ = XYZ( x, y, z ) + mesh: vtkUnstructuredGrid = build_rectilinear_blocks_mesh( [ xyzs ] ) + assert mesh.GetCells().GetNumberOfCells() == 2 + border_faces: tuple[ FaceNodesCoords ] = find_borders_faces_rectilinear_grid( mesh ) + for face in border_faces: + add_quad( mesh, face ) + assert mesh.GetCells().GetNumberOfCells() == 8 + # Create a quad cell to represent the fracture surface. + fracture: FaceNodesCoords = ( ( 1.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 1.0, 1.0, 1.0 ), ( 1.0, 0.0, 1.0 ) ) + add_quad( mesh, fracture ) + assert mesh.GetCells().GetNumberOfCells() == 9 + # Add a "TestField" array + assert mesh.GetCellData().GetNumberOfArrays() == 0 + add_simplified_field_for_cells( mesh, "TestField", 1 ) + assert mesh.GetCellData().GetNumberOfArrays() == 1 + assert mesh.GetCellData().GetArrayName( 0 ) == "TestField" + testField_values: list[ int ] = vtk_to_numpy( mesh.GetCellData().GetArray( 0 ) ).tolist() + assert testField_values == [ 1, 2, 3, 4, 5, 6, 7, 8, 9 ] + # Split the mesh along the fracture surface which is number 9 on TestField + options = Options( policy=FracturePolicy.INTERNAL_SURFACES, + field="TestField", + field_values=frozenset( map( int, [ "9" ] ) ), + vtk_output=None, + vtk_fracture_output=None ) + main_mesh, fracture_mesh = __split_mesh_on_fracture( mesh, options ) + assert main_mesh.GetCellData().GetNumberOfArrays() == 1 + assert fracture_mesh.GetCellData().GetNumberOfArrays() == 1 + assert main_mesh.GetCellData().GetArrayName( 0 ) == "TestField" + assert fracture_mesh.GetCellData().GetArrayName( 0 ) == "TestField" + # Make sure that only 1 correct value is in "TestField" array for fracture_mesh, 9 values for main_mesh + fracture_mesh_values: list[ int ] = vtk_to_numpy( fracture_mesh.GetCellData().GetArray( 0 ) ).tolist() + main_mesh_values: list[ int ] = vtk_to_numpy( main_mesh.GetCellData().GetArray( 0 ) ).tolist() + assert fracture_mesh_values == [ 9 ] # The value for the fracture surface + assert main_mesh_values == [ 1, 2, 3, 4, 5, 6, 7, 8, 9 ] + # Test for invalid point field name + add_simplified_field_for_cells( mesh, "GLOBAL_IDS_POINTS", 1 ) + with pytest.raises( ValueError ) as pytest_wrapped_e: + main_mesh, fracture_mesh = __split_mesh_on_fracture( mesh, options ) + assert pytest_wrapped_e.type == ValueError + # Test for invalid cell field name + mesh: vtkUnstructuredGrid = build_rectilinear_blocks_mesh( [ xyzs ] ) + border_faces: tuple[ FaceNodesCoords ] = find_borders_faces_rectilinear_grid( mesh ) + for face in border_faces: + add_quad( mesh, face ) + add_quad( mesh, fracture ) + add_simplified_field_for_cells( mesh, "TestField", 1 ) + add_simplified_field_for_cells( mesh, "GLOBAL_IDS_CELLS", 1 ) + assert mesh.GetCellData().GetNumberOfArrays() == 2 + with pytest.raises( ValueError ) as pytest_wrapped_e: + main_mesh, fracture_mesh = __split_mesh_on_fracture( mesh, options ) + assert pytest_wrapped_e.type == ValueError