Create multiple fractures with generate_fractures (#46)

* First attempt at outputting all fracture meshes. Ok when no more than 1 intersection with 2 faults.

* Corrected first commit, now works for multiple fractures at a same time.

* docs update

* Allows for flexibility in data mode for fractures

* yapf formatting

* Better docs readability

* Corrected tests with new formalism

* TypeAlias cannot be used in prod machines because of old python version

* yapf formatting

Author: alexbenedicto

docs/geos-mesh.rst

@@ -185,27 +185,26 @@ The ``generate_fractures`` module will split the mesh and generate the multi-blo
 .. code-block::
 
       $ python src/geos/mesh/doctor/mesh_doctor.py generate_fractures --help
-      usage: mesh_doctor.py generate_fractures [-h] --policy field, internal_surfaces [--name NAME] [--values VALUES]
-                                             --output OUTPUT [--data-mode binary, ascii] --fracture-output
-                                             FRACTURE_OUTPUT [--fracture-data-mode binary, ascii]
+      usage: mesh_doctor.py generate_fractures [-h] --policy field, internal_surfaces [--name NAME] [--values VALUES] --output OUTPUT
+                                               [--data-mode binary, ascii] [--fractures_output_dir FRACTURES_OUTPUT_DIR]
 
       options:
       -h, --help              show this help message and exit
       --policy field, internal_surfaces
-                              [string]: The criterion to define the surfaces that will be changed into fracture zones.
-                              Possible values are "field, internal_surfaces"
-      --name NAME             [string]: If the "field" policy is selected, defines which field will be considered to
-                              define the fractures. If the "internal_surfaces" policy is selected, defines the name of
-                              the attribute will be considered to identify the fractures.
-      --values VALUES         [list of comma separated integers]: If the "field" policy is selected, which changes of    
-                              the field will be considered as a fracture. If the "internal_surfaces" policy is
-                              selected, list of the fracture attributes.
+                              [string]: The criterion to define the surfaces that will be changed into fracture zones. Possible values are "field, internal_surfaces"
+      --name NAME             [string]: If the "field" policy is selected, defines which field will be considered to define the fractures.
+                              If the "internal_surfaces" policy is selected, defines the name of the attribute will be considered to identify the fractures.
+      --values VALUES         [list of comma separated integers]: If the "field" policy is selected, which changes of the field will be considered as a fracture.
+                              If the "internal_surfaces" policy is selected, list of the fracture attributes.
+                              You can create multiple fractures by separating the values with ':' like shown in this example.
+                              --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).
+                              If no ':' is found, all values specified will be assumed to create only 1 single fracture.
       --output OUTPUT         [string]: The vtk output file destination.
       --data-mode binary, ascii
                               [string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary.
-      --fracture-output FRACTURE_OUTPUT
-                              [string]: The vtk output file destination.
-      --fracture-data-mode binary, ascii
+      --fractures_output_dir FRACTURES_OUTPUT_DIR
+                              [string]: The output directory for the fractures meshes that will be generated from the mesh.
+      --fractures_data_mode FRACTURES_DATA_MODE
                               [string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary.
 
 ``generate_global_ids``

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

@@ -5,7 +5,7 @@
 from dataclasses import dataclass
 from enum import Enum
 from tqdm import tqdm
-from typing import Collection, Iterable, Mapping, Optional, Sequence, TypeAlias
+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,
@@ -16,11 +16,12 @@
                                                 has_invalid_field )
 from geos.mesh.doctor.checks.vtk_polyhedron import FaceStream
 """
-TypeAliases used in this file
+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: TypeAlias = Mapping[ int, int ]
-CellsPointsCoords: TypeAlias = dict[ int, list[ tuple[ float ] ] ]
-Coordinates3D: TypeAlias = tuple[ float ]
+IDMapping = Mapping[ int, int ]
+CellsPointsCoords = dict[ int, list[ tuple[ float ] ] ]
+Coordinates3D = tuple[ float ]
 
 
 class FracturePolicy( Enum ):
@@ -32,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 )
@@ -127,9 +129,15 @@ def __build_fracture_info_from_internal_surfaces( mesh: vtkUnstructuredGrid, f:
     return FractureInfo( node_to_cells=node_to_cells, face_nodes=face_nodes )
 
 
-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 ) )
@@ -538,21 +546,30 @@ def __generate_fracture_mesh( old_mesh: vtkUnstructuredGrid, fracture_info: Frac
     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 )
+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,
-                                                                        fracture.node_to_cells )
+                                                                        combined_fractures.node_to_cells )
     output_mesh: vtkUnstructuredGrid = __perform_split( mesh, cell_to_node_mapping )
-    fractured_mesh: vtkUnstructuredGrid = __generate_fracture_mesh( mesh, 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 )
-    write_mesh( output_mesh, options.vtk_output )
-    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" )
 

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,8 @@
 __FIELD_NAME = "name"
 __FIELD_VALUES = "values"
 
-__FRACTURE_PREFIX = "fracture"
+__FRACTURES_OUTPUT_DIR = "fractures_output_dir"
+__FRACTURES_DATA_MODE = "fractures_data_mode"
 
 
 def convert_to_fracture_policy( s: str ) -> FracturePolicy:
@@ -30,8 +31,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 +43,86 @@ 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,
+        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 ]
+    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: str, 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 does not exist." )
+
+    if not os.access( fracture_output_dir, os.W_OK ):
+        raise ValueError( f"The --{__FRACTURES_OUTPUT_DIR} given directory 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