diff --git a/docs/requirements.txt b/docs/requirements.txt index 3e8d02d..af812eb 100644 --- a/docs/requirements.txt +++ b/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 diff --git a/geos-ats/src/geos/ats/helpers/curve_check.py b/geos-ats/src/geos/ats/helpers/curve_check.py index a35b27f..31250c7 100644 --- a/geos-ats/src/geos/ats/helpers/curve_check.py +++ b/geos-ats/src/geos/ats/helpers/curve_check.py @@ -33,16 +33,18 @@ def interpolate_values_time( ta, xa, tb ): """ N = list( np.shape( xa ) ) M = len( tb ) + ta = np.ascontiguousarray( np.squeeze( ta ) ) + tb = np.ascontiguousarray( np.squeeze( tb ) ) if ( len( N ) == 1 ): return interp1d( ta, xa )( tb ) else: # Reshape the input array so that we can work on the non-time axes - S = np.product( N[ 1: ] ) + S = np.prod( N[ 1: ] ) xc = np.reshape( xa, ( N[ 0 ], S ) ) xd = np.zeros( ( len( tb ), S ) ) for ii in range( S ): - xd[ :, ii ] = interp1d( ta, xc[ :, ii ] )( tb ) + xd[ :, ii ] = interp1d( ta, xc[ :, ii ], bounds_error=False, fill_value='extrapolate' )( tb ) # Return the array to it's expected shape N[ 0 ] = M diff --git a/geos-ats/src/geos/ats/helpers/restart_check.py b/geos-ats/src/geos/ats/helpers/restart_check.py index 7aa282e..f491107 100644 --- a/geos-ats/src/geos/ats/helpers/restart_check.py +++ b/geos-ats/src/geos/ats/helpers/restart_check.py @@ -7,6 +7,7 @@ import argparse import logging import time +import string from pathlib import Path try: from geos.ats.helpers.permute_array import permuteArray # type: ignore[import] @@ -375,35 +376,100 @@ def compareIntArrays( self, path, arr, base_arr ): ARR [in]: The hdf5 Dataset to compare. BASE_ARR [in]: The hdf5 Dataset to compare against. """ - # If the shapes are different they can't be compared. + message = "" if arr.shape != base_arr.shape: - msg = "Datasets have different shapes and therefore can't be compared: %s, %s.\n" % ( arr.shape, - base_arr.shape ) - self.errorMsg( path, msg, True ) - return + message = "Datasets have different shapes and therefore can't be compared statistically: %s, %s.\n" % ( + arr.shape, base_arr.shape ) + else: + # Calculate the absolute difference. + difference = np.subtract( arr, base_arr ) + np.abs( difference, out=difference ) - # Create a copy of the arrays. + offenders = difference != 0.0 + n_offenders = np.sum( offenders ) - # Calculate the absolute difference. - difference = np.subtract( arr, base_arr ) - np.abs( difference, out=difference ) + if n_offenders != 0: + max_index = np.unravel_index( np.argmax( difference ), difference.shape ) + max_difference = difference[ max_index ] + offenders_mean = np.mean( difference[ offenders ] ) + offenders_std = np.std( difference[ offenders ] ) - offenders = difference != 0.0 - n_offenders = np.sum( offenders ) + message = "Arrays of types %s and %s have %s values of which %d have differing values.\n" % ( + arr.dtype, base_arr.dtype, offenders.size, n_offenders ) + message += "Statistics of the differences greater than 0:\n" + message += "\tmax_index = %s, max = %s, mean = %s, std = %s\n" % ( max_index, max_difference, + offenders_mean, offenders_std ) - if n_offenders != 0: - max_index = np.unravel_index( np.argmax( difference ), difference.shape ) - max_difference = difference[ max_index ] - offenders_mean = np.mean( difference[ offenders ] ) - offenders_std = np.std( difference[ offenders ] ) + # actually, int8 arrays are almost always char arrays, so we sould add a character comparison. + if arr.dtype == np.int8 and base_arr.dtype == np.int8: + message += self.compareCharArrays( arr, base_arr ) - message = "Arrays of types %s and %s have %s values of which %d have differing values.\n" % ( - arr.dtype, base_arr.dtype, offenders.size, n_offenders ) - message += "Statistics of the differences greater than 0:\n" - message += "\tmax_index = %s, max = %s, mean = %s, std = %s\n" % ( max_index, max_difference, - offenders_mean, offenders_std ) + if message != "": self.errorMsg( path, message, True ) + def compareCharArrays( self, comp_arr, base_arr ): + """ + Compare the valid characters of two arrays and return a formatted string showing differences. + + COMP_ARR [in]: The hdf5 Dataset to compare. + BASE_ARR [in]: The hdf5 Dataset to compare against. + + Returns a formatted string highlighting the differing characters. + """ + comp_ndarr = np.array( comp_arr ).flatten() + base_ndarr = np.array( base_arr ).flatten() + + # Replace invalid characters by group-separator characters ('\x1D') + valid_chars = set( string.printable ) + invalid_char = '\x1D' + comp_str = "".join( + [ chr( x ) if ( x >= 0 and chr( x ) in valid_chars ) else invalid_char for x in comp_ndarr ] ) + base_str = "".join( + [ chr( x ) if ( x >= 0 and chr( x ) in valid_chars ) else invalid_char for x in base_ndarr ] ) + + # replace whitespaces sequences by only one space (preventing indentation / spacing changes detection) + whitespace_pattern = r"[ \t\n\r\v\f]+" + comp_str = re.sub( whitespace_pattern, " ", comp_str ) + base_str = re.sub( whitespace_pattern, " ", base_str ) + # replace invalid characters sequences by a double space (for clear display) + invalid_char_pattern = r"\x1D+" + comp_str_display = re.sub( invalid_char_pattern, " ", comp_str ) + base_str_display = re.sub( invalid_char_pattern, " ", base_str ) + + message = "" + + def limited_display( n, string ): + return string[ :n ] + f"... ({len(string)-n} omitted chars)" if len( string ) > n else string + + if len( comp_str ) != len( base_str ): + max_display = 250 + message = f"Character arrays have different sizes: {len( comp_str )}, {len( base_str )}.\n" + message += f" {limited_display( max_display, comp_str_display )}\n" + message += f" {limited_display( max_display, base_str_display )}\n" + else: + # We need to trim arrays to the length of the shortest one for the comparisons + min_length = min( len( comp_str_display ), len( base_str_display ) ) + comp_str_trim = comp_str_display[ :min_length ] + base_str_trim = base_str_display[ :min_length ] + + differing_indices = np.where( np.array( list( comp_str_trim ) ) != np.array( list( base_str_trim ) ) )[ 0 ] + if differing_indices.size != 0: + # check for reordering + arr_set = sorted( set( comp_str.split( invalid_char ) ) ) + base_arr_set = sorted( set( base_str.split( invalid_char ) ) ) + reordering_detected = arr_set == base_arr_set + + max_display = 110 if reordering_detected else 250 + message = "Differing valid characters" + message += " (substrings reordering detected):\n" if reordering_detected else ":\n" + + message += f" {limited_display( max_display, comp_str_display )}\n" + message += f" {limited_display( max_display, base_str_display )}\n" + message += " " + "".join( + [ "^" if i in differing_indices else " " for i in range( min( max_display, min_length ) ) ] ) + "\n" + + return message + def compareStringArrays( self, path, arr, base_arr ): """ Compare two string datasets. Exact equality is used as the acceptance criteria. diff --git a/geos-ats/src/geos/ats/test_case.py b/geos-ats/src/geos/ats/test_case.py index 9c3c27f..7affe6c 100644 --- a/geos-ats/src/geos/ats/test_case.py +++ b/geos-ats/src/geos/ats/test_case.py @@ -1,4 +1,4 @@ -import ats # type: ignore[import] +import ats # type: ignore[import] import os import shutil import logging @@ -319,7 +319,7 @@ def testCreate( self ): priority = 1 # Setup a new test group - atsTest = None + parent_test = None ats.tests.AtsTest.newGroup( priority=priority, path=self.path ) for stepnum, step in enumerate( self.steps ): np = getattr( step.p, "np", 1 ) @@ -329,10 +329,10 @@ def testCreate( self ): label = "%s_%d_%s" % ( self.name, stepnum + 1, step.label() ) # call either 'test' or 'testif' - if atsTest is None: + if parent_test is None: func = lambda *a, **k: test( *a, **k ) else: - func = lambda *a, **k: testif( atsTest, *a, **k ) + func = lambda *a, **k: testif( parent_test, *a, **k ) # Set the time limit kw = {} @@ -359,6 +359,11 @@ def testCreate( self ): if self.last_status == PASSED: atsTest.status = SKIPPED + # Set the parent test + # Note: do not make tests dependent on the result of curvecheck + if step.label() != 'curvecheck': + parent_test = atsTest + # End the group ats.tests.AtsTest.endGroup() diff --git a/geos-mesh/src/geos/mesh/conversion/main.py b/geos-mesh/src/geos/mesh/conversion/main.py index 90b048f..9579f32 100644 --- a/geos-mesh/src/geos/mesh/conversion/main.py +++ b/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() diff --git a/geos-mesh/src/geos/mesh/doctor/__init__.py b/geos-mesh/src/geos/mesh/doctor/__init__.py new file mode 100644 index 0000000..b1cfe26 --- /dev/null +++ b/geos-mesh/src/geos/mesh/doctor/__init__.py @@ -0,0 +1 @@ +# Empty \ No newline at end of file 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..a9a5d7e 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,27 @@ +import logging +import networkx 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 numpy import empty, ones, zeros 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 +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 cannot be used with Python 3.9. A simple assignment like described there will be used: +https://docs.python.org/3/library/typing.html#typing.TypeAlias:~:text=through%20simple%20assignment%3A-,Vector%20%3D%20list%5Bfloat%5D,-Or%20marked%20with +""" +IDMapping = Mapping[ int, int ] +CellsPointsCoords = dict[ int, list[ tuple[ float ] ] ] +Coordinates3D = tuple[ float ] class FracturePolicy( Enum ): @@ -56,9 +33,10 @@ class FracturePolicy( Enum ): class Options: policy: FracturePolicy field: str - field_values: FrozenSet[ int ] - vtk_output: VtkOutput - vtk_fracture_output: VtkOutput + field_values_combined: frozenset[ int ] + field_values_per_fracture: list[ frozenset[ int ] ] + mesh_VtkOutput: VtkOutput + all_fractures_VtkOutput: list[ VtkOutput ] @dataclass( frozen=True ) @@ -70,20 +48,21 @@ class Result: class FractureInfo: node_to_cells: Mapping[ int, Iterable[ int ] ] # For each _fracture_ node, gives all the cells that use this node. face_nodes: Iterable[ Collection[ int ] ] # For each fracture face, returns the nodes of this face + face_cell_id: Iterable[ int ] # For each fracture face, returns the corresponding id of the cell in the mesh 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 +71,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 +88,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,25 +100,28 @@ 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 ) + face_cell_id: list = list() # no cell of the mesh corresponds to that face when fracture policy is 'field' - return FractureInfo( node_to_cells=node_to_cells, face_nodes=face_nodes ) + return FractureInfo( node_to_cells=node_to_cells, face_nodes=face_nodes, face_cell_id=face_cell_id ) 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() + face_cell_id: list[ 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 ) ) + face_cell_id.append( cell_id ) for cell_id in tqdm( range( mesh.GetNumberOfCells() ), desc="Computing the node to cells mapping" ): cell = mesh.GetCell( cell_id ) @@ -148,12 +130,18 @@ def __build_fracture_info_from_internal_surfaces( mesh: vtkUnstructuredGrid, f: if cell.GetPointId( v ) in node_to_cells: node_to_cells[ cell.GetPointId( v ) ].append( cell_id ) - return FractureInfo( node_to_cells=node_to_cells, face_nodes=face_nodes ) + return FractureInfo( node_to_cells=node_to_cells, face_nodes=face_nodes, face_cell_id=face_cell_id ) -def build_fracture_info( mesh: vtkUnstructuredGrid, options: Options ) -> FractureInfo: +def build_fracture_info( mesh: vtkUnstructuredGrid, + options: Options, + combined_fractures: bool, + fracture_id: int = 0 ) -> FractureInfo: field = options.field - field_values = options.field_values + if combined_fractures: + field_values = options.field_values_combined + else: + field_values = options.field_values_per_fracture[ fracture_id ] cell_data = mesh.GetCellData() if cell_data.HasArray( field ): f = vtk_to_numpy( cell_data.GetArray( field ) ) @@ -177,24 +165,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 +196,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 +216,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 +227,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,65 +239,134 @@ def __call__( self, index: int ) -> int: return result -def __copy_fields( old_mesh: vtkUnstructuredGrid, new_mesh: vtkUnstructuredGrid, - collocated_nodes: Sequence[ int ] ) -> None: +def __copy_fields_splitted_mesh( old_mesh: vtkUnstructuredGrid, splitted_mesh: vtkUnstructuredGrid, + added_points_with_old_id: list[ tuple[ int ] ] ) -> 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. :return: None """ - # Copying the cell data. - # The cells are the same, just their nodes support have changed. + # Copying the cell data. The cells are the same, just their nodes support have changed. 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 ) + tmp = input_array.NewInstance() + tmp.DeepCopy( input_array ) + splitted_mesh.GetCellData().AddArray( input_array ) # 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 ) + tmp = input_array.NewInstance() + tmp.DeepCopy( input_array ) + splitted_mesh.GetFieldData().AddArray( input_array ) - # Copying the point data. + # Copying copy data. Need to take into account the new points. input_point_data = old_mesh.GetPointData() + new_number_points: int = splitted_mesh.GetNumberOfPoints() for i in range( input_point_data.GetNumberOfArrays() ): - input_array = input_point_data.GetArray( i ) - 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 ] ) ) - new_mesh.GetPointData().AddArray( tmp ) + old_points_array = vtk_to_numpy( input_point_data.GetArray( i ) ) + name: str = input_point_data.GetArrayName( i ) + logging.info( f"Copying point data \"{name}\"." ) + old_nrows: int = old_points_array.shape[ 0 ] + old_ncols: int = 1 if len( old_points_array.shape ) == 1 else old_points_array.shape[ 1 ] + # Reshape old_points_array if it is 1-dimensional + if len( old_points_array.shape ) == 1: + old_points_array = old_points_array.reshape( ( old_nrows, 1 ) ) + new_points_array = empty( ( new_number_points, old_ncols ) ) + new_points_array[ :old_nrows, : ] = old_points_array + for new_and_old_id in added_points_with_old_id: + new_points_array[ new_and_old_id[ 0 ], : ] = old_points_array[ new_and_old_id[ 1 ], : ] + # Reshape the VTK array to match the original dimensions + if old_ncols > 1: + vtk_array = numpy_to_vtk( new_points_array.flatten() ) + vtk_array.SetNumberOfComponents( old_ncols ) + vtk_array.SetNumberOfTuples( new_number_points ) + else: + vtk_array = numpy_to_vtk( new_points_array ) + vtk_array.SetName( name ) + splitted_mesh.GetPointData().AddArray( vtk_array ) + + +def __copy_fields_fracture_mesh( old_mesh: vtkUnstructuredGrid, fracture_mesh: vtkUnstructuredGrid, + face_cell_id: list[ int ], node_3d_to_node_2d: IDMapping ) -> None: + """ + Copies the fields from the old mesh to the new fracture when using internal_surfaces policy. + :param old_mesh: The mesh before the split. + :param fracture: The fracture mesh generated from the fracture_info. + :return: None + """ + # No copy of field data will be done with the fracture mesh because may lose its relevance compared to the splitted. + # Copying the cell data. The interesting cells are the ones stored in face_cell_id. + new_number_cells: int = fracture_mesh.GetNumberOfCells() + input_cell_data = old_mesh.GetCellData() + for i in range( input_cell_data.GetNumberOfArrays() ): + old_cells_array = vtk_to_numpy( input_cell_data.GetArray( i ) ) + old_nrows: int = old_cells_array.shape[ 0 ] + if len( old_cells_array.shape ) == 1: + old_cells_array = old_cells_array.reshape( ( old_nrows, 1 ) ) + name: str = input_cell_data.GetArrayName( i ) + logging.info( f"Copying cell data \"{name}\"." ) + new_array = old_cells_array[ face_cell_id, : ] + # Reshape the VTK array to match the original dimensions + old_ncols: int = 1 if len( old_cells_array.shape ) == 1 else old_cells_array.shape[ 1 ] + if old_ncols > 1: + vtk_array = numpy_to_vtk( new_array.flatten() ) + vtk_array.SetNumberOfComponents( old_ncols ) + vtk_array.SetNumberOfTuples( new_number_cells ) + else: + vtk_array = numpy_to_vtk( new_array ) + vtk_array.SetName( name ) + fracture_mesh.GetCellData().AddArray( vtk_array ) + new_number_points: int = fracture_mesh.GetNumberOfPoints() + input_point_data = old_mesh.GetPointData() + for i in range( input_point_data.GetNumberOfArrays() ): + old_points_array = vtk_to_numpy( input_point_data.GetArray( i ) ) + old_nrows = old_points_array.shape[ 0 ] + if len( old_points_array.shape ) == 1: + old_points_array = old_points_array.reshape( ( old_nrows, 1 ) ) + name = input_point_data.GetArrayName( i ) + logging.info( f"Copying point data \"{name}\"." ) + new_array = old_points_array[ list( node_3d_to_node_2d.keys() ), : ] + old_ncols = 1 if len( old_points_array.shape ) == 1 else old_points_array.shape[ 1 ] + if old_ncols > 1: + vtk_array = numpy_to_vtk( new_array.flatten() ) + vtk_array.SetNumberOfComponents( old_ncols ) + vtk_array.SetNumberOfTuples( new_number_points ) + else: + vtk_array = numpy_to_vtk( new_array ) + vtk_array.SetName( name ) + fracture_mesh.GetPointData().AddArray( vtk_array ) -def __perform_split( old_mesh: vtkUnstructuredGrid, - cell_to_node_mapping: Mapping[ int, Mapping[ int, int ] ] ) -> vtkUnstructuredGrid: + +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() + added_points_with_old_id: list[ tuple[ int ] ] = list() for node_mapping in cell_to_node_mapping.values(): for i, o in node_mapping.items(): if i != o: added_points.add( o ) + added_points_with_old_id.append( ( o, i ) ) num_new_points: int = old_mesh.GetNumberOfPoints() + len( added_points ) # Creating the new points for the new mesh. old_points: vtkPoints = old_mesh.GetPoints() new_points = vtkPoints() new_points.SetNumberOfPoints( num_new_points ) - collocated_nodes = numpy.ones( num_new_points, dtype=int ) * -1 + collocated_nodes = ones( num_new_points, dtype=int ) * -1 # Copying old points into the new container. for p in range( old_points.GetNumberOfPoints() ): new_points.SetPoint( p, old_points.GetPoint( p ) ) @@ -336,16 +392,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 +417,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_splitted_mesh( old_mesh, new_mesh, added_points_with_old_id ) 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,35 +433,45 @@ def __generate_fracture_mesh( mesh_points: vtkPoints, fracture_info: FractureInf """ logging.info( "Generating the meshes" ) - is_node_duplicated = numpy.zeros( mesh_points.GetNumberOfPoints(), dtype=bool ) # defaults to False + mesh_points: vtkPoints = old_mesh.GetPoints() + is_node_duplicated = 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(): if not is_node_duplicated[ i ]: is_node_duplicated[ i ] = i != o # 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() - for ns in fracture_info.face_nodes: - if any( map( is_node_duplicated.__getitem__, ns ) ): - face_nodes.append( ns ) - else: - discarded_face_nodes.add( ns ) + # 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 ] ] = list() + discarded_face_nodes: set[ Iterable[ int ] ] = set() + if fracture_info.face_cell_id != list(): # The fracture policy is 'internal_surfaces' + face_cell_id: list[ int ] = list() + for ns, f_id in zip( fracture_info.face_nodes, fracture_info.face_cell_id ): + if any( map( is_node_duplicated.__getitem__, ns ) ): + face_nodes.append( ns ) + face_cell_id.append( f_id ) + else: + discarded_face_nodes.add( ns ) + else: # The fracture policy is 'field' + for ns in fracture_info.face_nodes: + if any( map( is_node_duplicated.__getitem__, ns ) ): + face_nodes.append( ns ) + else: + 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." - ) - - fracture_nodes_tmp = numpy.ones( mesh_points.GetNumberOfPoints(), dtype=int ) * -1 + # 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 = ones( mesh_points.GetNumberOfPoints(), dtype=int ) * -1 for ns in face_nodes: for n in ns: fracture_nodes_tmp[ n ] = n @@ -413,12 +479,14 @@ 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 = dict() # 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 + # The polygons are constructed in the same order as the faces defined in the fracture_info. Therefore, + # fracture_info.face_cell_id can be used to link old cells to fracture cells for copy with internal_surfaces. polygons = vtkCellArray() for ns in face_nodes: polygon = vtkPolygon() @@ -427,7 +495,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 ) ) @@ -436,7 +504,7 @@ def __generate_fracture_mesh( mesh_points: vtkPoints, fracture_info: FractureInf assert set( buckets.keys() ) == set( range( num_points ) ) max_collocated_nodes: int = max( map( len, buckets.values() ) ) if buckets.values() else 0 - collocated_nodes = numpy.ones( ( num_points, max_collocated_nodes ), dtype=int ) * -1 + collocated_nodes = ones( ( num_points, max_collocated_nodes ), dtype=int ) * -1 for i, bucket in buckets.items(): for j, val in enumerate( bucket ): collocated_nodes[ i, j ] = val @@ -448,31 +516,52 @@ 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 ) + + # The copy of fields from the old mesh to the fracture is only available when using the internal_surfaces policy + # because the FractureInfo is linked to 2D elements from the old_mesh + if fracture_info.face_cell_id != list(): + __copy_fields_fracture_mesh( old_mesh, fracture_mesh, face_cell_id, node_3d_to_node_2d ) + return fracture_mesh -def __split_mesh_on_fracture( mesh: 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 ) +def __split_mesh_on_fractures( mesh: vtkUnstructuredGrid, + options: Options ) -> tuple[ vtkUnstructuredGrid, list[ vtkUnstructuredGrid ] ]: + all_fracture_infos: list[ FractureInfo ] = list() + for fracture_id in range( len( options.field_values_per_fracture ) ): + fracture_info: FractureInfo = build_fracture_info( mesh, options, False, fracture_id ) + all_fracture_infos.append( fracture_info ) + combined_fractures: FractureInfo = build_fracture_info( mesh, options, True ) + cell_to_cell: networkx.Graph = build_cell_to_cell_graph( mesh, combined_fractures ) + cell_to_node_mapping: Mapping[ int, IDMapping ] = __identify_split( mesh.GetNumberOfPoints(), cell_to_cell, + combined_fractures.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 ) - return output_mesh, fractured_mesh + fracture_meshes: list[ vtkUnstructuredGrid ] = list() + for fracture_info_separated in all_fracture_infos: + fracture_mesh: vtkUnstructuredGrid = __generate_fracture_mesh( mesh, fracture_info_separated, + cell_to_node_mapping ) + fracture_meshes.append( fracture_mesh ) + return ( output_mesh, fracture_meshes ) 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 ) + output_mesh, fracture_meshes = __split_mesh_on_fractures( mesh, options ) + write_mesh( output_mesh, options.mesh_VtkOutput ) + for i, fracture_mesh in enumerate( fracture_meshes ): + write_mesh( fracture_mesh, options.all_fractures_VtkOutput[ i ] ) # 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 ) + # Mesh cannot contain global ids before splitting. + if has_invalid_field( 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." ) + logging.error( err_msg ) + raise ValueError( err_msg ) return __check( mesh, options ) except BaseException as e: logging.error( e ) diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/generate_fractures_parsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/generate_fractures_parsing.py index 9d94b89..0a549aa 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/generate_fractures_parsing.py +++ b/geos-mesh/src/geos/mesh/doctor/parsing/generate_fractures_parsing.py @@ -1,7 +1,7 @@ import logging - +import os from geos.mesh.doctor.checks.generate_fractures import Options, Result, FracturePolicy - +from geos.mesh.doctor.checks.vtk_utils import VtkOutput from . import vtk_output_parsing, GENERATE_FRACTURES __POLICY = "policy" @@ -12,7 +12,10 @@ __FIELD_NAME = "name" __FIELD_VALUES = "values" -__FRACTURE_PREFIX = "fracture" +__FRACTURES_OUTPUT_DIR = "fractures_output_dir" +__FRACTURES_DATA_MODE = "fractures_data_mode" +__FRACTURES_DATA_MODE_VALUES = "binary", "ascii" +__FRACTURES_DATA_MODE_DEFAULT = __FRACTURES_DATA_MODE_VALUES[ 0 ] def convert_to_fracture_policy( s: str ) -> FracturePolicy: @@ -30,8 +33,7 @@ def convert_to_fracture_policy( s: str ) -> FracturePolicy: def fill_subparser( subparsers ) -> None: - p = subparsers.add_parser( GENERATE_FRACTURES, - help="Splits the mesh to generate the faults and fractures. [EXPERIMENTAL]" ) + p = subparsers.add_parser( GENERATE_FRACTURES, help="Splits the mesh to generate the faults and fractures." ) p.add_argument( '--' + __POLICY, type=convert_to_fracture_policy, metavar=", ".join( __POLICIES ), @@ -43,30 +45,88 @@ def fill_subparser( subparsers ) -> None: type=str, help= f"[string]: If the \"{__FIELD_POLICY}\" {__POLICY} is selected, defines which field will be considered to define the fractures. " - f"If the \"{__INTERNAL_SURFACES_POLICY}\" {__POLICY} is selected, defines the name of the attribute will be considered to identify the fractures. " + f"If the \"{__INTERNAL_SURFACES_POLICY}\" {__POLICY} is selected, defines the name of the attribute will be considered to identify the fractures." ) p.add_argument( '--' + __FIELD_VALUES, type=str, help= - f"[list of comma separated integers]: If the \"{__FIELD_POLICY}\" {__POLICY} is selected, which changes of the field will be considered as a fracture. If the \"{__INTERNAL_SURFACES_POLICY}\" {__POLICY} is selected, list of the fracture attributes." - ) + f"[list of comma separated integers]: If the \"{__FIELD_POLICY}\" {__POLICY} is selected, which changes of the field will be considered " + f"as a fracture. If the \"{__INTERNAL_SURFACES_POLICY}\" {__POLICY} is selected, list of the fracture attributes. " + f"You can create multiple fractures by separating the values with ':' like shown in this example. " + f"--{__FIELD_VALUES} 10,12:13,14,16,18:22 will create 3 fractures identified respectively with the values (10,12), (13,14,16,18) and (22). " + f"If no ':' is found, all values specified will be assumed to create only 1 single fracture." ) vtk_output_parsing.fill_vtk_output_subparser( p ) - vtk_output_parsing.fill_vtk_output_subparser( p, prefix=__FRACTURE_PREFIX ) + p.add_argument( + '--' + __FRACTURES_OUTPUT_DIR, + type=str, + help=f"[string]: The output directory for the fractures meshes that will be generated from the mesh." ) + p.add_argument( + '--' + __FRACTURES_DATA_MODE, + type=str, + metavar=", ".join( __FRACTURES_DATA_MODE_VALUES ), + default=__FRACTURES_DATA_MODE_DEFAULT, + help=f'[string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary.' ) def convert( parsed_options ) -> Options: - policy = parsed_options[ __POLICY ] - field = parsed_options[ __FIELD_NAME ] - field_values = frozenset( map( int, parsed_options[ __FIELD_VALUES ].split( "," ) ) ) - vtk_output = vtk_output_parsing.convert( parsed_options ) - vtk_fracture_output = vtk_output_parsing.convert( parsed_options, prefix=__FRACTURE_PREFIX ) + policy: str = parsed_options[ __POLICY ] + field: str = parsed_options[ __FIELD_NAME ] + all_values: str = parsed_options[ __FIELD_VALUES ] + if not are_values_parsable( all_values ): + raise ValueError( + f"When entering --{__FIELD_VALUES}, respect this given format example:\n--{__FIELD_VALUES} " + + "10,12:13,14,16,18:22 to create 3 fractures identified with respectively the values (10,12), (13,14,16,18) and (22)." + ) + all_values_no_separator: str = all_values.replace( ":", "," ) + field_values_combined: frozenset[ int ] = frozenset( map( int, all_values_no_separator.split( "," ) ) ) + mesh_vtk_output = vtk_output_parsing.convert( parsed_options ) + # create the different fractures + per_fracture: list[ str ] = all_values.split( ":" ) + field_values_per_fracture: list[ frozenset[ int ] ] = [ + frozenset( map( int, fracture.split( "," ) ) ) for fracture in per_fracture + ] + fracture_names: list[ str ] = [ "fracture_" + frac.replace( ",", "_" ) + ".vtu" for frac in per_fracture ] + fractures_output_dir: str = parsed_options[ __FRACTURES_OUTPUT_DIR ] + fractures_data_mode: str = parsed_options[ __FRACTURES_DATA_MODE ] == __FRACTURES_DATA_MODE_DEFAULT + all_fractures_VtkOutput: list[ VtkOutput ] = build_all_fractures_VtkOutput( fractures_output_dir, + fractures_data_mode, mesh_vtk_output, + fracture_names ) return Options( policy=policy, field=field, - field_values=field_values, - vtk_output=vtk_output, - vtk_fracture_output=vtk_fracture_output ) + field_values_combined=field_values_combined, + field_values_per_fracture=field_values_per_fracture, + mesh_VtkOutput=mesh_vtk_output, + all_fractures_VtkOutput=all_fractures_VtkOutput ) def display_results( options: Options, result: Result ): pass + + +def are_values_parsable( values: str ) -> bool: + if not all( character.isdigit() or character in { ':', ',' } for character in values ): + return False + if values.startswith( ":" ) or values.startswith( "," ): + return False + if values.endswith( ":" ) or values.endswith( "," ): + return False + return True + + +def build_all_fractures_VtkOutput( fracture_output_dir: str, fractures_data_mode: bool, mesh_vtk_output: VtkOutput, + fracture_names: list[ str ] ) -> list[ VtkOutput ]: + if not os.path.exists( fracture_output_dir ): + raise ValueError( f"The --{__FRACTURES_OUTPUT_DIR} given directory '{fracture_output_dir}' does not exist." ) + + if not os.access( fracture_output_dir, os.W_OK ): + raise ValueError( f"The --{__FRACTURES_OUTPUT_DIR} given directory '{fracture_output_dir}' is not writable." ) + + output_name = os.path.basename( mesh_vtk_output.output ) + splitted_name_without_expension: list[ str ] = output_name.split( "." )[ :-1 ] + name_without_extension: str = '_'.join( splitted_name_without_expension ) + "_" + all_fractures_VtkOuput: list[ VtkOutput ] = list() + for fracture_name in fracture_names: + fracture_path = os.path.join( fracture_output_dir, name_without_extension + fracture_name ) + all_fractures_VtkOuput.append( VtkOutput( fracture_path, fractures_data_mode ) ) + return all_fractures_VtkOuput diff --git a/geos-mesh/tests/test_generate_fractures.py b/geos-mesh/tests/test_generate_fractures.py index 077d9d7..2592c17 100644 --- a/geos-mesh/tests/test_generate_fractures.py +++ b/geos-mesh/tests/test_generate_fractures.py @@ -1,31 +1,18 @@ -from dataclasses import dataclass - -from typing import ( - Tuple, - Iterable, - Iterator, - Sequence, -) - +import logging import numpy - import pytest +from dataclasses import dataclass +from typing import Iterable, Iterator, Sequence +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_fractures, Options, FracturePolicy, + Coordinates3D, IDMapping ) -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 +FaceNodesCoords = tuple[ tuple[ float ] ] +IDMatrix = Sequence[ Sequence[ int ] ] @dataclass( frozen=True ) @@ -42,11 +29,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 ): @@ -66,7 +53,12 @@ def __build_test_case( xs: Tuple[ numpy.ndarray, numpy.ndarray, numpy.ndarray ], else: fv = frozenset( field_values ) - options = Options( policy=policy, field="attribute", field_values=fv, vtk_output=None, vtk_fracture_output=None ) + options = Options( policy=policy, + field="attribute", + field_values_combined=fv, + field_values_per_fracture=[ fv ], + mesh_VtkOutput=None, + all_fractures_VtkOutput=None ) return mesh, options @@ -95,20 +87,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 +99,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 +114,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 +125,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 +136,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 +146,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 +156,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 +169,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 +188,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, ), @@ -277,7 +203,8 @@ def __generate_test_data() -> Iterator[ TestCase ]: @pytest.mark.parametrize( "test_case", __generate_test_data() ) def test_generate_fracture( test_case: TestCase ): - main_mesh, fracture_mesh = __split_mesh_on_fracture( test_case.input_mesh, test_case.options ) + main_mesh, fracture_meshes = __split_mesh_on_fractures( test_case.input_mesh, test_case.options ) + fracture_mesh: vtkUnstructuredGrid = fracture_meshes[ 0 ] assert main_mesh.GetNumberOfPoints() == test_case.result.main_mesh_num_points assert main_mesh.GetNumberOfCells() == test_case.result.main_mesh_num_cells assert fracture_mesh.GetNumberOfPoints() == test_case.result.fracture_mesh_num_points @@ -286,3 +213,148 @@ 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_fractures 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_combined=frozenset( map( int, [ "9" ] ) ), + field_values_per_fracture=[ frozenset( map( int, [ "9" ] ) ) ], + mesh_VtkOutput=None, + all_fractures_VtkOutput=None ) + main_mesh, fracture_meshes = __split_mesh_on_fractures( mesh, options ) + fracture_mesh: vtkUnstructuredGrid = fracture_meshes[ 0 ] + 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_meshes = __split_mesh_on_fractures( 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_meshes = __split_mesh_on_fractures( mesh, options ) + assert pytest_wrapped_e.type == ValueError