diff --git a/.github/workflows/README.md b/.github/workflows/README.md index 3f823457..9ad0a906 100644 --- a/.github/workflows/README.md +++ b/.github/workflows/README.md @@ -29,6 +29,7 @@ Tests each Python package independently to ensure: - `geos-xml-tools` - XML preprocessing and formatting - `geos-xml-viewer` - XML viewing tools - `hdf5-wrapper` - HDF5 file handling wrapper +- `mesh-doctor` - Tools to perform checks on vtkUnstructuredGrids - `pygeos-tools` - GEOS Python tools ### Jobs @@ -240,9 +241,10 @@ GEOS integration tests are **automatically triggered** when changes affect: #### GEOS-Integrated Packages - `geos-utils/` - Core utilities used of goesPythonPackages -- `geos-mesh/` - Mesh conversion (`convert_abaqus`, `mesh-doctor`) +- `geos-mesh/` - Mesh conversion (`convert_abaqus`) - `geos-xml-tools/` - XML preprocessing (`preprocess_xml`, `format_xml`) - `hdf5-wrapper/` - HDF5 I/O wrapper +- `mesh-doctor/` - Checks of vtkUnstructuredGrid before using in GEOS - `pygeos-tools/` - Python tools for GEOS - `geos-ats/` - Automated testing framework @@ -305,7 +307,6 @@ The CI uses the following decision matrix: ✅ **Tests Will Run (Required + Label)** ``` Changes: - - geos-mesh/src/mesh_converter.py - geos-xml-tools/src/preprocessor.py Labels: test-geos-integration Result: GEOS integration will run (changes affect integrated packages) @@ -374,6 +375,7 @@ Result: GEOS integration forced (tests will run regardless of changes) │ pip install geos-mesh │ │ pip install geos-xml-tools │ │ pip install hdf5-wrapper │ +│ pip install mesh-doctor │ │ pip install pygeos-tools │ │ pip install geos-ats │ └─────────────────────────────────────────────┘ @@ -402,9 +404,10 @@ Result: GEOS integration forced (tests will run regardless of changes) | Package | Tools | Purpose | |--------------------|-----------------------------------------------------|---------------------------------------------------------------------| | **geos-xml-tools** | `preprocess_xml`;`format_xml` | Preprocess XML input files;Format XML for readability | -| **geos-mesh** | `convert_abaqus`;`mesh-doctor` | Convert Abaqus meshes to VTU/GMSH;Validate and fix mesh quality | +| **geos-mesh** | `convert_abaqus`;`mesh-doctor` | Convert Abaqus meshes to VTU/GMSH | | **geos-ats** | `run_geos_ats`;`setup_ats_environment`;`geos_ats_*` | Run automated test suite;Setup test environment;Various test checks | | **hdf5-wrapper** | Python API | HDF5 file I/O operations | +| **mesh-doctor** | `mesh-doctor` | Validate and fix mesh quality | | **pygeos-tools** | Python API | GEOS workflow utilities | | **geos-utils** | Python API | Common utility functions | diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index 9c668d1b..48f10cce 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -68,6 +68,7 @@ jobs: - geos-trame - geos-xml-tools - geos-xml-viewer + - mesh-doctor - hdf5-wrapper - pygeos-tools include: @@ -77,6 +78,8 @@ jobs: dependencies: "geos-utils geos-geomechanics" - package-name: geos-processing dependencies: "geos-utils geos-mesh geos-geomechanics" + - package-name: mesh-doctor + dependencies: "geos-utils geos-mesh" - package-name: pygeos-tools dependencies: "geos-utils geos-mesh" - package-name: geos-timehistory @@ -201,6 +204,7 @@ jobs: "geos-mesh" "geos-xml-tools" "hdf5-wrapper" + "mesh-doctor" "pygeos-tools" "geos-ats" ) diff --git a/.github/workflows/typing-check.yml b/.github/workflows/typing-check.yml index 2cd6a7c9..6d22bf3b 100644 --- a/.github/workflows/typing-check.yml +++ b/.github/workflows/typing-check.yml @@ -16,7 +16,7 @@ jobs: max-parallel: 3 matrix: # add packages to check typing - package-name: ["geos-geomechanics", "geos-processing", "geos-timehistory", "geos-utils", "geos-trame", "geos-xml-tools", "hdf5-wrapper"] + package-name: ["geos-geomechanics", "geos-processing", "geos-timehistory", "geos-utils", "geos-trame", "geos-xml-tools", "hdf5-wrapper", "mesh-doctor"] steps: - uses: actions/checkout@v4 diff --git a/README.md b/README.md index b0b3f908..5736f9e9 100644 --- a/README.md +++ b/README.md @@ -73,8 +73,8 @@ GEOS Python packages dependency tree (inter-dependency and main external depende │ └── hdf5-wrapper | ├── mesh-doctor -│ ├── geos-prep -│ └── pyvista +│ ├── geos-utils +│ └── geos-mesh | ├── geos-trame │ ├── geos-xml-tools diff --git a/docs/conf.py b/docs/conf.py index 79fc5152..14c630da 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -14,20 +14,43 @@ # import os import sys +import subprocess # Add python modules to be documented python_root = '..' python_modules = ( 'geos-ats', 'geos-geomechanics', 'geos-mesh', 'geos-processing', 'geos-pv', 'geos-timehistory', - 'geos-utils', 'geos-xml-tools', 'geos-xml-viewer', 'hdf5-wrapper', 'pygeos-tools' ) + 'geos-utils', 'geos-xml-tools', 'geos-xml-viewer', 'hdf5-wrapper', 'mesh-doctor', 'pygeos-tools' ) for m in python_modules: sys.path.insert( 0, os.path.abspath( os.path.join( python_root, m, 'src' ) ) ) +# Install mesh-doctor in editable mode if not already available +# This ensures mesh-doctor command is available for sphinxcontrib.programoutput +try: + subprocess.run( [ sys.executable, '-m', 'geos.mesh_doctor.cli', '--help' ], + capture_output=True, check=True, timeout=5 ) +except ( subprocess.CalledProcessError, subprocess.TimeoutExpired, FileNotFoundError ): + # mesh-doctor not available, install it and its dependencies in editable mode + print( "Installing mesh-doctor dependencies for documentation generation..." ) + packages_to_install = [ 'geos-utils', 'geos-mesh', 'mesh-doctor' ] + for pkg in packages_to_install: + pkg_path = os.path.abspath( os.path.join( python_root, pkg ) ) + try: + print( f" Installing {pkg}..." ) + subprocess.run( [ sys.executable, '-m', 'pip', 'install', '-e', pkg_path ], + check=True, capture_output=True ) + except subprocess.CalledProcessError as e: + print( f"Warning: Could not install {pkg}: {e}" ) + print( "Documentation may be incomplete." ) + break + else: + print( "mesh-doctor installed successfully." ) + # -- Project information ----------------------------------------------------- project = u'GEOS Python Packages' -copyright = u'2018-2024 Lawrence Livermore National Security, The Board of Trustees of the Leland Stanford Junior University, TotalEnergies, and GEOSX Contributors.' +copyright = u'2018-2025 Lawrence Livermore National Security, The Board of Trustees of the Leland Stanford Junior University, TotalEnergies, and GEOS Contributors.' author = u'GEOS Contributors' # The short X.Y version diff --git a/docs/geos-mesh.rst b/docs/geos-mesh.rst index cb8a5264..ba604f8f 100644 --- a/docs/geos-mesh.rst +++ b/docs/geos-mesh.rst @@ -7,8 +7,6 @@ GEOS Mesh tools :maxdepth: 1 :caption: Contents: - ./geos_mesh_docs/doctor - ./geos_mesh_docs/converter ./geos_mesh_docs/io diff --git a/docs/geos_mesh_docs/doctor.rst b/docs/geos_mesh_docs/doctor.rst deleted file mode 100644 index 123d673a..00000000 --- a/docs/geos_mesh_docs/doctor.rst +++ /dev/null @@ -1,317 +0,0 @@ -Mesh Doctor ---------------- - -``mesh-doctor`` is a ``python`` executable that can be used through the command line to perform various checks, validations, and tiny fixes to the ``vtk`` mesh that are meant to be used in ``geos``. -``mesh-doctor`` is organized as a collection of modules with their dedicated sets of options. -The current page will introduce those modules, but the details and all the arguments can be retrieved by using the ``--help`` option for each module. - -Prerequisites -^^^^^^^^^^^^^ - -To use mesh-doctor, you first need to have installed the ``geos-mesh`` package using the following command: - -.. code-block:: bash - - python -m pip install --upgrade ./geos-mesh - -Once done, you can call ``mesh-doctor`` or ``meshDoctor`` in your command line as presented in the rest of this documentation. - -Modules -^^^^^^^ - -To list all the modules available through ``mesh-doctor``, you can simply use the ``--help`` option, which will list all available modules as well as a quick summary. - -.. code-block:: - - $ mesh-doctor --help - usage: meshDoctor.py [-h] [-v] [-q] -i VTK_MESH_FILE - {collocatedNodes,elementVolumes,fixElementsOrderings,generateCube,generateFractures,generateGlobalIds,nonConformal,selfIntersectingElements,supportedElements} - ... - Inspects meshes for GEOSX. - positional arguments: - {collocatedNodes,elementVolumes,fixElementsOrderings,generateCube,generateFractures,generateGlobalIds,nonConformal,selfIntersectingElements,supportedElements} - Modules - collocatedNodes - Checks if nodes are collocated. - elementVolumes - Checks if the volumes of the elements are greater than "min". - fixElementsOrderings - Reorders the support nodes for the given cell types. - generateCube - Generate a cube and its fields. - generateFractures - Splits the mesh to generate the faults and fractures. [EXPERIMENTAL] - generateGlobalIds - Adds globals ids for points and cells. - nonConformal - Detects non conformal elements. [EXPERIMENTAL] - selfIntersectingElements - Checks if the faces of the elements are self intersecting. - supportedElements - Check that all the elements of the mesh are supported by GEOSX. - options: - -h, --help - show this help message and exit - -v Use -v 'INFO', -vv for 'DEBUG'. Defaults to 'WARNING'. - -q Use -q to reduce the verbosity of the output. - -i VTK_MESH_FILE, --vtk-input-file VTK_MESH_FILE - Note that checks are dynamically loaded. - An option may be missing because of an unloaded module. - Increase verbosity (-v, -vv) to get full information. - -Then, if you are interested in a specific module, you can ask for its documentation using the ``mesh-doctor module_name --help`` pattern. -For example - -.. code-block:: - - $ mesh-doctor collocatedNodes --help - usage: meshDoctor.py collocatedNodes [-h] --tolerance TOLERANCE - options: - -h, --help show this help message and exit - --tolerance TOLERANCE [float]: The absolute distance between two nodes for them to be considered collocated. - -``mesh-doctor`` loads its module dynamically. -If a module can't be loaded, ``mesh-doctor`` will proceed and try to load other modules. -If you see a message like - -.. code-block:: bash - - [1970-04-14 03:07:15,625][WARNING] Could not load module "collocatedNodes": No module named 'vtkmodules' - -then most likely ``mesh-doctor`` could not load the ``collocatedNodes`` module, because the ``vtk`` python package was not found. -Thereafter, the documentation for module ``collocatedNodes`` will not be displayed. -You can solve this issue by installing the dependencies of ``mesh-doctor`` defined in its ``requirements.txt`` file (``python -m pip install -r requirements.txt``). - -Here is a list and brief description of all the modules available. - -``allChecks`` and ``mainChecks`` -"""""""""""""""""""""""""""""""" - -``mesh-doctor`` modules are called ``actions`` and they can be split into 2 different categories: -``check actions`` that will give you a feedback on a .vtu mesh that you would like to use in GEOS. -``operate actions`` that will either create a new mesh or modify an existing mesh. - -``allChecks`` aims at applying every single ``check`` action in one single command. The available list is of check is: -``collocatedNodes``, ``elementVolumes``, ``nonConformal``, ``selfIntersectingElements``, ``supportedElements``. - -``mainChecks`` does only the fastest checks ``collocatedNodes``, ``elementVolumes`` and ``selfIntersectingElements`` -that can quickly highlight some issues to deal with before investigating the other checks. - -Both ``allChecks`` and ``mainChecks`` have the same keywords and can be operated in the same way. The example below shows -the case of ``allChecks``, but it can be swapped for ``mainChecks``. - -.. code-block:: - - $ mesh-doctor allChecks --help - usage: mesh-doctor allChecks [-h] [--checksToPerform checksToPerform] [--set_parameters SET_PARAMETERS] - - options: - -h, --help show this help message and exit - --checksToPerform checksToPerform - Comma-separated list of mesh-doctor checks to perform. - If no input was given, all of the following checks will be executed by default: ['collocatedNodes', 'elementVolumes', 'selfIntersectingElements']. - The available choices for checks are ['collocatedNodes', 'elementVolumes', 'nonConformal', 'selfIntersectingElements', 'supportedElements']. - If you want to choose only certain of them, you can name them individually. - Example: --checksToPerform collocatedNodes,elementVolumes (default: ) - --setParameters setParameters - Comma-separated list of parameters to set for the checks (e.g., 'param_name:value'). These parameters override the defaults. - Default parameters are: For collocatedNodes: tolerance:0.0. For elementVolumes: minVolume:0.0. - For nonConformal: angleTolerance:10.0, pointTolerance:0.0, faceTolerance:0.0. - For selfIntersectingElements: minDistance:2.220446049250313e-16. For supportedElements: chunkSize:1, nproc:8. - Example: --setParameters parameter_name:10.5,other_param:25 (default: ) - -``collocatedNodes`` -"""""""""""""""""""" - -Displays the neighboring nodes that are closer to each other than a prescribed threshold. -It is not uncommon to define multiple nodes for the exact same position, which will typically be an issue for ``geos`` and should be fixed. - -.. code-block:: - - $ mesh-doctor collocatedNodes --help - usage: meshDoctor.py collocatedNodes [-h] --tolerance TOLERANCE - options: - -h, --help show this help message and exit - --tolerance TOLERANCE [float]: The absolute distance between two nodes for them to be considered collocated. - -``elementVolumes`` -"""""""""""""""""" - -Computes the volumes of all the cells and displays the ones that are below a prescribed threshold. -Cells with negative volumes will typically be an issue for ``geos`` and should be fixed. - -.. code-block:: - - $ mesh-doctor elementVolumes --help - usage: meshDoctor.py elementVolumes [-h] --minVolume 0.0 - options: - -h, --help show this help message and exit - --minVolume 0.0 [float]: The minimum acceptable volume. Defaults to 0.0. - -``fixElementsOrderings`` -"""""""""""""""""""""""" - -It sometimes happens that an exported mesh does not abide by the ``vtk`` orderings. -The ``fixElementsOrderings`` module can rearrange the nodes of given types of elements. -This can be convenient if you cannot regenerate the mesh. - -.. code-block:: - - $ mesh-doctor fixElementsOrderings --help - usage: meshDoctor.py fixElementsOrderings [-h] [--Hexahedron 1,6,5,4,7,0,2,3] [--Prism5 8,2,0,7,6,9,5,1,4,3] - [--Prism6 11,2,8,10,5,0,9,7,6,1,4,3] [--Pyramid 3,4,0,2,1] - [--Tetrahedron 2,0,3,1] [--Voxel 1,6,5,4,7,0,2,3] - [--Wedge 3,5,4,0,2,1] --output OUTPUT [--data-mode binary, ascii] - options: - -h, --help show this help message and exit - --Hexahedron 1,6,5,4,7,0,2,3 - [list of integers]: node permutation for "Hexahedron". - --Prism5 8,2,0,7,6,9,5,1,4,3 - [list of integers]: node permutation for "Prism5". - --Prism6 11,2,8,10,5,0,9,7,6,1,4,3 - [list of integers]: node permutation for "Prism6". - --Pyramid 3,4,0,2,1 [list of integers]: node permutation for "Pyramid". - --Tetrahedron 2,0,3,1 [list of integers]: node permutation for "Tetrahedron". - --Voxel 1,6,5,4,7,0,2,3 [list of integers]: node permutation for "Voxel". - --Wedge 3,5,4,0,2,1 [list of integers]: node permutation for "Wedge". - --output OUTPUT [string]: The vtk output file destination. - --data-mode binary, ascii - [string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary. - -``generateCube`` -"""""""""""""""" - -This module conveniently generates cubic meshes in ``vtk``. -It can also generate fields with simple values. -This tool can also be useful to generate a trial mesh that will later be refined or customized. - -.. code-block:: - - $ mesh-doctor generateCube --help - usage: meshDoctor.py generateCube [-h] [--x 0:1.5:3] [--y 0:5:10] [--z 0:1] [--nx 2:2] [--ny 1:1] [--nz 4] - [--fields name:support:dim [name:support:dim ...]] [--cells] [--no-cells] - [--points] [--no-points] --output OUTPUT [--data-mode binary, ascii] - options: - -h, --help show this help message and exit - --x 0:1.5:3 [list of floats]: X coordinates of the points. - --y 0:5:10 [list of floats]: Y coordinates of the points. - --z 0:1 [list of floats]: Z coordinates of the points. - --nx 2:2 [list of integers]: Number of elements in the X direction. - --ny 1:1 [list of integers]: Number of elements in the Y direction. - --nz 4 [list of integers]: Number of elements in the Z direction. - --fields name:support:dim - [name:support:dim ...]: Create fields on CELLS or POINTS, with given dimension (typically 1 or 3). - --cells [bool]: Generate global ids for cells. Defaults to true. - --no-cells [bool]: Don't generate global ids for cells. - --points [bool]: Generate global ids for points. Defaults to true. - --no-points [bool]: Don't generate global ids for points. - --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. - -``generateFractures`` -""""""""""""""""""""" - -For a conformal fracture to be defined in a mesh, ``geos`` requires the mesh to be split at the faces where the fracture gets across the mesh. -The ``generateFractures`` module will split the mesh and generate the multi-block ``vtk`` files. - -.. code-block:: - - $ mesh-doctor generateFractures --help - usage: meshDoctor.py generateFractures [-h] --policy field, internalSurfaces [--name NAME] [--values VALUES] --output OUTPUT - [--data-mode binary, ascii] [--fracturesOutputDir FRACTURES_OUTPUT_DIR] - options: - -h, --help show this help message and exit - --policy field, internalSurfaces - [string]: The criterion to define the surfaces that will be changed into fracture zones. Possible values are "field, internalSurfaces" - --name NAME [string]: If the "field" policy is selected, defines which field will be considered to define the fractures. - If the "internalSurfaces" 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 "internalSurfaces" 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. - --fracturesOutputDir FRACTURES_OUTPUT_DIR - [string]: The output directory for the fractures meshes that will be generated from the mesh. - --fracturesDataMode FRACTURES_DATA_MODE - [string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary. - -``generateGlobalIds`` -""""""""""""""""""""" - -When running ``geos`` in parallel, `global ids` can be used to refer to data across multiple ranks. -The ``generateGlobalIds`` can generate `global ids` for the imported ``vtk`` mesh. - -.. code-block:: - - $ mesh-doctor generateGlobalIds --help - usage: meshDoctor.py generateGlobalIds [-h] [--cells] [--no-cells] [--points] [--no-points] --output OUTPUT - [--data-mode binary, ascii] - options: - -h, --help show this help message and exit - --cells [bool]: Generate global ids for cells. Defaults to true. - --no-cells [bool]: Don't generate global ids for cells. - --points [bool]: Generate global ids for points. Defaults to true. - --no-points [bool]: Don't generate global ids for points. - --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. - -``nonConformal`` -"""""""""""""""" - -This module will detect elements which are close enough (there's a user defined threshold) but which are not in front of each other (another threshold can be defined). -`Close enough` can be defined in terms or proximity of the nodes and faces of the elements. -The angle between two faces can also be precribed. -This module can be a bit time consuming. - -.. code-block:: - - $ mesh-doctor nonConformal --help - usage: meshDoctor.py nonConformal [-h] [--angleTolerance 10.0] [--pointTolerance POINT_TOLERANCE] - [--faceTolerance FACE_TOLERANCE] - options: - -h, --help show this help message and exit - --angleTolerance 10.0 [float]: angle tolerance in degrees. Defaults to 10.0 - --pointTolerance POINT_TOLERANCE - [float]: tolerance for two points to be considered collocated. - --faceTolerance FACE_TOLERANCE - [float]: tolerance for two faces to be considered "touching". - -``selfIntersectingElements`` -"""""""""""""""""""""""""""" - -Some meshes can have cells that auto-intersect. -This module will display the elements that have faces intersecting. - -.. code-block:: - - $ mesh-doctor selfIntersectingElements --help - usage: meshDoctor.py selfIntersectingElements [-h] [--minDistance 2.220446049250313e-16] - options: - -h, --help show this help message and exit - --minDistance 2.220446049250313e-16 - [float]: The tolerance in the computation. Defaults to your machine precision 2.220446049250313e-16. - -``supportedElements`` -""""""""""""""""""""" - -``geos`` supports a specific set of elements. -Let's cite the standard elements like `tetrahedra`, `wedges`, `pyramids` or `hexahedra`. -But also prismes up to 11 faces. -``geos`` also supports the generic ``VTK_POLYHEDRON``/``42`` elements, which are converted on the fly into one of the elements just described. - -The ``supportedElements`` check will validate that no unsupported element is included in the input mesh. -It will also verify that the ``VTK_POLYHEDRON`` cells can effectively get converted into a supported type of element. - -.. code-block:: - - $ mesh-doctor supportedElements --help - usage: meshDoctor.py supportedElements [-h] [--chunkSize 1] [--nproc 8] - options: - -h, --help show this help message and exit - --chunkSize 1 [int]: Defaults chunk size for parallel processing to 1 - --nproc 8 [int]: Number of threads used for parallel processing. Defaults to your CPU count 8. \ No newline at end of file diff --git a/docs/index.rst b/docs/index.rst index 3f7b66d5..91c8c8a8 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -98,4 +98,6 @@ Packages geos-xml-viewer + mesh-doctor + pygeos-tools \ No newline at end of file diff --git a/docs/mesh-doctor.rst b/docs/mesh-doctor.rst new file mode 100644 index 00000000..2ff03f0e --- /dev/null +++ b/docs/mesh-doctor.rst @@ -0,0 +1,216 @@ +Mesh Doctor +----------- + +| ``mesh-doctor`` is a ``python`` executable that can be used through the command line to perform various checks, validations, and tiny fixes to the ``vtkUnstructuredGrid`` mesh that are meant to be used in ``geos``. + ``mesh-doctor`` is organized as a collection of modules with their dedicated sets of options. +| The current page will introduce those modules, but the details and all the arguments can be retrieved by using the ``--help`` option for each module. + +Prerequisites +^^^^^^^^^^^^^ + +To use mesh-doctor, you need to have installed the mandatory packages using the following commands (in this order): + +.. code-block:: bash + + python -m pip install --upgrade ./geos-utils + python -m pip install --upgrade ./geos-mesh + python -m pip install --upgrade ./mesh-doctor + +Once done, you can call ``mesh-doctor`` in your command line as presented in the rest of this documentation. + +Modules +^^^^^^^ + +To list all the modules available through ``mesh-doctor``, you can simply use the ``--help`` option, which will list all available modules as well as a quick summary. + +.. command-output:: mesh-doctor --help + :shell: + +Then, if you are interested in a specific module, you can ask for its documentation using the ``mesh-doctor module_name --help`` pattern. +For example + +.. command-output:: mesh-doctor collocatedNodes --help + :shell: + +Finally, to execute a specific module, you can use the following pattern: + +.. code-block:: bash + + mesh-doctor module_name -i /path/to/input.vtu --output /path/to/output.vtu [module_options] + +``mesh-doctor`` loads its module dynamically. +If a module can't be loaded, ``mesh-doctor`` will proceed and try to load other modules. +If you see a message like + +.. code-block:: bash + + [1970-04-14 03:07:15,625][WARNING] Could not load module "collocatedNodes": No module named 'vtkmodules' + +then most likely ``mesh-doctor`` could not load the ``collocatedNodes`` module, because the ``vtk`` python package was not found. +Thereafter, the documentation for module ``collocatedNodes`` will not be displayed. +You can solve this issue by installing the dependencies of ``mesh-doctor`` defined in its ``requirements.txt`` file (``python -m pip install -r requirements.txt``). + +Here is a list and brief description of all the modules available. + +``allChecks`` and ``mainChecks`` +"""""""""""""""""""""""""""""""" + +``mesh-doctor`` modules are called ``actions`` and they can be split into 2 different categories: +``check actions`` that will give you a feedback on a .vtu mesh that you would like to use in GEOS. +``operate actions`` that will either create a new mesh or modify an existing mesh. + +``allChecks`` aims at applying every single ``check`` action in one single command. The available list is of check is: +``collocatedNodes``, ``elementVolumes``, ``nonConformal``, ``selfIntersectingElements``, ``supportedElements``. + +``mainChecks`` does only the fastest checks ``collocatedNodes``, ``elementVolumes`` and ``selfIntersectingElements`` +that can quickly highlight some issues to deal with before investigating the other checks. + +Both ``allChecks`` and ``mainChecks`` have the same keywords and can be operated in the same way. The example below shows +the case of ``allChecks``, but it can be swapped for ``mainChecks``. + +.. command-output:: mesh-doctor allChecks --help + :shell: + +``collocatedNodes`` +""""""""""""""""""" + +Displays the neighboring nodes that are closer to each other than a prescribed threshold. +It is not uncommon to define multiple nodes for the exact same position, which will typically be an issue for ``geos`` and should be fixed. + +.. command-output:: mesh-doctor collocatedNodes --help + :shell: + +``elementVolumes`` +"""""""""""""""""" + +Computes the volumes of all the cells and displays the ones that are below a prescribed threshold. +Cells with negative volumes will typically be an issue for ``geos`` and should be fixed. + +.. command-output:: mesh-doctor elementVolumes --help + :shell: + +``fixElementsOrderings`` +"""""""""""""""""""""""" + +It sometimes happens that an exported mesh does not abide by the ``vtk`` orderings. +The ``fixElementsOrderings`` module can rearrange the nodes of given types of elements. +This can be convenient if you cannot regenerate the mesh. + +.. command-output:: mesh-doctor fixElementsOrderings --help + :shell: + +``generateCube`` +"""""""""""""""" + +This module conveniently generates cubic meshes in ``vtk``. +It can also generate fields with simple values. +This tool can also be useful to generate a trial mesh that will later be refined or customized. + +.. command-output:: mesh-doctor generateCube --help + :shell: + +Exceptionally, this module does not require an input ``vtk`` mesh because its purpose is to generate a new one. +The command to use would be something like: + +.. code-block:: bash + + mesh-doctor generateCube --output cube.vtu --x 0:10 --y 0:10 --z 0:10 --nx 10 --ny 15 --nz 20 --cells --nopoints + +``generateFractures`` +""""""""""""""""""""" + +For a conformal fracture to be defined in a mesh, ``geos`` requires the mesh to be split at the faces where the fracture gets across the mesh. +The ``generateFractures`` module will split the mesh and generate the multi-block ``vtk`` files. + +.. command-output:: mesh-doctor generateFractures --help + :shell: + +``generateGlobalIds`` +""""""""""""""""""""" + +When running ``geos`` in parallel, `global ids` can be used to refer to data across multiple ranks. +The ``generateGlobalIds`` can generate `global ids` for the imported ``vtk`` mesh. + +.. command-output:: mesh-doctor generateGlobalIds --help + :shell: + +``nonConformal`` +"""""""""""""""" + +This module will detect elements which are close enough (there's a user defined threshold) but which are not in front of each other (another threshold can be defined). +`Close enough` can be defined in terms or proximity of the nodes and faces of the elements. +The angle between two faces can also be precribed. +This module can be a bit time consuming. + +.. command-output:: mesh-doctor nonConformal --help + :shell: + +``selfIntersectingElements`` +"""""""""""""""""""""""""""" + +Some meshes can have cells that auto-intersect. +This module will display the elements that have faces intersecting. + +.. command-output:: mesh-doctor selfIntersectingElements --help + :shell: + +``supportedElements`` +""""""""""""""""""""" + +``geos`` supports a specific set of elements. +Let's cite the standard elements like `tetrahedra`, `wedges`, `pyramids` or `hexahedra`. +But also prismes up to 11 faces. +``geos`` also supports the generic ``VTK_POLYHEDRON``/``42`` elements, which are converted on the fly into one of the elements just described. + +The ``supportedElements`` check will validate that no unsupported element is included in the input mesh. +It will also verify that the ``VTK_POLYHEDRON`` cells can effectively get converted into a supported type of element. + +.. command-output:: mesh-doctor supportedElements --help + :shell: + + +Why only use vtkUnstructuredGrid? +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +| The mesh doctor is designed specifically for unstructured meshes used in GEOS. +| All input files are expected to be ``.vtu`` (VTK Unstructured Grid) format. +| What about other formats? + +VTK Hierarchy +""""""""""""" + +Supposedly, other grid types that are part of the following VTK hierarchy could be used:: + + vtkDataObject + └── vtkDataSet + └── vtkCartesianGrid + └── vtkRectilinearGrid + └── vtkImageData + └── vtkStructuredPoints + └── vtkUniformGrid + └── vtkPointSet + └── vtkExplicitStructuredGrid + └── vtkPolyData + └── vtkStructuredGrid + └── vtkUnstructuredGrid + +And when looking at specific methods used in mesh-doctor, it could suggest that other formats could be used: + +* Points access: ``mesh.GetPoints()`` - Available in all vtkPointSet subclasses ✓ +* Cell iteration: ``mesh.GetNumberOfCells()``, ``mesh.GetCell()`` - Available in all vtkDataSet subclasses ✓ +* Cell types: ``mesh.GetCellType()`` - Available in all vtkDataSet subclasses ✓ +* Cell/Point data: ``mesh.GetCellData()``, ``mesh.GetPointData()`` - Available in all vtkDataSet subclasses ✓ + +VTK Filter Compatibility +"""""""""""""""""""""""" + +| ``vtkCellSizeFilter``, ``vtkMeshQuality``, and other VTK filters used in the actions expect ``vtkDataSet`` or its subclasses + ``vtkUnstructuredGrid`` is compatible with all VTK filters used. +| ``vtkPolyData`` has a different data structure, not suitable for 3D volumetric meshes. + +Specific Operations Require vtkUnstructuredGrid +""""""""""""""""""""""""""""""""""""""""""""""" + +* ``GetCellNeighbors()`` - Only available in vtkUnstructuredGrid +* ``GetFaceStream()`` - Only available in vtkUnstructuredGrid (for polyhedron support) +* ``GetDistinctCellTypesArray()`` - Only available in vtkUnstructuredGrid \ No newline at end of file diff --git a/geos-mesh/pyproject.toml b/geos-mesh/pyproject.toml index 2daf4301..7cf4ef76 100644 --- a/geos-mesh/pyproject.toml +++ b/geos-mesh/pyproject.toml @@ -39,8 +39,6 @@ dependencies = [ ] [project.scripts] - mesh-doctor = "geos.mesh.doctor.meshDoctor:main" - meshDoctor = "geos.mesh.doctor.meshDoctor:main" convert_abaqus = "geos.mesh.conversion.main:main" [project.urls] diff --git a/geos-mesh/src/geos/mesh/doctor/__init__.py b/geos-mesh/src/geos/mesh/doctor/__init__.py deleted file mode 100644 index b1cfe267..00000000 --- a/geos-mesh/src/geos/mesh/doctor/__init__.py +++ /dev/null @@ -1 +0,0 @@ -# Empty \ No newline at end of file diff --git a/geos-mesh/src/geos/mesh/doctor/actions/allChecks.py b/geos-mesh/src/geos/mesh/doctor/actions/allChecks.py deleted file mode 100644 index 9b7f7c09..00000000 --- a/geos-mesh/src/geos/mesh/doctor/actions/allChecks.py +++ /dev/null @@ -1,26 +0,0 @@ -from dataclasses import dataclass -from geos.mesh.doctor.register import __loadModuleAction -from geos.mesh.doctor.parsing.cliParsing import setupLogger - - -@dataclass( frozen=True ) -class Options: - checksToPerform: list[ str ] - checksOptions: dict[ str, any ] - checkDisplays: dict[ str, any ] - - -@dataclass( frozen=True ) -class Result: - checkResults: dict[ str, any ] - - -def action( vtkInputFile: str, options: Options ) -> list[ Result ]: - checkResults: dict[ str, any ] = dict() - for checkName in options.checksToPerform: - checkAction = __loadModuleAction( checkName ) - setupLogger.info( f"Performing check '{checkName}'." ) - option = options.checksOptions[ checkName ] - checkResult = checkAction( vtkInputFile, option ) - checkResults[ checkName ] = checkResult - return Result( checkResults=checkResults ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/collocatedNodes.py b/geos-mesh/src/geos/mesh/doctor/actions/collocatedNodes.py deleted file mode 100644 index 992d23eb..00000000 --- a/geos-mesh/src/geos/mesh/doctor/actions/collocatedNodes.py +++ /dev/null @@ -1,67 +0,0 @@ -from collections import defaultdict -from dataclasses import dataclass -import numpy -from typing import Collection, Iterable -from vtkmodules.vtkCommonCore import reference, vtkPoints -from vtkmodules.vtkCommonDataModel import vtkIncrementalOctreePointLocator, vtkUnstructuredGrid -from geos.mesh.doctor.parsing.cliParsing import setupLogger -from geos.mesh.io.vtkIO import readUnstructuredGrid - - -@dataclass( frozen=True ) -class Options: - tolerance: float - - -@dataclass( frozen=True ) -class Result: - nodesBuckets: Iterable[ Iterable[ int ] ] # Each bucket contains the duplicated node indices. - wrongSupportElements: Collection[ int ] # Element indices with support node indices appearing more than once. - - -def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: - points = mesh.GetPoints() - - locator = vtkIncrementalOctreePointLocator() - locator.SetTolerance( options.tolerance ) - output = vtkPoints() - locator.InitPointInsertion( output, points.GetBounds() ) - - # original ids to/from filtered ids. - filteredToOriginal = numpy.ones( points.GetNumberOfPoints(), dtype=int ) * -1 - - rejectedPoints = defaultdict( list ) - pointId = reference( 0 ) - for i in range( points.GetNumberOfPoints() ): - isInserted = locator.InsertUniquePoint( points.GetPoint( i ), pointId ) - if not isInserted: - # If it's not inserted, `pointId` contains the node that was already at that location. - # But in that case, `pointId` is the new numbering in the destination points array. - # It's more useful for the user to get the old index in the original mesh so he can look for it in his data. - setupLogger.debug( f"Point {i} at {points.GetPoint(i)} has been rejected, " - f"point {filteredToOriginal[pointId.get()]} is already inserted." ) - rejectedPoints[ pointId.get() ].append( i ) - else: - # If it's inserted, `pointId` contains the new index in the destination array. - # We store this information to be able to connect the source and destination arrays. - # originalToFiltered[i] = pointId.get() - filteredToOriginal[ pointId.get() ] = i - - tmp = [] - for n, ns in rejectedPoints.items(): - tmp.append( ( n, *ns ) ) - - # Checking that the support node indices appear only once per element. - wrongSupportElements = [] - for c in range( mesh.GetNumberOfCells() ): - cell = mesh.GetCell( c ) - numPointsPerCell = cell.GetNumberOfPoints() - if len( { cell.GetPointId( i ) for i in range( numPointsPerCell ) } ) != numPointsPerCell: - wrongSupportElements.append( c ) - - return Result( nodesBuckets=tmp, wrongSupportElements=wrongSupportElements ) - - -def action( vtkInputFile: str, options: Options ) -> Result: - mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtkInputFile ) - return __action( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/elementVolumes.py b/geos-mesh/src/geos/mesh/doctor/actions/elementVolumes.py deleted file mode 100644 index 03430c66..00000000 --- a/geos-mesh/src/geos/mesh/doctor/actions/elementVolumes.py +++ /dev/null @@ -1,70 +0,0 @@ -from dataclasses import dataclass -import uuid -from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, VTK_HEXAHEDRON, VTK_PYRAMID, VTK_TETRA, VTK_WEDGE -from vtkmodules.vtkFiltersVerdict import vtkCellSizeFilter, vtkMeshQuality -from vtkmodules.util.numpy_support import vtk_to_numpy -from geos.mesh.doctor.parsing.cliParsing import setupLogger -from geos.mesh.io.vtkIO import readUnstructuredGrid - - -@dataclass( frozen=True ) -class Options: - minVolume: float - - -@dataclass( frozen=True ) -class Result: - elementVolumes: list[ tuple[ int, float ] ] - - -def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: - cs = vtkCellSizeFilter() - - cs.ComputeAreaOff() - cs.ComputeLengthOff() - cs.ComputeSumOff() - cs.ComputeVertexCountOff() - cs.ComputeVolumeOn() - volumeArrayName = "__MESH_DOCTOR_VOLUME-" + str( uuid.uuid4() ) # Making the name unique - cs.SetVolumeArrayName( volumeArrayName ) - - cs.SetInputData( mesh ) - cs.Update() - - mq = vtkMeshQuality() - SUPPORTED_TYPES = [ VTK_HEXAHEDRON, VTK_TETRA ] - - mq.SetTetQualityMeasureToVolume() - mq.SetHexQualityMeasureToVolume() - if hasattr( mq, "SetPyramidQualityMeasureToVolume" ): # This feature is quite recent - mq.SetPyramidQualityMeasureToVolume() - SUPPORTED_TYPES.append( VTK_PYRAMID ) - mq.SetWedgeQualityMeasureToVolume() - SUPPORTED_TYPES.append( VTK_WEDGE ) - else: - setupLogger.warning( - "Your \"pyvtk\" version does not bring pyramid nor wedge support with vtkMeshQuality. Using the fallback solution." - ) - - mq.SetInputData( mesh ) - mq.Update() - - volume = cs.GetOutput().GetCellData().GetArray( volumeArrayName ) - quality = mq.GetOutput().GetCellData().GetArray( "Quality" ) # Name is imposed by vtk. - - assert volume is not None - assert quality is not None - volume = vtk_to_numpy( volume ) - quality = vtk_to_numpy( quality ) - smallVolumes: list[ tuple[ int, float ] ] = [] - for i, pack in enumerate( zip( volume, quality ) ): - v, q = pack - vol = q if mesh.GetCellType( i ) in SUPPORTED_TYPES else v - if vol < options.minVolume: - smallVolumes.append( ( i, float( vol ) ) ) - return Result( elementVolumes=smallVolumes ) - - -def action( vtkInputFile: str, options: Options ) -> Result: - mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtkInputFile ) - return __action( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/generateCube.py b/geos-mesh/src/geos/mesh/doctor/actions/generateCube.py deleted file mode 100644 index 4247eddf..00000000 --- a/geos-mesh/src/geos/mesh/doctor/actions/generateCube.py +++ /dev/null @@ -1,149 +0,0 @@ -from dataclasses import dataclass -import numpy -from typing import Iterable, Sequence -from vtkmodules.util.numpy_support import numpy_to_vtk -from vtkmodules.vtkCommonCore import vtkPoints -from vtkmodules.vtkCommonDataModel import ( vtkCellArray, vtkHexahedron, vtkRectilinearGrid, vtkUnstructuredGrid, - VTK_HEXAHEDRON ) -from geos.mesh.doctor.actions.generateGlobalIds import __buildGlobalIds -from geos.mesh.doctor.parsing.cliParsing import setupLogger -from geos.mesh.io.vtkIO import VtkOutput, writeMesh - - -@dataclass( frozen=True ) -class Result: - info: str - - -@dataclass( frozen=True ) -class FieldInfo: - name: str - dimension: int - support: str - - -@dataclass( frozen=True ) -class Options: - vtkOutput: VtkOutput - generateCellsGlobalIds: bool - generatePointsGlobalIds: bool - xs: Sequence[ float ] - ys: Sequence[ float ] - zs: Sequence[ float ] - nxs: Sequence[ int ] - nys: Sequence[ int ] - nzs: Sequence[ int ] - fields: Iterable[ FieldInfo ] - - -@dataclass( frozen=True ) -class XYZ: - x: numpy.ndarray - y: numpy.ndarray - z: numpy.ndarray - - -def buildRectilinearBlocksMesh( xyzs: Iterable[ XYZ ] ) -> vtkUnstructuredGrid: - """Builds an unstructured vtk grid from the `xyzs` blocks. Kind of InternalMeshGenerator. - - Args: - xyzs (Iterable[ XYZ ]): The blocks. - - Returns: - vtkUnstructuredGrid: The unstructured mesh, even if it's topologically structured. - """ - rgs = [] - for xyz in xyzs: - rg = vtkRectilinearGrid() - rg.SetDimensions( len( xyz.x ), len( xyz.y ), len( xyz.z ) ) - rg.SetXCoordinates( numpy_to_vtk( xyz.x ) ) - rg.SetYCoordinates( numpy_to_vtk( xyz.y ) ) - rg.SetZCoordinates( numpy_to_vtk( xyz.z ) ) - rgs.append( rg ) - - numPoints = sum( map( lambda r: r.GetNumberOfPoints(), rgs ) ) - numCells = sum( map( lambda r: r.GetNumberOfCells(), rgs ) ) - - points = vtkPoints() - points.Allocate( numPoints ) - for rg in rgs: - for i in range( rg.GetNumberOfPoints() ): - points.InsertNextPoint( rg.GetPoint( i ) ) - - cellTypes = [ VTK_HEXAHEDRON ] * numCells - cells = vtkCellArray() - cells.AllocateExact( numCells, numCells * 8 ) - - m = ( 0, 1, 3, 2, 4, 5, 7, 6 ) # VTK_VOXEL and VTK_HEXAHEDRON do not share the same ordering. - offset = 0 - for rg in rgs: - for i in range( rg.GetNumberOfCells() ): - c = rg.GetCell( i ) - newCell = vtkHexahedron() - for j in range( 8 ): - newCell.GetPointIds().SetId( j, offset + c.GetPointId( m[ j ] ) ) - cells.InsertNextCell( newCell ) - offset += rg.GetNumberOfPoints() - - mesh = vtkUnstructuredGrid() - mesh.SetPoints( points ) - mesh.SetCells( cellTypes, cells ) - - return mesh - - -def __addFields( mesh: vtkUnstructuredGrid, fields: Iterable[ FieldInfo ] ) -> vtkUnstructuredGrid: - for fieldInfo in fields: - if fieldInfo.support == "CELLS": - data = mesh.GetCellData() - n = mesh.GetNumberOfCells() - elif fieldInfo.support == "POINTS": - data = mesh.GetPointData() - n = mesh.GetNumberOfPoints() - array = numpy.ones( ( n, fieldInfo.dimension ), dtype=float ) - vtkArray = numpy_to_vtk( array ) - vtkArray.SetName( fieldInfo.name ) - data.AddArray( vtkArray ) - return mesh - - -def __build( options: Options ): - - def buildCoordinates( positions, numElements ): - result = [] - it = zip( zip( positions, positions[ 1: ] ), numElements ) - try: - coords, n = next( it ) - while True: - start, stop = coords - endPoint = False - tmp = numpy.linspace( start=start, stop=stop, num=n + endPoint, endpoint=endPoint ) - coords, n = next( it ) - result.append( tmp ) - except StopIteration: - endPoint = True - tmp = numpy.linspace( start=start, stop=stop, num=n + endPoint, endpoint=endPoint ) - result.append( tmp ) - return numpy.concatenate( result ) - - x = buildCoordinates( options.xs, options.nxs ) - y = buildCoordinates( options.ys, options.nys ) - z = buildCoordinates( options.zs, options.nzs ) - cube = buildRectilinearBlocksMesh( ( XYZ( x, y, z ), ) ) - cube = __addFields( cube, options.fields ) - __buildGlobalIds( cube, options.generateCellsGlobalIds, options.generatePointsGlobalIds ) - return cube - - -def __action( options: Options ) -> Result: - outputMesh = __build( options ) - writeMesh( outputMesh, options.vtkOutput ) - return Result( info=f"Mesh was written to {options.vtkOutput.output}" ) - - -def action( vtkInputFile: str, options: Options ) -> Result: - try: - return __action( options ) - except BaseException as e: - setupLogger.error( e ) - return Result( info="Something went wrong." ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/mainChecks.py b/geos-mesh/src/geos/mesh/doctor/actions/mainChecks.py deleted file mode 100644 index f10cd49d..00000000 --- a/geos-mesh/src/geos/mesh/doctor/actions/mainChecks.py +++ /dev/null @@ -1 +0,0 @@ -from geos.mesh.doctor.actions.allChecks import action diff --git a/geos-mesh/src/geos/mesh/doctor/actions/selfIntersectingElements.py b/geos-mesh/src/geos/mesh/doctor/actions/selfIntersectingElements.py deleted file mode 100644 index 55956524..00000000 --- a/geos-mesh/src/geos/mesh/doctor/actions/selfIntersectingElements.py +++ /dev/null @@ -1,80 +0,0 @@ -from dataclasses import dataclass -from typing import Collection -from vtkmodules.util.numpy_support import vtk_to_numpy -from vtkmodules.vtkFiltersGeneral import vtkCellValidator -from vtkmodules.vtkCommonCore import vtkOutputWindow, vtkFileOutputWindow -from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid -from geos.mesh.io.vtkIO import readUnstructuredGrid - - -@dataclass( frozen=True ) -class Options: - minDistance: float - - -@dataclass( frozen=True ) -class Result: - wrongNumberOfPointsElements: Collection[ int ] - intersectingEdgesElements: Collection[ int ] - intersectingFacesElements: Collection[ int ] - nonContiguousEdgesElements: Collection[ int ] - nonConvexElements: Collection[ int ] - facesAreOrientedIncorrectlyElements: Collection[ int ] - - -def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: - errOut = vtkFileOutputWindow() - errOut.SetFileName( "/dev/null" ) # vtkCellValidator outputs loads for each cell... - vtkStdErrOut = vtkOutputWindow() - vtkStdErrOut.SetInstance( errOut ) - - valid = 0x0 - wrongNumberOfPoints = 0x01 - intersectingEdges = 0x02 - intersectingFaces = 0x04 - nonContiguousEdges = 0x08 - nonConvex = 0x10 - facesAreOrientedIncorrectly = 0x20 - - wrongNumberOfPointsElements: list[ int ] = [] - intersectingEdgesElements: list[ int ] = [] - intersectingFacesElements: list[ int ] = [] - nonContiguousEdgesElements: list[ int ] = [] - nonConvexElements: list[ int ] = [] - facesAreOrientedIncorrectlyElements: list[ int ] = [] - - f = vtkCellValidator() - f.SetTolerance( options.minDistance ) - - f.SetInputData( mesh ) - f.Update() - output = f.GetOutput() - - validity = output.GetCellData().GetArray( "ValidityState" ) # Could not change name using the vtk interface. - assert validity is not None - validity = vtk_to_numpy( validity ) - for i, v in enumerate( validity ): - if not v & valid: - if v & wrongNumberOfPoints: - wrongNumberOfPointsElements.append( i ) - if v & intersectingEdges: - intersectingEdgesElements.append( i ) - if v & intersectingFaces: - intersectingFacesElements.append( i ) - if v & nonContiguousEdges: - nonContiguousEdgesElements.append( i ) - if v & nonConvex: - nonConvexElements.append( i ) - if v & facesAreOrientedIncorrectly: - facesAreOrientedIncorrectlyElements.append( i ) - return Result( wrongNumberOfPointsElements=wrongNumberOfPointsElements, - intersectingEdgesElements=intersectingEdgesElements, - intersectingFacesElements=intersectingFacesElements, - nonContiguousEdgesElements=nonContiguousEdgesElements, - nonConvexElements=nonConvexElements, - facesAreOrientedIncorrectlyElements=facesAreOrientedIncorrectlyElements ) - - -def action( vtkInputFile: str, options: Options ) -> Result: - mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtkInputFile ) - return __action( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/supportedElements.py b/geos-mesh/src/geos/mesh/doctor/actions/supportedElements.py deleted file mode 100644 index e760c717..00000000 --- a/geos-mesh/src/geos/mesh/doctor/actions/supportedElements.py +++ /dev/null @@ -1,177 +0,0 @@ -from dataclasses import dataclass -import multiprocessing -import networkx -from tqdm import tqdm -from typing import Iterable, Mapping, Optional -from vtkmodules.util.numpy_support import vtk_to_numpy -from vtkmodules.vtkCommonCore import vtkIdList -from vtkmodules.vtkCommonDataModel import ( vtkCellTypes, vtkUnstructuredGrid, VTK_HEXAGONAL_PRISM, VTK_HEXAHEDRON, - VTK_PENTAGONAL_PRISM, VTK_POLYHEDRON, VTK_PYRAMID, VTK_TETRA, VTK_VOXEL, - VTK_WEDGE ) -from geos.mesh.doctor.actions.vtkPolyhedron import buildFaceToFaceConnectivityThroughEdges, FaceStream -from geos.mesh.doctor.parsing.cliParsing import setupLogger -from geos.mesh.io.vtkIO import readUnstructuredGrid -from geos.mesh.utils.genericHelpers import vtkIter - - -@dataclass( frozen=True ) -class Options: - nproc: int - chunkSize: int - - -@dataclass( frozen=True ) -class Result: - unsupportedStdElementsTypes: frozenset[ int ] # list of unsupported types - unsupportedPolyhedronElements: frozenset[ - int ] # list of polyhedron elements that could not be converted to supported std elements - - -# for multiprocessing, vtkUnstructuredGrid cannot be pickled. Let's use a global variable instead. -MESH: Optional[ vtkUnstructuredGrid ] = None - - -def initWorkerMesh( inputFileForWorker: str ): - """Initializer for multiprocessing.Pool to set the global MESH variable in each worker process. - - Args: - inputFileForWorker (str): Filepath to vtk grid - """ - global MESH - setupLogger.debug( - f"Worker process (PID: {multiprocessing.current_process().pid}) initializing MESH from file: {inputFileForWorker}" - ) - mesh: vtkUnstructuredGrid = readUnstructuredGrid( inputFileForWorker ) - if MESH is None: - setupLogger.error( - f"Worker process (PID: {multiprocessing.current_process().pid}) failed to load mesh from {inputFileForWorker}" - ) - # You might want to raise an error here or ensure MESH being None is handled downstream - # For now, the assert MESH is not None in __call__ will catch this. - - -class IsPolyhedronConvertible: - - def __init__( self ): - - def buildPrismGraph( n: int, name: str ) -> networkx.Graph: - """Builds the face to face connectivities (through edges) for prism graphs. - - Args: - n (int): The number of nodes of the basis (i.e. the pentagonal prims gets n = 5) - name (str): A human-readable name for logging purpose. - - Returns: - networkx.Graph: A graph instance. - """ - tmp = networkx.cycle_graph( n ) - for node in range( n ): - tmp.add_edge( node, n ) - tmp.add_edge( node, n + 1 ) - tmp.name = name - return tmp - - # Building the reference graphs - tetGraph = networkx.complete_graph( 4 ) - tetGraph.name = "Tetrahedron" - pyrGraph = buildPrismGraph( 4, "Pyramid" ) - pyrGraph.remove_node( 5 ) # Removing a node also removes its associated edges. - self.__referenceGraphs: Mapping[ int, Iterable[ networkx.Graph ] ] = { - 4: ( tetGraph, ), - 5: ( pyrGraph, buildPrismGraph( 3, "Wedge" ) ), - 6: ( buildPrismGraph( 4, "Hexahedron" ), ), - 7: ( buildPrismGraph( 5, "Prism5" ), ), - 8: ( buildPrismGraph( 6, "Prism6" ), ), - 9: ( buildPrismGraph( 7, "Prism7" ), ), - 10: ( buildPrismGraph( 8, "Prism8" ), ), - 11: ( buildPrismGraph( 9, "Prism9" ), ), - 12: ( buildPrismGraph( 10, "Prism10" ), ), - 13: ( buildPrismGraph( 11, "Prism11" ), ), - } - - def __isPolyhedronSupported( self, faceStream ) -> str: - """Checks if a polyhedron can be converted into a supported cell. - If so, returns the name of the type. If not, the returned name will be empty. - - Args: - faceStream (_type_): The polyhedron. - - Returns: - str: The name of the supported type or an empty string. - """ - cellGraph = buildFaceToFaceConnectivityThroughEdges( faceStream, add_compatibility=True ) - if cellGraph.order() not in self.__referenceGraphs: - return "" - for referenceGraph in self.__referenceGraphs[ cellGraph.order() ]: - if networkx.is_isomorphic( referenceGraph, cellGraph ): - return str( referenceGraph.name ) - return "" - - def __call__( self, ic: int ) -> int: - """Checks if a vtk polyhedron cell can be converted into a supported GEOS element. - - Args: - ic (int): The index element. - - Returns: - int: -1 if the polyhedron vtk element can be converted into a supported element type. The index otherwise. - """ - global MESH - assert MESH is not None, f"MESH global variable not initialized in worker process (PID: {multiprocessing.current_process().pid}). This should have been set by initWorkerMesh." - if MESH.GetCellType( ic ) != VTK_POLYHEDRON: - return -1 - ptIds = vtkIdList() - MESH.GetFaceStream( ic, ptIds ) - faceStream = FaceStream.buildFromVtkIdList( ptIds ) - convertedTypeName = self.__isPolyhedronSupported( faceStream ) - if convertedTypeName: - setupLogger.debug( f"Polyhedron cell {ic} can be converted into \"{convertedTypeName}\"" ) - return -1 - else: - setupLogger.debug( - f"Polyhedron cell {ic} (in PID {multiprocessing.current_process().pid}) cannot be converted into any supported element." - ) - return ic - - -def __action( vtkInputFile: str, options: Options ) -> Result: - # Main process loads the mesh for its own use - mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtkInputFile ) - if mesh is None: - setupLogger.error( f"Main process failed to load mesh from {vtkInputFile}. Aborting." ) - # Return an empty/error result or raise an exception - return Result( unsupportedStdElementsTypes=frozenset(), unsupportedPolyhedronElements=frozenset() ) - - if hasattr( mesh, "GetDistinctCellTypesArray" ): - cellTypesNumpy = vtk_to_numpy( mesh.GetDistinctCellTypesArray() ) - cellTypes = set( cellTypesNumpy.tolist() ) - else: - vtkCellTypesObj = vtkCellTypes() - mesh.GetCellTypes( vtkCellTypesObj ) - cellTypes = set( vtkIter( vtkCellTypesObj ) ) - - supportedCellTypes = { - VTK_HEXAGONAL_PRISM, VTK_HEXAHEDRON, VTK_PENTAGONAL_PRISM, VTK_POLYHEDRON, VTK_PYRAMID, VTK_TETRA, VTK_VOXEL, - VTK_WEDGE - } - unsupportedStdElementsTypes = cellTypes - supportedCellTypes - - # Dealing with polyhedron elements. - numCells = mesh.GetNumberOfCells() - polyhedronConverter = IsPolyhedronConvertible() - - unsupportedPolyhedronIndices = [] - # Pass the vtkInputFile to the initializer - with multiprocessing.Pool( processes=options.nproc, initializer=initWorkerMesh, - initargs=( vtkInputFile, ) ) as pool: # Comma makes it a tuple - generator = pool.imap_unordered( polyhedronConverter, range( numCells ), chunksize=options.chunkSize ) - for cellIndexOrNegOne in tqdm( generator, total=numCells, desc="Testing support for elements" ): - if cellIndexOrNegOne != -1: - unsupportedPolyhedronIndices.append( cellIndexOrNegOne ) - - return Result( unsupportedStdElementsTypes=frozenset( unsupportedStdElementsTypes ), - unsupportedPolyhedronElements=frozenset( unsupportedPolyhedronIndices ) ) - - -def action( vtkInputFile: str, options: Options ) -> Result: - return __action( vtkInputFile, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/meshDoctor.py b/geos-mesh/src/geos/mesh/doctor/meshDoctor.py deleted file mode 100644 index ab05af4b..00000000 --- a/geos-mesh/src/geos/mesh/doctor/meshDoctor.py +++ /dev/null @@ -1,24 +0,0 @@ -import sys -from geos.mesh.doctor.parsing import ActionHelper -from geos.mesh.doctor.parsing.cliParsing import parseAndSetVerbosity, setupLogger -from geos.mesh.doctor.register import registerParsingActions - - -def main(): - parseAndSetVerbosity( sys.argv ) - mainParser, allActions, allActionsHelpers = registerParsingActions() - args = mainParser.parse_args( sys.argv[ 1: ] ) - setupLogger.info( f"Working on mesh \"{args.vtkInputFile}\"." ) - actionOptions = allActionsHelpers[ args.subparsers ].convert( vars( args ) ) - try: - action = allActions[ args.subparsers ] - except KeyError: - setupLogger.error( f"Action {args.subparsers} is not a valid action." ) - sys.exit( 1 ) - helper: ActionHelper = allActionsHelpers[ args.subparsers ] - result = action( args.vtkInputFile, actionOptions ) - helper.displayResults( actionOptions, result ) - - -if __name__ == '__main__': - main() diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/checkFracturesParsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/checkFracturesParsing.py deleted file mode 100644 index 3f43bca6..00000000 --- a/geos-mesh/src/geos/mesh/doctor/parsing/checkFracturesParsing.py +++ /dev/null @@ -1 +0,0 @@ -# empty: the check is not available yet! \ No newline at end of file diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/collocatedNodesParsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/collocatedNodesParsing.py deleted file mode 100644 index 1b48c6f3..00000000 --- a/geos-mesh/src/geos/mesh/doctor/parsing/collocatedNodesParsing.py +++ /dev/null @@ -1,48 +0,0 @@ -from geos.mesh.doctor.actions.collocatedNodes import Options, Result -from geos.mesh.doctor.parsing import COLLOCATES_NODES -from geos.mesh.doctor.parsing._sharedChecksParsingLogic import getOptionsUsedMessage -from geos.mesh.doctor.parsing.cliParsing import setupLogger - -__TOLERANCE = "tolerance" -__TOLERANCE_DEFAULT = 0. - -__COLLOCATED_NODES_DEFAULT = { __TOLERANCE: __TOLERANCE_DEFAULT } - - -def convert( parsedOptions ) -> Options: - return Options( parsedOptions[ __TOLERANCE ] ) - - -def fillSubparser( subparsers ) -> None: - p = subparsers.add_parser( COLLOCATES_NODES, help="Checks if nodes are collocated." ) - p.add_argument( '--' + __TOLERANCE, - type=float, - metavar=__TOLERANCE_DEFAULT, - default=__TOLERANCE_DEFAULT, - required=True, - help="[float]: The absolute distance between two nodes for them to be considered collocated." ) - - -def displayResults( options: Options, result: Result ): - setupLogger.results( getOptionsUsedMessage( options ) ) - allCollocatedNodes: list[ int ] = [] - for bucket in result.nodesBuckets: - for node in bucket: - allCollocatedNodes.append( node ) - allCollocatedNodes: frozenset[ int ] = frozenset( allCollocatedNodes ) # Surely useless - if allCollocatedNodes: - setupLogger.results( f"You have {len( allCollocatedNodes )} collocated nodes." ) - setupLogger.results( "Here are all the buckets of collocated nodes." ) - tmp: list[ str ] = [] - for bucket in result.nodesBuckets: - tmp.append( f"({', '.join(map(str, bucket))})" ) - setupLogger.results( f"({', '.join(tmp)})" ) - else: - setupLogger.results( "You have no collocated node." ) - - if result.wrongSupportElements: - tmp: str = ", ".join( map( str, result.wrongSupportElements ) ) - setupLogger.results( f"You have {len(result.wrongSupportElements)} elements with duplicated support nodes.\n" + - tmp ) - else: - setupLogger.results( "You have no element with duplicated support nodes." ) diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/elementVolumesParsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/elementVolumesParsing.py deleted file mode 100644 index 102fbebd..00000000 --- a/geos-mesh/src/geos/mesh/doctor/parsing/elementVolumesParsing.py +++ /dev/null @@ -1,44 +0,0 @@ -from geos.mesh.doctor.actions.elementVolumes import Options, Result -from geos.mesh.doctor.parsing import ELEMENT_VOLUMES -from geos.mesh.doctor.parsing._sharedChecksParsingLogic import getOptionsUsedMessage -from geos.mesh.doctor.parsing.cliParsing import setupLogger - -__MIN_VOLUME = "minVolume" -__MIN_VOLUME_DEFAULT = 0. - -__ELEMENT_VOLUMES_DEFAULT = { __MIN_VOLUME: __MIN_VOLUME_DEFAULT } - - -def fillSubparser( subparsers ) -> None: - p = subparsers.add_parser( ELEMENT_VOLUMES, - help=f"Checks if the volumes of the elements are greater than \"{__MIN_VOLUME}\"." ) - p.add_argument( '--' + __MIN_VOLUME, - type=float, - metavar=__MIN_VOLUME_DEFAULT, - default=__MIN_VOLUME_DEFAULT, - required=True, - help=f"[float]: The minimum acceptable volume. Defaults to {__MIN_VOLUME_DEFAULT}." ) - - -def convert( parsedOptions ) -> Options: - """From the parsed cli options, return the converted options for elements volumes check. - - Args: - parsedOptions: Parsed cli options. - - Returns: - Options: The converted options for elements volumes check. - """ - return Options( minVolume=parsedOptions[ __MIN_VOLUME ] ) - - -def displayResults( options: Options, result: Result ): - setupLogger.results( getOptionsUsedMessage( options ) ) - setupLogger.results( - f"You have {len(result.elementVolumes)} elements with volumes smaller than {options.minVolume}." ) - if result.elementVolumes: - setupLogger.results( "Elements index | Volumes calculated" ) - setupLogger.results( "-----------------------------------" ) - maxLength: int = len( "Elements index " ) - for ( ind, volume ) in result.elementVolumes: - setupLogger.results( f"{ind:<{maxLength}}" + "| " + str( volume ) ) diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/generateGlobalIdsParsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/generateGlobalIdsParsing.py deleted file mode 100644 index e7e934a5..00000000 --- a/geos-mesh/src/geos/mesh/doctor/parsing/generateGlobalIdsParsing.py +++ /dev/null @@ -1,52 +0,0 @@ -from dataclasses import dataclass -from geos.mesh.doctor.actions.generateGlobalIds import Options, Result -from geos.mesh.doctor.parsing import vtkOutputParsing, GENERATE_GLOBAL_IDS -from geos.mesh.doctor.parsing.cliParsing import setupLogger - -__CELLS, __POINTS = "cells", "points" - - -@dataclass( frozen=True ) -class GlobalIdsInfo: - cells: bool - points: bool - - -def convertGlobalIds( parsedOptions ) -> GlobalIdsInfo: - return GlobalIdsInfo( cells=parsedOptions[ __CELLS ], points=parsedOptions[ __POINTS ] ) - - -def convert( parsedOptions ) -> Options: - gids: GlobalIdsInfo = convertGlobalIds( parsedOptions ) - return Options( vtkOutput=vtkOutputParsing.convert( parsedOptions ), - generateCellsGlobalIds=gids.cells, - generatePointsGlobalIds=gids.points ) - - -def fillGenerateGlobalIdsSubparser( p ): - p.add_argument( '--' + __CELLS, - action="store_true", - help="[bool]: Generate global ids for cells. Defaults to true." ) - p.add_argument( '--no-' + __CELLS, - action="store_false", - dest=__CELLS, - help="[bool]: Don't generate global ids for cells." ) - p.set_defaults( **{ __CELLS: True } ) - p.add_argument( '--' + __POINTS, - action="store_true", - help="[bool]: Generate global ids for points. Defaults to true." ) - p.add_argument( '--no-' + __POINTS, - action="store_false", - dest=__POINTS, - help="[bool]: Don't generate global ids for points." ) - p.set_defaults( **{ __POINTS: True } ) - - -def fillSubparser( subparsers ) -> None: - p = subparsers.add_parser( GENERATE_GLOBAL_IDS, help="Adds globals ids for points and cells." ) - fillGenerateGlobalIdsSubparser( p ) - vtkOutputParsing.fillVtkOutputSubparser( p ) - - -def displayResults( options: Options, result: Result ): - setupLogger.info( result.info ) diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/nonConformalParsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/nonConformalParsing.py deleted file mode 100644 index b57ca504..00000000 --- a/geos-mesh/src/geos/mesh/doctor/parsing/nonConformalParsing.py +++ /dev/null @@ -1,55 +0,0 @@ -from geos.mesh.doctor.actions.nonConformal import Options, Result -from geos.mesh.doctor.parsing import NON_CONFORMAL -from geos.mesh.doctor.parsing._sharedChecksParsingLogic import getOptionsUsedMessage -from geos.mesh.doctor.parsing.cliParsing import setupLogger - -__ANGLE_TOLERANCE = "angleTolerance" -__POINT_TOLERANCE = "pointTolerance" -__FACE_TOLERANCE = "faceTolerance" - -__ANGLE_TOLERANCE_DEFAULT = 10. -__POINT_TOLERANCE_DEFAULT = 0. -__FACE_TOLERANCE_DEFAULT = 0. - -__NON_CONFORMAL_DEFAULT = { - __ANGLE_TOLERANCE: __ANGLE_TOLERANCE_DEFAULT, - __POINT_TOLERANCE: __POINT_TOLERANCE_DEFAULT, - __FACE_TOLERANCE: __FACE_TOLERANCE_DEFAULT -} - - -def convert( parsedOptions ) -> Options: - return Options( angleTolerance=parsedOptions[ __ANGLE_TOLERANCE ], - pointTolerance=parsedOptions[ __POINT_TOLERANCE ], - faceTolerance=parsedOptions[ __FACE_TOLERANCE ] ) - - -def fillSubparser( subparsers ) -> None: - p = subparsers.add_parser( NON_CONFORMAL, help="Detects non conformal elements. [EXPERIMENTAL]" ) - p.add_argument( '--' + __ANGLE_TOLERANCE, - type=float, - metavar=__ANGLE_TOLERANCE_DEFAULT, - default=__ANGLE_TOLERANCE_DEFAULT, - help=f"[float]: angle tolerance in degrees. Defaults to {__ANGLE_TOLERANCE_DEFAULT}" ) - p.add_argument( - '--' + __POINT_TOLERANCE, - type=float, - metavar=__POINT_TOLERANCE_DEFAULT, - default=__POINT_TOLERANCE_DEFAULT, - help=f"[float]: tolerance for two points to be considered collocated. Defaults to {__POINT_TOLERANCE_DEFAULT}" ) - p.add_argument( - '--' + __FACE_TOLERANCE, - type=float, - metavar=__FACE_TOLERANCE_DEFAULT, - default=__FACE_TOLERANCE_DEFAULT, - help=f"[float]: tolerance for two faces to be considered \"touching\". Defaults to {__FACE_TOLERANCE_DEFAULT}" ) - - -def displayResults( options: Options, result: Result ): - setupLogger.results( getOptionsUsedMessage( options ) ) - nonConformalCells: list[ int ] = [] - for i, j in result.nonConformalCells: - nonConformalCells += i, j - nonConformalCells: frozenset[ int ] = frozenset( nonConformalCells ) - setupLogger.results( f"You have {len( nonConformalCells )} non conformal cells." ) - setupLogger.results( f"{', '.join( map( str, sorted( nonConformalCells ) ) )}" ) diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/selfIntersectingElementsParsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/selfIntersectingElementsParsing.py deleted file mode 100644 index a61e1fbd..00000000 --- a/geos-mesh/src/geos/mesh/doctor/parsing/selfIntersectingElementsParsing.py +++ /dev/null @@ -1,43 +0,0 @@ -import numpy -from geos.mesh.doctor.actions.selfIntersectingElements import Options, Result -from geos.mesh.doctor.parsing import SELF_INTERSECTING_ELEMENTS -from geos.mesh.doctor.parsing._sharedChecksParsingLogic import getOptionsUsedMessage -from geos.mesh.doctor.parsing.cliParsing import setupLogger - -__MIN_DISTANCE = "minDistance" -__MIN_DISTANCE_DEFAULT = numpy.finfo( float ).eps - -__SELF_INTERSECTING_ELEMENTS_DEFAULT = { __MIN_DISTANCE: __MIN_DISTANCE_DEFAULT } - - -def convert( parsedOptions ) -> Options: - minDistance = parsedOptions[ __MIN_DISTANCE ] - if minDistance == 0: - setupLogger.warning( - "Having minimum distance set to 0 can induce lots of false positive results (adjacent faces may be considered intersecting)." - ) - elif minDistance < 0: - raise ValueError( - f"Negative minimum distance ({minDistance}) in the {SELF_INTERSECTING_ELEMENTS} check is not allowed." ) - return Options( minDistance=minDistance ) - - -def fillSubparser( subparsers ) -> None: - p = subparsers.add_parser( SELF_INTERSECTING_ELEMENTS, - help="Checks if the faces of the elements are self intersecting." ) - p.add_argument( - '--' + __MIN_DISTANCE, - type=float, - required=False, - metavar=__MIN_DISTANCE_DEFAULT, - default=__MIN_DISTANCE_DEFAULT, - help= - f"[float]: The minimum distance in the computation. Defaults to your machine precision {__MIN_DISTANCE_DEFAULT}." - ) - - -def displayResults( options: Options, result: Result ): - setupLogger.results( getOptionsUsedMessage( options ) ) - setupLogger.results( f"You have {len(result.intersectingFacesElements)} elements with self intersecting faces." ) - if result.intersectingFacesElements: - setupLogger.results( "The elements indices are:\n" + ", ".join( map( str, result.intersectingFacesElements ) ) ) diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/supportedElementsParsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/supportedElementsParsing.py deleted file mode 100644 index 4ed3b1ef..00000000 --- a/geos-mesh/src/geos/mesh/doctor/parsing/supportedElementsParsing.py +++ /dev/null @@ -1,52 +0,0 @@ -import multiprocessing -from geos.mesh.doctor.actions.supportedElements import Options, Result -from geos.mesh.doctor.parsing import SUPPORTED_ELEMENTS -from geos.mesh.doctor.parsing._sharedChecksParsingLogic import getOptionsUsedMessage -from geos.mesh.doctor.parsing.cliParsing import setupLogger - -__CHUNK_SIZE = "chunkSize" -__NUM_PROC = "nproc" - -__CHUNK_SIZE_DEFAULT = 1 -__NUM_PROC_DEFAULT = multiprocessing.cpu_count() - -__SUPPORTED_ELEMENTS_DEFAULT = { __CHUNK_SIZE: __CHUNK_SIZE_DEFAULT, __NUM_PROC: __NUM_PROC_DEFAULT } - - -def convert( parsedOptions ) -> Options: - return Options( chunkSize=parsedOptions[ __CHUNK_SIZE ], nproc=parsedOptions[ __NUM_PROC ] ) - - -def fillSubparser( subparsers ) -> None: - p = subparsers.add_parser( SUPPORTED_ELEMENTS, - help="Check that all the elements of the mesh are supported by GEOSX." ) - p.add_argument( '--' + __CHUNK_SIZE, - type=int, - required=False, - metavar=__CHUNK_SIZE_DEFAULT, - default=__CHUNK_SIZE_DEFAULT, - help=f"[int]: Defaults chunk size for parallel processing to {__CHUNK_SIZE_DEFAULT}" ) - p.add_argument( - '--' + __NUM_PROC, - type=int, - required=False, - metavar=__NUM_PROC_DEFAULT, - default=__NUM_PROC_DEFAULT, - help=f"[int]: Number of threads used for parallel processing. Defaults to your CPU count {__NUM_PROC_DEFAULT}." - ) - - -def displayResults( options: Options, result: Result ): - setupLogger.results( getOptionsUsedMessage( options ) ) - if result.unsupportedPolyhedronElements: - setupLogger.results( f"There is/are {len(result.unsupportedPolyhedronElements)} polyhedra that may not be " - f"converted to supported elements." ) - setupLogger.results( - f"The list of the unsupported polyhedra is\n{tuple(sorted(result.unsupportedPolyhedronElements))}." ) - else: - setupLogger.results( "All the polyhedra (if any) can be converted to supported elements." ) - if result.unsupportedStdElementsTypes: - setupLogger.results( f"There are unsupported vtk standard element types. The list of those vtk types is " - f"{tuple(sorted(result.unsupportedStdElementsTypes))}." ) - else: - setupLogger.results( "All the standard vtk element types (if any) are supported." ) diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/vtkOutputParsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/vtkOutputParsing.py deleted file mode 100644 index 1431ce4e..00000000 --- a/geos-mesh/src/geos/mesh/doctor/parsing/vtkOutputParsing.py +++ /dev/null @@ -1,45 +0,0 @@ -import os.path -import textwrap -from geos.mesh.doctor.parsing.cliParsing import setupLogger -from geos.mesh.io.vtkIO import VtkOutput - -__OUTPUT_FILE = "output" -__OUTPUT_BINARY_MODE = "data-mode" -__OUTPUT_BINARY_MODE_VALUES = "binary", "ascii" -__OUTPUT_BINARY_MODE_DEFAULT = __OUTPUT_BINARY_MODE_VALUES[ 0 ] - - -def getVtkOutputHelp(): - msg = \ - f"""{__OUTPUT_FILE} [string]: The vtk output file destination. - {__OUTPUT_BINARY_MODE} [string]: For ".vtu" output format, the data mode can be {" or ".join(__OUTPUT_BINARY_MODE_VALUES)}. Defaults to {__OUTPUT_BINARY_MODE_DEFAULT}.""" - return textwrap.dedent( msg ) - - -def __buildArg( prefix, main ): - return "-".join( filter( None, ( prefix, main ) ) ) - - -def fillVtkOutputSubparser( parser, prefix="" ) -> None: - parser.add_argument( '--' + __buildArg( prefix, __OUTPUT_FILE ), - type=str, - required=True, - help=f"[string]: The vtk output file destination." ) - parser.add_argument( - '--' + __buildArg( prefix, __OUTPUT_BINARY_MODE ), - type=str, - metavar=", ".join( __OUTPUT_BINARY_MODE_VALUES ), - default=__OUTPUT_BINARY_MODE_DEFAULT, - help= - f"""[string]: For ".vtu" output format, the data mode can be {" or ".join(__OUTPUT_BINARY_MODE_VALUES)}. Defaults to {__OUTPUT_BINARY_MODE_DEFAULT}.""" - ) - - -def convert( parsedOptions, prefix="" ) -> VtkOutput: - outputKey = __buildArg( prefix, __OUTPUT_FILE ).replace( "-", "_" ) - binaryModeKey = __buildArg( prefix, __OUTPUT_BINARY_MODE ).replace( "-", "_" ) - output = parsedOptions[ outputKey ] - if parsedOptions[ binaryModeKey ] and os.path.splitext( output )[ -1 ] == ".vtk": - setupLogger.info( "VTK data mode will be ignored for legacy file format \"vtk\"." ) - isDataModeBinary: bool = parsedOptions[ binaryModeKey ] == __OUTPUT_BINARY_MODE_DEFAULT - return VtkOutput( output=output, isDataModeBinary=isDataModeBinary ) diff --git a/install_packages.sh b/install_packages.sh index bc6b6a13..f0384e5d 100755 --- a/install_packages.sh +++ b/install_packages.sh @@ -7,6 +7,7 @@ python -m pip install --upgrade ./geos-xml-tools python -m pip install --upgrade ./geos-xml-viewer python -m pip install --upgrade ./hdf5-wrapper python -m pip install --upgrade ./geos-timehistory +python -m pip install --upgrade ./mesh-doctor python -m pip install --upgrade ./pygeos-tools python -m pip install --upgrade ./geos-pv #! trame install requires npm diff --git a/mesh-doctor/pyproject.toml b/mesh-doctor/pyproject.toml new file mode 100644 index 00000000..06b94e7c --- /dev/null +++ b/mesh-doctor/pyproject.toml @@ -0,0 +1,67 @@ +[build-system] +requires = ["setuptools>=61.2", "wheel >= 0.37.1"] +build-backend = "setuptools.build_meta" + +[tool.setuptools] +include-package-data = true + +[tool.setuptools.packages.find] +where = ["src"] +include = ["geos.mesh_doctor*"] +exclude = ['tests*'] + +[project] +name = "mesh-doctor" +version = "0.1.0" +description = "GEOS mesh-doctor" +authors = [{name = "GEOS Contributors" }] +maintainers = [{name = "Alexandre Benedicto", email = "alexandre.benedicto@external.totalenergies.com" }, + {name = "Romain Baville", email = "romain.baville@external.totalenergies.com" }, + {name = "Paloma Martinez", email = "paloma.martinez@external.totalenergies.com" }] +license = {text = "Apache-2.0"} +classifiers = [ + "Development Status :: 4 - Beta", + "Programming Language :: Python" +] + +requires-python = ">=3.10" + +dependencies = [ + "geos-utils", + "geos-mesh", +] + +[project.scripts] + mesh-doctor = "geos.mesh_doctor.meshDoctor:main" + meshDoctor = "geos.mesh_doctor.meshDoctor:main" + +[project.urls] +Homepage = "https://github.com/GEOS-DEV/geosPythonPackages" +Documentation = "https://geosx-geosx.readthedocs-hosted.com/projects/geosx-geospythonpackages/en/latest/" +Repository = "https://github.com/GEOS-DEV/geosPythonPackages.git" +"Bug Tracker" = "https://github.com/GEOS-DEV/geosPythonPackages/issues" + +[project.optional-dependencies] +build = [ + "build ~= 1.2" +] +dev = [ + "mypy", + "ruff", + "yapf", +] +test = [ + "pytest-cov", + "pytest" +] + +[tool.pytest.ini_options] +addopts = "--import-mode=importlib" +console_output_style = "count" +pythonpath = ["src"] +python_classes = "Test" +python_files = "test*.py" +python_functions = "test*" +testpaths = ["tests"] +norecursedirs = "bin" +filterwarnings = [] diff --git a/geos-mesh/src/geos/mesh/doctor/actions/__init__.py b/mesh-doctor/src/geos/mesh_doctor/__init__.py similarity index 100% rename from geos-mesh/src/geos/mesh/doctor/actions/__init__.py rename to mesh-doctor/src/geos/mesh_doctor/__init__.py diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/__init__.py b/mesh-doctor/src/geos/mesh_doctor/actions/__init__.py new file mode 100644 index 00000000..b7db2541 --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/actions/__init__.py @@ -0,0 +1 @@ +# Empty diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/allChecks.py b/mesh-doctor/src/geos/mesh_doctor/actions/allChecks.py new file mode 100644 index 00000000..c193f63e --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/actions/allChecks.py @@ -0,0 +1,55 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from dataclasses import dataclass +from typing import Any +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from geos.mesh_doctor.parsing.cliParsing import setupLogger +from geos.mesh_doctor.register import __loadModuleAction +from geos.mesh.io.vtkIO import readUnstructuredGrid + + +@dataclass( frozen=True ) +class Options: + checksToPerform: list[ str ] + checksOptions: dict[ str, Any ] + checkDisplays: dict[ str, Any ] + + +@dataclass( frozen=True ) +class Result: + checkResults: dict[ str, Any ] + + +def meshAction( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + """Performs all checks available on a loaded mesh. + + Args: + mesh (vtkUnstructuredGrid): The loaded mesh to analyze. + options (Options): The options for processing. + + Returns: + Result: The result of all checks performed. + """ + checkResults: dict[ str, Any ] = {} + for checkName in options.checksToPerform: + checkMeshAction = __loadModuleAction( checkName, "meshAction" ) + setupLogger.info( f"Performing check '{checkName}'." ) + option = options.checksOptions[ checkName ] + checkResult = checkMeshAction( mesh, option ) + checkResults[ checkName ] = checkResult + return Result( checkResults=checkResults ) + + +def action( vtuInputFile: str, options: Options ) -> Result: + """Reads a vtu file and performs all checks available on it. + + Args: + vtuInputFile (str): The path to the input VTU file. + options (Options): The options for processing. + + Returns: + Result: The result of all checks performed. + """ + mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtuInputFile ) + return meshAction( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/checkFractures.py b/mesh-doctor/src/geos/mesh_doctor/actions/checkFractures.py similarity index 52% rename from geos-mesh/src/geos/mesh/doctor/actions/checkFractures.py rename to mesh-doctor/src/geos/mesh_doctor/actions/checkFractures.py index a7f303fc..50799f4d 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/checkFractures.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/checkFractures.py @@ -1,14 +1,18 @@ -import numpy +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto from dataclasses import dataclass +import numpy as np +import numpy.typing as npt from tqdm import tqdm from typing import Collection, Iterable, Sequence -from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkCell +from vtkmodules.vtkCommonDataModel import vtkCell, vtkMultiBlockDataSet, vtkUnstructuredGrid from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkIOXML import vtkXMLMultiBlockDataReader from vtkmodules.util.numpy_support import vtk_to_numpy -from geos.mesh.doctor.actions.generateFractures import Coordinates3D -from geos.mesh.doctor.parsing.cliParsing import setupLogger from geos.mesh.utils.genericHelpers import vtkIter +from geos.mesh_doctor.actions.generateFractures import Coordinates3D +from geos.mesh_doctor.parsing.cliParsing import setupLogger @dataclass( frozen=True ) @@ -27,18 +31,28 @@ class Result: errors: Sequence[ tuple[ int, int, int ] ] -def __readMultiblock( vtkInputFile: str, matrixName: str, +def __readMultiblock( vtuInputFile: str, matrixName: str, fractureName: str ) -> tuple[ vtkUnstructuredGrid, vtkUnstructuredGrid ]: + """Reads a multiblock VTU file and extracts the matrix and fracture meshes. + + Args: + vtuInputFile (str): The input VTU file path. + matrixName (str): The name of the matrix mesh block. + fractureName (str): The name of the fracture mesh block. + + Returns: + tuple[ vtkUnstructuredGrid, vtkUnstructuredGrid ]: The matrix and fracture meshes. + """ reader = vtkXMLMultiBlockDataReader() - reader.SetFileName( vtkInputFile ) + reader.SetFileName( vtuInputFile ) reader.Update() - multiBlock = reader.GetOutput() + multiBlock: vtkMultiBlockDataSet = reader.GetOutput() for b in range( multiBlock.GetNumberOfBlocks() ): blockName: str = multiBlock.GetMetaData( b ).Get( multiBlock.NAME() ) if blockName == matrixName: - matrix: vtkUnstructuredGrid = multiBlock.GetBlock( b ) + matrix = vtkUnstructuredGrid.SafeDownCast( multiBlock.GetBlock( b ) ) if blockName == fractureName: - fracture: vtkUnstructuredGrid = multiBlock.GetBlock( b ) + fracture = vtkUnstructuredGrid.SafeDownCast( multiBlock.GetBlock( b ) ) assert matrix and fracture return matrix, fracture @@ -52,28 +66,47 @@ def formatCollocatedNodes( fractureMesh: vtkUnstructuredGrid ) -> Sequence[ Iter Returns: Sequence[ Iterable[ int ] ]: An iterable over all the buckets of collocated nodes. """ - collocatedNodes: numpy.ndarray = vtk_to_numpy( fractureMesh.GetPointData().GetArray( "collocated_nodes" ) ) + collocatedNodes: npt.NDArray = vtk_to_numpy( fractureMesh.GetPointData().GetArray( "collocated_nodes" ) ) if len( collocatedNodes.shape ) == 1: - collocatedNodes: numpy.ndarray = collocatedNodes.reshape( ( collocatedNodes.shape[ 0 ], 1 ) ) + collocatedNodes = collocatedNodes.reshape( ( collocatedNodes.shape[ 0 ], 1 ) ) generator = ( tuple( sorted( bucket[ bucket > -1 ] ) ) for bucket in collocatedNodes ) return tuple( generator ) def __checkCollocatedNodesPositions( - matrixPoints: Sequence[ Coordinates3D ], fracturePoints: Sequence[ Coordinates3D ], g2l: Sequence[ int ], + matrixPoints: npt.NDArray, fracturePoints: npt.NDArray, g2l: npt.NDArray[ np.int64 ], collocatedNodes: Iterable[ Iterable[ int ] ] ) -> Collection[ tuple[ int, Iterable[ int ], Iterable[ Coordinates3D ] ] ]: + """Check that the collocated nodes are indeed collocated in space. + + Args: + matrixPoints (npt.NDArray): The points of the matrix mesh. + fracturePoints (npt.NDArray): The points of the fracture mesh. + g2l (npt.NDArray[ np.int64 ]): Mapping from global to local indices in the matrix mesh. + collocatedNodes (Iterable[ Iterable[ int ] ]): The collocated nodes information. + + Returns: + Collection[ tuple[ int, Iterable[ int ], Iterable[ Coordinates3D ] ] ]: A collection of issues found. + """ issues = [] for li, bucket in enumerate( collocatedNodes ): - matrix_nodes = ( fracturePoints[ li ], ) + tuple( map( lambda gi: matrixPoints[ g2l[ gi ] ], bucket ) ) - m = numpy.array( matrix_nodes ) - rank: int = numpy.linalg.matrix_rank( m ) + matrix_nodes = ( fracturePoints[ li ], ) + tuple( matrixPoints[ g2l[ gi ] ] for gi in bucket ) + m = np.array( matrix_nodes ) + rank: int = np.linalg.matrix_rank( m ) if rank > 1: - issues.append( ( li, bucket, tuple( map( lambda gi: matrixPoints[ g2l[ gi ] ], bucket ) ) ) ) + issues.append( ( li, bucket, tuple( matrixPoints[ g2l[ gi ] ] for gi in bucket ) ) ) return issues -def myIter( ccc ): +def myIter( ccc: tuple ) -> Iterable[ tuple ]: + """Recursively generates all combinations from nested sequences. + + Args: + ccc (tuple): A tuple of sequences to generate combinations from. + + Yields: + Iterable[ tuple ]: All possible combinations, one element from each sequence. + """ car, cdr = ccc[ 0 ], ccc[ 1: ] for i in car: if cdr: @@ -83,8 +116,16 @@ def myIter( ccc ): yield ( i, ) -def __checkNeighbors( matrix: vtkUnstructuredGrid, fracture: vtkUnstructuredGrid, g2l: Sequence[ int ], - collocatedNodes: Sequence[ Iterable[ int ] ] ): +def __checkNeighbors( matrix: vtkUnstructuredGrid, fracture: vtkUnstructuredGrid, g2l: npt.NDArray[ np.int64 ], + collocatedNodes: Sequence[ Iterable[ int ] ] ) -> None: + """Check that for each pair of collocated nodes, the corresponding fracture faces have exactly two neighbor cells. + + Args: + matrix (vtkUnstructuredGrid): The matrix mesh. + fracture (vtkUnstructuredGrid): The fracture mesh. + g2l (npt.NDArray[ np.int64 ]): Mapping from global to local indices in the matrix mesh. + collocatedNodes (Sequence[ Iterable[ int ] ]): The collocated nodes information. + """ fractureNodes: set[ int ] = set() for bucket in collocatedNodes: for gi in bucket: @@ -102,8 +143,8 @@ def __checkNeighbors( matrix: vtkUnstructuredGrid, fracture: vtkUnstructuredGrid fractureFaces.add( pointIds ) # 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 collocatedNodes + cell = fracture.GetCell( c ) + cns: set[ frozenset[ npt.NDArray[ np.int64 ] ] ] = set() # subset of collocatedNodes pointIds = frozenset( vtkIter( cell.GetPointIds() ) ) for pointId in pointIds: bucket = collocatedNodes[ pointId ] @@ -112,16 +153,25 @@ def __checkNeighbors( matrix: vtkUnstructuredGrid, fracture: vtkUnstructuredGrid found = 0 tmp = tuple( map( tuple, cns ) ) for nodeCombinations in myIter( tmp ): - f = frozenset( nodeCombinations ) - if f in fractureFaces: + faceCombination = frozenset( nodeCombinations ) + if faceCombination in fractureFaces: found += 1 if found != 2: setupLogger.warning( "Something went wrong since we should have found 2 fractures faces (we found" + f" {found}) for collocated nodes {cns}." ) -def __action( vtkInputFile: str, options: Options ) -> Result: - matrix, fracture = __readMultiblock( vtkInputFile, options.matrixName, options.fractureName ) +def meshAction( matrix: vtkUnstructuredGrid, fracture: vtkUnstructuredGrid, options: Options ) -> Result: + """Performs the check of the fractures and collocated nodes generated. + + Args: + matrix (vtkUnstructuredGrid): The matrix mesh. + fracture (vtkUnstructuredGrid): The fracture mesh. + options (Options): The options for processing. + + Returns: + Result: The result of the self intersecting elements check. + """ matrixPoints: vtkPoints = matrix.GetPoints() fracturePoints: vtkPoints = fracture.GetPoints() @@ -130,7 +180,7 @@ def __action( vtkInputFile: str, options: Options ) -> Result: fracture.GetPointData().GetGlobalIds() and fracture.GetCellData().GetGlobalIds() pointIds = vtk_to_numpy( matrix.GetPointData().GetGlobalIds() ) - g2l = numpy.ones( len( pointIds ), dtype=int ) * -1 + g2l: npt.NDArray[ np.int64 ] = np.ones( len( pointIds ), dtype=np.int64 ) * -1 for loc, glo in enumerate( pointIds ): g2l[ glo ] = loc g2l.flags.writeable = False @@ -146,14 +196,20 @@ def __action( vtkInputFile: str, options: Options ) -> Result: for duplicate in filter( lambda i: i > -1, duplicates ): p0 = matrixPoints.GetPoint( g2l[ duplicate ] ) p1 = fracturePoints.GetPoint( i ) - if numpy.linalg.norm( numpy.array( p1 ) - numpy.array( p0 ) ) > options.tolerance: + if np.linalg.norm( np.array( p1 ) - np.array( p0 ) ) > options.tolerance: errors.append( ( i, g2l[ duplicate ], duplicate ) ) return Result( errors=errors ) -def action( vtkInputFile: str, options: Options ) -> Result: - try: - return __action( vtkInputFile, options ) - except BaseException as e: - setupLogger.error( e ) - return Result( errors=() ) +def action( vtkMultiBlockInputFile: str, options: Options ) -> Result: + """Reads a vtkMultiblock and performs the check of the fractures and collocated nodes generated. + + Args: + vtkMultiBlockInputFile (str): The path to the input vtkMultiBlock file. + options (Options): The options for processing. + + Returns: + Result: The result of the check of the fractures and collocated nodes generated. + """ + matrix, fracture = __readMultiblock( vtkMultiBlockInputFile, options.matrixName, options.fractureName ) + return meshAction( matrix, fracture, options ) diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/collocatedNodes.py b/mesh-doctor/src/geos/mesh_doctor/actions/collocatedNodes.py new file mode 100644 index 00000000..e25df370 --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/actions/collocatedNodes.py @@ -0,0 +1,113 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from collections import defaultdict +from dataclasses import dataclass +import numpy +from typing import Collection, Iterable +from vtkmodules.vtkCommonCore import reference, vtkPoints +from vtkmodules.vtkCommonDataModel import vtkCell, vtkIncrementalOctreePointLocator, vtkUnstructuredGrid +from geos.mesh_doctor.parsing.cliParsing import setupLogger +from geos.mesh.io.vtkIO import readUnstructuredGrid + + +@dataclass( frozen=True ) +class Options: + tolerance: float + + +@dataclass( frozen=True ) +class Result: + nodesBuckets: Iterable[ Iterable[ int ] ] # Each bucket contains the duplicated node indices. + wrongSupportElements: Collection[ int ] # Element indices with support node indices appearing more than once. + + +def findCollocatedNodesBuckets( mesh: vtkUnstructuredGrid, tolerance: float ) -> list[ tuple[ int, ...] ]: + """Check all the nodes of a mesh and returns every bucket of nodes that are collocated within a tolerance. + + Args: + mesh (vtkUnstructuredGrid): The input mesh to analyze. + tolerance (float): The distance tolerance within which nodes are considered collocated. + + Returns: + list[ tuple[ int, ... ] ]: A list of tuples, each containing indices of nodes that are collocated. + """ + points: vtkPoints = mesh.GetPoints() + locator = vtkIncrementalOctreePointLocator() + locator.SetTolerance( tolerance ) + output = vtkPoints() + locator.InitPointInsertion( output, points.GetBounds() ) + + # original ids to/from filtered ids. + filteredToOriginal = numpy.ones( points.GetNumberOfPoints(), dtype=int ) * -1 + + rejectedPoints: defaultdict[ int, list[ int ] ] = defaultdict( list ) + pointId = reference( 0 ) + for i in range( points.GetNumberOfPoints() ): + isInserted = locator.InsertUniquePoint( points.GetPoint( i ), pointId ) # type: ignore[arg-type] + if not isInserted: + # If it's not inserted, `pointId` contains the node that was already at that location. + # But in that case, `pointId` is the new numbering in the destination points array. + # It's more useful for the user to get the old index in the original mesh so he can look for it in his data. + setupLogger.debug( + f"Point {i} at {points.GetPoint(i)} has been rejected, " + + f"point {filteredToOriginal[pointId.get()]} is already inserted.", # type: ignore[misc, call-overload] + ) + rejectedPoints[ pointId.get() ].append( i ) # type: ignore[misc, index] + else: + # If it's inserted, `pointId` contains the new index in the destination array. + # We store this information to be able to connect the source and destination arrays. + # originalToFiltered[i] = pointId.get() + filteredToOriginal[ pointId.get() ] = i # type: ignore[misc, call-overload] + + collocatedNodesBuckets: list[ tuple[ int, ...] ] = [] + for n, ns in rejectedPoints.items(): + collocatedNodesBuckets.append( ( n, *ns ) ) + return collocatedNodesBuckets + + +def findWrongSupportElements( mesh: vtkUnstructuredGrid ) -> list[ int ]: + """Checking that the support node indices appear only once per element. + + Args: + mesh (vtkUnstructuredGrid): The input mesh to analyze. + + Returns: + list[ int ]: A list of cell indices with support node indices appearing more than once. + """ + wrongSupportElements: list[ int ] = [] + for c in range( mesh.GetNumberOfCells() ): + cell: vtkCell = mesh.GetCell( c ) + numPointsPerCell: int = cell.GetNumberOfPoints() + if len( { cell.GetPointId( i ) for i in range( numPointsPerCell ) } ) != numPointsPerCell: + wrongSupportElements.append( c ) + return wrongSupportElements + + +def meshAction( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + """Performs the collocated nodes check on a vtkUnstructuredGrid. + + Args: + mesh (vtkUnstructuredGrid): The input mesh to analyze. + options (Options): The options for processing. + + Returns: + Result: The result of the collocated nodes check. + """ + nodesBuckets = findCollocatedNodesBuckets( mesh, options.tolerance ) + wrongSupportElements = findWrongSupportElements( mesh ) + return Result( nodesBuckets=nodesBuckets, wrongSupportElements=wrongSupportElements ) # type: ignore[arg-type] + + +def action( vtuInputFile: str, options: Options ) -> Result: + """Reads a vtu file and performs the element volumes check on it. + + Args: + vtuInputFile (str): The path to the input VTU file. + options (Options): The options for processing. + + Returns: + Result: The result of the element volumes check. + """ + mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtuInputFile ) + return meshAction( mesh, options ) diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/elementVolumes.py b/mesh-doctor/src/geos/mesh_doctor/actions/elementVolumes.py new file mode 100644 index 00000000..87ea9f6d --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/actions/elementVolumes.py @@ -0,0 +1,121 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from dataclasses import dataclass +import uuid +import numpy as np +from numpy.typing import NDArray +from vtkmodules.vtkCommonCore import vtkDoubleArray +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, VTK_HEXAHEDRON, VTK_PYRAMID, VTK_TETRA, VTK_WEDGE +from vtkmodules.vtkFiltersVerdict import vtkCellSizeFilter, vtkMeshQuality +from vtkmodules.util.numpy_support import vtk_to_numpy +from geos.mesh_doctor.parsing.cliParsing import setupLogger +from geos.mesh.io.vtkIO import readUnstructuredGrid + + +@dataclass( frozen=True ) +class Options: + minVolume: float + + +@dataclass( frozen=True ) +class Result: + elementVolumes: list[ tuple[ int, float ] ] + + +SUPPORTED_TYPES = [ VTK_HEXAHEDRON, VTK_TETRA ] + + +def getMeshQuality( mesh: vtkUnstructuredGrid ) -> vtkDoubleArray: + """Get the quality of the mesh. + + Args: + mesh (vtkUnstructuredGrid): The input mesh. + + Returns: + vtkDoubleArray: The array containing the quality of each cell. + """ + mq = vtkMeshQuality() + mq.SetTetQualityMeasureToVolume() + mq.SetHexQualityMeasureToVolume() + if hasattr( mq, "SetPyramidQualityMeasureToVolume" ): # This feature is quite recent + mq.SetPyramidQualityMeasureToVolume() + SUPPORTED_TYPES.append( VTK_PYRAMID ) + mq.SetWedgeQualityMeasureToVolume() + SUPPORTED_TYPES.append( VTK_WEDGE ) + else: + setupLogger.warning( "Your \"pyvtk\" version does not bring pyramid nor wedge support with vtkMeshQuality." + " Using the fallback solution." ) + mq.SetInputData( mesh ) + mq.Update() + + return mq.GetOutput().GetCellData().GetArray( "Quality" ) # Name is imposed by vtk. + + +def getMeshVolume( mesh: vtkUnstructuredGrid ) -> vtkDoubleArray: + """Get the volume of the mesh. + + Args: + mesh (vtkUnstructuredGrid): The input mesh. + + Returns: + vtkDataArray: The array containing the volume of each cell. + """ + cs = vtkCellSizeFilter() + cs.ComputeAreaOff() + cs.ComputeLengthOff() + cs.ComputeSumOff() + cs.ComputeVertexCountOff() + cs.ComputeVolumeOn() + + volumeArrayName: str = "__MESH_DOCTOR_VOLUME-" + str( uuid.uuid4() ) # Making the name unique + cs.SetVolumeArrayName( volumeArrayName ) + cs.SetInputData( mesh ) + cs.Update() + + return cs.GetOutput().GetCellData().GetArray( volumeArrayName ) + + +def meshAction( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + """Performs the supported elements check on a vtkUnstructuredGrid. + + Args: + mesh (vtkUnstructuredGrid): The input mesh to analyze. + options (Options): The options for processing. + + Returns: + Result: The result of the supported elements check. + """ + volume: vtkDoubleArray = getMeshVolume( mesh ) + if not volume: + setupLogger.error( "Volume computation failed." ) + raise RuntimeError( "Volume computation failed." ) + + quality: vtkDoubleArray = getMeshQuality( mesh ) + if not quality: + setupLogger.error( "Quality computation failed." ) + raise RuntimeError( "Quality computation failed." ) + + volumeArray: NDArray[ np.floating ] = vtk_to_numpy( volume ) + qualityArray: NDArray[ np.floating ] = vtk_to_numpy( quality ) + smallVolumes: list[ tuple[ int, float ] ] = [] + for i, pack in enumerate( zip( volumeArray, qualityArray ) ): + v, q = pack + vol = q if mesh.GetCellType( i ) in SUPPORTED_TYPES else v + if vol < options.minVolume: + smallVolumes.append( ( i, float( vol ) ) ) + return Result( elementVolumes=smallVolumes ) + + +def action( vtuInputFile: str, options: Options ) -> Result: + """Reads a vtu file and performs the element volumes check on it. + + Args: + vtuInputFile (str): The path to the input VTU file. + options (Options): The options for processing. + + Returns: + Result: The result of the element volumes check. + """ + mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtuInputFile ) + return meshAction( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/fixElementsOrderings.py b/mesh-doctor/src/geos/mesh_doctor/actions/fixElementsOrderings.py similarity index 65% rename from geos-mesh/src/geos/mesh/doctor/actions/fixElementsOrderings.py rename to mesh-doctor/src/geos/mesh_doctor/actions/fixElementsOrderings.py index e560551e..e29ea0fb 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/fixElementsOrderings.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/fixElementsOrderings.py @@ -1,8 +1,11 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto from dataclasses import dataclass from vtkmodules.vtkCommonCore import vtkIdList from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid -from geos.mesh.utils.genericHelpers import toVtkIdList from geos.mesh.io.vtkIO import VtkOutput, readUnstructuredGrid, writeMesh +from geos.mesh.utils.genericHelpers import toVtkIdList @dataclass( frozen=True ) @@ -17,7 +20,16 @@ class Result: unchangedCellTypes: frozenset[ int ] -def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: +def meshAction( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + """Performs the fix elements orderings on a vtkUnstructuredGrid. + + Args: + mesh (vtkUnstructuredGrid): The input mesh to reorder. + options (Options): The options for processing. + + Returns: + Result: The result of the fix elements orderings. + """ # The vtk cell type is an int and will be the key of the following mapping, # that will point to the relevant permutation. cellTypeToOrdering: dict[ int, list[ int ] ] = options.cellTypeToOrdering @@ -38,8 +50,8 @@ def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: supportPointIds = vtkIdList() cells.GetCellAtId( cellIdx, supportPointIds ) newSupportPointIds = [] - for i, v in enumerate( newOrdering ): - newSupportPointIds.append( supportPointIds.GetId( newOrdering[ i ] ) ) + for orderingId in newOrdering: + newSupportPointIds.append( supportPointIds.GetId( orderingId ) ) cells.ReplaceCellAtId( cellIdx, toVtkIdList( newSupportPointIds ) ) else: unchangedCellTypes.add( cellType ) @@ -48,6 +60,15 @@ def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: unchangedCellTypes=frozenset( unchangedCellTypes ) ) -def action( vtkInputFile: str, options: Options ) -> Result: - mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtkInputFile ) - return __action( mesh, options ) +def action( vtuInputFile: str, options: Options ) -> Result: + """Reads a vtu file and performs the fix elements orderings on it. + + Args: + vtuInputFile (str): The path to the input VTU file. + options (Options): The options for processing. + + Returns: + Result: The result of the fix elements orderings. + """ + mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtuInputFile ) + return meshAction( mesh, options ) diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/generateCube.py b/mesh-doctor/src/geos/mesh_doctor/actions/generateCube.py new file mode 100644 index 00000000..9af53ddb --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/actions/generateCube.py @@ -0,0 +1,239 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from dataclasses import dataclass +import numpy as np +import numpy.typing as npt +from typing import Iterable, Sequence +from vtkmodules.util.numpy_support import numpy_to_vtk +from vtkmodules.vtkCommonCore import vtkPoints +from vtkmodules.vtkCommonDataModel import ( vtkCellArray, vtkHexahedron, vtkRectilinearGrid, vtkUnstructuredGrid, + VTK_HEXAHEDRON ) +from geos.mesh.io.vtkIO import VtkOutput, writeMesh +from geos.mesh.utils.arrayModifiers import createConstantAttributeDataSet +from geos.mesh_doctor.actions.generateGlobalIds import buildGlobalIds +from geos.mesh_doctor.parsing.cliParsing import setupLogger + + +@dataclass( frozen=True ) +class Result: + info: str + + +@dataclass( frozen=True ) +class FieldInfo: + name: str + dimension: int + support: str + + +@dataclass( frozen=True ) +class Options: + vtkOutput: VtkOutput + generateCellsGlobalIds: bool + generatePointsGlobalIds: bool + xs: Sequence[ float ] + ys: Sequence[ float ] + zs: Sequence[ float ] + nxs: Sequence[ int ] + nys: Sequence[ int ] + nzs: Sequence[ int ] + fields: Iterable[ FieldInfo ] + + +@dataclass( frozen=True ) +class XYZ: + x: npt.NDArray + y: npt.NDArray + z: npt.NDArray + + +def buildCoordinates( positions: Sequence[ float ], numElements: Sequence[ int ] ) -> npt.NDArray: + """Builds the coordinates array from the positions of each mesh block vertex and number of elements in a given direction within each mesh block. + + Args: + positions (Sequence[ float ]): The positions of each mesh block vertex. + numElements (Sequence[ int ]): The number of elements in a given direction within each mesh block. + + Returns: + npt.NDArray: The array of coordinates. + """ + result: list[ npt.NDArray ] = [] + it = zip( zip( positions, positions[ 1: ] ), numElements ) + try: + coords, n = next( it ) + while True: + start, stop = coords + endPoint: bool = False + tmp = np.linspace( start=start, stop=stop, num=n + endPoint, endpoint=endPoint ) + coords, n = next( it ) + result.append( tmp ) + except StopIteration: + endPoint = True + tmp = np.linspace( start=start, stop=stop, num=n + endPoint, endpoint=endPoint ) + result.append( tmp ) + return np.concatenate( result ) + + +def buildRectilinearGrid( x: npt.NDArray, y: npt.NDArray, z: npt.NDArray ) -> vtkUnstructuredGrid: + """Builds an unstructured vtk grid from the x,y,z coordinates. + + Args: + x (npt.NDArray): The x coordinates. + y (npt.NDArray): The y coordinates. + z (npt.NDArray): The z coordinates. + + Returns: + The unstructured mesh, even if it's topologically structured. + """ + rg = vtkRectilinearGrid() + rg.SetDimensions( len( x ), len( y ), len( z ) ) + rg.SetXCoordinates( numpy_to_vtk( x ) ) + rg.SetYCoordinates( numpy_to_vtk( y ) ) + rg.SetZCoordinates( numpy_to_vtk( z ) ) + + numPoints = rg.GetNumberOfPoints() + numCells = rg.GetNumberOfCells() + + points = vtkPoints() + points.Allocate( numPoints ) + for i in range( rg.GetNumberOfPoints() ): + points.InsertNextPoint( rg.GetPoint( i ) ) + + cellTypes = [ VTK_HEXAHEDRON ] * numCells + cells = vtkCellArray() + cells.AllocateExact( numCells, numCells * 8 ) + + m = ( 0, 1, 3, 2, 4, 5, 7, 6 ) # VTK_VOXEL and VTK_HEXAHEDRON do not share the same ordering. + for i in range( rg.GetNumberOfCells() ): + c = rg.GetCell( i ) + newCell = vtkHexahedron() + for j in range( 8 ): + newCell.GetPointIds().SetId( j, c.GetPointId( m[ j ] ) ) + cells.InsertNextCell( newCell ) + + mesh = vtkUnstructuredGrid() + mesh.SetPoints( points ) + mesh.SetCells( cellTypes, cells ) + return mesh + + +def buildRectilinearBlocksMesh( xyzs: Iterable[ XYZ ] ) -> vtkUnstructuredGrid: + """Builds an unstructured vtk grid from the `xyzs` blocks. Kind of InternalMeshGenerator. + + Args: + xyzs (Iterable[ XYZ ]): The blocks. + + Returns: + The unstructured mesh, even if it's topologically structured. + """ + rgs: list[ vtkRectilinearGrid ] = [] + for xyz in xyzs: + rg = vtkRectilinearGrid() + rg.SetDimensions( len( xyz.x ), len( xyz.y ), len( xyz.z ) ) + rg.SetXCoordinates( numpy_to_vtk( xyz.x ) ) + rg.SetYCoordinates( numpy_to_vtk( xyz.y ) ) + rg.SetZCoordinates( numpy_to_vtk( xyz.z ) ) + rgs.append( rg ) + + numPoints: int = sum( r.GetNumberOfPoints() for r in rgs ) + numCells: int = sum( r.GetNumberOfCells() for r in rgs ) + + points = vtkPoints() + points.Allocate( numPoints ) + for rg in rgs: + for i in range( rg.GetNumberOfPoints() ): + points.InsertNextPoint( rg.GetPoint( i ) ) + + cellTypes: list[ int ] = [ VTK_HEXAHEDRON ] * numCells + cells = vtkCellArray() + cells.AllocateExact( numCells, numCells * 8 ) + + m = ( 0, 1, 3, 2, 4, 5, 7, 6 ) # VTK_VOXEL and VTK_HEXAHEDRON do not share the same ordering. + offset: int = 0 + for rg in rgs: + for i in range( rg.GetNumberOfCells() ): + c = rg.GetCell( i ) + newCell = vtkHexahedron() + for j in range( 8 ): + newCell.GetPointIds().SetId( j, offset + c.GetPointId( m[ j ] ) ) + cells.InsertNextCell( newCell ) + offset += rg.GetNumberOfPoints() + + mesh = vtkUnstructuredGrid() + mesh.SetPoints( points ) + mesh.SetCells( cellTypes, cells ) + + return mesh + + +def addFields( mesh: vtkUnstructuredGrid, fields: Iterable[ FieldInfo ] ) -> vtkUnstructuredGrid: + """Add constant fields to the mesh using arrayModifiers utilities. + + Each field is filled with ones (1.0) for all components. + + Args: + mesh (vtkUnstructuredGrid): The mesh to which fields will be added. + fields (Iterable[ FieldInfo ]): The fields to add. + + Returns: + vtkUnstructuredGrid: The mesh with added fields. + """ + for fieldInfo in fields: + onPoints = fieldInfo.support == "POINTS" + # Create list of values (all 1.0) for each component + listValues = [ 1.0 ] * fieldInfo.dimension + # Use the robust createConstantAttributeDataSet function + success = createConstantAttributeDataSet( dataSet=mesh, + listValues=listValues, + attributeName=fieldInfo.name, + onPoints=onPoints, + logger=setupLogger ) + if not success: + setupLogger.warning( f"Failed to create field {fieldInfo.name}" ) + return mesh + + +def buildCube( options: Options ) -> vtkUnstructuredGrid: + """Creates a cube vtkUnstructuredGrid from options given. + + Args: + options (Options): The options for processing. + + Returns: + vtkUnstructuredGrid: The created mesh. + """ + x = buildCoordinates( options.xs, options.nxs ) + y = buildCoordinates( options.ys, options.nys ) + z = buildCoordinates( options.zs, options.nzs ) + mesh = buildRectilinearBlocksMesh( ( XYZ( x, y, z ), ) ) + mesh = addFields( mesh, options.fields ) + buildGlobalIds( mesh, options.generateCellsGlobalIds, options.generatePointsGlobalIds ) + return mesh + + +def meshAction( options: Options ) -> Result: + """Creates a vtkUnstructuredGrid from options given. + + Args: + options (Options): The options for processing. + + Returns: + Result: The result of creation of the mesh. + """ + outputMesh = buildCube( options ) + writeMesh( outputMesh, options.vtkOutput ) + return Result( info=f"Mesh was written to {options.vtkOutput.output}" ) + + +def action( vtuInputFile: str, options: Options ) -> Result: + """Creates a vtkUnstructuredGrid from options given. + + Args: + vtuInputFile (str): The path to the input VTU file. This is unused. + options (Options): The options for processing. + + Returns: + Result: The result of the fix elements orderings. + """ + return meshAction( options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/generateFractures.py b/mesh-doctor/src/geos/mesh_doctor/actions/generateFractures.py similarity index 79% rename from geos-mesh/src/geos/mesh/doctor/actions/generateFractures.py rename to mesh-doctor/src/geos/mesh_doctor/actions/generateFractures.py index 089354d9..8e38f8fe 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/generateFractures.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/generateFractures.py @@ -1,29 +1,28 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto from collections import defaultdict from dataclasses import dataclass from enum import Enum import networkx from numpy import empty, ones, zeros from tqdm import tqdm -from typing import Collection, Iterable, Mapping, Optional, Sequence +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 numpy_to_vtk, vtk_to_numpy from vtkmodules.util.vtkConstants import VTK_ID_TYPE -from geos.mesh.doctor.actions.vtkPolyhedron import FaceStream -from geos.mesh.doctor.parsing.cliParsing import setupLogger +from geos.mesh_doctor.actions.vtkPolyhedron import FaceStream +from geos.mesh_doctor.parsing.cliParsing import setupLogger +from geos.mesh.io.vtkIO import VtkOutput, readUnstructuredGrid, writeMesh from geos.mesh.utils.arrayHelpers import hasArray from geos.mesh.utils.genericHelpers import toVtkIdList, vtkIter -from geos.mesh.io.vtkIO import VtkOutput, readUnstructuredGrid, writeMesh -""" -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 ] +IDMapping: TypeAlias = Mapping[ int, int ] +CellsPointsCoords: TypeAlias = dict[ int, list[ tuple[ float ] ] ] +Coordinates3D: TypeAlias = tuple[ float, float, float ] class FracturePolicy( Enum ): @@ -54,7 +53,16 @@ class FractureInfo: def buildNodeToCells( mesh: vtkUnstructuredGrid, - faceNodes: Iterable[ Iterable[ int ] ] ) -> dict[ int, Iterable[ int ] ]: + faceNodes: Iterable[ Iterable[ int ] ] ) -> Mapping[ int, Iterable[ int ] ]: + """Builds the mapping from each fracture node to the cells that use this node. + + Args: + mesh (vtkUnstructuredGrid): The input mesh. + faceNodes (Iterable[ Iterable[ int ] ]): The nodes of each face of the fracture. + + Returns: + Mapping[ int, Iterable[ int ] ]: A mapping from each fracture node to the cells that use this node. + """ # TODO normally, just a list and not a set should be enough. nodeToCells: dict[ int, set[ int ] ] = defaultdict( set ) @@ -74,6 +82,16 @@ def buildNodeToCells( mesh: vtkUnstructuredGrid, def __buildFractureInfoFromFields( mesh: vtkUnstructuredGrid, f: Sequence[ int ], fieldValues: frozenset[ int ] ) -> FractureInfo: + """Build the fracture info from field values. + + Args: + mesh (vtkUnstructuredGrid): The input mesh. + f (Sequence[ int ]): The field values for each cell. + fieldValues (frozenset[ int ]): The set of field values to consider. + + Returns: + FractureInfo: The fracture information. + """ cellsToFaces: 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, @@ -93,7 +111,7 @@ def __buildFractureInfoFromFields( mesh: vtkUnstructuredGrid, f: Sequence[ int ] if f[ neighborCellId ] != f[ cellId ] and f[ neighborCellId ] in fieldValues: # TODO add this (cellIds, faceId) information to the fractureInfo? cellsToFaces[ cellId ].append( i ) - faceNodes: list[ Collection[ int ] ] = list() + faceNodes: list[ Collection[ int ] ] = [] faceNodesHashes: set[ frozenset[ int ] ] = set() # A temporary not to add multiple times the same face. for cellId, facesIds in tqdm( cellsToFaces.items(), desc="Extracting the faces of the fractures" ): cell = mesh.GetCell( cellId ) @@ -103,28 +121,37 @@ def __buildFractureInfoFromFields( mesh: vtkUnstructuredGrid, f: Sequence[ int ] if fnh not in faceNodesHashes: faceNodesHashes.add( fnh ) faceNodes.append( fn ) - nodeToCells: dict[ int, Iterable[ int ] ] = buildNodeToCells( mesh, faceNodes ) - faceCellId: list = list() # no cell of the mesh corresponds to that face when fracture policy is 'field' + nodeToCells: Mapping[ int, Iterable[ int ] ] = buildNodeToCells( mesh, faceNodes ) + faceCellId: list = [] # no cell of the mesh corresponds to that face when fracture policy is 'field' return FractureInfo( nodeToCells=nodeToCells, faceNodes=faceNodes, faceCellId=faceCellId ) def __buildFractureInfoFromInternalSurfaces( mesh: vtkUnstructuredGrid, f: Sequence[ int ], fieldValues: frozenset[ int ] ) -> FractureInfo: + """Build the fracture info from internal surfaces defined by the field values. + + Args: + mesh (vtkUnstructuredGrid): The input mesh. + f (Sequence[ int ]): The field values for each cell. + fieldValues (frozenset[ int ]): The set of field values to consider. + + Returns: + FractureInfo: The fracture information. + """ nodeToCells: dict[ int, list[ int ] ] = defaultdict( list ) - faceNodes: list[ Collection[ int ] ] = list() - faceCellId: list[ int ] = list() + faceNodes: list[ Collection[ int ] ] = [] + faceCellId: list[ int ] = [] for cellId in tqdm( range( mesh.GetNumberOfCells() ), desc="Computing the face to nodes mapping" ): cell = mesh.GetCell( cellId ) - if cell.GetCellDimension() == 2: - if f[ cellId ] in fieldValues: - nodes = list() - for v in range( cell.GetNumberOfPoints() ): - pointId: int = cell.GetPointId( v ) - nodeToCells[ pointId ] = list() - nodes.append( pointId ) - faceNodes.append( tuple( nodes ) ) - faceCellId.append( cellId ) + if cell.GetCellDimension() == 2 and f[ cellId ] in fieldValues: + nodes = [] + for v in range( cell.GetNumberOfPoints() ): + pointId: int = cell.GetPointId( v ) + nodeToCells[ pointId ] = [] + nodes.append( pointId ) + faceNodes.append( tuple( nodes ) ) + faceCellId.append( cellId ) for cellId in tqdm( range( mesh.GetNumberOfCells() ), desc="Computing the node to cells mapping" ): cell = mesh.GetCell( cellId ) @@ -140,11 +167,22 @@ def buildFractureInfo( mesh: vtkUnstructuredGrid, options: Options, combinedFractures: bool, fractureId: int = 0 ) -> FractureInfo: + """For a given mesh and options, builds the fracture information. + + Args: + mesh (vtkUnstructuredGrid): The input mesh. + options (Options): Options for processing. + combinedFractures (bool): If True, considers all fractures together. Else, only the fracture with id. + fractureId (int, optional): The fracture id to consider if combinedFractures is False. Defaults to 0. + + Raises: + ValueError: If the field does not exist in the mesh. + + Returns: + FractureInfo: The fracture information. + """ field = options.field - if combinedFractures: - fieldValues = options.fieldValuesCombined - else: - fieldValues = options.fieldValuesPerFracture[ fractureId ] + fieldValues = options.fieldValuesCombined if combinedFractures else options.fieldValuesPerFracture[ fractureId ] cellData = mesh.GetCellData() if cellData.HasArray( field ): f = vtk_to_numpy( cellData.GetArray( field ) ) @@ -159,6 +197,7 @@ def buildFractureInfo( mesh: vtkUnstructuredGrid, def buildCellToCellGraph( mesh: vtkUnstructuredGrid, fracture: FractureInfo ) -> networkx.Graph: """Connects all the cells that touch the fracture by at least one node. + Two cells are connected when they share at least a face which is not a face of the fracture. Args: @@ -171,7 +210,7 @@ def buildCellToCellGraph( 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 ] ] = list() + tmp: list[ frozenset[ int ] ] = [] for fn in fracture.faceNodes: tmp.append( frozenset( fn ) ) faceHashes: frozenset[ frozenset[ int ] ] = frozenset( tmp ) @@ -202,28 +241,28 @@ def buildCellToCellGraph( mesh: vtkUnstructuredGrid, fracture: FractureInfo ) -> def _identifySplit( numPoints: int, cellToCell: networkx.Graph, - nodeToCells: dict[ int, Iterable[ int ] ] ) -> dict[ int, IDMapping ]: + nodeToCells: Mapping[ int, Iterable[ int ] ] ) -> Mapping[ int, IDMapping ]: """For each cell, compute the node indices replacements. Args: numPoints (int): The number of points in the whole mesh (not the fracture). cellToCell (networkx.Graph): The cell to cell graph (connection through common faces). - nodeToCells (dict[ int, Iterable[ int ] ]): Maps the nodes of the fracture to the cells relying on this node. + nodeToCells (Mapping[ int, Iterable[ int ] ]): Maps the nodes of the fracture to the cells relying on this node. Returns: - dict[ int, IDMapping ]: For each cell (first key), returns a mapping from the current index + Mapping[ int, IDMapping ]: For each cell (first key), returns a mapping from the current index and the new index that should replace the current index. Note that the current index and the new index can be identical: no replacement should be done then. """ class NewIndex: - """ - Returns the next available index. + """Returns the next available index. + Note that the first time an index is met, the index itself is returned: we do not want to change an index if we do not have to. """ - def __init__( self, numNodes: int ): + def __init__( self, numNodes: int ) -> None: self.__currentLastIndex = numNodes - 1 self.__seen: set[ int ] = set() @@ -236,7 +275,7 @@ def __call__( self, index: int ) -> int: return index buildNewIndex = NewIndex( numPoints ) - result: dict[ int, IDMapping ] = defaultdict( dict ) + result: dict[ int, dict[ int, int ] ] = defaultdict( dict ) # Iteration over `sorted` nodes to have a predictable result for tests. for node, cells in tqdm( sorted( nodeToCells.items() ), desc="Identifying the node splits" ): for connectedCells in networkx.connected_components( cellToCell.subgraph( cells ) ): @@ -249,14 +288,13 @@ def __call__( self, index: int ) -> int: def __copyFieldsSplitMesh( oldMesh: vtkUnstructuredGrid, splitMesh: vtkUnstructuredGrid, - addedPointsWithOldId: list[ tuple[ int ] ] ) -> None: - """Copies the fields from the old mesh to the new one. - Point data will be duplicated for collocated nodes. + addedPointsWithOldId: list[ tuple[ int, int ] ] ) -> None: + """Copies the fields from the old mesh to the new one. Point data will be duplicated for collocated nodes. Args: - oldMesh (vtkUnstructuredGrid): The mesh before the split._ + oldMesh (vtkUnstructuredGrid): The mesh before the split. splitMesh (vtkUnstructuredGrid): The mesh after the split. Will receive the fields in place. - addedPointsWithOldId (list[ tuple[ int ] ]): _description_ + addedPointsWithOldId (list[ tuple[ int, int ] ]): List of tuples (newId, oldId) for each duplicated point. """ # Copying the cell data. The cells are the same, just their nodes support have changed. inputCellData = oldMesh.GetCellData() @@ -358,7 +396,7 @@ def __copyFieldsFractureMesh( oldMesh: vtkUnstructuredGrid, fractureMesh: vtkUns def __performSplit( oldMesh: vtkUnstructuredGrid, cellToNodeMapping: Mapping[ int, IDMapping ] ) -> vtkUnstructuredGrid: - """Split the main 3d mesh based on the node duplication information contained in @p cellToNodeMapping + """Split the main 3d mesh based on the node duplication information contained in @p cellToNodeMapping. Args: oldMesh (vtkUnstructuredGrid): The main 3d mesh. @@ -369,7 +407,7 @@ def __performSplit( oldMesh: vtkUnstructuredGrid, cellToNodeMapping: Mapping[ in vtkUnstructuredGrid: The main 3d mesh split at the fracture location. """ addedPoints: set[ int ] = set() - addedPointsWithOldId: list[ tuple[ int ] ] = list() + addedPointsWithOldId: list[ tuple[ int, int ] ] = [] for nodeMapping in cellToNodeMapping.values(): for i, o in nodeMapping.items(): if i != o: @@ -387,8 +425,8 @@ def __performSplit( oldMesh: vtkUnstructuredGrid, cellToNodeMapping: Mapping[ in newPoints.SetPoint( p, oldPoints.GetPoint( p ) ) collocatedNodes[ p ] = p # Creating the new collocated/duplicated points based on the old points positions. - for nodeMapping in cellToNodeMapping.values(): - for i, o in nodeMapping.items(): + for nodeMappingDup in cellToNodeMapping.values(): + for i, o in nodeMappingDup.items(): if i != o: newPoints.SetPoint( o, oldPoints.GetPoint( i ) ) collocatedNodes[ o ] = i @@ -407,18 +445,18 @@ def __performSplit( oldMesh: vtkUnstructuredGrid, cellToNodeMapping: Mapping[ in newMesh.Allocate( oldMesh.GetNumberOfCells() ) for c in tqdm( range( oldMesh.GetNumberOfCells() ), desc="Performing the mesh split" ): - nodeMapping: IDMapping = cellToNodeMapping.get( c, {} ) + cellNodeMapping: IDMapping = cellToNodeMapping.get( c, {} ) cell: vtkCell = oldMesh.GetCell( c ) cellType: int = cell.GetCellType() # For polyhedron, we'll manipulate the face stream directly. if cellType == VTK_POLYHEDRON: faceStream = vtkIdList() oldMesh.GetFaceStream( c, faceStream ) - newFaceNodes: list[ list[ int ] ] = list() + newFaceNodes: list[ list[ int ] ] = [] for faceNodes in FaceStream.buildFromVtkIdList( faceStream ).faceNodes: - newPointIds = list() + newPointIds = [] for currentPointId in faceNodes: - newPointId: int = nodeMapping.get( currentPointId, currentPointId ) + newPointId: int = cellNodeMapping.get( currentPointId, currentPointId ) newPointIds.append( newPointId ) newFaceNodes.append( newPointIds ) newMesh.InsertNextCell( cellType, toVtkIdList( FaceStream( newFaceNodes ).dump() ) ) @@ -427,9 +465,9 @@ def __performSplit( oldMesh: vtkUnstructuredGrid, cellToNodeMapping: Mapping[ in # Then the values will be (potentially) overwritten in place, before being sent back into the cell. cellPointIds: vtkIdList = cell.GetPointIds() for i in range( cellPointIds.GetNumberOfIds() ): - currentPointId: int = cellPointIds.GetId( i ) - newPointId: int = nodeMapping.get( currentPointId, currentPointId ) - cellPointIds.SetId( i, newPointId ) + currentPtId: int = cellPointIds.GetId( i ) + newPtId: int = cellNodeMapping.get( currentPtId, currentPtId ) + cellPointIds.SetId( i, newPtId ) newMesh.InsertNextCell( cellType, cellPointIds ) __copyFieldsSplitMesh( oldMesh, newMesh, addedPointsWithOldId ) @@ -461,10 +499,10 @@ def __generateFractureMesh( oldMesh: vtkUnstructuredGrid, fractureInfo: Fracture # 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. - faceNodes: list[ Collection[ int ] ] = list() + faceNodes: list[ Collection[ int ] ] = [] discardedFaceNodes: set[ Iterable[ int ] ] = set() - if fractureInfo.faceCellId != list(): # The fracture policy is 'internalSurfaces' - faceCellId: list[ int ] = list() + if fractureInfo.faceCellId != []: # The fracture policy is 'internalSurfaces' + faceCellId: list[ int ] = [] for ns, fId in zip( fractureInfo.faceNodes, fractureInfo.faceCellId ): if any( map( isNodeDuplicated.__getitem__, ns ) ): faceNodes.append( ns ) @@ -479,7 +517,7 @@ def __generateFractureMesh( oldMesh: vtkUnstructuredGrid, fractureInfo: Fracture discardedFaceNodes.add( ns ) if discardedFaceNodes: - msg: str = "(" + '), ('.join( map( lambda dfns: ", ".join( map( str, dfns ) ), discardedFaceNodes ) ) + ")" + msg: str = "(" + '), ('.join( ", ".join( map( str, dfns ) ) for dfns in discardedFaceNodes ) + ")" setupLogger.info( f"The faces made of nodes [{msg}] were/was discarded" " from the fracture mesh because none of their/its nodes were duplicated." ) @@ -491,7 +529,7 @@ def __generateFractureMesh( oldMesh: vtkUnstructuredGrid, fractureInfo: Fracture numPoints: int = len( fractureNodes ) points = vtkPoints() points.SetNumberOfPoints( numPoints ) - node3dToNode2d: IDMapping = dict() # Building the node mapping, from 3d mesh nodes to 2d fracture nodes. + node3dToNode2d: dict[ int, int ] = {} # Building the node mapping, from 3d mesh nodes to 2d fracture nodes. for i, n in enumerate( fractureNodes ): coords: Coordinates3D = meshPoints.GetPoint( n ) points.SetPoint( i, coords ) @@ -508,8 +546,8 @@ def __generateFractureMesh( oldMesh: vtkUnstructuredGrid, fractureInfo: Fracture polygons.InsertNextCell( polygon ) buckets: dict[ int, set[ int ] ] = defaultdict( set ) - for nodeMapping in cellToNodeMapping.values(): - for i, o in nodeMapping.items(): + for nodeMappingBucket in cellToNodeMapping.values(): + for i, o in nodeMappingBucket.items(): k: Optional[ int ] = node3dToNode2d.get( min( i, o ) ) if k is not None: buckets[ k ].update( ( i, o ) ) @@ -532,7 +570,7 @@ def __generateFractureMesh( oldMesh: vtkUnstructuredGrid, fractureInfo: Fracture # The copy of fields from the old mesh to the fracture is only available when using the internalSurfaces policy # because the FractureInfo is linked to 2D elements from the oldMesh - if fractureInfo.faceCellId != list(): + if fractureInfo.faceCellId != []: __copyFieldsFractureMesh( oldMesh, fractureMesh, faceCellId, node3dToNode2d ) return fractureMesh @@ -540,7 +578,16 @@ def __generateFractureMesh( oldMesh: vtkUnstructuredGrid, fractureInfo: Fracture def __splitMeshOnFractures( mesh: vtkUnstructuredGrid, options: Options ) -> tuple[ vtkUnstructuredGrid, list[ vtkUnstructuredGrid ] ]: - allFractureInfos: list[ FractureInfo ] = list() + """Splits the input mesh based on the fracture information contained in options. + + Args: + mesh (vtkUnstructuredGrid): The input mesh to be split. + options (Options): The options for processing. + + Returns: + tuple[ vtkUnstructuredGrid, list[ vtkUnstructuredGrid ] ]: The output mesh and a list of fracture meshes. + """ + allFractureInfos: list[ FractureInfo ] = [] for fractureId in range( len( options.fieldValuesPerFracture ) ): fractureInfo: FractureInfo = buildFractureInfo( mesh, options, False, fractureId ) allFractureInfos.append( fractureInfo ) @@ -549,14 +596,23 @@ def __splitMeshOnFractures( mesh: vtkUnstructuredGrid, cellToNodeMapping: Mapping[ int, IDMapping ] = _identifySplit( mesh.GetNumberOfPoints(), cellToCell, combinedFractures.nodeToCells ) outputMesh: vtkUnstructuredGrid = __performSplit( mesh, cellToNodeMapping ) - fractureMeshes: list[ vtkUnstructuredGrid ] = list() + fractureMeshes: list[ vtkUnstructuredGrid ] = [] for fractureInfoSeparated in allFractureInfos: fractureMesh: vtkUnstructuredGrid = __generateFractureMesh( mesh, fractureInfoSeparated, cellToNodeMapping ) fractureMeshes.append( fractureMesh ) return ( outputMesh, fractureMeshes ) -def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: +def meshAction( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + """Performs the generation of fractures on a vtkUnstructuredGrid if options given to do so. + + Args: + mesh (vtkUnstructuredGrid): The input mesh to generate fractures on. + options (Options): The options for processing. + + Returns: + Result: The result of the fracture generation. + """ outputMesh, fractureMeshes = __splitMeshOnFractures( mesh, options ) writeMesh( outputMesh, options.meshVtkOutput ) for i, fractureMesh in enumerate( fractureMeshes ): @@ -565,16 +621,21 @@ def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: return Result( info="OK" ) -def action( vtkInputFile: str, options: Options ) -> Result: - try: - mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtkInputFile ) - # Mesh cannot contain global ids before splitting. - if hasArray( mesh, [ "GLOBAL_IDS_POINTS", "GLOBAL_IDS_CELLS" ] ): - errMsg: 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." ) - setupLogger.error( errMsg ) - raise ValueError( errMsg ) - return __action( mesh, options ) - except BaseException as e: - setupLogger.error( e ) - return Result( info="Something went wrong" ) +def action( vtuInputFile: str, options: Options ) -> Result: + """Reads a vtkUnstructuredGrid and generates fractures based on the provided options. + + Args: + vtuInputFile (str): The path to the input VTU file. + options (Options): The options for processing. + + Returns: + Result: The result of the global ids addition. + """ + mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtuInputFile ) + # Mesh cannot contain global ids before splitting. + if hasArray( mesh, [ "GLOBAL_IDS_POINTS", "GLOBAL_IDS_CELLS" ] ): + errMsg: 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." ) + setupLogger.error( errMsg ) + raise ValueError( errMsg ) + return meshAction( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/generateGlobalIds.py b/mesh-doctor/src/geos/mesh_doctor/actions/generateGlobalIds.py similarity index 61% rename from geos-mesh/src/geos/mesh/doctor/actions/generateGlobalIds.py rename to mesh-doctor/src/geos/mesh_doctor/actions/generateGlobalIds.py index 34948350..70ce5a56 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/generateGlobalIds.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/generateGlobalIds.py @@ -1,8 +1,11 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto from dataclasses import dataclass from vtkmodules.vtkCommonCore import vtkIdTypeArray from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid -from geos.mesh.doctor.parsing.cliParsing import setupLogger from geos.mesh.io.vtkIO import VtkOutput, readUnstructuredGrid, writeMesh +from geos.mesh_doctor.parsing.cliParsing import setupLogger @dataclass( frozen=True ) @@ -17,7 +20,7 @@ class Result: info: str -def __buildGlobalIds( mesh: vtkUnstructuredGrid, generateCellsGlobalIds: bool, generatePointsGlobalIds: bool ) -> None: +def buildGlobalIds( mesh: vtkUnstructuredGrid, generateCellsGlobalIds: bool, generatePointsGlobalIds: bool ) -> None: """Adds the global ids for cells and points in place into the mesh instance. Args: @@ -48,16 +51,30 @@ def __buildGlobalIds( mesh: vtkUnstructuredGrid, generateCellsGlobalIds: bool, g mesh.GetCellData().SetGlobalIds( cellsGlobalIds ) -def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: - __buildGlobalIds( mesh, options.generateCellsGlobalIds, options.generatePointsGlobalIds ) +def meshAction( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + """Performs the addition of global ids on a vtkUnstructuredGrid if options given to do so. + + Args: + mesh (vtkUnstructuredGrid): The input mesh to add global ids. + options (Options): The options for processing. + + Returns: + Result: The result of the global ids addition. + """ + buildGlobalIds( mesh, options.generateCellsGlobalIds, options.generatePointsGlobalIds ) writeMesh( mesh, options.vtkOutput ) return Result( info=f"Mesh was written to {options.vtkOutput.output}" ) -def action( vtkInputFile: str, options: Options ) -> Result: - try: - mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtkInputFile ) - return __action( mesh, options ) - except BaseException as e: - setupLogger.error( e ) - return Result( info="Something went wrong." ) +def action( vtuInputFile: str, options: Options ) -> Result: + """Reads a vtkUnstructuredGrid and adds global ids as per options. + + Args: + vtuInputFile (str): The path to the input VTU file. + options (Options): The options for processing. + + Returns: + Result: The result of the global ids addition. + """ + mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtuInputFile ) + return meshAction( mesh, options ) diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/mainChecks.py b/mesh-doctor/src/geos/mesh_doctor/actions/mainChecks.py new file mode 100644 index 00000000..7d20c225 --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/actions/mainChecks.py @@ -0,0 +1,4 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from geos.mesh_doctor.actions.allChecks import action # noqa: F401 diff --git a/geos-mesh/src/geos/mesh/doctor/actions/nonConformal.py b/mesh-doctor/src/geos/mesh_doctor/actions/nonConformal.py similarity index 66% rename from geos-mesh/src/geos/mesh/doctor/actions/nonConformal.py rename to mesh-doctor/src/geos/mesh_doctor/actions/nonConformal.py index b7683540..7d2ce18a 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/nonConformal.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/nonConformal.py @@ -1,9 +1,14 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto from dataclasses import dataclass import math -import numpy +import numpy as np +import numpy.typing as npt +from typing import Union from tqdm import tqdm from vtk import reference as vtkReference -from vtkmodules.vtkCommonCore import vtkDataArray, vtkIdList, vtkPoints +from vtkmodules.vtkCommonCore import vtkIntArray, vtkFloatArray, vtkIdList, vtkPoints from vtkmodules.vtkCommonDataModel import ( vtkBoundingBox, vtkCell, vtkCellArray, vtkPointSet, vtkPolyData, vtkStaticCellLocator, vtkStaticPointLocator, vtkUnstructuredGrid, VTK_POLYHEDRON ) @@ -11,9 +16,9 @@ from vtkmodules.vtkFiltersCore import vtkPolyDataNormals from vtkmodules.vtkFiltersGeometry import vtkDataSetSurfaceFilter from vtkmodules.vtkFiltersModeling import vtkCollisionDetectionFilter, vtkLinearExtrusionFilter -from geos.mesh.doctor.actions import reorientMesh, triangleDistance -from geos.mesh.utils.genericHelpers import vtkIter +from geos.mesh_doctor.actions import reorientMesh, triangleDistance from geos.mesh.io.vtkIO import readUnstructuredGrid +from geos.mesh.utils.genericHelpers import vtkIter @dataclass( frozen=True ) @@ -29,8 +34,8 @@ class Result: class BoundaryMesh: - """ - A BoundaryMesh is the envelope of the 3d mesh on which we want to perform the simulations. + """A BoundaryMesh is the envelope of the 3d mesh on which we want to perform the simulations. + It is computed by vtk. But we want to be sure that the normals of the envelope are directed outwards. The `vtkDataSetSurfaceFilter` does not have the same behavior for standard vtk cells (like tets or hexs), and for polyhedron meshes, for which the result is a bit brittle. @@ -38,7 +43,7 @@ class BoundaryMesh: And then we compute the boundary meshes for both meshes, given that the computing options are not identical. """ - def __init__( self, mesh: vtkUnstructuredGrid ): + def __init__( self, mesh: vtkUnstructuredGrid ) -> None: """Builds a boundary mesh. Args: @@ -46,20 +51,22 @@ def __init__( self, mesh: vtkUnstructuredGrid ): """ # Building the boundary meshes boundaryMesh, __normals, self.__originalCells = BoundaryMesh.__buildBoundaryMesh( mesh ) - cellsToReorient = filter( - lambda c: mesh.GetCell( c ).GetCellType() == VTK_POLYHEDRON, - map( self.__originalCells.GetValue, range( self.__originalCells.GetNumberOfValues() ) ) ) + cellsToReorient: filter[ int ] = filter( + lambda c: mesh.GetCell( c ).GetCellType() == VTK_POLYHEDRON, # type: ignore[arg-type] + map( + self.__originalCells.GetValue, # type: ignore[attr-defined] + range( self.__originalCells.GetNumberOfValues() ) ) ) reorientedMesh = reorientMesh.reorientMesh( mesh, cellsToReorient ) self.reBoundaryMesh, reNormals, _ = BoundaryMesh.__buildBoundaryMesh( reorientedMesh, consistency=False ) numCells = boundaryMesh.GetNumberOfCells() # Precomputing the underlying cell type - self.__isUnderlyingCellTypeAPolyhedron = numpy.zeros( numCells, dtype=bool ) + self.__isUnderlyingCellTypeAPolyhedron = np.zeros( numCells, dtype=bool ) for ic in range( numCells ): self.__isUnderlyingCellTypeAPolyhedron[ ic ] = mesh.GetCell( - self.__originalCells.GetValue( ic ) ).GetCellType() == VTK_POLYHEDRON + self.__originalCells.GetValue( ic ) ).GetCellType() == VTK_POLYHEDRON # type: ignore[attr-defined] # Precomputing the normals - self.__normals: numpy.ndarray = numpy.empty( ( numCells, 3 ), dtype=numpy.double, - order='C' ) # Do not modify the storage layout + self.__normals: np.ndarray = np.empty( ( numCells, 3 ), dtype=np.double, + order='C' ) # Do not modify the storage layout for ic in range( numCells ): if self.__isUnderlyingCellTypeAPolyhedron[ ic ]: self.__normals[ ic, : ] = reNormals.GetTuple3( ic ) @@ -68,7 +75,7 @@ def __init__( self, mesh: vtkUnstructuredGrid ): @staticmethod def __buildBoundaryMesh( mesh: vtkUnstructuredGrid, - consistency=True ) -> tuple[ vtkUnstructuredGrid, vtkDataArray, vtkDataArray ]: + consistency: bool = True ) -> tuple[ vtkPolyData, vtkFloatArray, vtkIntArray ]: """From a 3d mesh, build the envelope meshes. Args: @@ -76,8 +83,9 @@ def __buildBoundaryMesh( mesh: vtkUnstructuredGrid, consistency (bool, optional): The vtk option passed to the `vtkDataSetSurfaceFilter`. Defaults to True. Returns: - tuple[ vtkUnstructuredGrid, Any, Any ]: A tuple containing the boundary mesh, the normal vectors array, - an array that maps the id of the boundary element to the id of the 3d cell it touches. + tuple[ vtkPolyData, vtkDataArray, vtkDataArray ]: A tuple containing the boundary mesh, + the normal vectors array, and an array that maps the id of the boundary element + to the id of the 3d cell it touches. """ f = vtkDataSetSurfaceFilter() f.PassThroughCellIdsOn() @@ -98,11 +106,11 @@ def __buildBoundaryMesh( mesh: vtkUnstructuredGrid, n.ComputeCellNormalsOn() n.SetInputData( boundaryMesh ) n.Update() - normals: vtkDataArray = n.GetOutput().GetCellData().GetArray( "Normals" ) + normals: vtkFloatArray = n.GetOutput().GetCellData().GetArray( "Normals" ) assert normals assert normals.GetNumberOfComponents() == 3 assert normals.GetNumberOfTuples() == boundaryMesh.GetNumberOfCells() - originalCells: vtkDataArray = boundaryMesh.GetCellData().GetArray( originalCellsKey ) + originalCells: vtkIntArray = boundaryMesh.GetCellData().GetArray( originalCellsKey ) assert originalCells return boundaryMesh, normals, originalCells @@ -133,14 +141,14 @@ def bounds( self, i: int ) -> tuple[ float, float, float, float, float, float ]: """ return self.reBoundaryMesh.GetCell( i ).GetBounds() - def normals( self, i: int ) -> numpy.ndarray: - """The normal of cell `i`. This normal will be directed outwards + def normals( self, i: int ) -> np.ndarray: + """The normal of cell `i`. This normal will be directed outwards. Args: i (int): The boundary cell index. Returns: - numpy.ndarray: The normal as a length-3 numpy array. + np.ndarray: The normal as a length-3 numpy array. """ return self.__normals[ i ] @@ -157,6 +165,7 @@ def GetCell( self, i: int ) -> vtkCell: def GetPoint( self, i: int ) -> tuple[ float, float, float ]: """Point i of the boundary mesh. + Args: i (int): The boundary point index. @@ -166,44 +175,15 @@ def GetPoint( self, i: int ) -> tuple[ float, float, float ]: return self.reBoundaryMesh.GetPoint( i ) @property - def originalCells( self ) -> vtkDataArray: + def originalCells( self ) -> vtkIntArray: """Returns the 2d boundary cell to the 3d cell index of the original mesh. Returns: - vtkDataArray: A 1d array. + vtkIntArray: A 1d array. """ return self.__originalCells -def buildPolyDataForExtrusion( i: int, boundaryMesh: BoundaryMesh ) -> vtkPolyData: - """Creates a vtkPolyData containing the unique cell `i` of the boundary mesh. - - Args: - i (int): The boundary cell index that will eventually be extruded. - boundaryMesh (BoundaryMesh): The boundary mesh containing the cell. - - Returns: - vtkPolyData: The created vtkPolyData. - """ - cell = boundaryMesh.GetCell( i ) - copiedCell = cell.NewInstance() - copiedCell.DeepCopy( cell ) - pointsIdsMapping = [] - for i in range( copiedCell.GetNumberOfPoints() ): - copiedCell.GetPointIds().SetId( i, i ) - pointsIdsMapping.append( cell.GetPointId( i ) ) - polygons = vtkCellArray() - polygons.InsertNextCell( copiedCell ) - points = vtkPoints() - points.SetNumberOfPoints( len( pointsIdsMapping ) ) - for i, v in enumerate( pointsIdsMapping ): - points.SetPoint( i, boundaryMesh.GetPoint( v ) ) - polygonPolyData = vtkPolyData() - polygonPolyData.SetPoints( points ) - polygonPolyData.SetPolys( polygons ) - return polygonPolyData - - def arePointsConformal( pointTolerance: float, cellI: vtkCell, cellJ: vtkCell ) -> bool: """Checks if points of cell `i` matches, one by one, the points of cell `j`. @@ -233,32 +213,67 @@ def arePointsConformal( pointTolerance: float, cellI: vtkCell, cellJ: vtkCell ) return foundPoints == set( range( cellI.GetNumberOfPoints() ) ) -class Extruder: +def buildPolyDataForExtrusion( i: int, boundaryMesh: BoundaryMesh ) -> vtkPolyData: + """Creates a vtkPolyData containing the unique cell `i` of the boundary mesh. + + Args: + i (int): The boundary cell index that will eventually be extruded. + boundaryMesh (BoundaryMesh): The boundary mesh containing the cell. + + Returns: + vtkPolyData: The created vtkPolyData. """ - Computes and stores all the extrusions of the boundary faces. + cell = boundaryMesh.GetCell( i ) + copiedCell = cell.NewInstance() + copiedCell.DeepCopy( cell ) + pointsIdsMapping = [] + for i in range( copiedCell.GetNumberOfPoints() ): + copiedCell.GetPointIds().SetId( i, i ) + pointsIdsMapping.append( cell.GetPointId( i ) ) + polygons = vtkCellArray() + polygons.InsertNextCell( copiedCell ) + points = vtkPoints() + points.SetNumberOfPoints( len( pointsIdsMapping ) ) + for i, v in enumerate( pointsIdsMapping ): + points.SetPoint( i, boundaryMesh.GetPoint( v ) ) + polygonPolyData = vtkPolyData() + polygonPolyData.SetPoints( points ) + polygonPolyData.SetPolys( polygons ) + return polygonPolyData + + +class Extruder: + """Computes and stores all the extrusions of the boundary faces. + The main reason for this class is to be lazy and cache the extrusions. """ - def __init__( self, boundaryMesh: BoundaryMesh, faceTolerance: float ): - self.__extrusions: list[ vtkPolyData ] = [ + def __init__( self, boundaryMesh: BoundaryMesh, faceTolerance: float ) -> None: + """Initializes the extruder. + + Args: + boundaryMesh (BoundaryMesh): Boundary mesh. + faceTolerance (float): Tolerance value. + """ + self.__extrusions: list[ Union[ vtkPolyData, None ] ] = [ None, ] * boundaryMesh.GetNumberOfCells() self.__boundaryMesh = boundaryMesh self.__faceTolerance = faceTolerance - def __extrude( self, polygonPolyData: vtkPolyData, normal: numpy.ndarray ) -> vtkPolyData: + def __extrude( self, polygonPolyData: vtkPolyData, normal: np.ndarray ) -> vtkPolyData: """Extrude the polygon data to create a surface that will be used for intersection. Args: polygonPolyData (vtkPolyData): The data to extrude - normal (numpy.ndarray): The (uniform) direction of the extrusion. + normal (np.ndarray): The (uniform) direction of the extrusion. Returns: vtkPolyData: The extruded surface. """ extruder = vtkLinearExtrusionFilter() extruder.SetExtrusionTypeToVectorExtrusion() - extruder.SetVector( normal ) + extruder.SetVector( float( normal[ 0 ] ), float( normal[ 1 ] ), float( normal[ 2 ] ) ) extruder.SetScaleFactor( self.__faceTolerance / 2. ) extruder.SetInputData( polygonPolyData ) extruder.Update() @@ -282,16 +297,15 @@ def __getitem__( self, i: int ) -> vtkPolyData: return extrusion -def areFacesConformalUsingExtrusions( extrusions: Extruder, i: int, j: int, boundaryMesh: vtkUnstructuredGrid, +def areFacesConformalUsingExtrusions( extrusions: Extruder, i: int, j: int, boundaryMesh: BoundaryMesh, pointTolerance: float ) -> bool: - """ - Tests if two boundary faces are conformal, checking for intersection between their normal extruded volumes. + """Tests if two boundary faces are conformal, checking for intersection between their normal extruded volumes. Args: extrusions (Extruder): The extrusions cache. i (int): The cell index of the first cell. j (int): The cell index of the second cell. - boundaryMesh (vtkUnstructuredGrid): The boundary mesh. + boundaryMesh (BoundaryMesh): The boundary mesh. pointTolerance (float): The point tolerance to consider that two points match. Returns: @@ -332,17 +346,17 @@ def areFacesConformalUsingDistances( i: int, j: int, boundaryMesh: vtkUnstructur Returns: bool: True if the faces are conformal, False otherwise. """ - cpI = boundaryMesh.GetCell( i ).NewInstance() + cpI: vtkCell = boundaryMesh.GetCell( i ).NewInstance() cpI.DeepCopy( boundaryMesh.GetCell( i ) ) - cpJ = boundaryMesh.GetCell( j ).NewInstance() + cpJ: vtkCell = boundaryMesh.GetCell( j ).NewInstance() cpJ.DeepCopy( boundaryMesh.GetCell( j ) ) - def triangulate( cell ): + def triangulate( cell: vtkCell ) -> tuple[ tuple[ int, ...], vtkPoints ]: assert cell.GetCellDimension() == 2 - _pointsIds = vtkIdList() + _pointsIdsList = vtkIdList() _points = vtkPoints() - cell.Triangulate( 0, _pointsIds, _points ) - _pointsIds = tuple( vtkIter( _pointsIds ) ) + cell.Triangulate( 0, _pointsIdsList, _points ) + _pointsIds: tuple[ int, ...] = tuple( vtkIter( _pointsIdsList ) ) assert len( _pointsIds ) % 3 == 0 assert _points.GetNumberOfPoints() % 3 == 0 return _pointsIds, _points @@ -350,19 +364,19 @@ def triangulate( cell ): pointsIdsI, pointsI = triangulate( cpI ) pointsIdsJ, pointsJ = triangulate( cpJ ) - def buildNumpyTriangles( pointsIds ): + def buildNumpyTriangles( pointsIds: tuple[ int, ...] ) -> list[ npt.NDArray[ np.float64 ] ]: __triangles = [] for __i in range( 0, len( pointsIds ), 3 ): __t = [] for __pi in pointsIds[ __i:__i + 3 ]: __t.append( boundaryMesh.GetPoint( __pi ) ) - __triangles.append( numpy.array( __t, dtype=float ) ) + __triangles.append( np.array( __t, dtype=float ) ) return __triangles trianglesI = buildNumpyTriangles( pointsIdsI ) trianglesJ = buildNumpyTriangles( pointsIdsJ ) - minDist = numpy.inf + minDist = np.inf for ti, tj in [ ( ti, tj ) for ti in trianglesI for tj in trianglesJ ]: # Note that here, we compute the exact distance to compare with the threshold. # We could improve by exiting the iterative distance computation as soon as @@ -378,48 +392,93 @@ def buildNumpyTriangles( pointsIds ): return arePointsConformal( pointTolerance, cpI, cpJ ) -def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: - """Checks if the mesh is "conformal" (i.e. if some of its boundary faces may not be too close to each other - without matching nodes). +def computeBoundingBox( boundaryMesh: BoundaryMesh, faceTolerance: float ) -> npt.NDArray: + """Computes the bounding boxes of all boundary cells, inflated by 2 * faceTolerance. Args: - mesh (vtkUnstructuredGrid): The vtk mesh - options (Options): The check options. + boundaryMesh (BoundaryMesh): The boundary mesh. + faceTolerance (float): The tolerance for face proximity. Returns: - Result: The result of the conformity check. + npt.NDArray[ np.float64 ]: An array of bounding boxes for each cell. """ - boundaryMesh = BoundaryMesh( mesh ) - cosTheta = abs( math.cos( numpy.deg2rad( options.angleTolerance ) ) ) - numCells = boundaryMesh.GetNumberOfCells() + # The options are important to directly interact with memory in C++. + boundingBoxes = np.empty( ( boundaryMesh.GetNumberOfCells(), 6 ), dtype=np.double, order="C" ) + for i in range( boundaryMesh.GetNumberOfCells() ): + bb = vtkBoundingBox( boundaryMesh.bounds( i ) ) + bb.Inflate( 2 * faceTolerance ) + assert boundingBoxes[ + i, : ].data.contiguous # Do not modify the storage layout since vtk deals with raw memory here. + boundsList: list[ float ] = [ 0.0 ] * 6 + bb.GetBounds( boundsList ) + boundingBoxes[ i, : ] = boundsList + return boundingBoxes + + +def computeNumberCellsPerNode( boundaryMesh: BoundaryMesh ) -> npt.NDArray[ np.int64 ]: + """Computes the number of cells connected to each node in the boundary mesh. + + Args: + boundaryMesh (BoundaryMesh): The boundary mesh. + Returns: + npt.NDArray[ np.int64 ]: The number of cells per node. + """ # Computing the exact number of cells per node - numCellsPerNode = numpy.zeros( boundaryMesh.GetNumberOfPoints(), dtype=int ) + numCellsPerNode = np.zeros( boundaryMesh.GetNumberOfPoints(), dtype=int ) for ic in range( boundaryMesh.GetNumberOfCells() ): c = boundaryMesh.GetCell( ic ) pointIds = c.GetPointIds() for pointId in vtkIter( pointIds ): numCellsPerNode[ pointId ] += 1 + return numCellsPerNode + + +def buildCellLocator( reBoundaryMesh: vtkPolyData, numberMaxCellPerNode: int ) -> vtkStaticCellLocator: + """Builds a vtkStaticCellLocator for the boundary mesh. + + Args: + reBoundaryMesh (vtkPolyData): The reBoundary mesh. + numberMaxCellPerNode (int): The maximum number of cells per node. + Returns: + vtkStaticCellLocator: The built cell locator. + """ cellLocator = vtkStaticCellLocator() cellLocator.Initialize() - cellLocator.SetNumberOfCellsPerNode( numCellsPerNode.max() ) - cellLocator.SetDataSet( boundaryMesh.reBoundaryMesh ) + cellLocator.SetNumberOfCellsPerNode( numberMaxCellPerNode ) + cellLocator.SetDataSet( reBoundaryMesh ) cellLocator.BuildLocator() + return cellLocator - # Precomputing the bounding boxes. - # The options are important to directly interact with memory in C++. - boundingBoxes = numpy.empty( ( boundaryMesh.GetNumberOfCells(), 6 ), dtype=numpy.double, order="C" ) - for i in range( boundaryMesh.GetNumberOfCells() ): - bb = vtkBoundingBox( boundaryMesh.bounds( i ) ) - bb.Inflate( 2 * options.faceTolerance ) - assert boundingBoxes[ - i, : ].data.contiguous # Do not modify the storage layout since vtk deals with raw memory here. - bb.GetBounds( boundingBoxes[ i, : ] ) - nonConformalCells = [] +def findNonConformalCells( mesh: vtkUnstructuredGrid, options: Options ) -> list[ tuple[ int, int ] ]: + """Finds all pairs of non-conformal boundary faces in the mesh. + + Args: + mesh (vtkUnstructuredGrid): The input mesh. + options (Options): The options for the non-conformal check. + + Returns: + list[ tuple[ int, int ] ]: A list of pairs of non-conformal cell indices. + """ + # Extracts the outer surface of the 3D mesh. + # Ensures that face normals are consistently oriented outward. + boundaryMesh = BoundaryMesh( mesh ) + numCells: int = boundaryMesh.GetNumberOfCells() + + # Used to filter out face pairs that are not facing each other. + cosTheta: float = abs( math.cos( np.deg2rad( options.angleTolerance ) ) ) + + # Prepares extruded volumes of boundary faces for intersection testing. extrusions = Extruder( boundaryMesh, options.faceTolerance ) + + numCellsPerNode = computeNumberCellsPerNode( boundaryMesh ) + boundingBoxes = computeBoundingBox( boundaryMesh, options.faceTolerance ) + cellLocator = buildCellLocator( boundaryMesh.reBoundaryMesh, numCellsPerNode.max() ) + closeCells = vtkIdList() + nonConformalCellsBoundaryId: list[ tuple[ int, int ] ] = [] # Looping on all the pairs of boundary cells. We'll hopefully discard most of the pairs. for i in tqdm( range( numCells ), desc="Non conformal elements" ): cellLocator.FindCellsWithinBounds( boundingBoxes[ i ], closeCells ) @@ -428,19 +487,42 @@ def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: continue # Discarding pairs that are not facing each others (with a threshold). normalI, normalJ = boundaryMesh.normals( i ), boundaryMesh.normals( j ) - if numpy.dot( normalI, normalJ ) > -cosTheta: # opposite directions only (can be facing or not) + if np.dot( normalI, normalJ ) > -cosTheta: # opposite directions only (can be facing or not) continue # At this point, back-to-back and face-to-face pairs of elements are considered. if not areFacesConformalUsingExtrusions( extrusions, i, j, boundaryMesh, options.pointTolerance ): - nonConformalCells.append( ( i, j ) ) + nonConformalCellsBoundaryId.append( ( i, j ) ) # Extracting the original 3d element index (and not the index of the boundary mesh). - tmp = [] - for i, j in nonConformalCells: - tmp.append( ( boundaryMesh.originalCells.GetValue( i ), boundaryMesh.originalCells.GetValue( j ) ) ) + nonConformalCells: list[ tuple[ int, int ] ] = [] + for i, j in nonConformalCellsBoundaryId: + nonConformalCells.append( + ( boundaryMesh.originalCells.GetValue( i ), boundaryMesh.originalCells.GetValue( j ) ) ) + return nonConformalCells + + +def meshAction( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + """Checks if the mesh is "conformal" (i.e. if some of its boundary faces may not be too close to each other without matching nodes). + + Args: + mesh (vtkUnstructuredGrid): The input mesh to analyze. + options (Options): The check options. - return Result( nonConformalCells=tmp ) + Returns: + Result: The Result instance. + """ + nonConformalCells = findNonConformalCells( mesh, options ) + return Result( nonConformalCells=nonConformalCells ) -def action( vtkInputFile: str, options: Options ) -> Result: - mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtkInputFile ) - return __action( mesh, options ) +def action( vtuInputFile: str, options: Options ) -> Result: + """Reads a vtu file and performs the self intersecting elements check on it. + + Args: + vtuInputFile (str): The path to the input VTU file. + options (Options): The options for processing. + + Returns: + Result: The result of the self intersecting elements check. + """ + mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtuInputFile ) + return meshAction( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/reorientMesh.py b/mesh-doctor/src/geos/mesh_doctor/actions/reorientMesh.py similarity index 94% rename from geos-mesh/src/geos/mesh/doctor/actions/reorientMesh.py rename to mesh-doctor/src/geos/mesh_doctor/actions/reorientMesh.py index a9fd11a3..af341a43 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/reorientMesh.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/reorientMesh.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto import networkx import numpy from tqdm import tqdm @@ -6,8 +9,8 @@ from vtkmodules.vtkCommonDataModel import ( VTK_POLYHEDRON, VTK_TRIANGLE, vtkCellArray, vtkPolyData, vtkPolygon, vtkUnstructuredGrid, vtkTetra ) from vtkmodules.vtkFiltersCore import vtkTriangleFilter -from geos.mesh.doctor.actions.vtkPolyhedron import FaceStream, buildFaceToFaceConnectivityThroughEdges -from geos.mesh.doctor.parsing.cliParsing import setupLogger +from geos.mesh_doctor.actions.vtkPolyhedron import FaceStream, buildFaceToFaceConnectivityThroughEdges +from geos.mesh_doctor.parsing.cliParsing import setupLogger from geos.mesh.utils.genericHelpers import toVtkIdList @@ -54,7 +57,7 @@ def __computeVolume( meshPoints: vtkPoints, faceStream: FaceStream ) -> float: tmpBarycenter = numpy.empty( ( faceStream.numSupportPoints, 3 ), dtype=float ) for i, pointId in enumerate( faceStream.supportPointIds ): tmpBarycenter[ i, : ] = meshPoints.GetPoint( pointId ) - barycenter = tmpBarycenter[ :, 0 ].mean(), tmpBarycenter[ :, 1 ].mean(), tmpBarycenter[ :, 2 ].mean() + barycenter = [ tmpBarycenter[ :, 0 ].mean(), tmpBarycenter[ :, 1 ].mean(), tmpBarycenter[ :, 2 ].mean() ] # Looping on all the triangles of the envelope of the polyhedron, creating the matching tetra. # Then the volume of all the tetra are added to get the final polyhedron volume. cellVolume = 0. @@ -68,8 +71,7 @@ def __computeVolume( meshPoints: vtkPoints, faceStream: FaceStream ) -> float: def __selectAndFlipFaces( meshPoints: vtkPoints, colors: dict[ frozenset[ int ], int ], faceStream: FaceStream ) -> FaceStream: - """Given a polyhedra, given that we were able to paint the faces in two colors, - we now need to select which faces/color to flip such that the volume of the element is positive. + """Given a polyhedra, given that we were able to paint the faces in two colors, we now need to select which faces/color to flip such that the volume of the element is positive. Args: meshPoints (vtkPoints): The mesh points, needed to compute the volume. diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/selfIntersectingElements.py b/mesh-doctor/src/geos/mesh_doctor/actions/selfIntersectingElements.py new file mode 100644 index 00000000..1329005f --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/actions/selfIntersectingElements.py @@ -0,0 +1,116 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from dataclasses import dataclass +from vtkmodules.util.numpy_support import vtk_to_numpy +from vtkmodules.vtkFiltersGeneral import vtkCellValidator +from vtkmodules.vtkCommonCore import vtkOutputWindow, vtkFileOutputWindow +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from geos.mesh.io.vtkIO import readUnstructuredGrid + + +@dataclass( frozen=True ) +class Options: + minDistance: float + + +@dataclass( frozen=True ) +class Result: + invalidCellIds: dict[ str, list[ int ] ] + + +def getInvalidCellIds( mesh: vtkUnstructuredGrid, minDistance: float ) -> dict[ str, list[ int ] ]: + """For every cell element in a vtk mesh, check if the cell is invalid regarding 8 specific criteria. + + "wrongNumberOfPoints", "intersectingEdges", "intersectingFaces", "nonContiguousEdges","nonConvex", + "facesAreOrientedIncorrectly", "nonPlanarFacesElements" and "degenerateFacesElements". + + If any of this criteria was met, the cell index is added to a list corresponding to this specific criteria. + The dict with the complete list of cell indices by criteria is returned. + + Args: + mesh (vtkUnstructuredGrid): A vtk grid. + minDistance (float): Minimum distance in the computation. + + Returns: + dict[ str, list[ int ] ]: + { + "wrongNumberOfPointsElements": [ 10, 34, ... ], + "intersectingEdgesElements": [ ... ], + "intersectingFacesElements": [ ... ], + "nonContiguousEdgesElements": [ ... ], + "nonConvexElements": [ ... ], + "facesOrientedIncorrectlyElements": [ ... ], + "nonPlanarFacesElements": [ ... ], + "degenerateFacesElements": [ ... ] + } + """ + errOut = vtkFileOutputWindow() + errOut.SetFileName( "/dev/null" ) + vtkStdErrOut = vtkOutputWindow() + vtkStdErrOut.SetInstance( errOut ) + + # Different types of cell invalidity are defined as hexadecimal values, specific to vtkCellValidator + # Complete set of validity checks available in vtkCellValidator + errorMasks: dict[ str, int ] = { + "wrongNumberOfPointsElements": 0x01, # 0000 0001 + "intersectingEdgesElements": 0x02, # 0000 0010 + "intersectingFacesElements": 0x04, # 0000 0100 + "nonContiguousEdgesElements": 0x08, # 0000 1000 + "nonConvexElements": 0x10, # 0001 0000 + "facesOrientedIncorrectlyElements": 0x20, # 0010 0000 + "nonPlanarFacesElements": 0x40, # 0100 0000 + "degenerateFacesElements": 0x80, # 1000 0000 + } + + # The results can be stored in a dictionary where keys are the error names + # and values are the lists of cell indices with that error. + # We can initialize it directly from the keys of our errorMasks dictionary. + invalidCellIds: dict[ str, list[ int ] ] = { _errorName: [] for _errorName in errorMasks } + + f = vtkCellValidator() + f.SetTolerance( minDistance ) + f.SetInputData( mesh ) + f.Update() # executes the filter + output = f.GetOutput() + + validity = output.GetCellData().GetArray( "ValidityState" ) + assert validity is not None + # array of np.int16 that combines the flags using a bitwise OR operation for each cell index. + validity = vtk_to_numpy( validity ) + for cellIndex, validityFlag in enumerate( validity ): + if not validityFlag: # Skip valid cells ( validityFlag == 0 or 0000 0000 ) + continue + for errorName, errorMask in errorMasks.items(): # Check only invalid cells against all possible errors. + if validityFlag & errorMask: + invalidCellIds[ errorName ].append( cellIndex ) + + return invalidCellIds + + +def meshAction( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + """Performs the self intersecting elements check on a vtkUnstructuredGrid. + + Args: + mesh (vtkUnstructuredGrid): The input mesh to analyze. + options (Options): The options for processing. + + Returns: + Result: The result of the self intersecting elements check. + """ + invalidCellIds: dict[ str, list[ int ] ] = getInvalidCellIds( mesh, options.minDistance ) + return Result( invalidCellIds ) + + +def action( vtuInputFile: str, options: Options ) -> Result: + """Reads a vtu file and performs the self intersecting elements check on it. + + Args: + vtuInputFile (str): The path to the input VTU file. + options (Options): The options for processing. + + Returns: + Result: The result of the self intersecting elements check. + """ + mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtuInputFile ) + return meshAction( mesh, options ) diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/supportedElements.py b/mesh-doctor/src/geos/mesh_doctor/actions/supportedElements.py new file mode 100644 index 00000000..23c6a28c --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/actions/supportedElements.py @@ -0,0 +1,203 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from dataclasses import dataclass +import multiprocessing +import networkx +from numpy import ones +from tqdm import tqdm +from typing import Iterable, Mapping, Optional +from vtkmodules.util.numpy_support import vtk_to_numpy +from vtkmodules.vtkCommonCore import vtkIdList +from vtkmodules.vtkCommonDataModel import ( vtkCellTypes, vtkUnstructuredGrid, VTK_HEXAGONAL_PRISM, VTK_HEXAHEDRON, + VTK_LINE, VTK_PENTAGONAL_PRISM, VTK_POLYGON, VTK_POLYHEDRON, VTK_PYRAMID, + VTK_QUAD, VTK_TETRA, VTK_TRIANGLE, VTK_VERTEX, VTK_VOXEL, VTK_WEDGE ) +from geos.mesh_doctor.actions.vtkPolyhedron import buildFaceToFaceConnectivityThroughEdges, FaceStream +from geos.mesh_doctor.parsing.cliParsing import setupLogger +from geos.mesh.io.vtkIO import readUnstructuredGrid +from geos.mesh.utils.genericHelpers import vtkIter + + +@dataclass( frozen=True ) +class Options: + nproc: int + chunkSize: int + + +@dataclass( frozen=True ) +class Result: + unsupportedStdElementsTypes: list[ str ] # list of unsupported types + unsupportedPolyhedronElements: frozenset[ + int ] # list of polyhedron elements that could not be converted to supported std elements + + +# for multiprocessing, vtkUnstructuredGrid cannot be pickled. Let's use a global variable instead. +# Global variable to be set in each worker process +MESH: Optional[ vtkUnstructuredGrid ] = None + + +def initWorker( meshToInit: vtkUnstructuredGrid ) -> None: + """Initializer for each worker process to set the global mesh.""" + global MESH + MESH = meshToInit + + +supportedCellTypes: set[ int ] = { + VTK_HEXAGONAL_PRISM, VTK_HEXAHEDRON, VTK_LINE, VTK_PENTAGONAL_PRISM, VTK_POLYGON, VTK_POLYHEDRON, VTK_PYRAMID, + VTK_QUAD, VTK_TETRA, VTK_TRIANGLE, VTK_VERTEX, VTK_VOXEL, VTK_WEDGE +} + + +class IsPolyhedronConvertible: + + def __init__( self ) -> None: + """Initializer for the IsPolyhedronConvertible class.""" + + def buildPrismGraph( n: int, name: str ) -> networkx.Graph: + """Builds the face to face connectivities (through edges) for prism graphs. + + Args: + n (int): The number of nodes of the base (e.g., pentagonal prism gets n = 5). + name (str): A human-readable name for logging purposes. + + Returns: + networkx.Graph: A graph instance representing the prism. + """ + tmp = networkx.cycle_graph( n ) + for node in range( n ): + tmp.add_edge( node, n ) + tmp.add_edge( node, n + 1 ) + tmp.name = name + return tmp + + # Building the reference graphs + tetGraph = networkx.complete_graph( 4 ) + tetGraph.name = "Tetrahedron" + pyrGraph = buildPrismGraph( 4, "Pyramid" ) + pyrGraph.remove_node( 5 ) # Removing a node also removes its associated edges. + self.__referenceGraphs: Mapping[ int, Iterable[ networkx.Graph ] ] = { + 4: ( tetGraph, ), + 5: ( pyrGraph, buildPrismGraph( 3, "Wedge" ) ), + 6: ( buildPrismGraph( 4, "Hexahedron" ), ), + 7: ( buildPrismGraph( 5, "Prism5" ), ), + 8: ( buildPrismGraph( 6, "Prism6" ), ), + 9: ( buildPrismGraph( 7, "Prism7" ), ), + 10: ( buildPrismGraph( 8, "Prism8" ), ), + 11: ( buildPrismGraph( 9, "Prism9" ), ), + 12: ( buildPrismGraph( 10, "Prism10" ), ), + 13: ( buildPrismGraph( 11, "Prism11" ), ), + } + + def __isPolyhedronSupported( self, faceStream: FaceStream ) -> str: + """Checks if a polyhedron can be converted into a supported cell. + + Args: + faceStream: The polyhedron. + + Returns: + str: The name of the supported type or an empty string. + """ + cellGraph = buildFaceToFaceConnectivityThroughEdges( faceStream, addCompatibility=True ) + for referenceGraph in self.__referenceGraphs[ cellGraph.order() ]: + if networkx.is_isomorphic( referenceGraph, cellGraph ): + return str( referenceGraph.name ) + return "" + + def __call__( self, ic: int ) -> int: + """Check if a vtk polyhedron cell can be converted into a supported GEOS element. + + Args: + ic (int): The index of the element. + + Returns: + int: -1 if the polyhedron vtk element can be converted into a supported element type, the index otherwise. + """ + global MESH + assert MESH is not None + if MESH.GetCellType( ic ) != VTK_POLYHEDRON: + return -1 + ptIds = vtkIdList() + MESH.GetFaceStream( ic, ptIds ) + faceStream = FaceStream.buildFromVtkIdList( ptIds ) + convertedTypeName = self.__isPolyhedronSupported( faceStream ) + if convertedTypeName: + setupLogger.debug( f"Polyhedron cell {ic} can be converted into \"{convertedTypeName}\"" ) + return -1 + else: + setupLogger.debug( f"Polyhedron cell {ic} cannot be converted into any supported element." ) + return ic + + +def findUnsupportedStdElementsTypes( mesh: vtkUnstructuredGrid ) -> list[ str ]: + """Find unsupported standard element types in the mesh. + + Args: + mesh (vtkUnstructuredGrid): The input mesh to analyze. + + Returns: + list[ str ]: List of unsupported element type descriptions. + """ + if hasattr( mesh, "GetDistinctCellTypesArray" ): # For more recent versions of vtk. + uniqueCellTypes = set( vtk_to_numpy( mesh.GetDistinctCellTypesArray() ) ) + else: + _vtkCellTypes = vtkCellTypes() + mesh.GetCellTypes( _vtkCellTypes ) + uniqueCellTypes = set( vtkIter( _vtkCellTypes ) ) + resultValues: set[ int ] = uniqueCellTypes - supportedCellTypes + results = [ f"Type {i}: {vtkCellTypes.GetClassNameFromTypeId( i )}" for i in frozenset( resultValues ) ] + return results + + +def findUnsupportedPolyhedronElements( mesh: vtkUnstructuredGrid, options: Options ) -> list[ int ]: + """Find unsupported polyhedron elements in the mesh. + + Args: + mesh (vtkUnstructuredGrid): The input mesh to analyze. + options (Options): The options for processing. + + Returns: + list[ int ]: List of element indices for unsupported polyhedrons. + """ + # Dealing with polyhedron elements. + numCells: int = mesh.GetNumberOfCells() + result = ones( numCells, dtype=int ) * -1 + # Use the initializer to set up each worker process + # Pass the mesh to the initializer + with multiprocessing.Pool( processes=options.nproc, initializer=initWorker, initargs=( mesh, ) ) as pool: + # Pass a mesh-free instance of the class to the workers. + # The MESH global will already be set in each worker. + generator = pool.imap_unordered( IsPolyhedronConvertible(), range( numCells ), chunksize=options.chunkSize ) + for i, val in enumerate( tqdm( generator, total=numCells, desc="Testing support for elements" ) ): + result[ i ] = val + + return [ i for i in result if i > -1 ] + + +def meshAction( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + """Performs the supported elements check on a vtkUnstructuredGrid. + + Args: + mesh (vtkUnstructuredGrid): The input mesh to analyze. + options (Options): The options for processing. + + Returns: + Result: The result of the supported elements check. + """ + unsupportedStdElementsTypes: list[ str ] = findUnsupportedStdElementsTypes( mesh ) + unsupportedPolyhedronElements: list[ int ] = findUnsupportedPolyhedronElements( mesh, options ) + return Result( unsupportedStdElementsTypes=unsupportedStdElementsTypes, + unsupportedPolyhedronElements=frozenset( unsupportedPolyhedronElements ) ) + + +def action( vtuInputFile: str, options: Options ) -> Result: + """Reads a vtu file and performs the supported elements check on it. + + Args: + vtuInputFile (str): The path to the input VTU file. + options (Options): The options for processing. + + Returns: + Result: The result of the supported elements check. + """ + mesh: vtkUnstructuredGrid = readUnstructuredGrid( vtuInputFile ) + return meshAction( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/triangleDistance.py b/mesh-doctor/src/geos/mesh_doctor/actions/triangleDistance.py similarity index 88% rename from geos-mesh/src/geos/mesh/doctor/actions/triangleDistance.py rename to mesh-doctor/src/geos/mesh_doctor/actions/triangleDistance.py index a8e309fe..03c7fc1a 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/triangleDistance.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/triangleDistance.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto import itertools from math import sqrt import numpy @@ -7,6 +10,7 @@ def __divClamp( num: float, den: float ) -> float: """Computes the division `num / den`. and clamps the result between 0 and 1. + If `den` is zero, the result of the division is set to 0. Args: @@ -67,14 +71,14 @@ def distanceBetweenTwoSegments( x0: numpy.ndarray, d0: numpy.ndarray, x1: numpy. # Step 3: compute t1 for point on line 1 closest to point at t0. t1: float = __divClamp( t0 * R - S1, D1 ) # Eq (10, right) sol1: numpy.ndarray = x1 + t1 * d1 # Eq (3) - t0: float = __divClamp( t1 * R + S0, D0 ) # Eq (10, left) + t0 = __divClamp( t1 * R + S0, D0 ) # Eq (10, left) sol0: numpy.ndarray = x0 + t0 * d0 # Eq (4) return sol0, sol1 def __computeNodesToTriangleDistance( - tri0: numpy.ndarray, edges0, tri1: numpy.ndarray + tri0: numpy.ndarray, edges0: numpy.ndarray, tri1: numpy.ndarray ) -> tuple[ Union[ float, None ], Union[ numpy.ndarray, None ], Union[ numpy.ndarray, None ], bool ]: """Computes the distance from nodes of `tri1` points onto `tri0`. @@ -105,28 +109,29 @@ def __computeNodesToTriangleDistance( # If so, let's take the closest point. point: int = -1 if numpy.all( tri1Proj > 0 ): - point = numpy.argmin( tri1Proj ) + point = int( numpy.argmin( tri1Proj ) ) elif numpy.all( tri1Proj < 0 ): - point = numpy.argmax( tri1Proj ) + point = int( numpy.argmax( tri1Proj ) ) # So if `tri1` is actually "on one side", # point `tri1[point]` is candidate to be the closest point. if point > -1: areDisjoint = True # But we must check that its projection is inside `tri0`. - if numpy.dot( tri1[ point ] - tri0[ 0 ], numpy.cross( tri0Normal, edges0[ 0 ] ) ) > 0: - if numpy.dot( tri1[ point ] - tri0[ 1 ], numpy.cross( tri0Normal, edges0[ 1 ] ) ) > 0: - if numpy.dot( tri1[ point ] - tri0[ 2 ], numpy.cross( tri0Normal, edges0[ 2 ] ) ) > 0: - # It is! - sol0 = tri1[ point ] - sol1 = tri1[ point ] + ( tri1Proj[ point ] / tri0NormalNorm ) * tri0Normal - return norm( sol1 - sol0 ), sol0, sol1, areDisjoint + if numpy.dot( tri1[ point ] - tri0[ 0 ], numpy.cross( tri0Normal, edges0[ 0 ] ) ) > 0 and \ + numpy.dot( tri1[ point ] - tri0[ 1 ], numpy.cross( tri0Normal, edges0[ 1 ] ) ) > 0 and \ + numpy.dot( tri1[ point ] - tri0[ 2 ], numpy.cross( tri0Normal, edges0[ 2 ] ) ) > 0: + # It is! + sol0 = tri1[ point ] + sol1 = tri1[ point ] + ( tri1Proj[ point ] / tri0NormalNorm ) * tri0Normal + return float( norm( sol1 - sol0 ) ), sol0, sol1, areDisjoint return None, None, None, areDisjoint def distanceBetweenTwoTriangles( tri0: numpy.ndarray, tri1: numpy.ndarray ) -> tuple[ float, numpy.ndarray, numpy.ndarray ]: """Returns the minimum distance between two triangles, and the two points where this minimum occurs. + If the two triangles touch, then distance is exactly 0. But the two points are dummy and cannot be used as contact points (they are still though). @@ -183,12 +188,14 @@ def distanceBetweenTwoTriangles( tri0: numpy.ndarray, areDisjoint = True # No edge pair contained the closest points. # Checking the node/face situation. - distance, sol0, sol1, areDisjointTmp = __computeNodesToTriangleDistance( tri0, edges0, tri1 ) + # yapf: disable + distance, sol0, sol1, areDisjointTmp = __computeNodesToTriangleDistance( tri0, edges0, tri1 ) # type: ignore[ assignment ] if distance: return distance, sol0, sol1 areDisjoint = areDisjoint or areDisjointTmp - distance, sol0, sol1, areDisjointTmp = __computeNodesToTriangleDistance( tri1, edges1, tri0 ) + distance, sol0, sol1, areDisjointTmp = __computeNodesToTriangleDistance( tri1, edges1, tri0 ) # type: ignore[ assignment ] + # yapf: enable if distance: return distance, sol0, sol1 areDisjoint = areDisjoint or areDisjointTmp diff --git a/geos-mesh/src/geos/mesh/doctor/actions/vtkPolyhedron.py b/mesh-doctor/src/geos/mesh_doctor/actions/vtkPolyhedron.py similarity index 82% rename from geos-mesh/src/geos/mesh/doctor/actions/vtkPolyhedron.py rename to mesh-doctor/src/geos/mesh_doctor/actions/vtkPolyhedron.py index 3093d96a..40cb5e6e 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/vtkPolyhedron.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/vtkPolyhedron.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto from collections import defaultdict from dataclasses import dataclass import networkx @@ -17,8 +20,9 @@ class Result: def parseFaceStream( ids: vtkIdList ) -> Sequence[ Sequence[ int ] ]: - """Parses the face stream raw information and converts it into a tuple of tuple of integers, - each tuple of integer being the nodes of a face. + """Parses the face stream raw information and converts it into a tuple of tuple of integers. + + Each tuple of integer corresponds to the nodes of a face. Args: ids (vtkIdList): The raw vtk face stream. @@ -33,7 +37,7 @@ def parseFaceStream( ids: vtkIdList ) -> Sequence[ Sequence[ int ] ]: while True: numNodes = next( it ) tmp = [] - for i in range( numNodes ): + for _ in range( numNodes ): tmp.append( next( it ) ) result.append( tuple( tmp ) ) except StopIteration: @@ -45,18 +49,17 @@ def parseFaceStream( ids: vtkIdList ) -> Sequence[ Sequence[ int ] ]: class FaceStream: - """ - Helper class to manipulate the vtk face streams. - """ + """Helper class to manipulate the vtk face streams.""" - def __init__( self, data: Sequence[ Sequence[ int ] ] ): + def __init__( self, data: Sequence[ Sequence[ int ] ] ) -> None: + """Initializes the FaceStream instance.""" # self.__data contains the list of faces nodes, like it appears in vtk face streams. # Except that the additional size information is removed # in favor of the __len__ of the containers. self.__data: Sequence[ Sequence[ int ] ] = data @staticmethod - def buildFromVtkIdList( ids: vtkIdList ): + def buildFromVtkIdList( ids: vtkIdList ) -> "FaceStream": """Builds a FaceStream from the raw vtk face stream. Args: @@ -78,7 +81,7 @@ def faceNodes( self ) -> Iterable[ Sequence[ int ] ]: @property def numFaces( self ) -> int: - """Number of faces in the face stream + """Number of faces in the face stream. Returns: int: The number of faces. @@ -133,6 +136,7 @@ def flipFaces( self, faceIndices: Collection[ int ] ) -> "FaceStream": def dump( self ) -> Sequence[ int ]: """Returns the face stream awaited by vtk, but in a python container. + The content can be used, once converted to a vtkIdList, to define another polyhedron in vtk. Returns: @@ -144,7 +148,8 @@ def dump( self ) -> Sequence[ int ]: result += faceNodes return tuple( result ) - def __repr__( self ): + def __repr__( self ) -> str: + """String representation of the face stream.""" result = [ str( len( self.__data ) ) ] for faceNodes in self.__data: result.append( str( len( faceNodes ) ) ) @@ -152,8 +157,9 @@ def __repr__( self ): return ",\n".join( result ) -def buildFaceToFaceConnectivityThroughEdges( faceStream: FaceStream, addCompatibility=False ) -> networkx.Graph: +def buildFaceToFaceConnectivityThroughEdges( faceStream: FaceStream, addCompatibility: bool = False ) -> networkx.Graph: """Given a face stream/polyhedron, builds the connections between the faces. + Those connections happen when two faces share an edge. Args: @@ -170,26 +176,26 @@ def buildFaceToFaceConnectivityThroughEdges( faceStream: FaceStream, addCompatib edgesToFaceIndices: dict[ frozenset[ int ], list[ int ] ] = defaultdict( list ) for faceIndex, faceNodes in enumerate( faceStream.faceNodes ): # Each edge is defined by two nodes. We do a small trick to loop on consecutive points. - faceIndices: tuple[ int, int ] - for faceIndices in zip( faceNodes, faceNodes[ 1: ] + ( faceNodes[ 0 ], ) ): - edgesToFaceIndices[ frozenset( faceIndices ) ].append( faceIndex ) + faceNodesTuple = tuple( faceNodes ) + for edgeNodes in zip( faceNodesTuple, faceNodesTuple[ 1: ] + ( faceNodesTuple[ 0 ], ) ): + edgesToFaceIndices[ frozenset( edgeNodes ) ].append( faceIndex ) # We are doing here some small validations w.r.t. the connections of the faces # which may only make sense in the context of numerical simulations. # As such, an error will be thrown in case the polyhedron is not closed. # So there may be a lack of absolute genericity, and the code may evolve if needed. - for faceIndices in edgesToFaceIndices.values(): - assert len( faceIndices ) == 2 + for connectedFaces in edgesToFaceIndices.values(): + assert len( connectedFaces ) == 2 # Computing the graph degree for validation degrees: dict[ int, int ] = defaultdict( int ) - for faceIndices in edgesToFaceIndices.values(): - for faceIndex in faceIndices: + for connectedFaces in edgesToFaceIndices.values(): + for faceIndex in connectedFaces: degrees[ faceIndex ] += 1 for faceIndex, degree in degrees.items(): assert len( faceStream[ faceIndex ] ) == degree # Validation that there is one unique edge connecting two faces. faceIndicesToEdgeIndex = defaultdict( list ) - for edgeIndex, faceIndices in edgesToFaceIndices.items(): - faceIndicesToEdgeIndex[ frozenset( faceIndices ) ].append( edgeIndex ) + for edgeIndex, connectedFaces in edgesToFaceIndices.items(): + faceIndicesToEdgeIndex[ frozenset( connectedFaces ) ].append( edgeIndex ) for edgeIndices in faceIndicesToEdgeIndex.values(): assert len( edgeIndices ) == 1 # Connecting the faces. Neighbor faces with consistent normals (i.e. facing both inward or outward) @@ -198,10 +204,10 @@ def buildFaceToFaceConnectivityThroughEdges( faceStream: FaceStream, addCompatib # will consistently point outward. graph = networkx.Graph() graph.add_nodes_from( range( faceStream.numFaces ) ) - for edge, faceIndices in edgesToFaceIndices.items(): - faceIndex0, faceIndex1 = faceIndices - faceNodes0 = faceStream[ faceIndex0 ] + ( faceStream[ faceIndex0 ][ 0 ], ) - faceNodes1 = faceStream[ faceIndex1 ] + ( faceStream[ faceIndex1 ][ 0 ], ) + for edge, connectedFaces in edgesToFaceIndices.items(): + faceIndex0, faceIndex1 = connectedFaces + faceNodes0 = tuple( faceStream[ faceIndex0 ] ) + ( faceStream[ faceIndex0 ][ 0 ], ) + faceNodes1 = tuple( faceStream[ faceIndex1 ] ) + ( faceStream[ faceIndex1 ][ 0 ], ) node0, node1 = edge order0 = 1 if faceNodes0[ faceNodes0.index( node0 ) + 1 ] == node1 else -1 order1 = 1 if faceNodes1[ faceNodes1.index( node0 ) + 1 ] == node1 else -1 diff --git a/mesh-doctor/src/geos/mesh_doctor/meshDoctor.py b/mesh-doctor/src/geos/mesh_doctor/meshDoctor.py new file mode 100644 index 00000000..417d3cf0 --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/meshDoctor.py @@ -0,0 +1,33 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +import sys +from geos.mesh_doctor.parsing import ActionHelper +from geos.mesh_doctor.parsing.cliParsing import parseAndSetVerbosity, setupLogger +from geos.mesh_doctor.register import registerParsingActions + + +def main() -> None: + """Main function for mesh-doctor CLI.""" + parseAndSetVerbosity( sys.argv ) + mainParser, allActions, allActionsHelpers = registerParsingActions() + args = mainParser.parse_args( sys.argv[ 1: ] ) + + # Extract vtuInputFile from parsed arguments + vtuInputFile = getattr( args, 'vtuInputFile', None ) + if vtuInputFile: + setupLogger.info( f"Working on mesh \"{vtuInputFile}\"." ) + + actionOptions = allActionsHelpers[ args.subparsers ].convert( vars( args ) ) + try: + action = allActions[ args.subparsers ] + except KeyError: + setupLogger.error( f"Action {args.subparsers} is not a valid action." ) + sys.exit( 1 ) + helper: ActionHelper = allActionsHelpers[ args.subparsers ] + result = action( vtuInputFile, actionOptions ) + helper.displayResults( actionOptions, result ) + + +if __name__ == '__main__': + main() diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/__init__.py b/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py similarity index 81% rename from geos-mesh/src/geos/mesh/doctor/parsing/__init__.py rename to mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py index 9c747332..7169fad3 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/__init__.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto import argparse from dataclasses import dataclass from typing import Callable, Any diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/_sharedChecksParsingLogic.py b/mesh-doctor/src/geos/mesh_doctor/parsing/_sharedChecksParsingLogic.py similarity index 73% rename from geos-mesh/src/geos/mesh/doctor/parsing/_sharedChecksParsingLogic.py rename to mesh-doctor/src/geos/mesh_doctor/parsing/_sharedChecksParsingLogic.py index a8c6cfb1..21a1d090 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/_sharedChecksParsingLogic.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/_sharedChecksParsingLogic.py @@ -1,13 +1,16 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations import argparse from copy import deepcopy from dataclasses import dataclass -from typing import Type, Any -from geos.mesh.doctor.actions.allChecks import Options as AllChecksOptions -from geos.mesh.doctor.actions.allChecks import Result as AllChecksResult -from geos.mesh.doctor.parsing.cliParsing import parseCommaSeparatedString, setupLogger +from typing import Any, Type, Union +from geos.mesh_doctor.actions.allChecks import Options as AllChecksOptions +from geos.mesh_doctor.actions.allChecks import Result as AllChecksResult +from geos.mesh_doctor.parsing.cliParsing import parseCommaSeparatedString, setupLogger, addVtuInputFileArgument -# --- Data Structure for Check Features --- @dataclass( frozen=True ) class CheckFeature: """A container for a check's configuration and associated classes.""" @@ -18,7 +21,6 @@ class CheckFeature: display: Type[ Any ] -# --- Argument Parser Constants --- CHECKS_TO_DO_ARG = "checksToPerform" PARAMETERS_ARG = "setParameters" @@ -34,27 +36,26 @@ def _generateParametersHelp( orderedCheckNames: list[ str ], checkFeaturesConfig return helpText -def getOptionsUsedMessage( optionsUsed: dataclass ) -> str: +def getOptionsUsedMessage( optionsUsed: object ) -> str: """Dynamically generates the description of every parameter used when loaching a check. Args: - optionsUsed (dataclass) + optionsUsed (dataclass object): The options dataclass used for a specific check. Returns: str: A message like "Parameters used: ( param1:value1 param2:value2 )" for as many paramters found. """ optionsMsg: str = "Parameters used: (" - for attrName in optionsUsed.__dataclass_fields__: - attrValue = getattr( optionsUsed, attrName ) - optionsMsg += f" {attrName} = {attrValue}" + if hasattr( optionsUsed, "__dataclass_fields__" ): + for attrName in optionsUsed.__dataclass_fields__: + attrValue = getattr( optionsUsed, attrName ) + optionsMsg += f" {attrName} = {attrValue}" return optionsMsg + " )." -# --- Generic Argument Parser Setup --- def fillSubparser( subparsers: argparse._SubParsersAction, subparserName: str, helpMessage: str, orderedCheckNames: list[ str ], checkFeaturesConfig: dict[ str, CheckFeature ] ) -> None: - """ - Fills a subparser with arguments for performing a set of checks. + """Fills a subparser with arguments for performing a set of checks. Args: subparsers: The subparsers action from argparse. @@ -66,7 +67,7 @@ def fillSubparser( subparsers: argparse._SubParsersAction, subparserName: str, h parser = subparsers.add_parser( subparserName, help=helpMessage, formatter_class=argparse.ArgumentDefaultsHelpFormatter ) - + addVtuInputFileArgument( parser ) parametersHelp: str = _generateParametersHelp( orderedCheckNames, checkFeaturesConfig ) parser.add_argument( f"--{CHECKS_TO_DO_ARG}", @@ -86,24 +87,38 @@ def fillSubparser( subparsers: argparse._SubParsersAction, subparserName: str, h f"Example: --{PARAMETERS_ARG} parameter_name:10.5,other_param:25" ) ) -def convert( parsedArgs: argparse.Namespace, orderedCheckNames: list[ str ], +def convert( parsedArgs: Union[ dict[ str, Any ], argparse.Namespace ], orderedCheckNames: list[ str ], checkFeaturesConfig: dict[ str, CheckFeature ] ) -> AllChecksOptions: - """ - Converts parsed command-line arguments into an AllChecksOptions object based on the provided configuration. + """Converts parsed command-line arguments into an AllChecksOptions object based on the provided configuration. + + Args: + parsedArgs: Parsed command-line arguments (dict or Namespace). + orderedCheckNames (list[ str ]): Ordered list of check names. + checkFeaturesConfig (dict[ str, CheckFeature ]): Configuration dictionary for check features. + + Raises: + ValueError: If no valid checks are selected. + + Returns: + AllChecksOptions: The options for all checks to be performed. """ # 1. Determine which checks to perform - if not parsedArgs[ CHECKS_TO_DO_ARG ]: # handles default and if user explicitly provides --checksToPerform "" + # Support both dict and Namespace + if isinstance( parsedArgs, dict ): + checksToDo = parsedArgs.get( CHECKS_TO_DO_ARG, "" ) + else: + checksToDo = getattr( parsedArgs, CHECKS_TO_DO_ARG ) + if not checksToDo: finalSelectedCheckNames: list[ str ] = deepcopy( orderedCheckNames ) setupLogger.info( "All configured checks will be performed by default." ) else: - userChecks = parseCommaSeparatedString( parsedArgs[ CHECKS_TO_DO_ARG ] ) - finalSelectedCheckNames = list() + userChecks = parseCommaSeparatedString( checksToDo ) + finalSelectedCheckNames = [] for name in userChecks: if name not in checkFeaturesConfig: setupLogger.warning( f"Check '{name}' does not exist. Choose from: {orderedCheckNames}." ) elif name not in finalSelectedCheckNames: finalSelectedCheckNames.append( name ) - if not finalSelectedCheckNames: raise ValueError( "No valid checks were selected. No operations will be configured." ) @@ -111,10 +126,15 @@ def convert( parsedArgs: argparse.Namespace, orderedCheckNames: list[ str ], defaultParams = { name: feature.defaultParams.copy() for name, feature in checkFeaturesConfig.items() } finalCheckParams = { name: defaultParams[ name ] for name in finalSelectedCheckNames } - if not parsedArgs[ PARAMETERS_ARG ]: # handles default and if user explicitly provides --setParameters "" + # Support both dict and Namespace + if isinstance( parsedArgs, dict ): + parametersArg = parsedArgs.get( PARAMETERS_ARG, "" ) + else: + parametersArg = getattr( parsedArgs, PARAMETERS_ARG ) + if not parametersArg: setupLogger.info( "Default configuration of parameters adopted for every check to perform." ) else: - setParameters = parseCommaSeparatedString( parsedArgs[ PARAMETERS_ARG ] ) + setParameters = parseCommaSeparatedString( parametersArg ) for param in setParameters: if ':' not in param: setupLogger.warning( f"Parameter '{param}' is not in 'name:value' format. Skipping." ) @@ -140,8 +160,8 @@ def convert( parsedArgs: argparse.Namespace, orderedCheckNames: list[ str ], finalCheckParams[ checkNameKey ][ name ] = valueFloat # 3. Instantiate Options objects for the selected checks - individualCheckOptions: dict[ str, Any ] = dict() - individualCheckDisplay: dict[ str, Any ] = dict() + individualCheckOptions: dict[ str, Any ] = {} + individualCheckDisplay: dict[ str, Any ] = {} for checkName in list( finalCheckParams.keys() ): params = finalCheckParams[ checkName ] @@ -160,7 +180,12 @@ def convert( parsedArgs: argparse.Namespace, orderedCheckNames: list[ str ], # Generic display of Results def displayResults( options: AllChecksOptions, result: AllChecksResult ) -> None: - """Displays the results of all the checks that have been performed.""" + """Displays the results of all the checks that have been performed. + + Args: + options (AllChecksOptions): The options used for the checks. + result (AllChecksResult): The result of the checks. + """ if not options.checksToPerform: setupLogger.results( "No checks were performed or all failed during configuration." ) return diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/allChecksParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/allChecksParsing.py similarity index 77% rename from geos-mesh/src/geos/mesh/doctor/parsing/allChecksParsing.py rename to mesh-doctor/src/geos/mesh_doctor/parsing/allChecksParsing.py index 2411f140..a4d200da 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/allChecksParsing.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/allChecksParsing.py @@ -1,9 +1,14 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations import argparse from copy import deepcopy -from geos.mesh.doctor.parsing._sharedChecksParsingLogic import ( CheckFeature, convert as sharedConvert, fillSubparser - as sharedFillSubparser, displayResults ) -from geos.mesh.doctor.actions.allChecks import Options as AllChecksOptions -from geos.mesh.doctor.parsing import ( +from typing import Any +from geos.mesh_doctor.actions.allChecks import Options as AllChecksOptions +from geos.mesh_doctor.parsing._sharedChecksParsingLogic import ( CheckFeature, convert as sharedConvert, fillSubparser + as sharedFillSubparser, displayResults ) # noqa: F401 +from geos.mesh_doctor.parsing import ( ALL_CHECKS, COLLOCATES_NODES, ELEMENT_VOLUMES, @@ -11,11 +16,11 @@ SELF_INTERSECTING_ELEMENTS, SUPPORTED_ELEMENTS, ) -from geos.mesh.doctor.parsing import collocatedNodesParsing as cnParser -from geos.mesh.doctor.parsing import elementVolumesParsing as evParser -from geos.mesh.doctor.parsing import nonConformalParsing as ncParser -from geos.mesh.doctor.parsing import selfIntersectingElementsParsing as sieParser -from geos.mesh.doctor.parsing import supportedElementsParsing as seParser +from geos.mesh_doctor.parsing import collocatedNodesParsing as cnParser +from geos.mesh_doctor.parsing import elementVolumesParsing as evParser +from geos.mesh_doctor.parsing import nonConformalParsing as ncParser +from geos.mesh_doctor.parsing import selfIntersectingElementsParsing as sieParser +from geos.mesh_doctor.parsing import supportedElementsParsing as seParser # Ordered list of check names for this configuration ORDERED_CHECK_NAMES = [ @@ -70,7 +75,7 @@ def fillSubparser( subparsers: argparse._SubParsersAction ) -> None: checkFeaturesConfig=CHECK_FEATURES_CONFIG ) -def convert( parsedArgs: argparse.Namespace ) -> AllChecksOptions: +def convert( parsedArgs: dict[ str, Any ] ) -> AllChecksOptions: """Converts arguments by calling the shared logic with the 'allChecks' configuration.""" return sharedConvert( parsedArgs=parsedArgs, orderedCheckNames=ORDERED_CHECK_NAMES, diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/checkFracturesParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/checkFracturesParsing.py new file mode 100644 index 00000000..f0f2303d --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/checkFracturesParsing.py @@ -0,0 +1 @@ +# TODO create the parsing for checkFractures action diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/cliParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/cliParsing.py similarity index 76% rename from geos-mesh/src/geos/mesh/doctor/parsing/cliParsing.py rename to mesh-doctor/src/geos/mesh_doctor/parsing/cliParsing.py index d185d981..cef2a2bd 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/cliParsing.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/cliParsing.py @@ -1,3 +1,6 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto import argparse import logging import textwrap @@ -16,12 +19,35 @@ def parseCommaSeparatedString( value: str ) -> list[ str ]: """Helper to parse comma-separated strings, stripping whitespace and removing empty items.""" if not value or not value.strip(): - return list() + return [] return [ item.strip() for item in value.split( ',' ) if item.strip() ] +def addVtuInputFileArgument( parser: argparse.ArgumentParser, required: bool = True ) -> None: + """Add the VTU input file argument to a parser. + + Args: + parser: The argument parser to add the input file argument to. + required: Whether the input file argument is required. Defaults to True. + """ + parser.add_argument( '-i', + '--vtu-input-file', + metavar='VTU_MESH_FILE', + type=str, + required=required, + dest='vtuInputFile', + help="[string]: The VTU mesh file to process." + ( "" if required else " (optional)" ) ) + + def parseAndSetVerbosity( cliArgs: list[ str ] ) -> None: """Parse the verbosity flag only and set the root logger's level accordingly. + + The verbosity is controlled via -v and -q flags: + No flags: WARNING level + -v: INFO level + -vv or more: DEBUG level + -q: ERROR level + -qq or more: CRITICAL level Messages from loggers created with `get_custom_logger` will inherit this level if their own level is set to NOTSET. @@ -72,13 +98,16 @@ def parseAndSetVerbosity( cliArgs: list[ str ] ) -> None: def initParser() -> argparse.ArgumentParser: - vtkInputFileKey = "vtkInputFile" + """Initialize the main argument parser for mesh-doctor.""" epilogMsg = f"""\ Note that checks are dynamically loaded. An option may be missing because of an unloaded module. Increase verbosity (-{__VERBOSITY_FLAG}, -{__VERBOSITY_FLAG * 2}) to get full information. """ - formatter = lambda prog: argparse.RawTextHelpFormatter( prog, max_help_position=8 ) + + def formatter( prog: str ) -> argparse.RawTextHelpFormatter: + return argparse.RawTextHelpFormatter( prog, max_help_position=8 ) + parser = argparse.ArgumentParser( description='Inspects meshes for GEOS.', epilog=textwrap.dedent( epilogMsg ), formatter_class=formatter ) @@ -96,10 +125,4 @@ def initParser() -> argparse.ArgumentParser: default=0, dest=__QUIET_KEY, help=f"Use -{__QUIET_FLAG} to reduce the verbosity of the output." ) - parser.add_argument( '-i', - '--vtk-input-file', - metavar='VTK_MESH_FILE', - type=str, - required=True, - dest=vtkInputFileKey ) return parser diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/collocatedNodesParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/collocatedNodesParsing.py new file mode 100644 index 00000000..7abe6eaf --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/collocatedNodesParsing.py @@ -0,0 +1,89 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations +from argparse import _SubParsersAction +from typing import Any +from geos.mesh_doctor.actions.collocatedNodes import Options, Result +from geos.mesh_doctor.parsing import COLLOCATES_NODES +from geos.mesh_doctor.parsing._sharedChecksParsingLogic import getOptionsUsedMessage +from geos.mesh_doctor.parsing.cliParsing import setupLogger, addVtuInputFileArgument + +__TOLERANCE = "tolerance" +__TOLERANCE_DEFAULT = 0. + +__COLLOCATED_NODES_DEFAULT = { __TOLERANCE: __TOLERANCE_DEFAULT } + + +def convert( parsedOptions: dict[ str, Any ] ) -> Options: + """Convert parsed command-line options to Options object. + + Args: + parsedOptions: Dictionary of parsed command-line options. + + Returns: + Options: Configuration options for supported elements check. + """ + return Options( parsedOptions[ __TOLERANCE ] ) + + +def fillSubparser( subparsers: _SubParsersAction[ Any ] ) -> None: + """Add supported elements check subparser with its arguments. + + Args: + subparsers: The subparsers action to add the parser to. + """ + p = subparsers.add_parser( COLLOCATES_NODES, help="Checks if nodes are collocated." ) + addVtuInputFileArgument( p ) + p.add_argument( '--' + __TOLERANCE, + type=float, + metavar=__TOLERANCE_DEFAULT, + default=__TOLERANCE_DEFAULT, + required=True, + help="[float]: The absolute distance between two nodes for them to be considered collocated." ) + + +def displayResults( options: Options, result: Result ) -> None: + """Display the results of the collocated nodes check. + + Args: + options: The options used for the check. + result: The result of the collocated nodes check. + """ + setupLogger.results( getOptionsUsedMessage( options ) ) + loggerResults( setupLogger, result.nodesBuckets, result.wrongSupportElements ) + + +def loggerResults( logger: Any, nodesBuckets: list[ tuple[ int ] ], wrongSupportElements: list[ int ] ) -> None: + """Log the results of the collocated nodes check. + + Args: + logger: Logger instance for output. + nodesBuckets (list[ tuple[ int ] ]): List of collocated nodes buckets. + wrongSupportElements (list[ int ]): List of elements with wrong support nodes. + """ + # Accounts for external logging object that would not contain 'results' attribute + logMethod = logger.info + if hasattr( logger, 'results' ): + logMethod = logger.results + + allCollocatedNodes: list[ int ] = [] + for bucket in nodesBuckets: + for node in bucket: + allCollocatedNodes.append( node ) + allCollocatedNodesUnique = list( set( allCollocatedNodes ) ) + if allCollocatedNodesUnique: + logMethod( f"You have {len( allCollocatedNodesUnique )} collocated nodes." ) + logMethod( "Here are all the buckets of collocated nodes." ) + tmp: list[ str ] = [] + for bucket in nodesBuckets: + tmp.append( f"({', '.join( map( str, bucket ) )})" ) + logMethod( f"({', '.join( tmp )})" ) + else: + logMethod( "You have no collocated node." ) + + if wrongSupportElements: + wsElements: str = ", ".join( map( str, wrongSupportElements ) ) + logMethod( f"You have {len( wrongSupportElements )} elements with duplicated support nodes.\n" + wsElements ) + else: + logMethod( "You have no element with duplicated support nodes." ) diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/elementVolumesParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/elementVolumesParsing.py new file mode 100644 index 00000000..1ac55c1c --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/elementVolumesParsing.py @@ -0,0 +1,74 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations +from argparse import _SubParsersAction +from typing import Any +from geos.mesh_doctor.actions.elementVolumes import Options, Result +from geos.mesh_doctor.parsing import ELEMENT_VOLUMES +from geos.mesh_doctor.parsing._sharedChecksParsingLogic import getOptionsUsedMessage +from geos.mesh_doctor.parsing.cliParsing import setupLogger, addVtuInputFileArgument + +__MIN_VOLUME = "minVolume" +__MIN_VOLUME_DEFAULT = 0. + +__ELEMENT_VOLUMES_DEFAULT = { __MIN_VOLUME: __MIN_VOLUME_DEFAULT } + + +def fillSubparser( subparsers: _SubParsersAction[ Any ] ) -> None: + """Add supported elements check subparser with its arguments. + + Args: + subparsers: The subparsers action to add the parser to. + """ + p = subparsers.add_parser( ELEMENT_VOLUMES, + help=f"Checks if the volumes of the elements are greater than \"{__MIN_VOLUME}\"." ) + addVtuInputFileArgument( p ) + p.add_argument( '--' + __MIN_VOLUME, + type=float, + metavar=__MIN_VOLUME_DEFAULT, + default=__MIN_VOLUME_DEFAULT, + required=True, + help=f"[float]: The minimum acceptable volume. Defaults to {__MIN_VOLUME_DEFAULT}." ) + + +def convert( parsedOptions: dict[ str, Any ] ) -> Options: + """From the parsed cli options, return the converted options for elements volumes check. + + Args: + parsedOptions: Parsed cli options. + + Returns: + Options: The converted options for elements volumes check. + """ + return Options( minVolume=parsedOptions[ __MIN_VOLUME ] ) + + +def displayResults( options: Options, result: Result ) -> None: + """Display the results of the element volumes check. + + Args: + options: The options used for the check. + result: The result of the element volumes check. + """ + setupLogger.results( getOptionsUsedMessage( options ) ) + loggerResults( setupLogger, result.elementVolumes ) + + +def loggerResults( logger: Any, elementVolumes: list[ tuple[ int, float ] ] ) -> None: + """Show the results of the element volumes check. + + Args: + logger: Logger instance for output. + elementVolumes (list[ tuple[ int, float ] ]): List of element volumes. + """ + # Accounts for external logging object that would not contain 'results' attribute + logMethod = logger.info + if hasattr( logger, 'results' ): + logMethod = logger.results + + logMethod( "Elements index | Volumes calculated" ) + logMethod( "-----------------------------------" ) + max_length: int = len( "Elements index " ) + for ( ind, volume ) in elementVolumes: + logMethod( f"{ind:<{max_length}}" + "| " + str( volume ) ) diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/fixElementsOrderingsParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/fixElementsOrderingsParsing.py similarity index 69% rename from geos-mesh/src/geos/mesh/doctor/parsing/fixElementsOrderingsParsing.py rename to mesh-doctor/src/geos/mesh_doctor/parsing/fixElementsOrderingsParsing.py index 83fbbc34..f40ce53b 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/fixElementsOrderingsParsing.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/fixElementsOrderingsParsing.py @@ -1,3 +1,8 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations +from argparse import _SubParsersAction import random from vtkmodules.vtkCommonDataModel import ( VTK_HEXAGONAL_PRISM, @@ -8,9 +13,10 @@ VTK_VOXEL, VTK_WEDGE, ) -from geos.mesh.doctor.actions.fixElementsOrderings import Options, Result -from geos.mesh.doctor.parsing import vtkOutputParsing, FIX_ELEMENTS_ORDERINGS -from geos.mesh.doctor.parsing.cliParsing import setupLogger +from typing import Any +from geos.mesh_doctor.actions.fixElementsOrderings import Options, Result +from geos.mesh_doctor.parsing import vtkOutputParsing, FIX_ELEMENTS_ORDERINGS +from geos.mesh_doctor.parsing.cliParsing import setupLogger, addVtuInputFileArgument __CELL_TYPE_MAPPING = { "Hexahedron": VTK_HEXAHEDRON, @@ -33,8 +39,14 @@ } -def fillSubparser( subparsers ) -> None: +def fillSubparser( subparsers: _SubParsersAction[ Any ] ) -> None: + """Add supported elements check subparser with its arguments. + + Args: + subparsers: The subparsers action to add the parser to. + """ p = subparsers.add_parser( FIX_ELEMENTS_ORDERINGS, help="Reorders the support nodes for the given cell types." ) + addVtuInputFileArgument( p ) for key, vtkKey in __CELL_TYPE_MAPPING.items(): tmp = list( range( __CELL_TYPE_SUPPORT_SIZE[ vtkKey ] ) ) random.Random( 4 ).shuffle( tmp ) @@ -47,7 +59,7 @@ def fillSubparser( subparsers ) -> None: vtkOutputParsing.fillVtkOutputSubparser( p ) -def convert( parsedOptions ) -> Options: +def convert( parsedOptions: dict[ str, Any ] ) -> Options: """From the parsed cli options, return the converted options for self intersecting elements check. Args: @@ -64,7 +76,7 @@ def convert( parsedOptions ) -> Options: rawMapping = parsedOptions[ key ] if rawMapping: tmp = tuple( map( int, rawMapping.split( "," ) ) ) - if not set( tmp ) == set( range( __CELL_TYPE_SUPPORT_SIZE[ vtkKey ] ) ): + if set( tmp ) != set( range( __CELL_TYPE_SUPPORT_SIZE[ vtkKey ] ) ): errMsg = f"Permutation {rawMapping} for type {key} is not valid." setupLogger.error( errMsg ) raise ValueError( errMsg ) @@ -73,7 +85,13 @@ def convert( parsedOptions ) -> Options: return Options( vtkOutput=vtkOutput, cellTypeToOrdering=cellTypeToOrdering ) -def displayResults( options: Options, result: Result ): +def displayResults( options: Options, result: Result ) -> None: + """Display the results of the fix elements orderings feature. + + Args: + options: The options used for the fix. + result: The result of the fix elements orderings feature. + """ if result.output: setupLogger.info( f"New mesh was written to file '{result.output}'" ) if result.unchangedCellTypes: diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/generateCubeParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/generateCubeParsing.py similarity index 62% rename from geos-mesh/src/geos/mesh/doctor/parsing/generateCubeParsing.py rename to mesh-doctor/src/geos/mesh_doctor/parsing/generateCubeParsing.py index f0bc8a38..f4b1d8ad 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/generateCubeParsing.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/generateCubeParsing.py @@ -1,15 +1,29 @@ -from geos.mesh.doctor.actions.generateCube import Options, Result, FieldInfo -from geos.mesh.doctor.parsing import vtkOutputParsing, generateGlobalIdsParsing, GENERATE_CUBE -from geos.mesh.doctor.parsing.cliParsing import setupLogger -from geos.mesh.doctor.parsing.generateGlobalIdsParsing import GlobalIdsInfo +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations +from argparse import _SubParsersAction +from typing import Any +from geos.mesh_doctor.actions.generateCube import Options, Result, FieldInfo +from geos.mesh_doctor.parsing import vtkOutputParsing, generateGlobalIdsParsing, GENERATE_CUBE +from geos.mesh_doctor.parsing.cliParsing import setupLogger, addVtuInputFileArgument +from geos.mesh_doctor.parsing.generateGlobalIdsParsing import GlobalIdsInfo, convertToGlobalIdsInfo __X, __Y, __Z, __NX, __NY, __NZ = "x", "y", "z", "nx", "ny", "nz" __FIELDS = "fields" -def convert( parsedOptions ) -> Options: +def convert( parsedOptions: dict[ str, Any ] ) -> Options: + """Convert parsed command-line options to Options object. - def checkDiscretizations( x, nx, title ): + Args: + parsedOptions: Dictionary of parsed command-line options. + + Returns: + Options: Configuration options for supported elements check. + """ + + def checkDiscretizations( x: tuple[ float, ...], nx: tuple[ int, ...], title: str ) -> None: if len( x ) != len( nx ) + 1: raise ValueError( f"{title} information (\"{x}\" and \"{nx}\") does not have consistent size." ) @@ -17,20 +31,20 @@ def checkDiscretizations( x, nx, title ): checkDiscretizations( parsedOptions[ __Y ], parsedOptions[ __NY ], __Y ) checkDiscretizations( parsedOptions[ __Z ], parsedOptions[ __NZ ], __Z ) - def parseFields( s ): + def parseFields( s: str ) -> FieldInfo: name, support, dim = s.split( ":" ) if support not in ( "CELLS", "POINTS" ): raise ValueError( f"Support {support} for field \"{name}\" must be one of \"CELLS\" or \"POINTS\"." ) try: - dim = int( dim ) - assert dim > 0 - except ValueError: - raise ValueError( f"Dimension {dim} cannot be converted to an integer." ) - except AssertionError: - raise ValueError( f"Dimension {dim} must be a positive integer" ) - return FieldInfo( name=name, support=support, dimension=dim ) + dimension = int( dim ) + assert dimension > 0 + except ValueError as e: + raise ValueError( f"Dimension {dimension} cannot be converted to an integer." ) from e + except AssertionError as e: + raise ValueError( f"Dimension {dimension} must be a positive integer" ) from e + return FieldInfo( name=name, support=support, dimension=dimension ) - gids: GlobalIdsInfo = generateGlobalIdsParsing.convertGlobalIds( parsedOptions ) + gids: GlobalIdsInfo = convertToGlobalIdsInfo( parsedOptions ) return Options( vtkOutput=vtkOutputParsing.convert( parsedOptions ), generateCellsGlobalIds=gids.cells, @@ -44,8 +58,14 @@ def parseFields( s ): fields=tuple( map( parseFields, parsedOptions[ __FIELDS ] ) ) ) -def fillSubparser( subparsers ) -> None: +def fillSubparser( subparsers: _SubParsersAction[ Any ] ) -> None: + """Add supported elements check subparser with its arguments. + + Args: + subparsers: The subparsers action to add the parser to. + """ p = subparsers.add_parser( GENERATE_CUBE, help="Generate a cube and its fields." ) + addVtuInputFileArgument( p, required=False ) p.add_argument( '--' + __X, type=lambda s: tuple( map( float, s.split( ":" ) ) ), metavar="0:1.5:3", @@ -77,9 +97,15 @@ def fillSubparser( subparsers ) -> None: required=False, default=(), help="Create fields on CELLS or POINTS, with given dimension (typically 1 or 3)." ) - generateGlobalIdsParsing.fillGenerateGlobalIdsSubparser( p ) + generateGlobalIdsParsing.addArguments( p ) vtkOutputParsing.fillVtkOutputSubparser( p ) -def displayResults( options: Options, result: Result ): +def displayResults( options: Options, result: Result ) -> None: + """Display the results of the generate cube feature. + + Args: + options: The options used for the check. + result: The result of the generate cube feature. + """ setupLogger.info( result.info ) diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/generateFracturesParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/generateFracturesParsing.py similarity index 66% rename from geos-mesh/src/geos/mesh/doctor/parsing/generateFracturesParsing.py rename to mesh-doctor/src/geos/mesh_doctor/parsing/generateFracturesParsing.py index ffc3f17e..15d05db7 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/generateFracturesParsing.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/generateFracturesParsing.py @@ -1,6 +1,13 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations +from argparse import _SubParsersAction import os -from geos.mesh.doctor.actions.generateFractures import Options, Result, FracturePolicy -from geos.mesh.doctor.parsing import vtkOutputParsing, GENERATE_FRACTURES +from typing import Any +from geos.mesh_doctor.actions.generateFractures import Options, Result, FracturePolicy +from geos.mesh_doctor.parsing import vtkOutputParsing, GENERATE_FRACTURES +from geos.mesh_doctor.parsing.cliParsing import addVtuInputFileArgument from geos.mesh.io.vtkIO import VtkOutput __POLICY = "policy" @@ -19,6 +26,7 @@ def convertToFracturePolicy( s: str ) -> FracturePolicy: """Converts the user input to the proper enum chosen. + I do not want to use the auto conversion already available to force explicit conversion. Args: @@ -37,8 +45,14 @@ def convertToFracturePolicy( s: str ) -> FracturePolicy: raise ValueError( f"Policy {s} is not valid. Please use one of \"{', '.join(map(str, __POLICIES))}\"." ) -def fillSubparser( subparsers ) -> None: +def fillSubparser( subparsers: _SubParsersAction[ Any ] ) -> None: + """Add supported elements check subparser with its arguments. + + Args: + subparsers: The subparsers action to add the parser to. + """ p = subparsers.add_parser( GENERATE_FRACTURES, help="Splits the mesh to generate the faults and fractures." ) + addVtuInputFileArgument( p ) p.add_argument( '--' + __POLICY, type=convertToFracturePolicy, metavar=", ".join( __POLICIES ), @@ -58,23 +72,31 @@ def fillSubparser( subparsers ) -> None: help= 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. " + "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." ) + "If no ':' is found, all values specified will be assumed to create only 1 single fracture." ) vtkOutputParsing.fillVtkOutputSubparser( p ) 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." ) + help="[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.' ) + help='[string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary.' ) -def convert( parsedOptions ) -> Options: +def convert( parsedOptions: dict[ str, Any ] ) -> Options: + """Convert parsed command-line options to Options object. + + Args: + parsedOptions: Dictionary of parsed command-line options. + + Returns: + Options: Configuration options for supported elements check. + """ policy: str = parsedOptions[ __POLICY ] field: str = parsedOptions[ __FIELD_NAME ] allValues: str = parsedOptions[ __FIELD_VALUES ] @@ -93,8 +115,9 @@ def convert( parsedOptions ) -> Options: ] fractureNames: list[ str ] = [ "fracture_" + frac.replace( ",", "_" ) + ".vtu" for frac in perFracture ] fracturesOutputDir: str = parsedOptions[ __FRACTURES_OUTPUT_DIR ] - fracturesDataMode: str = parsedOptions[ __FRACTURES_DATA_MODE ] == __FRACTURES_DATA_MODE_DEFAULT - allFracturesVtkOutput: list[ VtkOutput ] = buildAllFracturesVtkOutput( fracturesOutputDir, fracturesDataMode, + fracturesDataMode: str = parsedOptions[ __FRACTURES_DATA_MODE ] + fracturesDataModeResult: bool = fracturesDataMode == __FRACTURES_DATA_MODE_DEFAULT + allFracturesVtkOutput: list[ VtkOutput ] = buildAllFracturesVtkOutput( fracturesOutputDir, fracturesDataModeResult, meshVtkOutput, fractureNames ) return Options( policy=policy, field=field, @@ -104,22 +127,49 @@ def convert( parsedOptions ) -> Options: allFracturesVtkOutput=allFracturesVtkOutput ) -def displayResults( options: Options, result: Result ): +def displayResults( options: Options, result: Result ) -> None: + """Display the results of the generate fractures feature. + + Args: + options: The options used for the check. + result: The result of the generate fractures feature. + """ pass def areValuesParsable( values: str ) -> bool: + """Check if a string contains parsable values. + + Args: + values (str): The string containing values to be checked. + + Returns: + bool: True if the string contains parsable values, False otherwise. + """ 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 + return not ( values.endswith( ":" ) or values.endswith( "," ) ) def buildAllFracturesVtkOutput( fractureOutputDir: str, fracturesDataMode: bool, meshVtkOutput: VtkOutput, fractureNames: list[ str ] ) -> list[ VtkOutput ]: + """Create all the VtkOutput objects for every fracture mesh that will be output. + + Args: + fractureOutputDir (str): The directory where fracture output files will be saved. + fracturesDataMode (bool): Whether the fractures data mode is enabled. + meshVtkOutput (VtkOutput): The VtkOutput object for the mesh. + fractureNames (list[ str ]): The list of fracture file names. + + Raises: + ValueError: If the given fracture output directory does not exist. + ValueError: If the given fracture output directory is not writable. + + Returns: + list[ VtkOutput ]: All the VtkOutput objects for every fracture mesh that will be output. + """ if not os.path.exists( fractureOutputDir ): raise ValueError( f"The --{__FRACTURES_OUTPUT_DIR} given directory '{fractureOutputDir}' does not exist." ) @@ -129,7 +179,7 @@ def buildAllFracturesVtkOutput( fractureOutputDir: str, fracturesDataMode: bool, outputName = os.path.basename( meshVtkOutput.output ) splittedNameWithoutExtension: list[ str ] = outputName.split( "." )[ :-1 ] nameWithoutExtension: str = '_'.join( splittedNameWithoutExtension ) + "_" - allFracturesVtkOutput: list[ VtkOutput ] = list() + allFracturesVtkOutput: list[ VtkOutput ] = [] for fractureName in fractureNames: fracturePath = os.path.join( fractureOutputDir, nameWithoutExtension + fractureName ) allFracturesVtkOutput.append( VtkOutput( fracturePath, fracturesDataMode ) ) diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/generateGlobalIdsParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/generateGlobalIdsParsing.py new file mode 100644 index 00000000..0bd8063c --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/generateGlobalIdsParsing.py @@ -0,0 +1,91 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations +from argparse import _SubParsersAction, ArgumentParser, _ArgumentGroup +from dataclasses import dataclass +from typing import Any +from geos.mesh_doctor.actions.generateGlobalIds import Options, Result +from geos.mesh_doctor.parsing import vtkOutputParsing, GENERATE_GLOBAL_IDS +from geos.mesh_doctor.parsing.cliParsing import setupLogger, addVtuInputFileArgument + +__CELLS, __POINTS = "cells", "points" + + +@dataclass( frozen=True ) +class GlobalIdsInfo: + cells: bool + points: bool + + +def convertToGlobalIdsInfo( parsedOptions: dict[ str, Any ] ) -> GlobalIdsInfo: + """Extract GlobalIdsInfo from parsed options. + + Args: + parsedOptions: Dictionary of parsed command-line options. + + Returns: + GlobalIdsInfo: The global IDs configuration. + """ + return GlobalIdsInfo( cells=parsedOptions[ __CELLS ], points=parsedOptions[ __POINTS ] ) + + +def convert( parsedOptions: dict[ str, Any ] ) -> Options: + """Convert parsed command-line options to Options object. + + Args: + parsedOptions: Dictionary of parsed command-line options. + + Returns: + Options: Configuration options for supported elements check. + """ + gids: GlobalIdsInfo = convertToGlobalIdsInfo( parsedOptions ) + return Options( vtkOutput=vtkOutputParsing.convert( parsedOptions ), + generateCellsGlobalIds=gids.cells, + generatePointsGlobalIds=gids.points ) + + +def addArguments( parser: ArgumentParser | _ArgumentGroup ) -> None: + """Add global IDs arguments to an existing parser. + + Args: + parser: The argument parser or argument group to add arguments to. + """ + parser.add_argument( '--' + __CELLS, + action="store_true", + help="[bool]: Generate global ids for cells. Defaults to true." ) + parser.add_argument( '--no-' + __CELLS, + action="store_false", + dest=__CELLS, + help="[bool]: Don't generate global ids for cells." ) + parser.set_defaults( **{ __CELLS: True } ) + parser.add_argument( '--' + __POINTS, + action="store_true", + help="[bool]: Generate global ids for points. Defaults to true." ) + parser.add_argument( '--no-' + __POINTS, + action="store_false", + dest=__POINTS, + help="[bool]: Don't generate global ids for points." ) + parser.set_defaults( **{ __POINTS: True } ) + + +def fillSubparser( subparsers: _SubParsersAction[ Any ] ) -> None: + """Add supported elements check subparser with its arguments. + + Args: + subparsers: The subparsers action to add the parser to. + """ + p = subparsers.add_parser( GENERATE_GLOBAL_IDS, help="Adds globals ids for points and cells." ) + addVtuInputFileArgument( p ) + addArguments( p ) + vtkOutputParsing.fillVtkOutputSubparser( p ) + + +def displayResults( options: Options, result: Result ) -> None: + """Display the results of the generate global ids feature. + + Args: + options: The options used for the check. + result: The result of the generate global ids feature. + """ + setupLogger.info( result.info ) diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/mainChecksParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/mainChecksParsing.py similarity index 76% rename from geos-mesh/src/geos/mesh/doctor/parsing/mainChecksParsing.py rename to mesh-doctor/src/geos/mesh_doctor/parsing/mainChecksParsing.py index d7d124ff..f45dac16 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/mainChecksParsing.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/mainChecksParsing.py @@ -1,17 +1,22 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations import argparse from copy import deepcopy -from geos.mesh.doctor.actions.allChecks import Options as AllChecksOptions -from geos.mesh.doctor.parsing._sharedChecksParsingLogic import ( CheckFeature, convert as sharedConvert, fillSubparser - as sharedFillSubparser, displayResults ) -from geos.mesh.doctor.parsing import ( +from typing import Any +from geos.mesh_doctor.actions.allChecks import Options as AllChecksOptions +from geos.mesh_doctor.parsing._sharedChecksParsingLogic import ( CheckFeature, convert as sharedConvert, fillSubparser + as sharedFillSubparser, displayResults ) # noqa: F401 +from geos.mesh_doctor.parsing import ( MAIN_CHECKS, COLLOCATES_NODES, ELEMENT_VOLUMES, SELF_INTERSECTING_ELEMENTS, ) -from geos.mesh.doctor.parsing import collocatedNodesParsing as cnParser -from geos.mesh.doctor.parsing import elementVolumesParsing as evParser -from geos.mesh.doctor.parsing import selfIntersectingElementsParsing as sieParser +from geos.mesh_doctor.parsing import collocatedNodesParsing as cnParser +from geos.mesh_doctor.parsing import elementVolumesParsing as evParser +from geos.mesh_doctor.parsing import selfIntersectingElementsParsing as sieParser # Ordered list of check names for this configuration ORDERED_CHECK_NAMES = [ @@ -52,7 +57,7 @@ def fillSubparser( subparsers: argparse._SubParsersAction ) -> None: checkFeaturesConfig=CHECK_FEATURES_CONFIG ) -def convert( parsedArgs: argparse.Namespace ) -> AllChecksOptions: +def convert( parsedArgs: dict[ str, Any ] ) -> AllChecksOptions: """Converts arguments by calling the shared logic with the 'main_checks' configuration.""" return sharedConvert( parsedArgs=parsedArgs, orderedCheckNames=ORDERED_CHECK_NAMES, diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/nonConformalParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/nonConformalParsing.py new file mode 100644 index 00000000..45f1facd --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/nonConformalParsing.py @@ -0,0 +1,96 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations +from argparse import _SubParsersAction +from typing import Any +from geos.mesh_doctor.actions.nonConformal import Options, Result +from geos.mesh_doctor.parsing import NON_CONFORMAL +from geos.mesh_doctor.parsing._sharedChecksParsingLogic import getOptionsUsedMessage +from geos.mesh_doctor.parsing.cliParsing import setupLogger, addVtuInputFileArgument + +__ANGLE_TOLERANCE = "angleTolerance" +__POINT_TOLERANCE = "pointTolerance" +__FACE_TOLERANCE = "faceTolerance" + +__ANGLE_TOLERANCE_DEFAULT = 10. +__POINT_TOLERANCE_DEFAULT = 0. +__FACE_TOLERANCE_DEFAULT = 0. + +__NON_CONFORMAL_DEFAULT = { + __ANGLE_TOLERANCE: __ANGLE_TOLERANCE_DEFAULT, + __POINT_TOLERANCE: __POINT_TOLERANCE_DEFAULT, + __FACE_TOLERANCE: __FACE_TOLERANCE_DEFAULT +} + + +def convert( parsedOptions: dict[ str, Any ] ) -> Options: + """Convert parsed command-line options to Options object. + + Args: + parsedOptions: Dictionary of parsed command-line options. + + Returns: + Options: Configuration options for supported elements check. + """ + return Options( angleTolerance=parsedOptions[ __ANGLE_TOLERANCE ], + pointTolerance=parsedOptions[ __POINT_TOLERANCE ], + faceTolerance=parsedOptions[ __FACE_TOLERANCE ] ) + + +def fillSubparser( subparsers: _SubParsersAction[ Any ] ) -> None: + """Add supported elements check subparser with its arguments. + + Args: + subparsers: The subparsers action to add the parser to. + """ + p = subparsers.add_parser( NON_CONFORMAL, help="Detects non conformal elements. [EXPERIMENTAL]" ) + addVtuInputFileArgument( p ) + p.add_argument( '--' + __ANGLE_TOLERANCE, + type=float, + metavar=__ANGLE_TOLERANCE_DEFAULT, + default=__ANGLE_TOLERANCE_DEFAULT, + help=f"[float]: angle tolerance in degrees. Defaults to {__ANGLE_TOLERANCE_DEFAULT}" ) + p.add_argument( + '--' + __POINT_TOLERANCE, + type=float, + metavar=__POINT_TOLERANCE_DEFAULT, + default=__POINT_TOLERANCE_DEFAULT, + help=f"[float]: tolerance for two points to be considered collocated. Defaults to {__POINT_TOLERANCE_DEFAULT}" ) + p.add_argument( + '--' + __FACE_TOLERANCE, + type=float, + metavar=__FACE_TOLERANCE_DEFAULT, + default=__FACE_TOLERANCE_DEFAULT, + help=f"[float]: tolerance for two faces to be considered \"touching\". Defaults to {__FACE_TOLERANCE_DEFAULT}" ) + + +def displayResults( options: Options, result: Result ) -> None: + """Display the results of the non conformal elements check. + + Args: + options: The options used for the check. + result: The result of the non conformal elements check. + """ + setupLogger.results( getOptionsUsedMessage( options ) ) + loggerResults( setupLogger, result.nonConformalCells ) + + +def loggerResults( logger: Any, nonConformalCells: list[ tuple[ int, int ] ] ) -> None: + """Log the results of the non-conformal cells check. + + Args: + logger: Logger instance for output. + nonConformalCells (list[ tuple[ int, int ] ]): List of non-conformal cells. + """ + # Accounts for external logging object that would not contain 'results' attribute + logMethod = logger.info + if hasattr( logger, 'results' ): + logMethod = logger.results + + cells: list[ int ] = [] + for i, j in nonConformalCells: + cells += i, j + uniqueCells: frozenset[ int ] = frozenset( cells ) + logMethod( f"You have {len( uniqueCells )} non conformal cells." ) + logMethod( f"{', '.join( map( str, sorted( uniqueCells ) ) )}" ) diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/selfIntersectingElementsParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/selfIntersectingElementsParsing.py new file mode 100644 index 00000000..5838394c --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/selfIntersectingElementsParsing.py @@ -0,0 +1,100 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations +from argparse import _SubParsersAction +import numpy +from typing import Any +from geos.mesh_doctor.actions.selfIntersectingElements import Options, Result +from geos.mesh_doctor.parsing import SELF_INTERSECTING_ELEMENTS +from geos.mesh_doctor.parsing._sharedChecksParsingLogic import getOptionsUsedMessage +from geos.mesh_doctor.parsing.cliParsing import setupLogger, addVtuInputFileArgument + +__MIN_DISTANCE = "minDistance" +__MIN_DISTANCE_DEFAULT = numpy.finfo( float ).eps + +__SELF_INTERSECTING_ELEMENTS_DEFAULT = { __MIN_DISTANCE: __MIN_DISTANCE_DEFAULT } + + +def convert( parsedOptions: dict[ str, Any ] ) -> Options: + """Convert parsed command-line options to Options object. + + Args: + parsedOptions: Dictionary of parsed command-line options. + + Returns: + Options: Configuration options for self intersecting elements check. + + Raises: + ValueError: If minimum distance is negative. + """ + minDistance = parsedOptions[ __MIN_DISTANCE ] + if minDistance == 0: + setupLogger.warning( "Having minimum distance set to 0 can induce lots of false positive results " + "(adjacent faces may be considered intersecting)." ) + elif minDistance < 0: + raise ValueError( + f"Negative minimum distance ({minDistance}) in the {SELF_INTERSECTING_ELEMENTS} check is not allowed." ) + return Options( minDistance=minDistance ) + + +def fillSubparser( subparsers: _SubParsersAction[ Any ] ) -> None: + """Add self intersecting elements check subparser with its arguments. + + Args: + subparsers: The subparsers action to add the parser to. + """ + p = subparsers.add_parser( SELF_INTERSECTING_ELEMENTS, + help="Checks if the faces of the elements are self intersecting." ) + addVtuInputFileArgument( p ) + help_text = ( "[float]: The minimum distance in the computation. " + f"Defaults to your machine precision {__MIN_DISTANCE_DEFAULT}." ) + p.add_argument( '--' + __MIN_DISTANCE, + type=float, + required=False, + metavar=__MIN_DISTANCE_DEFAULT, + default=__MIN_DISTANCE_DEFAULT, + help=help_text ) + + +def displayResults( options: Options, result: Result ) -> None: + """Display the results of the self intersecting elements check. + + Args: + options: The options used for the check. + result: The result of the self intersecting elements check. + """ + setupLogger.results( getOptionsUsedMessage( options ) ) + loggerResults( setupLogger, result.invalidCellIds ) + + +def loggerResults( logger: Any, invalidCellIds: dict[ str, list[ int ] ] ) -> None: + """Log the results of the self-intersecting elements check. + + Args: + logger: Logger instance for output. + invalidCellIds: Dictionary of invalid cell IDs by error type. + """ + # Accounts for external logging object that would not contain 'results' attribute + logMethod: Any = logger.info + if hasattr( logger, 'results' ): + logMethod = logger.results + + # Human-readable descriptions for each error type + error_descriptions = { + 'wrongNumberOfPointsElements': 'elements with wrong number of points', + 'intersectingEdgesElements': 'elements with intersecting edges', + 'intersectingFacesElements': 'elements with self intersecting faces', + 'nonContiguousEdgesElements': 'elements with non-contiguous edges', + 'nonConvexElements': 'non-convex elements', + 'facesOrientedIncorrectlyElements': 'elements with incorrectly oriented faces', + 'nonPlanarFacesElements': 'elements with non-planar faces', + 'degenerateFacesElements': 'elements with degenerate faces' + } + + # Log results for each error type that has invalid elements + for errorType, invalidIds in invalidCellIds.items(): + if invalidIds: + description = error_descriptions.get( errorType, f'elements with {errorType}' ) + logMethod( f"You have {len(invalidIds)} {description}." ) + logMethod( "The elements indices are:\n" + ", ".join( map( str, invalidIds ) ) ) diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/supportedElementsParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/supportedElementsParsing.py new file mode 100644 index 00000000..89f57e93 --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/supportedElementsParsing.py @@ -0,0 +1,94 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations +import multiprocessing +from argparse import _SubParsersAction +from typing import Any +from geos.mesh_doctor.actions.supportedElements import Options, Result +from geos.mesh_doctor.parsing import SUPPORTED_ELEMENTS +from geos.mesh_doctor.parsing._sharedChecksParsingLogic import getOptionsUsedMessage +from geos.mesh_doctor.parsing.cliParsing import setupLogger, addVtuInputFileArgument + +__CHUNK_SIZE = "chunkSize" +__NUM_PROC = "nproc" + +__CHUNK_SIZE_DEFAULT = 1 +__NUM_PROC_DEFAULT = multiprocessing.cpu_count() + +__SUPPORTED_ELEMENTS_DEFAULT = { __CHUNK_SIZE: __CHUNK_SIZE_DEFAULT, __NUM_PROC: __NUM_PROC_DEFAULT } + + +def convert( parsedOptions: dict[ str, Any ] ) -> Options: + """Convert parsed command-line options to Options object. + + Args: + parsedOptions: Dictionary of parsed command-line options. + + Returns: + Options: Configuration options for supported elements check. + """ + return Options( chunkSize=parsedOptions[ __CHUNK_SIZE ], nproc=parsedOptions[ __NUM_PROC ] ) + + +def fillSubparser( subparsers: _SubParsersAction[ Any ] ) -> None: + """Add supported elements check subparser with its arguments. + + Args: + subparsers: The subparsers action to add the parser to. + """ + p = subparsers.add_parser( SUPPORTED_ELEMENTS, + help="Check that all the elements of the mesh are supported by GEOSX." ) + addVtuInputFileArgument( p ) + p.add_argument( '--' + __CHUNK_SIZE, + type=int, + required=False, + metavar=__CHUNK_SIZE_DEFAULT, + default=__CHUNK_SIZE_DEFAULT, + help=f"[int]: Defaults chunk size for parallel processing to {__CHUNK_SIZE_DEFAULT}" ) + p.add_argument( + '--' + __NUM_PROC, + type=int, + required=False, + metavar=__NUM_PROC_DEFAULT, + default=__NUM_PROC_DEFAULT, + help=f"[int]: Number of threads used for parallel processing. Defaults to your CPU count {__NUM_PROC_DEFAULT}." + ) + + +def displayResults( options: Options, result: Result ) -> None: + """Display the results of the supported elements check. + + Args: + options: The options used for the check. + result: The result of the supported elements check. + """ + setupLogger.results( getOptionsUsedMessage( options ) ) + loggerResults( setupLogger, result.unsupportedPolyhedronElements, result.unsupportedStdElementsTypes ) + + +def loggerResults( logger: Any, unsupportedPolyhedronElements: frozenset[ int ], + unsupportedStdElementsTypes: list[ str ] ) -> None: + """Log the results of the supported elements check. + + Args: + logger: Logger instance for output. + unsupportedPolyhedronElements (frozenset[ int ]): List of unsupported polyhedron elements. + unsupportedStdElementsTypes (list[ str ]): List of unsupported standard element types. + """ + # Accounts for external logging object that would not contain 'results' attribute + logMethod: Any = logger.info + if hasattr( logger, 'results' ): + logMethod = logger.results + + if unsupportedPolyhedronElements: + logMethod( f"There is/are {len(unsupportedPolyhedronElements)} polyhedra that may not be converted to" + " supported elements." ) + logMethod( f"The list of the unsupported polyhedra is\n{tuple( sorted( unsupportedPolyhedronElements ) )}." ) + else: + logMethod( "All the polyhedra (if any) can be converted to supported elements." ) + if unsupportedStdElementsTypes: + logMethod( "There are unsupported vtk standard element types. The list of those vtk types is" + f" {tuple( sorted( unsupportedStdElementsTypes ) )}." ) + else: + logMethod( "All the standard vtk element types (if any) are supported." ) diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/vtkOutputParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/vtkOutputParsing.py new file mode 100644 index 00000000..b382e8e8 --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/vtkOutputParsing.py @@ -0,0 +1,76 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations +import textwrap +from argparse import ArgumentParser, _ArgumentGroup +from typing import Any +from geos.mesh.io.vtkIO import VtkOutput + +__OUTPUT_FILE = "output" +__OUTPUT_BINARY_MODE = "data-mode" +__OUTPUT_BINARY_MODE_VALUES = "binary", "ascii" +__OUTPUT_BINARY_MODE_DEFAULT = __OUTPUT_BINARY_MODE_VALUES[ 0 ] + + +def getVtkOutputHelp() -> str: + """Get help text for VTK output options. + + Returns: + str: Formatted help text describing output file and data mode options. + """ + msg = ( f"{__OUTPUT_FILE} [string]: The vtk output file destination.\n" + f" {__OUTPUT_BINARY_MODE} [string]: For \".vtu\" output format, " + f"the data mode can be {' or '.join(__OUTPUT_BINARY_MODE_VALUES)}. " + f"Defaults to {__OUTPUT_BINARY_MODE_DEFAULT}." ) + return textwrap.dedent( msg ) + + +def __buildArg( prefix: str, main: str ) -> str: + """Build an argument name by joining prefix and main with a hyphen. + + Args: + prefix: The prefix string (can be empty). + main: The main argument name. + + Returns: + str: The joined argument name. + """ + return "-".join( filter( None, ( prefix, main ) ) ) + + +def fillVtkOutputSubparser( parser: ArgumentParser | _ArgumentGroup, prefix: str = "" ) -> None: + """Add VTK output arguments to an argument parser. + + Args: + parser: The argument parser or argument group to add arguments to. + prefix: Optional prefix for argument names. + """ + parser.add_argument( '--' + __buildArg( prefix, __OUTPUT_FILE ), + type=str, + required=True, + help="[string]: The vtk output file destination." ) + help_text = ( f"[string]: For \".vtu\" output format, the data mode can be " + f"{' or '.join(__OUTPUT_BINARY_MODE_VALUES)}. Defaults to {__OUTPUT_BINARY_MODE_DEFAULT}." ) + parser.add_argument( '--' + __buildArg( prefix, __OUTPUT_BINARY_MODE ), + type=str, + metavar=", ".join( __OUTPUT_BINARY_MODE_VALUES ), + default=__OUTPUT_BINARY_MODE_DEFAULT, + help=help_text ) + + +def convert( parsedOptions: dict[ str, Any ], prefix: str = "" ) -> VtkOutput: + """Convert parsed command-line options to a VtkOutput object. + + Args: + parsedOptions: Dictionary of parsed command-line options. + prefix: Optional prefix used when parsing arguments. + + Returns: + VtkOutput: Configured VTK output object. + """ + outputKey = __buildArg( prefix, __OUTPUT_FILE ).replace( "-", "_" ) + binaryModeKey = __buildArg( prefix, __OUTPUT_BINARY_MODE ).replace( "-", "_" ) + output = parsedOptions[ outputKey ] + isDataModeBinary: bool = parsedOptions[ binaryModeKey ] == __OUTPUT_BINARY_MODE_DEFAULT + return VtkOutput( output=output, isDataModeBinary=isDataModeBinary ) diff --git a/geos-mesh/src/geos/mesh/doctor/register.py b/mesh-doctor/src/geos/mesh_doctor/register.py similarity index 58% rename from geos-mesh/src/geos/mesh/doctor/register.py rename to mesh-doctor/src/geos/mesh_doctor/register.py index e265c0dd..bb677f4e 100644 --- a/geos-mesh/src/geos/mesh/doctor/register.py +++ b/mesh-doctor/src/geos/mesh_doctor/register.py @@ -1,21 +1,24 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto import argparse import importlib from typing import Callable, Any -import geos.mesh.doctor.parsing as parsing -from geos.mesh.doctor.parsing import ActionHelper, cliParsing -from geos.mesh.doctor.parsing.cliParsing import setupLogger +import geos.mesh_doctor.parsing as parsing +from geos.mesh_doctor.parsing import ActionHelper, cliParsing +from geos.mesh_doctor.parsing.cliParsing import setupLogger -__HELPERS: dict[ str, Callable[ [ None ], ActionHelper ] ] = dict() -__ACTIONS: dict[ str, Callable[ [ None ], Any ] ] = dict() +__HELPERS: dict[ str, str ] = {} +__ACTIONS: dict[ str, str ] = {} -def __loadModuleAction( moduleName: str, actionFct="action" ): - module = importlib.import_module( "geos.mesh.doctor.actions." + moduleName ) +def __loadModuleAction( moduleName: str, actionFct: str = "action" ) -> Callable[ [ str, Any ], Any ]: + module = importlib.import_module( "geos.mesh_doctor.actions." + moduleName ) return getattr( module, actionFct ) -def __loadModuleActionHelper( moduleName: str, parsingFctSuffix="Parsing" ): - module = importlib.import_module( "geos.mesh.doctor.parsing." + moduleName + parsingFctSuffix ) +def __loadModuleActionHelper( moduleName: str, parsingFctSuffix: str = "Parsing" ) -> ActionHelper: + module = importlib.import_module( "geos.mesh_doctor.parsing." + moduleName + parsingFctSuffix ) return ActionHelper( fillSubparser=module.fillSubparser, convert=module.convert, displayResults=module.displayResults ) @@ -23,16 +26,17 @@ def __loadModuleActionHelper( moduleName: str, parsingFctSuffix="Parsing" ): def __loadActions() -> dict[ str, Callable[ [ str, Any ], Any ] ]: """Loads all the actions. + This function acts like a protection layer if a module fails to load. A action that fails to load won't stop the process. Returns: dict[ str, Callable[ [ str, Any ], Any ] ]: The actions. """ - loadedActions: dict[ str, Callable[ [ str, Any ], Any ] ] = dict() - for actionName, actionProvider in __ACTIONS.items(): + loadedActions: dict[ str, Callable[ [ str, Any ], Any ] ] = {} + for actionName, moduleName in __ACTIONS.items(): try: - loadedActions[ actionName ] = actionProvider() + loadedActions[ actionName ] = __loadModuleAction( moduleName ) setupLogger.debug( f"Action \"{actionName}\" is loaded." ) except Exception as e: setupLogger.warning( f"Could not load module \"{actionName}\": {e}" ) @@ -50,21 +54,23 @@ def registerParsingActions( parser = cliParsing.initParser() subparsers = parser.add_subparsers( help="Modules", dest="subparsers" ) - def closureTrick( cn: str ): - __HELPERS[ actionName ] = lambda: __loadModuleActionHelper( cn ) - __ACTIONS[ actionName ] = lambda: __loadModuleAction( cn ) - # Register the modules to load here. for actionName in ( parsing.ALL_CHECKS, parsing.COLLOCATES_NODES, parsing.ELEMENT_VOLUMES, parsing.FIX_ELEMENTS_ORDERINGS, parsing.GENERATE_CUBE, parsing.GENERATE_FRACTURES, parsing.GENERATE_GLOBAL_IDS, parsing.MAIN_CHECKS, parsing.NON_CONFORMAL, parsing.SELF_INTERSECTING_ELEMENTS, parsing.SUPPORTED_ELEMENTS ): - closureTrick( actionName ) + __HELPERS[ actionName ] = actionName + __ACTIONS[ actionName ] = actionName + loadedActions: dict[ str, Callable[ [ str, Any ], Any ] ] = __loadActions() - loadedActionsHelpers: dict[ str, ActionHelper ] = dict() - for actionName in loadedActions.keys(): - h = __HELPERS[ actionName ]() - h.fillSubparser( subparsers ) - loadedActionsHelpers[ actionName ] = h - setupLogger.debug( f"Parsing for action \"{actionName}\" is loaded." ) + loadedActionsHelpers: dict[ str, ActionHelper ] = {} + for actionName in loadedActions: + moduleName = __HELPERS[ actionName ] + try: + h = __loadModuleActionHelper( moduleName ) + h.fillSubparser( subparsers ) + loadedActionsHelpers[ actionName ] = h + setupLogger.debug( f"Parsing for action \"{actionName}\" is loaded." ) + except Exception as e: + setupLogger.warning( f"Could not load parsing for action \"{actionName}\": {e}" ) return parser, loadedActions, loadedActionsHelpers diff --git a/mesh-doctor/tests/data/supportedElements.vtu b/mesh-doctor/tests/data/supportedElements.vtu new file mode 100644 index 00000000..830e0c81 --- /dev/null +++ b/mesh-doctor/tests/data/supportedElements.vtu @@ -0,0 +1,47 @@ + + + + + + + + + + + AQAAAACAAADAAwAAHwEAAA==eJxlk7FuwjAURT2xlS9g6NCFL4AB6lepW8f+AFs3ZtgshjYDQ8REhwp/CEMkFsbsXSLxBZUqZmI79/kqWIp04tx3/JzExvBwlm6EkNgRLyTXeJqveswZXeOlvSIf//5lelxHfh1eZDt7U64f3p+Rua7Os/v+og/rWeq15Ueh50J7IdZ8j01XKz0WcoIdOT2xI6cnRj7UemLkMQ+uyNkQ1+SsyNMQ1+TRtQ31bagPQ3sIo50vae9lnF9dxvI1OUQuRgPZfXwrL/efmlmf4Uets7nW2VybONWmTKpFX2Xn2cv896fjgthRppDT00bzmeEJTngCF8SOMvCkfOb4P3fOhji8Z08MTxg4D/otDX0nk10YODO6ns15OBPfAMsrmQg= + + + 0 + + + 8.2006097334 + + + + + + + AQAAAACAAACAAgAAjgAAAA==eJwtxddCAQAAAECbStJQGdGm0paRkS0R//83Hty9XCCwFXTIYUccdcxxJ7zjXe856X2nfOC0D33kY58441Of+dxZ55x3wRcuuuRLX/naN771ne9ddsUPfvSTq372i1/95nd/+NM1f7nuhptu+dttd9x1zz/ue+ChRx574qln/vXcf1546X+vvPYGbVYMWQ== + + + AQAAAACAAABwAAAALQAAAA==eJxjZIAAZijNBqW5oDQ/lBaG0tJQWhlKa0JpPShtAaVdoLQHlA6A0gBDeAHg + + + AQAAAACAAAAOAAAAFgAAAA==eJxjZGblZOfi4ebl4xfQ0gIAA60AyQ== + + + AQAAAACAAAAgAQAARwAAAA==eJx1jjkOwDAIBPM4ErCd4/+/SeHZZiVoRhwDxLHjhBcMq2czl5aH1QsOOOGCN3zgC7/G15w83ydPe6rpy/e/5OnOD9yzClc= + + + AQAAAACAAABQAAAAIwAAAA==eJxjZoAANijNCaV5oLQAlBaB0hJQWgZKK0BpFSgNABawALs= + + + AQAAAACAAABQAAAAIgAAAA==eJwtxbcBACAIADAsiP7/sAPJkog2PL28nT4uXz9/BXgALg== + + + AQAAAACAAABwAAAAEQAAAA==eJxjYKAtYIHSXFAaAAEAAA8= + + + + + diff --git a/geos-mesh/tests/test_cliParsing.py b/mesh-doctor/tests/test_cliParsing.py similarity index 75% rename from geos-mesh/tests/test_cliParsing.py rename to mesh-doctor/tests/test_cliParsing.py index fdaff897..5fb0b3a9 100644 --- a/geos-mesh/tests/test_cliParsing.py +++ b/mesh-doctor/tests/test_cliParsing.py @@ -1,9 +1,12 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto import argparse from dataclasses import dataclass import pytest from typing import Iterator, Sequence -from geos.mesh.doctor.actions.generateFractures import FracturePolicy, Options -from geos.mesh.doctor.parsing.generateFracturesParsing import convert, displayResults, fillSubparser +from geos.mesh_doctor.actions.generateFractures import FracturePolicy, Options +from geos.mesh_doctor.parsing.generateFracturesParsing import convert, displayResults, fillSubparser from geos.mesh.io.vtkIO import VtkOutput @@ -17,10 +20,13 @@ class TestCase: def __generateGenerateFracturesParsingTestData() -> Iterator[ TestCase ]: field: str = "attribute" + inputMesh: str = "input.vtu" mainMesh: str = "output.vtu" fractureMesh: str = "fracture.vtu" - cliGen: str = f"generateFractures --policy {{}} --name {field} --values 0,1 --output {mainMesh} --fracturesOutputDir ." + cliBase: str = f"generateFractures -i {inputMesh} --policy {{}} --name {field} --values 0,1" + cliEnd: str = f"--output {mainMesh} --fracturesOutputDir ." + cliGen: str = f"{cliBase} {cliEnd}" allCliArgs = cliGen.format( "field" ).split(), cliGen.format( "internalSurfaces" ).split(), cliGen.format( "dummy" ).split() policies = FracturePolicy.FIELD, FracturePolicy.INTERNAL_SURFACES, FracturePolicy.FIELD @@ -35,9 +41,8 @@ def __generateGenerateFracturesParsingTestData() -> Iterator[ TestCase ]: yield TestCase( cliArgs, options, exception ) -def __parseAndValidateOptions( testCase: TestCase ): - """ - Parse CLI arguments and validate that the resulting options match expected values. +def __parseAndValidateOptions( testCase: TestCase ) -> None: + """Parse CLI arguments and validate that the resulting options match expected values. This helper function simulates the CLI parsing process by: 1. Creating an argument parser with the generateFractures subparser @@ -61,17 +66,17 @@ def __parseAndValidateOptions( testCase: TestCase ): assert options.fieldValuesCombined == testCase.options.fieldValuesCombined -def test_displayResults(): +def test_displayResults() -> None: + """Test displayResults function for code coverage.""" # Dummy test for code coverage only. Shame on me! displayResults( None, None ) @pytest.mark.parametrize( "testCase", __generateGenerateFracturesParsingTestData() ) -def test( testCase: TestCase ): +def test( testCase: TestCase ) -> None: + """Test CLI parsing for generateFractures action.""" if testCase.exception: with pytest.raises( SystemExit ): - pytest.skip( "Test to be fixed" ) __parseAndValidateOptions( testCase ) else: - pytest.skip( "Test to be fixed" ) __parseAndValidateOptions( testCase ) diff --git a/geos-mesh/tests/test_collocatedNodes.py b/mesh-doctor/tests/test_collocatedNodes.py similarity index 69% rename from geos-mesh/tests/test_collocatedNodes.py rename to mesh-doctor/tests/test_collocatedNodes.py index 7232a2a8..281f7df2 100644 --- a/geos-mesh/tests/test_collocatedNodes.py +++ b/mesh-doctor/tests/test_collocatedNodes.py @@ -1,15 +1,20 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto import pytest from typing import Iterator, Tuple from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkTetra, vtkUnstructuredGrid, VTK_TETRA -from geos.mesh.doctor.actions.collocatedNodes import Options, __action +from geos.mesh_doctor.actions.collocatedNodes import Options, meshAction def getPoints() -> Iterator[ Tuple[ vtkPoints, int ] ]: """Generates the data for the cases. - One case has two nodes at the exact same position. - The other has two differente nodes - :return: Generator to (vtk points, number of expected duplicated locations) + + One case has two nodes at the exact same position. The other has two different nodes. + + Returns: + Iterator[ Tuple[ vtkPoints, int ] ]: The points and the expected number of buckets of collocated nodes. """ for p0, p1 in ( ( 0, 0, 0 ), ( 1, 1, 1 ) ), ( ( 0, 0, 0 ), ( 0, 0, 0 ) ): points = vtkPoints() @@ -21,13 +26,14 @@ def getPoints() -> Iterator[ Tuple[ vtkPoints, int ] ]: @pytest.mark.parametrize( "data", getPoints() ) -def test_simpleCollocatedPoints( data: Tuple[ vtkPoints, int ] ): +def test_simpleCollocatedPoints( data: Tuple[ vtkPoints, int ] ) -> None: + """Tests the detection of collocated nodes for a simple case.""" points, numNodesBucket = data mesh = vtkUnstructuredGrid() mesh.SetPoints( points ) - result = __action( mesh, Options( tolerance=1.e-12 ) ) + result = meshAction( mesh, Options( tolerance=1.e-12 ) ) assert len( result.wrongSupportElements ) == 0 assert len( result.nodesBuckets ) == numNodesBucket @@ -35,7 +41,8 @@ def test_simpleCollocatedPoints( data: Tuple[ vtkPoints, int ] ): assert len( result.nodesBuckets[ 0 ] ) == points.GetNumberOfPoints() -def test_wrongSupportElements(): +def test_wrongSupportElements() -> None: + """Tests that a tetrahedron with two nodes at the same location is detected.""" points = vtkPoints() points.SetNumberOfPoints( 4 ) points.SetPoint( 0, ( 0, 0, 0 ) ) @@ -58,7 +65,7 @@ def test_wrongSupportElements(): mesh.SetPoints( points ) mesh.SetCells( cellTypes, cells ) - result = __action( mesh, Options( tolerance=1.e-12 ) ) + result = meshAction( mesh, Options( tolerance=1.e-12 ) ) assert len( result.nodesBuckets ) == 0 assert len( result.wrongSupportElements ) == 1 diff --git a/geos-mesh/tests/test_elementVolumes.py b/mesh-doctor/tests/test_elementVolumes.py similarity index 69% rename from geos-mesh/tests/test_elementVolumes.py rename to mesh-doctor/tests/test_elementVolumes.py index f62c74ec..669a5e92 100644 --- a/geos-mesh/tests/test_elementVolumes.py +++ b/mesh-doctor/tests/test_elementVolumes.py @@ -1,10 +1,14 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto import numpy from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkCommonDataModel import VTK_TETRA, vtkCellArray, vtkTetra, vtkUnstructuredGrid -from geos.mesh.doctor.actions.elementVolumes import Options, __action +from geos.mesh_doctor.actions.elementVolumes import Options, meshAction -def test_simpleTet(): +def test_simpleTet() -> None: + """Tests the calculation of element volumes for a simple tetrahedron.""" # creating a simple tetrahedron points = vtkPoints() points.SetNumberOfPoints( 4 ) @@ -28,12 +32,12 @@ def test_simpleTet(): mesh.SetPoints( points ) mesh.SetCells( cellTypes, cells ) - result = __action( mesh, Options( minVolume=1. ) ) + result = meshAction( mesh, Options( minVolume=1. ) ) assert len( result.elementVolumes ) == 1 assert result.elementVolumes[ 0 ][ 0 ] == 0 assert abs( result.elementVolumes[ 0 ][ 1 ] - 1. / 6. ) < 10 * numpy.finfo( float ).eps - result = __action( mesh, Options( minVolume=0. ) ) + result = meshAction( mesh, Options( minVolume=0. ) ) assert len( result.elementVolumes ) == 0 diff --git a/geos-mesh/tests/test_generateCube.py b/mesh-doctor/tests/test_generateCube.py similarity index 66% rename from geos-mesh/tests/test_generateCube.py rename to mesh-doctor/tests/test_generateCube.py index 069aba5e..c087bc46 100644 --- a/geos-mesh/tests/test_generateCube.py +++ b/mesh-doctor/tests/test_generateCube.py @@ -1,7 +1,11 @@ -from geos.mesh.doctor.actions.generateCube import __build, Options, FieldInfo +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from geos.mesh_doctor.actions.generateCube import Options, FieldInfo, buildCube -def test_generateCube(): +def test_generateCube() -> None: + """Tests the generation of a cube mesh with specific options.""" options = Options( vtkOutput=None, generateCellsGlobalIds=True, generatePointsGlobalIds=False, @@ -12,7 +16,7 @@ def test_generateCube(): nys=( 1, 1 ), nzs=( 1, ), fields=( FieldInfo( name="test", dimension=2, support="CELLS" ), ) ) - output = __build( options ) + output = buildCube( options ) assert output.GetNumberOfCells() == 14 assert output.GetNumberOfPoints() == 48 assert output.GetCellData().GetArray( "test" ).GetNumberOfComponents() == 2 diff --git a/geos-mesh/tests/test_generateFractures.py b/mesh-doctor/tests/test_generateFractures.py similarity index 59% rename from geos-mesh/tests/test_generateFractures.py rename to mesh-doctor/tests/test_generateFractures.py index cef3f325..9848b9c3 100644 --- a/geos-mesh/tests/test_generateFractures.py +++ b/mesh-doctor/tests/test_generateFractures.py @@ -1,17 +1,22 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto from dataclasses import dataclass -import numpy +import numpy as np +import numpy.typing as npt import pytest -from typing import Iterable, Iterator, Sequence +from typing import Iterable, Iterator, Optional, 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.actions.checkFractures import formatCollocatedNodes -from geos.mesh.doctor.actions.generateCube import buildRectilinearBlocksMesh, XYZ -from geos.mesh.doctor.actions.generateFractures import ( __splitMeshOnFractures, Options, FracturePolicy, Coordinates3D, - IDMapping ) from geos.mesh.utils.genericHelpers import toVtkIdList +from geos.mesh_doctor.actions.checkFractures import formatCollocatedNodes +from geos.mesh_doctor.actions.generateCube import buildRectilinearBlocksMesh, XYZ +from geos.mesh_doctor.actions.generateFractures import ( __splitMeshOnFractures, Options, FracturePolicy, Coordinates3D, + IDMapping ) -FaceNodesCoords = tuple[ tuple[ float ] ] -IDMatrix = Sequence[ Sequence[ int ] ] +BorderFacesNodesCoords: TypeAlias = tuple[ tuple[ Coordinates3D, ...], ...] +FaceNodesCoords: TypeAlias = tuple[ Coordinates3D, ...] +IDMatrix: TypeAlias = Sequence[ Sequence[ int ] ] @dataclass( frozen=True ) @@ -32,25 +37,22 @@ class TestCase: result: TestResult -def __buildTestCase( xs: tuple[ numpy.ndarray, numpy.ndarray, numpy.ndarray ], +def __buildTestCase( xs: tuple[ npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ], npt.NDArray[ np.float64 ] ], attribute: Iterable[ int ], - fieldValues: Iterable[ int ] = None, - policy: FracturePolicy = FracturePolicy.FIELD ): + fieldValues: Optional[ Iterable[ int ] ] = None, + policy: FracturePolicy = FracturePolicy.FIELD ) -> tuple[ vtkUnstructuredGrid, Options ]: + """Builds a test case mesh and options for fracture generation testing.""" xyz = XYZ( *xs ) - mesh: vtkUnstructuredGrid = buildRectilinearBlocksMesh( ( xyz, ) ) - ref = numpy.array( attribute, dtype=int ) + ref = np.array( attribute, dtype=int ) if policy == FracturePolicy.FIELD: assert len( ref ) == mesh.GetNumberOfCells() attr = numpy_to_vtk( ref ) attr.SetName( "attribute" ) mesh.GetCellData().AddArray( attr ) - if fieldValues is None: - fv = frozenset( attribute ) - else: - fv = frozenset( fieldValues ) + fv = frozenset( attribute ) if fieldValues is None else frozenset( fieldValues ) options = Options( policy=policy, field="attribute", @@ -64,112 +66,115 @@ def __buildTestCase( xs: tuple[ numpy.ndarray, numpy.ndarray, numpy.ndarray ], # Utility class to generate the new indices of the newly created collocated nodes. class Incrementor: - def __init__( self, start ): + def __init__( self, start: int ) -> None: + """Initializes the incrementor with a starting value.""" self.__val = start def next( self, num: int ) -> Iterable[ int ]: + """Generates the next 'num' values in the incrementor sequence.""" self.__val += num return range( self.__val - num, self.__val ) def __generateTestData() -> Iterator[ TestCase ]: - twoNodes = numpy.arange( 2, dtype=float ) - threeNodes = numpy.arange( 3, dtype=float ) - fourNodes = numpy.arange( 4, dtype=float ) + """Generates test data for fracture generation tests.""" + twoNodes = np.arange( 2, dtype=float ) + threeNodes = np.arange( 3, dtype=float ) + fourNodes = np.arange( 4, dtype=float ) # Split in 2 mesh, options = __buildTestCase( ( threeNodes, threeNodes, threeNodes ), ( 0, 1, 0, 1, 0, 1, 0, 1 ) ) yield TestCase( inputMesh=mesh, options=options, - collocatedNodes=tuple( map( lambda i: ( 1 + 3 * i, 27 + i ), range( 9 ) ) ), + collocatedNodes=tuple( ( 1 + 3 * i, 27 + i ) for i in range( 9 ) ), result=TestResult( 9 * 4, 8, 9, 4 ) ) # Split in 3 inc = Incrementor( 27 ) - collocatedNodes: 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 ) ) ) + collocatedNodes0: 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 = __buildTestCase( ( threeNodes, threeNodes, threeNodes ), ( 0, 1, 2, 1, 0, 1, 2, 1 ) ) yield TestCase( inputMesh=mesh, options=options, - collocatedNodes=collocatedNodes, + collocatedNodes=collocatedNodes0, result=TestResult( 9 * 4 + 6, 8, 12, 6 ) ) # Split in 8 inc = Incrementor( 27 ) - collocatedNodes: 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 ) ) ) + collocatedNodes2: 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 = __buildTestCase( ( threeNodes, threeNodes, threeNodes ), range( 8 ) ) yield TestCase( inputMesh=mesh, options=options, - collocatedNodes=collocatedNodes, + collocatedNodes=collocatedNodes2, result=TestResult( 8 * 8, 8, 3 * 3 * 3 - 8, 12 ) ) # Straight notch inc = Incrementor( 27 ) - collocatedNodes: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 4, ), ( 1 + 9, *inc.next( 1 ) ), ( 4 + 9, ), - ( 1 + 18, *inc.next( 1 ) ), ( 4 + 18, ) ) + collocatedNodes3: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 4, ), ( 1 + 9, *inc.next( 1 ) ), ( 4 + 9, ), + ( 1 + 18, *inc.next( 1 ) ), ( 4 + 18, ) ) mesh, options = __buildTestCase( ( threeNodes, threeNodes, threeNodes ), ( 0, 1, 2, 2, 0, 1, 2, 2 ), fieldValues=( 0, 1 ) ) yield TestCase( inputMesh=mesh, options=options, - collocatedNodes=collocatedNodes, + collocatedNodes=collocatedNodes3, result=TestResult( 3 * 3 * 3 + 3, 8, 6, 2 ) ) # L-shaped notch inc = Incrementor( 27 ) - collocatedNodes: 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, ) ) + collocatedNodes4: 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 = __buildTestCase( ( threeNodes, threeNodes, threeNodes ), ( 0, 1, 0, 1, 0, 1, 2, 2 ), fieldValues=( 0, 1 ) ) yield TestCase( inputMesh=mesh, options=options, - collocatedNodes=collocatedNodes, + collocatedNodes=collocatedNodes4, result=TestResult( 3 * 3 * 3 + 5, 8, 8, 3 ) ) # 3x1x1 split inc = Incrementor( 2 * 2 * 4 ) - collocatedNodes: 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 ) ) ) + collocatedNodes5: 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 = __buildTestCase( ( fourNodes, twoNodes, twoNodes ), ( 0, 1, 2 ) ) yield TestCase( inputMesh=mesh, options=options, - collocatedNodes=collocatedNodes, + collocatedNodes=collocatedNodes5, result=TestResult( 6 * 4, 3, 2 * 4, 2 ) ) # Discarded fracture element if no node duplication. - collocatedNodes: IDMatrix = tuple() + collocatedNodes6: IDMatrix = () mesh, options = __buildTestCase( ( threeNodes, fourNodes, fourNodes ), ( 0, ) * 8 + ( 1, 2 ) + ( 0, ) * 8, fieldValues=( 1, 2 ) ) yield TestCase( inputMesh=mesh, options=options, - collocatedNodes=collocatedNodes, + collocatedNodes=collocatedNodes6, result=TestResult( 3 * 4 * 4, 2 * 3 * 3, 0, 0 ) ) # Fracture on a corner inc = Incrementor( 3 * 4 * 4 ) - collocatedNodes: 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, ) ) + collocatedNodes7: 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 = __buildTestCase( ( threeNodes, fourNodes, fourNodes ), ( 0, ) * 6 + ( 1, 2, 1, 2, 0, 0, 1, 2, 1, 2, 0, 0 ), fieldValues=( 1, 2 ) ) yield TestCase( inputMesh=mesh, options=options, - collocatedNodes=collocatedNodes, + collocatedNodes=collocatedNodes7, result=TestResult( 3 * 4 * 4 + 4, 2 * 3 * 3, 9, 4 ) ) # Generate mesh with 2 hexs, one being a standard hex, the other a 42 hex. inc = Incrementor( 3 * 2 * 2 ) - collocatedNodes: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 1 + 3, *inc.next( 1 ) ), ( 1 + 6, *inc.next( 1 ) ), - ( 1 + 9, *inc.next( 1 ) ) ) + collocatedNodes8: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 1 + 3, *inc.next( 1 ) ), ( 1 + 6, *inc.next( 1 ) ), + ( 1 + 9, *inc.next( 1 ) ) ) mesh, options = __buildTestCase( ( threeNodes, twoNodes, twoNodes ), ( 0, 1 ) ) polyhedronMesh = vtkUnstructuredGrid() polyhedronMesh.SetPoints( mesh.GetPoints() ) @@ -182,13 +187,13 @@ def __generateTestData() -> Iterator[ TestCase ]: yield TestCase( inputMesh=polyhedronMesh, options=options, - collocatedNodes=collocatedNodes, + collocatedNodes=collocatedNodes8, result=TestResult( 4 * 4, 2, 4, 1 ) ) # Split in 2 using the internal fracture description inc = Incrementor( 3 * 2 * 2 ) - collocatedNodes: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 1 + 3, *inc.next( 1 ) ), ( 1 + 6, *inc.next( 1 ) ), - ( 1 + 9, *inc.next( 1 ) ) ) + collocatedNodes9: IDMatrix = ( ( 1, *inc.next( 1 ) ), ( 1 + 3, *inc.next( 1 ) ), ( 1 + 6, *inc.next( 1 ) ), + ( 1 + 9, *inc.next( 1 ) ) ) mesh, options = __buildTestCase( ( threeNodes, twoNodes, twoNodes ), attribute=( 0, 0, 0 ), fieldValues=( 0, ), @@ -196,12 +201,13 @@ def __generateTestData() -> Iterator[ TestCase ]: mesh.InsertNextCell( VTK_QUAD, toVtkIdList( ( 1, 4, 7, 10 ) ) ) # Add a fracture on the fly yield TestCase( inputMesh=mesh, options=options, - collocatedNodes=collocatedNodes, + collocatedNodes=collocatedNodes9, result=TestResult( 4 * 4, 3, 4, 1 ) ) @pytest.mark.parametrize( "TestCase", __generateTestData() ) -def test_generateFracture( TestCase: TestCase ): +def test_generateFracture( TestCase: TestCase ) -> None: + """Tests the generation of fractures on a mesh according to various test cases.""" mainMesh, fractureMeshes = __splitMeshOnFractures( TestCase.inputMesh, TestCase.options ) fractureMesh: vtkUnstructuredGrid = fractureMeshes[ 0 ] assert mainMesh.GetNumberOfPoints() == TestCase.result.mainMeshNumPoints @@ -214,9 +220,10 @@ def test_generateFracture( TestCase: TestCase ): assert len( res ) == TestCase.result.fractureMeshNumPoints -def addSimplifiedFieldForCells( mesh: vtkUnstructuredGrid, field_name: str, fieldDimension: int ): - """Reduce functionality obtained from src.geos.mesh.doctor.actions.generateFractures.__add_fields - where the goal is to add a cell data array with incrementing values. +def addSimplifiedFieldForCells( mesh: vtkUnstructuredGrid, field_name: str, fieldDimension: int ) -> None: + """Reduce functionality obtained from src.geos.mesh_doctor.actions.generateFractures.__add_fields. + + The goal is to add a cell data array with incrementing values. Args: mesh (vtkUnstructuredGrid): Unstructured mesh. @@ -225,15 +232,16 @@ def addSimplifiedFieldForCells( mesh: vtkUnstructuredGrid, field_name: str, fiel """ data = mesh.GetCellData() n = mesh.GetNumberOfCells() - array = numpy.ones( ( n, fieldDimension ), dtype=float ) - array = numpy.arange( 1, n * fieldDimension + 1 ).reshape( n, fieldDimension ) + array = np.ones( ( n, fieldDimension ), dtype=float ) + array = np.arange( 1, n * fieldDimension + 1 ).reshape( n, fieldDimension ) vtkArray = numpy_to_vtk( array ) vtkArray.SetName( field_name ) data.AddArray( vtkArray ) -def findBordersFacesRectilinearGrid( mesh: vtkUnstructuredGrid ) -> tuple[ FaceNodesCoords ]: - """ +def findBordersFacesRectilinearGrid( mesh: vtkUnstructuredGrid ) -> BorderFacesNodesCoords: + """For a vtk rectilinear grid, gives the coordinates of each of its borders face nodes. + 6+--------+7 / /| / / | @@ -244,20 +252,18 @@ def findBordersFacesRectilinearGrid( mesh: vtkUnstructuredGrid ) -> tuple[ FaceN | |/ 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. + BorderFacesNodesCoords: For a rectilinear grid, returns a tuple of 6 faces nodeset. """ - meshBounds: tuple[ float ] = mesh.GetBounds() - minBound: Coordinates3D = [ meshBounds[ i ] for i in range( len( meshBounds ) ) if i % 2 == 0 ] - maxBound: Coordinates3D = [ meshBounds[ i ] for i in range( len( meshBounds ) ) if i % 2 == 1 ] + meshBounds: tuple[ float, float, float, float, float, float ] = mesh.GetBounds() + minBound: tuple[ float, ...] = tuple( [ meshBounds[ i ] for i in range( len( meshBounds ) ) if i % 2 == 0 ] ) + maxBound: tuple[ float, ...] = tuple( [ meshBounds[ i ] for i in range( len( meshBounds ) ) if i % 2 == 1 ] ) center: Coordinates3D = mesh.GetCenter() - faceDiag: tuple[ float ] = ( ( maxBound[ 0 ] - minBound[ 0 ] ) / 2, ( maxBound[ 1 ] - minBound[ 1 ] ) / 2, - ( maxBound[ 2 ] - minBound[ 2 ] ) / 2 ) + faceDiag: Coordinates3D = ( ( maxBound[ 0 ] - minBound[ 0 ] ) / 2, ( maxBound[ 1 ] - minBound[ 1 ] ) / 2, + ( maxBound[ 2 ] - minBound[ 2 ] ) / 2 ) node0: Coordinates3D = ( center[ 0 ] - faceDiag[ 0 ], center[ 1 ] - faceDiag[ 1 ], center[ 2 ] - faceDiag[ 2 ] ) node1: Coordinates3D = ( center[ 0 ] + faceDiag[ 0 ], center[ 1 ] - faceDiag[ 1 ], center[ 2 ] - faceDiag[ 2 ] ) node2: Coordinates3D = ( center[ 0 ] - faceDiag[ 0 ], center[ 1 ] + faceDiag[ 1 ], center[ 2 ] - faceDiag[ 2 ] ) @@ -266,17 +272,18 @@ def findBordersFacesRectilinearGrid( mesh: vtkUnstructuredGrid ) -> tuple[ FaceN node5: Coordinates3D = ( center[ 0 ] + faceDiag[ 0 ], center[ 1 ] - faceDiag[ 1 ], center[ 2 ] + faceDiag[ 2 ] ) node6: Coordinates3D = ( center[ 0 ] - faceDiag[ 0 ], center[ 1 ] + faceDiag[ 1 ], center[ 2 ] + faceDiag[ 2 ] ) node7: Coordinates3D = ( center[ 0 ] + faceDiag[ 0 ], center[ 1 ] + faceDiag[ 1 ], center[ 2 ] + faceDiag[ 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 ) ) + faces: BorderFacesNodesCoords = ( ( 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 addQuad( mesh: vtkUnstructuredGrid, face: FaceNodesCoords ): +def addQuad( mesh: vtkUnstructuredGrid, face: FaceNodesCoords ) -> None: """Adds a quad cell to each border of an unstructured mesh. Args: mesh (vtkUnstructuredGrid): Unstructured mesh. + face (FaceNodesCoords): Coordinates of the quad face to add. """ pointsCoords = mesh.GetPoints().GetData() quad: vtkQuad = vtkQuad() @@ -296,32 +303,29 @@ def addQuad( mesh: vtkUnstructuredGrid, face: FaceNodesCoords ): mesh.InsertNextCell( quad.GetCellType(), quad.GetPointIds() ) -@pytest.mark.skip( "Test to be fixed" ) -def test_copyFieldsWhenSplittingMesh(): - """This test is designed to check the __copyFields method from generateFractures, - that will be called when using __splitMeshOnFractures method from generateFractures. - """ +def test_copyFieldsWhenSplittingMesh() -> None: + """This test is designed to check the __copyFields method from generateFractures, that will be called when using __splitMeshOnFractures method from generateFractures.""" # 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 ] ) + x: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 1.0, 2.0 ] ) + y: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 1.0 ] ) + z: npt.NDArray[ np.float64 ] = np.array( [ 0.0, 1.0 ] ) xyzs: XYZ = XYZ( x, y, z ) - mesh: vtkUnstructuredGrid = buildRectilinearBlocksMesh( [ xyzs ] ) - assert mesh.GetCells().GetNumberOfCells() == 2 - borderFaces: tuple[ FaceNodesCoords ] = findBordersFacesRectilinearGrid( mesh ) - for face in borderFaces: - addQuad( mesh, face ) - assert mesh.GetCells().GetNumberOfCells() == 8 + mesh0: vtkUnstructuredGrid = buildRectilinearBlocksMesh( [ xyzs ] ) + assert mesh0.GetCells().GetNumberOfCells() == 2 + borderFaces0: BorderFacesNodesCoords = findBordersFacesRectilinearGrid( mesh0 ) + for face in borderFaces0: + addQuad( mesh0, face ) + assert mesh0.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 ) ) - addQuad( mesh, fracture ) - assert mesh.GetCells().GetNumberOfCells() == 9 + addQuad( mesh0, fracture ) + assert mesh0.GetCells().GetNumberOfCells() == 9 # Add a "TestField" array - assert mesh.GetCellData().GetNumberOfArrays() == 0 - addSimplifiedFieldForCells( mesh, "TestField", 1 ) - assert mesh.GetCellData().GetNumberOfArrays() == 1 - assert mesh.GetCellData().GetArrayName( 0 ) == "TestField" - testFieldValues: list[ int ] = vtk_to_numpy( mesh.GetCellData().GetArray( 0 ) ).tolist() + assert mesh0.GetCellData().GetNumberOfArrays() == 0 + addSimplifiedFieldForCells( mesh0, "TestField", 1 ) + assert mesh0.GetCellData().GetNumberOfArrays() == 1 + assert mesh0.GetCellData().GetArrayName( 0 ) == "TestField" + testFieldValues: list[ int ] = vtk_to_numpy( mesh0.GetCellData().GetArray( 0 ) ).tolist() assert testFieldValues == [ 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, @@ -330,7 +334,7 @@ def test_copyFieldsWhenSplittingMesh(): fieldValuesPerFracture=[ frozenset( map( int, [ "9" ] ) ) ], meshVtkOutput=None, allFracturesVtkOutput=None ) - mainMesh, fractureMeshes = __splitMeshOnFractures( mesh, options ) + mainMesh, fractureMeshes = __splitMeshOnFractures( mesh0, options ) fractureMesh: vtkUnstructuredGrid = fractureMeshes[ 0 ] assert mainMesh.GetCellData().GetNumberOfArrays() == 1 assert fractureMesh.GetCellData().GetNumberOfArrays() == 1 @@ -341,20 +345,3 @@ def test_copyFieldsWhenSplittingMesh(): mainMeshValues: list[ int ] = vtk_to_numpy( mainMesh.GetCellData().GetArray( 0 ) ).tolist() assert fractureMeshValues == [ 9 ] # The value for the fracture surface assert mainMeshValues == [ 1, 2, 3, 4, 5, 6, 7, 8, 9 ] - # Test for invalid point field name - addSimplifiedFieldForCells( mesh, "GLOBAL_IDS_POINTS", 1 ) - with pytest.raises( ValueError ) as pytestWrappedError: - mainMesh, fractureMeshes = __splitMeshOnFractures( mesh, options ) - assert pytestWrappedError.type == ValueError - # Test for invalid cell field name - mesh: vtkUnstructuredGrid = buildRectilinearBlocksMesh( [ xyzs ] ) - borderFaces: tuple[ FaceNodesCoords ] = findBordersFacesRectilinearGrid( mesh ) - for face in borderFaces: - addQuad( mesh, face ) - addQuad( mesh, fracture ) - addSimplifiedFieldForCells( mesh, "TestField", 1 ) - addSimplifiedFieldForCells( mesh, "GLOBAL_IDS_CELLS", 1 ) - assert mesh.GetCellData().GetNumberOfArrays() == 2 - with pytest.raises( ValueError ) as pytestWrappedError: - mainMesh, fractureMeshes = __splitMeshOnFractures( mesh, options ) - assert pytestWrappedError.type == ValueError diff --git a/geos-mesh/tests/test_generateGlobalIds.py b/mesh-doctor/tests/test_generateGlobalIds.py similarity index 63% rename from geos-mesh/tests/test_generateGlobalIds.py rename to mesh-doctor/tests/test_generateGlobalIds.py index f623ca27..a5549369 100644 --- a/geos-mesh/tests/test_generateGlobalIds.py +++ b/mesh-doctor/tests/test_generateGlobalIds.py @@ -1,9 +1,13 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkUnstructuredGrid, vtkVertex, VTK_VERTEX -from geos.mesh.doctor.actions.generateGlobalIds import __buildGlobalIds +from geos.mesh_doctor.actions.generateGlobalIds import buildGlobalIds -def test_generateGlobalIds(): +def test_generateGlobalIds() -> None: + """Tests the generation of global IDs for a simple mesh with one vertex.""" points = vtkPoints() points.InsertNextPoint( 0, 0, 0 ) @@ -17,7 +21,7 @@ def test_generateGlobalIds(): mesh.SetPoints( points ) mesh.SetCells( [ VTK_VERTEX ], vertices ) - __buildGlobalIds( mesh, True, True ) + buildGlobalIds( mesh, True, True ) globalCellIds = mesh.GetCellData().GetGlobalIds() globalPointIds = mesh.GetPointData().GetGlobalIds() diff --git a/geos-mesh/tests/test_nonConformal.py b/mesh-doctor/tests/test_nonConformal.py similarity index 66% rename from geos-mesh/tests/test_nonConformal.py rename to mesh-doctor/tests/test_nonConformal.py index 231d01ed..3d8df190 100644 --- a/geos-mesh/tests/test_nonConformal.py +++ b/mesh-doctor/tests/test_nonConformal.py @@ -1,9 +1,13 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto import numpy -from geos.mesh.doctor.actions.generateCube import buildRectilinearBlocksMesh, XYZ -from geos.mesh.doctor.actions.nonConformal import Options, __action +from geos.mesh_doctor.actions.generateCube import buildRectilinearBlocksMesh, XYZ +from geos.mesh_doctor.actions.nonConformal import Options, meshAction -def test_twoCloseHexs(): +def test_twoCloseHexs() -> None: + """Tests two close hexahedrons for non-conformality detection.""" delta = 1.e-6 tmp = numpy.arange( 2, dtype=float ) xyz0 = XYZ( tmp, tmp, tmp ) @@ -12,17 +16,18 @@ def test_twoCloseHexs(): # Close enough, but points tolerance is too strict to consider the faces matching. options = Options( angleTolerance=1., pointTolerance=delta / 2, faceTolerance=delta * 2 ) - results = __action( mesh, options ) + results = meshAction( mesh, options ) assert len( results.nonConformalCells ) == 1 assert set( results.nonConformalCells[ 0 ] ) == { 0, 1 } # Close enough, and points tolerance is loose enough to consider the faces matching. options = Options( angleTolerance=1., pointTolerance=delta * 2, faceTolerance=delta * 2 ) - results = __action( mesh, options ) + results = meshAction( mesh, options ) assert len( results.nonConformalCells ) == 0 -def test_twoDistantHexs(): +def test_twoDistantHexs() -> None: + """Tests two distant hexahedrons for non-conformality detection.""" delta = 1 tmp = numpy.arange( 2, dtype=float ) xyz0 = XYZ( tmp, tmp, tmp ) @@ -31,11 +36,12 @@ def test_twoDistantHexs(): options = Options( angleTolerance=1., pointTolerance=delta / 2., faceTolerance=delta / 2. ) - results = __action( mesh, options ) + results = meshAction( mesh, options ) assert len( results.nonConformalCells ) == 0 -def test_twoCloseShiftedHexs(): +def test_twoCloseShiftedHexs() -> None: + """Tests two close but shifted hexahedrons for non-conformality detection.""" deltaX, deltaY = 1.e-6, 0.5 tmp = numpy.arange( 2, dtype=float ) xyz0 = XYZ( tmp, tmp, tmp ) @@ -44,12 +50,13 @@ def test_twoCloseShiftedHexs(): options = Options( angleTolerance=1., pointTolerance=deltaX * 2, faceTolerance=deltaX * 2 ) - results = __action( mesh, options ) + results = meshAction( mesh, options ) assert len( results.nonConformalCells ) == 1 assert set( results.nonConformalCells[ 0 ] ) == { 0, 1 } -def test_bigElemNextToSmallElem(): +def test_bigElemNextToSmallElem() -> None: + """Tests a big element next to a small element for non-conformality detection.""" delta = 1.e-6 tmp = numpy.arange( 2, dtype=float ) xyz0 = XYZ( tmp, tmp + 1, tmp + 1 ) @@ -58,6 +65,6 @@ def test_bigElemNextToSmallElem(): options = Options( angleTolerance=1., pointTolerance=delta * 2, faceTolerance=delta * 2 ) - results = __action( mesh, options ) + results = meshAction( mesh, options ) assert len( results.nonConformalCells ) == 1 assert set( results.nonConformalCells[ 0 ] ) == { 0, 1 } diff --git a/geos-mesh/tests/test_reorientMesh.py b/mesh-doctor/tests/test_reorientMesh.py similarity index 86% rename from geos-mesh/tests/test_reorientMesh.py rename to mesh-doctor/tests/test_reorientMesh.py index 3f7277a8..dd3ffa8f 100644 --- a/geos-mesh/tests/test_reorientMesh.py +++ b/mesh-doctor/tests/test_reorientMesh.py @@ -1,11 +1,14 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto from dataclasses import dataclass import numpy import pytest from typing import Generator from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, VTK_POLYHEDRON -from geos.mesh.doctor.actions.reorientMesh import reorientMesh -from geos.mesh.doctor.actions.vtkPolyhedron import FaceStream +from geos.mesh_doctor.actions.reorientMesh import reorientMesh +from geos.mesh_doctor.actions.vtkPolyhedron import FaceStream from geos.mesh.utils.genericHelpers import toVtkIdList, vtkIter @@ -16,6 +19,7 @@ class Expected: def __buildTestMeshes() -> Generator[ Expected, None, None ]: + """Builds test meshes for reorientMesh testing.""" # Creating the support nodes for the polyhedron. # It has a C shape and is actually non-convex, non star-shaped. frontNodes = numpy.array( ( @@ -42,7 +46,7 @@ def __buildTestMeshes() -> Generator[ Expected, None, None ]: points.InsertNextPoint( coords ) # Creating the polyhedron with faces all directed outward. - faces = [] + faces: list[ tuple[ int, ...] ] = [] # Creating the side faces for i in range( n ): faces.append( ( i % n + n, ( i + 1 ) % n + n, ( i + 1 ) % n, i % n ) ) @@ -76,7 +80,8 @@ def __buildTestMeshes() -> Generator[ Expected, None, None ]: @pytest.mark.parametrize( "expected", __buildTestMeshes() ) -def test_reorientPolyhedron( expected: Expected ): +def test_reorientPolyhedron( expected: Expected ) -> None: + """Tests reorientMesh on polyhedron elements.""" outputMesh = reorientMesh( expected.mesh, range( expected.mesh.GetNumberOfCells() ) ) assert outputMesh.GetNumberOfCells() == 1 assert outputMesh.GetCell( 0 ).GetCellType() == VTK_POLYHEDRON diff --git a/geos-mesh/tests/test_selfIntersectingElements.py b/mesh-doctor/tests/test_selfIntersectingElements.py similarity index 67% rename from geos-mesh/tests/test_selfIntersectingElements.py rename to mesh-doctor/tests/test_selfIntersectingElements.py index a1783e59..a477e912 100644 --- a/geos-mesh/tests/test_selfIntersectingElements.py +++ b/mesh-doctor/tests/test_selfIntersectingElements.py @@ -1,9 +1,13 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkHexahedron, vtkUnstructuredGrid, VTK_HEXAHEDRON -from geos.mesh.doctor.actions.selfIntersectingElements import Options, __action +from geos.mesh_doctor.actions.selfIntersectingElements import Options, meshAction -def test_jumbledHex(): +def test_jumbledHex() -> None: + """Tests that a hexahedron with intersecting faces is detected.""" # creating a simple hexahedron points = vtkPoints() points.SetNumberOfPoints( 8 ) @@ -35,7 +39,7 @@ def test_jumbledHex(): mesh.SetPoints( points ) mesh.SetCells( cellTypes, cells ) - result = __action( mesh, Options( minDistance=0. ) ) + result = meshAction( mesh, Options( minDistance=0. ) ) - assert len( result.intersectingFacesElements ) == 1 - assert result.intersectingFacesElements[ 0 ] == 0 + assert len( result.invalidCellIds[ "intersectingFacesElements" ] ) == 1 + assert result.invalidCellIds[ "intersectingFacesElements" ][ 0 ] == 0 diff --git a/geos-mesh/tests/test_sharedChecksParsingLogic.py b/mesh-doctor/tests/test_sharedChecksParsingLogic.py similarity index 59% rename from geos-mesh/tests/test_sharedChecksParsingLogic.py rename to mesh-doctor/tests/test_sharedChecksParsingLogic.py index ec6420f5..91ca47d5 100644 --- a/geos-mesh/tests/test_sharedChecksParsingLogic.py +++ b/mesh-doctor/tests/test_sharedChecksParsingLogic.py @@ -1,11 +1,14 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto import argparse from dataclasses import dataclass import pytest -from unittest.mock import patch +from unittest.mock import MagicMock, patch # Import the module to test -from geos.mesh.doctor.actions.allChecks import Options as AllChecksOptions -from geos.mesh.doctor.actions.allChecks import Result as AllChecksResult -from geos.mesh.doctor.parsing._sharedChecksParsingLogic import ( CheckFeature, _generateParametersHelp, +from geos.mesh_doctor.actions.allChecks import Options as AllChecksOptions +from geos.mesh_doctor.actions.allChecks import Result as AllChecksResult +from geos.mesh_doctor.parsing._sharedChecksParsingLogic import ( CheckFeature, _generateParametersHelp, getOptionsUsedMessage, fillSubparser, convert, displayResults, CHECKS_TO_DO_ARG, PARAMETERS_ARG ) @@ -22,12 +25,14 @@ class MockResult: value: str = "testResult" -def mockDisplayFunc( options, result ): +def mockDisplayFunc( options: MockOptions, result: MockResult ) -> None: + """Mock display function for testing.""" pass @pytest.fixture -def checkFeaturesConfig(): +def checkFeaturesConfig() -> dict[ str, CheckFeature ]: + """Provides a mock check features configuration for testing.""" return { "check1": CheckFeature( name="check1", @@ -51,17 +56,21 @@ def checkFeaturesConfig(): @pytest.fixture -def orderedCheckNames(): +def orderedCheckNames() -> list[ str ]: + """Provides an ordered list of check names for testing.""" return [ "check1", "check2" ] -def test_generateParametersHelp( checkFeaturesConfig, orderedCheckNames ): +def test_generateParametersHelp( checkFeaturesConfig: dict[ str, CheckFeature ], + orderedCheckNames: list[ str ] ) -> None: + """Tests _generateParametersHelp functionality.""" helpText = _generateParametersHelp( orderedCheckNames, checkFeaturesConfig ) assert "For check1: param1:1.0, param2:2.0" in helpText assert "For check2: param1:3.0, param2:4.0" in helpText -def test_getOptionsUsedMessage(): +def test_getOptionsUsedMessage() -> None: + """Tests getOptionsUsedMessage functionality.""" options = MockOptions( param1=10.0, param2=20.0 ) message = getOptionsUsedMessage( options ) assert "Parameters used: (" in message @@ -70,23 +79,28 @@ def test_getOptionsUsedMessage(): assert ")." in message -def test_fillSubparser( checkFeaturesConfig, orderedCheckNames ): +def test_fillSubparser( checkFeaturesConfig: dict[ str, CheckFeature ], orderedCheckNames: list[ str ] ) -> None: + """Tests fillSubparser functionality.""" parser = argparse.ArgumentParser() subparsers = parser.add_subparsers( dest="command" ) fillSubparser( subparsers, "test-command", "Test help message", orderedCheckNames, checkFeaturesConfig ) - # Parse with no args should use defaults - args = parser.parse_args( [ "test-command" ] ) + # Parse with no args (except required vtu-input-file) should use defaults + args = parser.parse_args( [ "test-command", "-i", "test.vtu" ] ) assert getattr( args, CHECKS_TO_DO_ARG ) == "" assert getattr( args, PARAMETERS_ARG ) == "" + assert args.vtuInputFile == "test.vtu" # Parse with specified args args = parser.parse_args( - [ "test-command", f"--{CHECKS_TO_DO_ARG}", "check1", f"--{PARAMETERS_ARG}", "param1:10.5" ] ) + [ "test-command", "-i", "test.vtu", f"--{CHECKS_TO_DO_ARG}", "check1", f"--{PARAMETERS_ARG}", "param1:10.5" ] ) assert getattr( args, CHECKS_TO_DO_ARG ) == "check1" assert getattr( args, PARAMETERS_ARG ) == "param1:10.5" + assert args.vtuInputFile == "test.vtu" -@patch( 'geos.mesh.doctor.parsing._sharedChecksParsingLogic.setupLogger' ) -def test_convertDefaultChecks( mockLogger, checkFeaturesConfig, orderedCheckNames ): +@patch( 'geos.mesh_doctor.parsing._sharedChecksParsingLogic.setupLogger' ) +def test_convertDefaultChecks( mockLogger: MagicMock, checkFeaturesConfig: dict[ str, CheckFeature ], + orderedCheckNames: list[ str ] ) -> None: + """Tests convert when no specific checks or parameters are specified.""" parsedArgs = { CHECKS_TO_DO_ARG: "", PARAMETERS_ARG: "" } options = convert( parsedArgs, orderedCheckNames, checkFeaturesConfig ) assert options.checksToPerform == orderedCheckNames @@ -95,8 +109,10 @@ def test_convertDefaultChecks( mockLogger, checkFeaturesConfig, orderedCheckName assert options.checksOptions[ "check2" ].param2 == 4.0 -@patch( 'geos.mesh.doctor.parsing._sharedChecksParsingLogic.setupLogger' ) -def test_convertSpecificChecks( mockLogger, checkFeaturesConfig, orderedCheckNames ): +@patch( 'geos.mesh_doctor.parsing._sharedChecksParsingLogic.setupLogger' ) +def test_convertSpecificChecks( mockLogger: MagicMock, checkFeaturesConfig: dict[ str, CheckFeature ], + orderedCheckNames: list[ str ] ) -> None: + """Tests convert when specific checks are specified.""" parsedArgs = { CHECKS_TO_DO_ARG: "check1", PARAMETERS_ARG: "" } options = convert( parsedArgs, orderedCheckNames, checkFeaturesConfig ) assert options.checksToPerform == [ "check1" ] @@ -105,8 +121,10 @@ def test_convertSpecificChecks( mockLogger, checkFeaturesConfig, orderedCheckNam assert "check2" not in options.checksOptions -@patch( 'geos.mesh.doctor.parsing._sharedChecksParsingLogic.setupLogger' ) -def test_convertWithParameters( mockLogger, checkFeaturesConfig, orderedCheckNames ): +@patch( 'geos.mesh_doctor.parsing._sharedChecksParsingLogic.setupLogger' ) +def test_convertWithParameters( mockLogger: MagicMock, checkFeaturesConfig: dict[ str, CheckFeature ], + orderedCheckNames: list[ str ] ) -> None: + """Tests convert when parameters are specified.""" parsedArgs = { CHECKS_TO_DO_ARG: "", PARAMETERS_ARG: "param1:10.5,param2:20.5" } options = convert( parsedArgs, orderedCheckNames, checkFeaturesConfig ) assert options.checksToPerform == orderedCheckNames @@ -116,8 +134,10 @@ def test_convertWithParameters( mockLogger, checkFeaturesConfig, orderedCheckNam assert options.checksOptions[ "check2" ].param2 == 20.5 -@patch( 'geos.mesh.doctor.parsing._sharedChecksParsingLogic.setupLogger' ) -def test_convertWithInvalidParameters( mockLogger, checkFeaturesConfig, orderedCheckNames ): +@patch( 'geos.mesh_doctor.parsing._sharedChecksParsingLogic.setupLogger' ) +def test_convertWithInvalidParameters( mockLogger: MagicMock, checkFeaturesConfig: dict[ str, CheckFeature ], + orderedCheckNames: list[ str ] ) -> None: + """Tests convert when some invalid parameters are specified.""" parsedArgs = { CHECKS_TO_DO_ARG: "", PARAMETERS_ARG: "param1:invalid,param2:20.5" } options = convert( parsedArgs, orderedCheckNames, checkFeaturesConfig ) # The invalid parameter should be skipped, but the valid one applied @@ -125,8 +145,10 @@ def test_convertWithInvalidParameters( mockLogger, checkFeaturesConfig, orderedC assert options.checksOptions[ "check1" ].param2 == 20.5 # Updated -@patch( 'geos.mesh.doctor.parsing._sharedChecksParsingLogic.setupLogger' ) -def test_convertWithInvalidCheck( mockLogger, checkFeaturesConfig, orderedCheckNames ): +@patch( 'geos.mesh_doctor.parsing._sharedChecksParsingLogic.setupLogger' ) +def test_convertWithInvalidCheck( mockLogger: MagicMock, checkFeaturesConfig: dict[ str, CheckFeature ], + orderedCheckNames: list[ str ] ) -> None: + """Tests convert when an invalid check is specified.""" parsedArgs = { CHECKS_TO_DO_ARG: "invalid_check,check1", PARAMETERS_ARG: "" } options = convert( parsedArgs, orderedCheckNames, checkFeaturesConfig ) # The invalid check should be skipped @@ -135,16 +157,20 @@ def test_convertWithInvalidCheck( mockLogger, checkFeaturesConfig, orderedCheckN assert "invalid_check" not in options.checksOptions -@patch( 'geos.mesh.doctor.parsing._sharedChecksParsingLogic.setupLogger' ) -def test_convertWithAllInvalidChecks( mockLogger, checkFeaturesConfig, orderedCheckNames ): +@patch( 'geos.mesh_doctor.parsing._sharedChecksParsingLogic.setupLogger' ) +def test_convertWithAllInvalidChecks( mockLogger: MagicMock, checkFeaturesConfig: dict[ str, CheckFeature ], + orderedCheckNames: list[ str ] ) -> None: + """Tests convert when all checks are invalid.""" parsedArgs = { CHECKS_TO_DO_ARG: "invalid_check1,invalid_check2", PARAMETERS_ARG: "" } # Should raise ValueError since no valid checks were selected with pytest.raises( ValueError, match="No valid checks were selected" ): convert( parsedArgs, orderedCheckNames, checkFeaturesConfig ) -@patch( 'geos.mesh.doctor.parsing._sharedChecksParsingLogic.setupLogger' ) -def test_displayResultsWithChecks( mockLogger, checkFeaturesConfig, orderedCheckNames ): +@patch( 'geos.mesh_doctor.parsing._sharedChecksParsingLogic.setupLogger' ) +def test_displayResultsWithChecks( mockLogger: MagicMock, checkFeaturesConfig: dict[ str, CheckFeature ], + orderedCheckNames: list[ str ] ) -> None: + """Tests displayResults when checks were performed.""" options = AllChecksOptions( checksToPerform=[ "check1", "check2" ], checksOptions={ "check1": MockOptions(), @@ -163,8 +189,9 @@ def test_displayResultsWithChecks( mockLogger, checkFeaturesConfig, orderedCheck assert mockLogger.results.call_count >= 2 -@patch( 'geos.mesh.doctor.parsing._sharedChecksParsingLogic.setupLogger' ) -def test_displayResultsNoChecks( mockLogger ): +@patch( 'geos.mesh_doctor.parsing._sharedChecksParsingLogic.setupLogger' ) +def test_displayResultsNoChecks( mockLogger: MagicMock ) -> None: + """Tests displayResults when no checks were performed.""" options = AllChecksOptions( checksToPerform=[], checksOptions={}, checkDisplays={} ) result = AllChecksResult( checkResults={} ) displayResults( options, result ) diff --git a/geos-mesh/tests/test_supportedElements.py b/mesh-doctor/tests/test_supportedElements.py similarity index 63% rename from geos-mesh/tests/test_supportedElements.py rename to mesh-doctor/tests/test_supportedElements.py index 4e8ecc28..f55f9e80 100644 --- a/geos-mesh/tests/test_supportedElements.py +++ b/mesh-doctor/tests/test_supportedElements.py @@ -1,33 +1,36 @@ -# import os -import pytest +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +import os from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, VTK_POLYHEDRON -# from geos.mesh.doctor.actions.supportedElements import Options, action, __action -from geos.mesh.doctor.actions.vtkPolyhedron import parseFaceStream, FaceStream from geos.mesh.utils.genericHelpers import toVtkIdList +from geos.mesh_doctor.actions.supportedElements import Options, action, meshAction +from geos.mesh_doctor.actions.vtkPolyhedron import FaceStream, parseFaceStream +dataRoot: str = os.path.join( os.path.dirname( os.path.abspath( __file__ ) ), "data" ) +supportElementsFile: str = os.path.join( dataRoot, "supportedElements.vtu" ) -# TODO Update this test to have access to another meshTests file -@pytest.mark.parametrize( "baseName", ( "supportedElements.vtk", "supportedElementsAsVTKPolyhedra.vtk" ) ) -def test_supportedElements( baseName: str ) -> None: + +def test_supportedElements() -> None: """Testing that the supported elements are properly detected as supported! Args: baseName (str): Supported elements are provided as standard elements or polyhedron elements. """ - ... - # directory = os.path.dirname( os.path.realpath( __file__ ) ) - # supportedElementsFileName = os.path.join( directory, "../../../../unitTests/meshTests", baseName ) - # options = Options( chunkSize=1, numProc=4 ) - # result = check( supportedElementsFileName, options ) - # assert not result.unsupportedStdElementsTypes - # assert not result.unsupportedPolyhedronElements + options = Options( chunkSize=1, nproc=4 ) + result = action( supportElementsFile, options ) + assert not result.unsupportedStdElementsTypes + assert not result.unsupportedPolyhedronElements def makeDodecahedron() -> tuple[ vtkPoints, vtkIdList ]: """Returns the points and faces for a dodecahedron. + This code was adapted from an official vtk example. - :return: The tuple of points and faces (as vtk instances). + + Returns: + The tuple of points and faces (as vtk instances). """ # yapf: disable points = ( @@ -78,9 +81,11 @@ def makeDodecahedron() -> tuple[ vtkPoints, vtkIdList ]: return p, f -# TODO make this test work def test_dodecahedron() -> None: - """Tests whether a dodecahedron is support by GEOS or not. + """Tests whether a dodecahedron is supported by GEOS or not. + + A dodecahedron has 12 pentagonal faces and is not supported by GEOS, + which only supports hexahedra, tetrahedra, pyramids, wedges, and polygons. """ points, faces = makeDodecahedron() mesh = vtkUnstructuredGrid() @@ -88,13 +93,16 @@ def test_dodecahedron() -> None: mesh.SetPoints( points ) mesh.InsertNextCell( VTK_POLYHEDRON, faces ) - # TODO Why does __check triggers an assertion error with 'assert MESH is not None' ? - # result = __check( mesh, Options( num_proc=1, chunk_size=1 ) ) - # assert set( result.unsupported_polyhedron_elements ) == { 0 } - # assert not result.unsupported_std_elements_types + # Test using meshAction directly instead of the internal __check function + options = Options( nproc=1, chunkSize=1 ) + result = meshAction( mesh, options ) + # Dodecahedron (12 pentagonal faces) is not supported by GEOS + assert set( result.unsupportedPolyhedronElements ) == { 0 } + assert not result.unsupportedStdElementsTypes def test_parseFaceStream() -> None: + """Tests the parsing of a face stream for a dodecahedron.""" _, faces = makeDodecahedron() result = parseFaceStream( faces ) # yapf: disable @@ -114,6 +122,6 @@ def test_parseFaceStream() -> None: ) # yapf: enable assert result == expected - face_stream = FaceStream.buildFromVtkIdList( faces ) - assert face_stream.numFaces == 12 - assert face_stream.numSupportPoints == 20 + faceStream = FaceStream.buildFromVtkIdList( faces ) + assert faceStream.numFaces == 12 + assert faceStream.numSupportPoints == 20 diff --git a/geos-mesh/tests/test_triangleDistance.py b/mesh-doctor/tests/test_triangleDistance.py similarity index 73% rename from geos-mesh/tests/test_triangleDistance.py rename to mesh-doctor/tests/test_triangleDistance.py index 326585b6..0938f04f 100644 --- a/geos-mesh/tests/test_triangleDistance.py +++ b/mesh-doctor/tests/test_triangleDistance.py @@ -1,26 +1,35 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto from dataclasses import dataclass import numpy +import numpy.typing as npt from numpy.linalg import norm import pytest -from geos.mesh.doctor.actions.triangleDistance import distanceBetweenTwoSegments, distanceBetweenTwoTriangles +from typing import Iterable +from geos.mesh_doctor.actions.triangleDistance import distanceBetweenTwoSegments, distanceBetweenTwoTriangles @dataclass( frozen=True ) class ExpectedSeg: - p0: numpy.array - u0: numpy.array - p1: numpy.array - u1: numpy.array - x: numpy.array - y: numpy.array + p0: npt.NDArray[ numpy.float64 ] + u0: npt.NDArray[ numpy.float64 ] + p1: npt.NDArray[ numpy.float64 ] + u1: npt.NDArray[ numpy.float64 ] + x: npt.NDArray[ numpy.float64 ] + y: npt.NDArray[ numpy.float64 ] @classmethod - def fromTuples( cls, p0, u0, p1, u1, x, y ): + def fromTuples( cls, p0: tuple[ float, float, float ], u0: tuple[ float, float, float ], + p1: tuple[ float, float, float ], u1: tuple[ float, float, float ], x: tuple[ float, float, float ], + y: tuple[ float, float, float ] ) -> "ExpectedSeg": + """Creates an ExpectedSeg from tuples.""" return cls( numpy.array( p0 ), numpy.array( u0 ), numpy.array( p1 ), numpy.array( u1 ), numpy.array( x ), numpy.array( y ) ) -def __getSegmentsReferences(): +def __getSegmentsReferences() -> Iterable[ ExpectedSeg ]: + """Provides reference segments for testing.""" # Node to node configuration. yield ExpectedSeg.fromTuples( p0=( 0., 0., 0. ), @@ -80,7 +89,12 @@ def __getSegmentsReferences(): @pytest.mark.parametrize( "expected", __getSegmentsReferences() ) -def test_segments( expected: ExpectedSeg ): +def test_segments( expected: ExpectedSeg ) -> None: + """Tests the distance between two segments. + + Args: + expected (ExpectedSeg): The expected segment data. + """ eps = numpy.finfo( float ).eps x, y = distanceBetweenTwoSegments( expected.p0, expected.u0, expected.p1, expected.u1 ) if norm( expected.x - expected.y ) == 0: @@ -92,18 +106,21 @@ def test_segments( expected: ExpectedSeg ): @dataclass( frozen=True ) class ExpectedTri: - t0: numpy.array - t1: numpy.array + t0: npt.NDArray[ numpy.float64 ] + t1: npt.NDArray[ numpy.float64 ] d: float - p0: numpy.array - p1: numpy.array + p0: npt.NDArray[ numpy.float64 ] + p1: npt.NDArray[ numpy.float64 ] @classmethod - def fromTuples( cls, t0, t1, d, p0, p1 ): + def fromTuples( cls, t0: tuple[ tuple[ float, float, float ], ...], t1: tuple[ tuple[ float, float, float ], ...], + d: float, p0: tuple[ float, float, float ], p1: tuple[ float, float, float ] ) -> "ExpectedTri": + """Creates an ExpectedTri from tuples.""" return cls( numpy.array( t0 ), numpy.array( t1 ), float( d ), numpy.array( p0 ), numpy.array( p1 ) ) -def __getTrianglesReferences(): +def __getTrianglesReferences() -> Iterable[ ExpectedTri ]: + """Provides reference triangles for testing.""" # Node to node configuration. yield ExpectedTri.fromTuples( t0=( ( 0., 0., 0. ), ( 1., 0., 0. ), ( 0., 1., 1. ) ), t1=( ( 2., 0., 0. ), ( 3., 0., 0. ), ( 2., 1., 1. ) ), @@ -143,7 +160,12 @@ def __getTrianglesReferences(): @pytest.mark.parametrize( "expected", __getTrianglesReferences() ) -def test_triangles( expected: ExpectedTri ): +def test_triangles( expected: ExpectedTri ) -> None: + """Tests the distance between two triangles. + + Args: + expected (ExpectedTri): The expected triangle data. + """ eps = numpy.finfo( float ).eps d, p0, p1 = distanceBetweenTwoTriangles( expected.t0, expected.t1 ) assert abs( d - expected.d ) < eps