From c654c23b964d594aa2eb2910165697362b031663 Mon Sep 17 00:00:00 2001 From: alexbenedicto Date: Mon, 16 Sep 2024 17:49:38 -0700 Subject: [PATCH 1/5] New branch up to date with previous PR 39 before having merge errors --- docs/geos-mesh.rst | 220 ++++- .../doctor/checks/fix_elements_orderings.py | 200 ++++- .../src/geos/mesh/doctor/checks/vtk_utils.py | 54 +- geos-mesh/src/geos/mesh/doctor/mesh_doctor.py | 12 +- .../parsing/fix_elements_orderings_parsing.py | 90 +- .../tests/test_fix_elements_orderings.py | 814 ++++++++++++++++++ 6 files changed, 1257 insertions(+), 133 deletions(-) create mode 100644 geos-mesh/tests/test_fix_elements_orderings.py diff --git a/docs/geos-mesh.rst b/docs/geos-mesh.rst index 8582f106..739d4936 100644 --- a/docs/geos-mesh.rst +++ b/docs/geos-mesh.rst @@ -15,14 +15,59 @@ 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:: python src/geos/mesh/doctor/mesh_doctor.py --help - :cwd: ../geos-mesh +.. code-block:: + + $ python src/geos/mesh/doctor/mesh_doctor.py --help + usage: mesh_doctor.py [-h] [-v] [-q] -i VTK_MESH_FILE + {collocated_nodes,element_volumes,fix_elements_orderings,generate_cube,generate_fractures,generate_global_ids,non_conformal,self_intersecting_elements,supported_elements} + ... + + Inspects meshes for GEOSX. + + positional arguments: + {collocated_nodes,element_volumes,fix_elements_orderings,generate_cube,generate_fractures,generate_global_ids,non_conformal,self_intersecting_elements,supported_elements} + Modules + collocated_nodes + Checks if nodes are collocated. + element_volumes + Checks if the volumes of the elements are greater than "min". + fix_elements_orderings + Reorders the support nodes for the given cell types. + generate_cube + Generate a cube and its fields. + generate_fractures + Splits the mesh to generate the faults and fractures. [EXPERIMENTAL] + generate_global_ids + Adds globals ids for points and cells. + non_conformal + Detects non conformal elements. [EXPERIMENTAL] + self_intersecting_elements + Checks if the faces of the elements are self intersecting. + supported_elements + 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 -.. command-output:: python src/geos/mesh/doctor/mesh_doctor.py collocated_nodes --help - :cwd: ../geos-mesh +.. code-block:: + + $ python src/geos/mesh/doctor/mesh_doctor.py collocated_nodes --help + usage: mesh_doctor.py collocated_nodes [-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. @@ -44,8 +89,14 @@ Here is a list and brief description of all the modules available. 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:: python src/geos/mesh/doctor/mesh_doctor.py collocated_nodes --help - :cwd: ../geos-mesh +.. code-block:: + + $ python src/geos/mesh/doctor/mesh_doctor.py collocated_nodes --help + usage: mesh_doctor.py collocated_nodes [-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. ``element_volumes`` """"""""""""""""""" @@ -53,8 +104,14 @@ It is not uncommon to define multiple nodes for the exact same position, which w 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:: python src/geos/mesh/doctor/mesh_doctor.py element_volumes --help - :cwd: ../geos-mesh +.. code-block:: + + $ python src/geos/mesh/doctor/mesh_doctor.py element_volumes --help + usage: mesh_doctor.py element_volumes [-h] --min 0.0 + + options: + -h, --help show this help message and exit + --min 0.0 [float]: The minimum acceptable volume. Defaults to 0.0. ``fix_elements_orderings`` """""""""""""""""""""""""" @@ -63,8 +120,38 @@ It sometimes happens that an exported mesh does not abide by the ``vtk`` orderin The ``fix_elements_orderings`` module can rearrange the nodes of given types of elements. This can be convenient if you cannot regenerate the mesh. -.. command-output:: python src/geos/mesh/doctor/mesh_doctor.py fix_elements_orderings --help - :cwd: ../geos-mesh +.. code-block:: + + $ python src/geos/mesh/doctor/mesh_doctor.py fix_elements_orderings --help + usage: mesh_doctor.py fix_elements_orderings [-h] [--Tetrahedron 2,0,3,1] [--Pyramid 3,4,0,2,1] [--Wedge 3,5,4,0,2,1] + [--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] + --volume_to_reorder all,positive,negative --output OUTPUT [--data-mode binary, ascii] + + options: + -h, --help show this help message and exit + --Tetrahedron 2,0,3,1 [list of integers]: node permutation for "Tetrahedron". + --Pyramid 3,4,0,2,1 [list of integers]: node permutation for "Pyramid". + --Wedge 3,5,4,0,2,1 [list of integers]: node permutation for "Wedge". + --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". + --volume_to_reorder all,positive,negative + [str]: Select which element volume is invalid and needs reordering. 'all' will allow reordering of nodes for every element, regarding of their volume. + 'positive' or 'negative' will only reorder the element with the corresponding volume. + --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. + +For example, assume that you have a mesh that contains 3 different cell types like Hexahedrons, Pyramids and Tetrahedrons. +After checking ``element_volumes``, you found that all your Hexahedrons and half of your Tetrahedrons have a negative volume. +To correct that, assuming you know the correct node ordering to correct these volumes, you can use this command: + +.. code-block:: + + $ python src/geos/mesh/doctor/mesh_doctor.py -i /path/to/your/mesh/file fix_elements_orderings --Hexahedron 1,6,5,4,7,0,2,3 + --Tetrahedron 2,0,3,1 --volume_to_reorder negative --output /new/path/for/your/new/mesh/reordered/file --data-mode binary ``generate_cube`` """"""""""""""""" @@ -73,8 +160,30 @@ 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:: python src/geos/mesh/doctor/mesh_doctor.py generate_cube --help - :cwd: ../geos-mesh +.. code-block:: + + $ python src/geos/mesh/doctor/mesh_doctor.py generate_cube --help + usage: mesh_doctor.py generate_cube [-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. ``generate_fractures`` """""""""""""""""""""" @@ -82,8 +191,31 @@ This tool can also be useful to generate a trial mesh that will later be refined 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 ``generate_fractures`` module will split the mesh and generate the multi-block ``vtk`` files. -.. command-output:: python src/geos/mesh/doctor/mesh_doctor.py generate_fractures --help - :cwd: ../geos-mesh +.. code-block:: + + $ python src/geos/mesh/doctor/mesh_doctor.py generate_fractures --help + usage: mesh_doctor.py generate_fractures [-h] --policy field, internal_surfaces [--name NAME] [--values VALUES] + --output OUTPUT [--data-mode binary, ascii] --fracture-output + FRACTURE_OUTPUT [--fracture-data-mode binary, ascii] + + options: + -h, --help show this help message and exit + --policy field, internal_surfaces + [string]: The criterion to define the surfaces that will be changed into fracture zones. + Possible values are "field, internal_surfaces" + --name NAME [string]: If the "field" policy is selected, defines which field will be considered to + define the fractures. If the "internal_surfaces" policy is selected, defines the name of + the attribute will be considered to identify the fractures. + --values VALUES [list of comma separated integers]: If the "field" policy is selected, which changes of + the field will be considered as a fracture. If the "internal_surfaces" policy is + selected, list of the fracture attributes. + --output OUTPUT [string]: The vtk output file destination. + --data-mode binary, ascii + [string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary. + --fracture-output FRACTURE_OUTPUT + [string]: The vtk output file destination. + --fracture-data-mode binary, ascii + [string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary. ``generate_global_ids`` """"""""""""""""""""""" @@ -91,8 +223,21 @@ The ``generate_fractures`` module will split the mesh and generate the multi-blo When running ``geos`` in parallel, `global ids` can be used to refer to data across multiple ranks. The ``generate_global_ids`` can generate `global ids` for the imported ``vtk`` mesh. -.. command-output:: python src/geos/mesh/doctor/mesh_doctor.py generate_global_ids --help - :cwd: ../geos-mesh +.. code-block:: + + $ python src/geos/mesh/doctor/mesh_doctor.py generate_global_ids --help + usage: mesh_doctor.py generate_global_ids [-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. ``non_conformal`` """"""""""""""""" @@ -102,8 +247,19 @@ This module will detect elements which are close enough (there's a user defined The angle between two faces can also be precribed. This module can be a bit time consuming. -.. command-output:: python src/geos/mesh/doctor/mesh_doctor.py non_conformal --help - :cwd: ../geos-mesh +.. code-block:: + + $ python src/geos/mesh/doctor/mesh_doctor.py non_conformal --help + usage: mesh_doctor.py non_conformal [-h] [--angle_tolerance 10.0] [--point_tolerance POINT_TOLERANCE] + [--face_tolerance FACE_TOLERANCE] + + options: + -h, --help show this help message and exit + --angle_tolerance 10.0 [float]: angle tolerance in degrees. Defaults to 10.0 + --point_tolerance POINT_TOLERANCE + [float]: tolerance for two points to be considered collocated. + --face_tolerance FACE_TOLERANCE + [float]: tolerance for two faces to be considered "touching". ``self_intersecting_elements`` """""""""""""""""""""""""""""" @@ -111,8 +267,15 @@ This module can be a bit time consuming. Some meshes can have cells that auto-intersect. This module will display the elements that have faces intersecting. -.. command-output:: python src/geos/mesh/doctor/mesh_doctor.py self_intersecting_elements --help - :cwd: ../geos-mesh +.. code-block:: + + $ python src/geos/mesh/doctor/mesh_doctor.py self_intersecting_elements --help + usage: mesh_doctor.py self_intersecting_elements [-h] [--min 2.220446049250313e-16] + + options: + -h, --help show this help message and exit + --min 2.220446049250313e-16 + [float]: The tolerance in the computation. Defaults to your machine precision 2.220446049250313e-16. ``supported_elements`` """""""""""""""""""""" @@ -125,8 +288,15 @@ But also prismes up to 11 faces. The ``supported_elements`` 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:: python src/geos/mesh/doctor/mesh_doctor.py supported_elements --help - :cwd: ../geos-mesh +.. code-block:: + + $ python src/geos/mesh/doctor/mesh_doctor.py supported_elements --help + usage: mesh_doctor.py supported_elements [-h] [--chunck_size 1] [--nproc 8] + + options: + -h, --help show this help message and exit + --chunck_size 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. @@ -179,8 +349,4 @@ API ^^^ .. automodule:: geos.mesh.conversion.abaqus_converter - :members: - - - - + :members: \ No newline at end of file diff --git a/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py b/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py index eb603b4e..bc693a99 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py @@ -1,65 +1,173 @@ from dataclasses import dataclass +import numpy as np import logging -from typing import ( - List, - Dict, - Set, - FrozenSet, -) - -from vtkmodules.vtkCommonCore import ( - vtkIdList, ) - -from . import vtk_utils -from .vtk_utils import ( - to_vtk_id_list, - VtkOutput, -) +from vtk import vtkCellSizeFilter +from vtkmodules.vtkCommonCore import vtkIdList +from vtkmodules.util.numpy_support import vtk_to_numpy +from vtkmodules.vtkCommonDataModel import ( vtkDataSet, VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, + VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM ) +from .vtk_utils import VtkOutput, to_vtk_id_list, write_mesh, read_mesh @dataclass( frozen=True ) class Options: vtk_output: VtkOutput - cell_type_to_ordering: Dict[ int, List[ int ] ] + cell_name_to_ordering: dict[ str, list[ int ] ] + volume_to_reorder: str @dataclass( frozen=True ) class Result: output: str - unchanged_cell_types: FrozenSet[ int ] + reordering_stats: dict[ str, list[ int ] ] + + +GEOS_ACCEPTED_TYPES = [ VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM ] +# the number of different nodes that needs to be entered in parsing when dealing with a specific vtk element +NAME_TO_VTK_TYPE = { + "Hexahedron": VTK_HEXAHEDRON, + "Tetrahedron": VTK_TETRA, + "Pyramid": VTK_PYRAMID, + "Wedge": VTK_WEDGE, + "Prism5": VTK_PENTAGONAL_PRISM, + "Prism6": VTK_HEXAGONAL_PRISM +} +VTK_TYPE_TO_NAME = { vtk_type: name for name, vtk_type in NAME_TO_VTK_TYPE.items() } + + +# Knowing the calculation of cell volumes in vtk was discussed there: https://github.com/GEOS-DEV/GEOS/issues/2253 +# Here, we do not need to have the exact volumes matching the simulation softwares results +# because we are mostly interested in knowing the sign of the volume for the rest of the workflow. +# Therefore, there is no need to use vtkMeshQuality for specific cell types when vtkCellSizeFilter works with all types. +def compute_mesh_cells_volume( mesh: vtkDataSet ) -> np.array: + """Generates a volume array that was calculated on all cells of a mesh. + + Args: + mesh (vtkDataSet): A vtk grid. + + Returns: + np.array: Volume for every cell of a mesh. + """ + cell_size_filter = vtkCellSizeFilter() + cell_size_filter.SetInputData( mesh ) + cell_size_filter.SetComputeVolume( True ) + cell_size_filter.Update() + return vtk_to_numpy( cell_size_filter.GetOutput().GetCellData().GetArray( "Volume" ) ) + + +def get_cell_types_and_number( mesh: vtkDataSet ) -> tuple[ list[ int ] ]: + """Gets the cell type for every cell of a mesh and the amount for each cell type. + + Args: + mesh (vtkDataSet): A vtk grid. + + Raises: + ValueError: "Invalid type '{cell_type}' for GEOS is in the mesh. Dying ..." + + Returns: + tuple[ list[ int ] ]: ( unique_cell_types, total_per_cell_types ) + """ + number_cells: int = mesh.GetNumberOfCells() + all_cells_type: np.array = np.ones( number_cells, dtype=int ) + for cell_id in range( number_cells ): + all_cells_type[ cell_id ] = mesh.GetCellType( cell_id ) + unique_cell_types, total_per_cell_types = np.unique( all_cells_type, return_counts=True ) + unique_cell_types, total_per_cell_types = unique_cell_types.tolist(), total_per_cell_types.tolist() + for cell_type in unique_cell_types: + if cell_type not in GEOS_ACCEPTED_TYPES: + err_msg: str = f"Invalid type '{cell_type}' for GEOS is in the mesh. Dying ..." + logging.error( err_msg ) + raise ValueError( err_msg ) + return ( unique_cell_types, total_per_cell_types ) + + +def is_cell_to_reorder( cell_volume: str, options: Options ) -> bool: + """Check if the volume of vtkCell qualifies the cell to be reordered. + + Args: + cell_volume (float): The volume of a vtkCell. + options (Options): Options defined by the user. + + Returns: + bool: True if cell needs to be reordered + """ + if options.volume_to_reorder == "all": + return True + if cell_volume == 0.0: + return True + sign_of_cell_volume: int = int( cell_volume / abs( cell_volume ) ) + if options.volume_to_reorder == "positive" and sign_of_cell_volume == 1: + return True + elif options.volume_to_reorder == "negative" and sign_of_cell_volume == -1: + return True + return False + + +def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple: + """From an old mesh, creates a new one where certain cell nodes are reordered to obtain a correct volume. + + Args: + old_mesh (vtkDataSet): A vtk grid needing nodes to be reordered. + options (Options): Options defined by the user. + + Returns: + tuple: ( vtkDataSet, reordering_stats ) + """ + unique_cell_types, total_per_cell_types = get_cell_types_and_number( old_mesh ) + unique_cell_names: list[ str ] = [ VTK_TYPE_TO_NAME[ vtk_type ] for vtk_type in unique_cell_types ] + names_with_totals: dict[ str, int ] = { n: v for n, v in zip( unique_cell_names, total_per_cell_types ) } + # sorted dict allow for sorted output of reordering_stats + names_with_totals_sorted: dict[ str, int ] = dict( sorted( names_with_totals.items() ) ) + useful_VTK_TYPEs: list[ int ] = [ NAME_TO_VTK_TYPE[ name ] for name in options.cell_name_to_ordering.keys() ] + all_cells_volume: np.array = compute_mesh_cells_volume( old_mesh ) + # a new mesh with the same data is created from the old mesh + new_mesh: vtkDataSet = old_mesh.NewInstance() + new_mesh.CopyStructure( old_mesh ) + new_mesh.CopyAttributes( old_mesh ) + cells = new_mesh.GetCells() + # Statistics on how many cells have or have not been reordered + reordering_stats: dict[ str, list[ any ] ] = { + "Types reordered": list(), + "Number of cells reordered": list(), + "Types non reordered": list( names_with_totals_sorted.keys() ), + "Number of cells non reordered": list( names_with_totals_sorted.values() ) + } + counter_cells_reordered: dict[ str, int ] = { name: 0 for name in options.cell_name_to_ordering.keys() } + # sorted dict allow for sorted output of reordering_stats + ounter_cells_reordered_sorted: dict[ str, int ] = dict( sorted( counter_cells_reordered.items() ) ) + # Reordering operations + for cell_id in range( new_mesh.GetNumberOfCells() ): + vtk_type: int = new_mesh.GetCellType( cell_id ) + if vtk_type in useful_VTK_TYPEs: + if is_cell_to_reorder( float( all_cells_volume[ cell_id ] ), options ): + cell_name: str = VTK_TYPE_TO_NAME[ vtk_type ] + support_point_ids = vtkIdList() + cells.GetCellAtId( cell_id, support_point_ids ) + new_support_point_ids: list[ int ] = list() + node_ordering: list[ int ] = options.cell_name_to_ordering[ cell_name ] + for i in range( len( node_ordering ) ): + new_support_point_ids.append( support_point_ids.GetId( node_ordering[ i ] ) ) + cells.ReplaceCellAtId( cell_id, to_vtk_id_list( new_support_point_ids ) ) + ounter_cells_reordered_sorted[ cell_name ] += 1 + # Calculation of stats + for cell_name_reordered, amount in ounter_cells_reordered_sorted.items(): + if amount > 0: + reordering_stats[ "Types reordered" ].append( cell_name_reordered ) + reordering_stats[ "Number of cells reordered" ].append( amount ) + index_non_reordered: int = reordering_stats[ "Types non reordered" ].index( cell_name_reordered ) + reordering_stats[ "Number of cells non reordered" ][ index_non_reordered ] -= amount + if reordering_stats[ "Number of cells non reordered" ][ index_non_reordered ] == 0: + reordering_stats[ "Types non reordered" ].pop( index_non_reordered ) + reordering_stats[ "Number of cells non reordered" ].pop( index_non_reordered ) + return ( new_mesh, reordering_stats ) def __check( mesh, options: Options ) -> Result: - # The vtk cell type is an int and will be the key of the following mapping, - # that will point to the relevant permutation. - cell_type_to_ordering: Dict[ int, List[ int ] ] = options.cell_type_to_ordering - unchanged_cell_types: Set[ int ] = set() # For logging purpose - - # Preparing the output mesh by first keeping the same instance type. - output_mesh = mesh.NewInstance() - output_mesh.CopyStructure( mesh ) - output_mesh.CopyAttributes( mesh ) - - # `output_mesh` now contains a full copy of the input mesh. - # We'll now modify the support nodes orderings in place if needed. - cells = output_mesh.GetCells() - for cell_idx in range( output_mesh.GetNumberOfCells() ): - cell_type: int = output_mesh.GetCell( cell_idx ).GetCellType() - new_ordering = cell_type_to_ordering.get( cell_type ) - if new_ordering: - support_point_ids = vtkIdList() - cells.GetCellAtId( cell_idx, support_point_ids ) - new_support_point_ids = [] - for i, v in enumerate( new_ordering ): - new_support_point_ids.append( support_point_ids.GetId( new_ordering[ i ] ) ) - cells.ReplaceCellAtId( cell_idx, to_vtk_id_list( new_support_point_ids ) ) - else: - unchanged_cell_types.add( cell_type ) - is_written_error = vtk_utils.write_mesh( output_mesh, options.vtk_output ) - return Result( output=options.vtk_output.output if not is_written_error else "", - unchanged_cell_types=frozenset( unchanged_cell_types ) ) + output_mesh, reordering_stats = reorder_nodes_to_new_mesh( mesh, options ) + write_mesh( output_mesh, options.vtk_output ) + return Result( output=options.vtk_output.output, reordering_stats=reordering_stats ) def check( vtk_input_file: str, options: Options ) -> Result: - mesh = vtk_utils.read_mesh( vtk_input_file ) + mesh = read_mesh( vtk_input_file ) return __check( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py index 9beb3758..c27872ac 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py @@ -1,7 +1,7 @@ from dataclasses import dataclass -import os.path +from os import path, access, R_OK, W_OK import logging -import sys +from sys import exit from typing import ( Any, Iterator, @@ -50,6 +50,39 @@ def vtk_iter( l ) -> Iterator[ Any ]: yield l.GetCellType( i ) +# Inspired from https://stackoverflow.com/a/78870363 +def is_filepath_valid( filepath: str, is_input=True ) -> bool: + """Checks if a filepath can be used to read or to create a new file. + + Args: + filepath (str): A filepath. + is_input (bool, optional): If the filepath is used to read a file, use True. + If the filepath is used to create a new file, use False. Defaults to True. + + Returns: + bool: False if invalid, True instead. + """ + dirname = path.dirname( filepath ) + if not path.isdir( dirname ): + logging.error( f"The directory '{dirname}' specified does not exist." ) + return False + if is_input: + if not access( dirname, R_OK ): + logging.error( f"You do not have rights to read in directory '{dirname}' specified for the file." ) + return False + if not path.exists( filepath ): + logging.error( f"The file specified does not exist." ) + return False + else: + if not access( dirname, W_OK ): + logging.error( f"You do not have rights to write in directory '{dirname}' specified for the file." ) + return False + if path.exists( filepath ): + logging.error( f"A file with the same name already exists. No over-writing possible." ) + return False + return True + + def __read_vtk( vtk_input_file: str ) -> Optional[ vtkUnstructuredGrid ]: reader = vtkUnstructuredGridReader() logging.info( f"Testing file format \"{vtk_input_file}\" using legacy format reader..." ) @@ -83,7 +116,10 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid: If first guess does not work, eventually all the others reader available will be tested. :return: A unstructured grid. """ - file_extension = os.path.splitext( vtk_input_file )[ -1 ] + if not is_filepath_valid( vtk_input_file, True ): + logging.error( f"Invalid filepath to read. Dying ..." ) + exit( 1 ) + file_extension = path.splitext( vtk_input_file )[ -1 ] extension_to_reader = { ".vtk": __read_vtk, ".vtu": __read_vtu } # Testing first the reader that should match if file_extension in extension_to_reader: @@ -97,7 +133,7 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid: return output_mesh # No reader did work. Dying. logging.critical( f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\". Dying..." ) - sys.exit( 1 ) + exit( 1 ) def __write_vtk( mesh: vtkUnstructuredGrid, output: str ) -> int: @@ -125,10 +161,10 @@ def write_mesh( mesh: vtkUnstructuredGrid, vtk_output: VtkOutput ) -> int: :param vtk_output: Where to write. The file extension will be used to select the VTK file format. :return: 0 in case of success. """ - if os.path.exists( vtk_output.output ): - logging.error( f"File \"{vtk_output.output}\" already exists, nothing done." ) - return 1 - file_extension = os.path.splitext( vtk_output.output )[ -1 ] + if not is_filepath_valid( vtk_output.output, False ): + logging.error( f"Invalid filepath to write. Dying ..." ) + exit( 1 ) + file_extension = path.splitext( vtk_output.output )[ -1 ] if file_extension == ".vtk": success_code = __write_vtk( mesh, vtk_output.output ) elif file_extension == ".vtu": @@ -136,5 +172,5 @@ def write_mesh( mesh: vtkUnstructuredGrid, vtk_output: VtkOutput ) -> int: else: # No writer found did work. Dying. logging.critical( f"Could not find the appropriate VTK writer for extension \"{file_extension}\". Dying..." ) - sys.exit( 1 ) + exit( 1 ) return 0 if success_code else 2 # the Write member function return 1 in case of success, 0 otherwise. diff --git a/geos-mesh/src/geos/mesh/doctor/mesh_doctor.py b/geos-mesh/src/geos/mesh/doctor/mesh_doctor.py index ea1bfe8a..1311145f 100644 --- a/geos-mesh/src/geos/mesh/doctor/mesh_doctor.py +++ b/geos-mesh/src/geos/mesh/doctor/mesh_doctor.py @@ -1,17 +1,17 @@ import sys +import logging +from geos.mesh.doctor.parsing import CheckHelper +from geos.mesh.doctor.parsing.cli_parsing import parse_and_set_verbosity +import geos.mesh.doctor.register as register +min_python_version = ( 3, 7 ) try: - min_python_version = ( 3, 7 ) assert sys.version_info >= min_python_version except AssertionError as e: print( f"Please update python to at least version {'.'.join(map(str, min_python_version))}." ) sys.exit( 1 ) -import logging - -from geos.mesh.doctor.parsing import CheckHelper -from geos.mesh.doctor.parsing.cli_parsing import parse_and_set_verbosity -import geos.mesh.doctor.register as register +MESH_DOCTOR_FILEPATH = __file__ def main(): diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py index 71fb3a51..641a9505 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py +++ b/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py @@ -1,52 +1,41 @@ import logging import random - -from vtkmodules.vtkCommonDataModel import ( - VTK_HEXAGONAL_PRISM, - VTK_HEXAHEDRON, - VTK_PENTAGONAL_PRISM, - VTK_PYRAMID, - VTK_TETRA, - VTK_VOXEL, - VTK_WEDGE, -) - from geos.mesh.doctor.checks.fix_elements_orderings import Options, Result - from . import vtk_output_parsing, FIX_ELEMENTS_ORDERINGS -__CELL_TYPE_MAPPING = { - "Hexahedron": VTK_HEXAHEDRON, - "Prism5": VTK_PENTAGONAL_PRISM, - "Prism6": VTK_HEXAGONAL_PRISM, - "Pyramid": VTK_PYRAMID, - "Tetrahedron": VTK_TETRA, - "Voxel": VTK_VOXEL, - "Wedge": VTK_WEDGE, +__CELL_NAME_WITH_NUMBER_NODES = { + "Tetrahedron": 4, + "Pyramid": 5, + "Wedge": 6, + "Hexahedron": 8, + "Prism5": 10, + "Prism6": 12 } -__CELL_TYPE_SUPPORT_SIZE = { - VTK_HEXAHEDRON: 8, - VTK_PENTAGONAL_PRISM: 10, - VTK_HEXAGONAL_PRISM: 12, - VTK_PYRAMID: 5, - VTK_TETRA: 4, - VTK_VOXEL: 8, - VTK_WEDGE: 6, -} +__VOLUME_TO_REORDER = "volume_to_reorder" +__VOLUME_TO_REORDER_DEFAULT = "all" +__VOLUME_TO_REORDER_CHOICES = [ "all", "positive", "negative" ] def fill_subparser( subparsers ) -> None: p = subparsers.add_parser( FIX_ELEMENTS_ORDERINGS, help="Reorders the support nodes for the given cell types." ) - for key, vtk_key in __CELL_TYPE_MAPPING.items(): - tmp = list( range( __CELL_TYPE_SUPPORT_SIZE[ vtk_key ] ) ) + for element_name, size in __CELL_NAME_WITH_NUMBER_NODES.items(): + tmp = list( range( size ) ) random.Random( 4 ).shuffle( tmp ) - p.add_argument( '--' + key, + p.add_argument( '--' + element_name, type=str, metavar=",".join( map( str, tmp ) ), default=None, required=False, - help=f"[list of integers]: node permutation for \"{key}\"." ) + help=f"[list of integers]: node permutation for \"{element_name}\"." ) + p.add_argument( '--' + __VOLUME_TO_REORDER, + type=str, + default=__VOLUME_TO_REORDER_DEFAULT, + metavar=",".join( map( str, __VOLUME_TO_REORDER_CHOICES ) ), + required=True, + help="[str]: Select which element volume is invalid and needs reordering." + + " 'all' will allow reordering of nodes for every element, regarding of their volume." + + " 'positive' or 'negative' will only reorder the element with the corresponding volume." ) vtk_output_parsing.fill_vtk_output_subparser( p ) @@ -56,26 +45,37 @@ def convert( parsed_options ) -> Options: :param options_str: Parsed cli options. :return: Options instance. """ - cell_type_to_ordering = {} - for key, vtk_key in __CELL_TYPE_MAPPING.items(): - raw_mapping = parsed_options[ key ] + cell_name_to_ordering: dict[ str, list[ int ] ] = {} + for element_name, size in __CELL_NAME_WITH_NUMBER_NODES.items(): + raw_mapping = parsed_options[ element_name ] if raw_mapping: - tmp = tuple( map( int, raw_mapping.split( "," ) ) ) - if not set( tmp ) == set( range( __CELL_TYPE_SUPPORT_SIZE[ vtk_key ] ) ): - err_msg = f"Permutation {raw_mapping} for type {key} is not valid." + nodes_ordering = tuple( map( int, raw_mapping.split( "," ) ) ) + if not set( nodes_ordering ) == set( range( size ) ): + err_msg: str = f"Permutation {raw_mapping} for type {element_name} is not valid." logging.error( err_msg ) raise ValueError( err_msg ) - cell_type_to_ordering[ vtk_key ] = tmp + cell_name_to_ordering[ element_name ] = nodes_ordering vtk_output = vtk_output_parsing.convert( parsed_options ) - return Options( vtk_output=vtk_output, cell_type_to_ordering=cell_type_to_ordering ) + volume_to_reorder: str = parsed_options[ __VOLUME_TO_REORDER ] + if volume_to_reorder.lower() not in __VOLUME_TO_REORDER_CHOICES: + raise ValueError( f"Please use one of these options for --volume_to_reorder: {__VOLUME_TO_REORDER_CHOICES}." ) + return Options( vtk_output=vtk_output, + cell_name_to_ordering=cell_name_to_ordering, + volume_to_reorder=volume_to_reorder ) def display_results( options: Options, result: Result ): if result.output: logging.info( f"New mesh was written to file '{result.output}'" ) - if result.unchanged_cell_types: - logging.info( f"Those vtk types were not reordered: [{', '.join(map(str, result.unchanged_cell_types))}]." ) - else: - logging.info( "All the cells of the mesh were reordered." ) else: logging.info( "No output file was written." ) + logging.info( f"Number of cells reordered:" ) + logging.info( f"\tCellType\tNumber" ) + for i in range( len( result.reordering_stats[ "Types reordered" ] ) ): + logging.info( f"\t{result.reordering_stats[ 'Types reordered' ][ i ]}" + + f"\t\t{result.reordering_stats[ 'Number of cells reordered' ][ i ]}" ) + logging.info( f"Number of cells non reordered:" ) + logging.info( f"\tCellType\tNumber" ) + for i in range( len( result.reordering_stats[ "Types non reordered" ] ) ): + logging.info( f"\t{result.reordering_stats[ 'Types non reordered' ][ i ]}" + + f"\t\t{result.reordering_stats[ 'Number of cells non reordered' ][ i ]}" ) diff --git a/geos-mesh/tests/test_fix_elements_orderings.py b/geos-mesh/tests/test_fix_elements_orderings.py new file mode 100644 index 00000000..0f1ec40d --- /dev/null +++ b/geos-mesh/tests/test_fix_elements_orderings.py @@ -0,0 +1,814 @@ +import os +import re +import pytest +import logging +import subprocess +import numpy as np +from geos.mesh.doctor.mesh_doctor import MESH_DOCTOR_FILEPATH +from geos.mesh.doctor.checks import fix_elements_orderings as feo +from geos.mesh.doctor.checks.generate_cube import Options, __build +from geos.mesh.doctor.checks.vtk_utils import ( VtkOutput, to_vtk_id_list, write_mesh ) +from geos.mesh.doctor.checks.fix_elements_orderings import Options as opt, VTK_TYPE_TO_NAME +from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints +from vtkmodules.vtkCommonDataModel import ( vtkDataSet, vtkUnstructuredGrid, vtkCellArray, vtkHexahedron, vtkTetra, + vtkPyramid, vtkVoxel, vtkWedge, vtkPentagonalPrism, vtkHexagonalPrism, + VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_VOXEL, + VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM ) + + +def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int ] ): + """Utility function to reorder the nodes of one cell for test purposes. + + Args: + mesh (vtkDataSet): A vtk grid. + cell_id (int): Cell id to identify the cell which will be modified. + node_ordering (list[ int ]): Nodes id ordering to construct a cell. + """ + if mesh.GetCell( cell_id ).GetNumberOfPoints() != len( node_ordering ): + raise ValueError( f"The cell to reorder needs to have '{mesh.GetCell( cell_id ).GetNumberOfPoints()}'" + + "nodes in reordering." ) + cells = mesh.GetCells() + support_point_ids = vtkIdList() + cells.GetCellAtId( cell_id, support_point_ids ) + new_support_point_ids = [] + node_ordering: list[ int ] = node_ordering + for i in range( len( node_ordering ) ): + new_support_point_ids.append( support_point_ids.GetId( node_ordering[ i ] ) ) + cells.ReplaceCellAtId( cell_id, to_vtk_id_list( new_support_point_ids ) ) + + +""" +For creation of output test meshes +""" +current_file_path: str = __file__ +dir_name: str = os.path.dirname( current_file_path ) +filepath_non_ordered_mesh: str = os.path.join( dir_name, "to_reorder_mesh.vtu" ) +filepath_reordered_mesh: str = os.path.join( dir_name, "reordered_mesh.vtu" ) +test_file: VtkOutput = VtkOutput( filepath_non_ordered_mesh, True ) +""" +Dict used to apply false nodes orderings for test purposes +""" +to_change_order: dict[ int, list[ int ] ] = { + "Hexahedron": [ 0, 3, 2, 1, 4, 5, 6, 7 ], + "Tetrahedron": [ 0, 2, 1, 3 ], + "Pyramid": [ 0, 3, 2, 1, 4 ], + "Wedge": [ 0, 2, 1, 3, 4, 5 ], + "Prism5": [ 0, 4, 3, 2, 1, 5, 6, 7, 8, 9 ], + "Prism6": [ 0, 1, 4, 2, 3, 5, 11, 10, 9, 8, 7, 6 ] +} +to_change_order = dict( sorted( to_change_order.items() ) ) +""" +1 Hexahedron: no invalid ordering +""" +out: VtkOutput = VtkOutput( "test", False ) +options_one_hex: Options = Options( vtk_output=out, + generate_cells_global_ids=False, + generate_points_global_ids=False, + xs=np.array( [ 0.0, 1.0 ] ), + ys=np.array( [ 0.0, 1.0 ] ), + zs=np.array( [ 0.0, 1.0 ] ), + nxs=[ 1 ], + nys=[ 1 ], + nzs=[ 1 ], + fields=[] ) +one_hex: vtkDataSet = __build( options_one_hex ) +""" +4 Hexahedrons: no invalid ordering +""" +out: VtkOutput = VtkOutput( "test", False ) +options_hexahedrons_grid: Options = Options( vtk_output=out, + generate_cells_global_ids=False, + generate_points_global_ids=False, + xs=np.array( [ 0.0, 1.0, 2.0 ] ), + ys=np.array( [ 0.0, 1.0, 2.0 ] ), + zs=np.array( [ 0.0, 1.0 ] ), + nxs=[ 1, 1 ], + nys=[ 1, 1 ], + nzs=[ 1 ], + fields=[] ) +hexahedrons_grid: vtkDataSet = __build( options_hexahedrons_grid ) +""" +4 Hexahedrons: 2 Hexahedrons with invalid ordering +""" +hexahedrons_grid_invalid: vtkDataSet = __build( options_hexahedrons_grid ) +for i in range( 2 ): + reorder_cell_nodes( hexahedrons_grid_invalid, i * 2 + 1, to_change_order[ "Hexahedron" ] ) +""" +4 tetrahedrons +""" +points_tetras: vtkPoints = vtkPoints() +points_tetras_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), + ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ) ] +for point_tetra in points_tetras_coords: + points_tetras.InsertNextPoint( point_tetra ) + +tetra1: vtkTetra = vtkTetra() +tetra1.GetPointIds().SetId( 0, 3 ) +tetra1.GetPointIds().SetId( 1, 0 ) +tetra1.GetPointIds().SetId( 2, 1 ) +tetra1.GetPointIds().SetId( 3, 4 ) + +tetra2: vtkTetra = vtkTetra() +tetra2.GetPointIds().SetId( 0, 6 ) +tetra2.GetPointIds().SetId( 1, 5 ) +tetra2.GetPointIds().SetId( 2, 4 ) +tetra2.GetPointIds().SetId( 3, 1 ) + +tetra3: vtkTetra = vtkTetra() +tetra3.GetPointIds().SetId( 0, 1 ) +tetra3.GetPointIds().SetId( 1, 2 ) +tetra3.GetPointIds().SetId( 2, 3 ) +tetra3.GetPointIds().SetId( 3, 6 ) + +tetra4: vtkTetra = vtkTetra() +tetra4.GetPointIds().SetId( 0, 4 ) +tetra4.GetPointIds().SetId( 1, 7 ) +tetra4.GetPointIds().SetId( 2, 6 ) +tetra4.GetPointIds().SetId( 3, 3 ) + +tetras_cells: vtkCellArray = vtkCellArray() +tetras_cells.InsertNextCell( tetra1 ) +tetras_cells.InsertNextCell( tetra2 ) +tetras_cells.InsertNextCell( tetra3 ) +tetras_cells.InsertNextCell( tetra4 ) + +tetras_grid: vtkUnstructuredGrid = vtkUnstructuredGrid() +tetras_grid.SetPoints( points_tetras ) +tetras_grid.SetCells( VTK_TETRA, tetras_cells ) + +# one of every other wedge has invalid ordering +tetras_grid_invalid = vtkUnstructuredGrid() +tetras_grid_invalid.DeepCopy( tetras_grid ) +for i in range( 2 ): + reorder_cell_nodes( tetras_grid_invalid, i * 2 + 1, to_change_order[ "Tetrahedron" ] ) +""" +4 pyramids +""" +points_pyramids: vtkPoints = vtkPoints() +points_pyramids_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), ( 0.5, 0.5, 1.0 ), ( 2.0, 0.0, 0.0 ), + ( 3.0, 0.0, 0.0 ), ( 3.0, 1.0, 0.0 ), ( 2.0, 1.0, 0.0 ), + ( 2.5, 0.5, 1.0 ), ( 0.0, 2.0, 0.0 ), ( 1.0, 2.0, 0.0 ), + ( 1.0, 3.0, 0.0 ), ( 0.0, 3.0, 0.0 ), ( 0.5, 2.5, 1.0 ), + ( 2.0, 2.0, 0.0 ), ( 3.0, 2.0, 0.0 ), ( 3.0, 3.0, 0.0 ), + ( 2.0, 3.0, 0.0 ), ( 2.5, 2.5, 1.0 ) ] +for point_pyramid in points_pyramids_coords: + points_pyramids.InsertNextPoint( point_pyramid ) + +pyramid1: vtkPyramid = vtkPyramid() +pyramid1.GetPointIds().SetId( 0, 0 ) +pyramid1.GetPointIds().SetId( 1, 1 ) +pyramid1.GetPointIds().SetId( 2, 2 ) +pyramid1.GetPointIds().SetId( 3, 3 ) +pyramid1.GetPointIds().SetId( 4, 4 ) + +pyramid2: vtkPyramid = vtkPyramid() +pyramid2.GetPointIds().SetId( 0, 5 ) +pyramid2.GetPointIds().SetId( 1, 6 ) +pyramid2.GetPointIds().SetId( 2, 7 ) +pyramid2.GetPointIds().SetId( 3, 8 ) +pyramid2.GetPointIds().SetId( 4, 9 ) + +pyramid3: vtkPyramid = vtkPyramid() +pyramid3.GetPointIds().SetId( 0, 10 ) +pyramid3.GetPointIds().SetId( 1, 11 ) +pyramid3.GetPointIds().SetId( 2, 12 ) +pyramid3.GetPointIds().SetId( 3, 13 ) +pyramid3.GetPointIds().SetId( 4, 14 ) + +pyramid4: vtkPyramid = vtkPyramid() +pyramid4.GetPointIds().SetId( 0, 15 ) +pyramid4.GetPointIds().SetId( 1, 16 ) +pyramid4.GetPointIds().SetId( 2, 17 ) +pyramid4.GetPointIds().SetId( 3, 18 ) +pyramid4.GetPointIds().SetId( 4, 19 ) + +pyramids_cells: vtkCellArray = vtkCellArray() +pyramids_cells.InsertNextCell( pyramid1 ) +pyramids_cells.InsertNextCell( pyramid2 ) +pyramids_cells.InsertNextCell( pyramid3 ) +pyramids_cells.InsertNextCell( pyramid4 ) + +pyramids_grid: vtkUnstructuredGrid = vtkUnstructuredGrid() +pyramids_grid.SetPoints( points_pyramids ) +pyramids_grid.SetCells( VTK_PYRAMID, pyramids_cells ) + +# one of every other wedge has invalid ordering +pyramids_grid_invalid = vtkUnstructuredGrid() +pyramids_grid_invalid.DeepCopy( pyramids_grid ) +for i in range( 2 ): + reorder_cell_nodes( pyramids_grid_invalid, i * 2 + 1, to_change_order[ "Pyramid" ] ) +""" +4 voxels: this type of element cannot be used in GEOS, we just test that the feature rejects them +""" +points_voxels: vtkPoints = vtkPoints() +points_voxels_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), + ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ), ( 2.0, 0.0, 0.0 ), + ( 3.0, 0.0, 0.0 ), ( 3.0, 1.0, 0.0 ), ( 2.0, 1.0, 0.0 ), + ( 2.0, 0.0, 1.0 ), ( 3.0, 0.0, 1.0 ), ( 3.0, 1.0, 1.0 ), + ( 2.0, 1.0, 1.0 ), ( 0.0, 2.0, 0.0 ), ( 1.0, 2.0, 0.0 ), + ( 1.0, 3.0, 0.0 ), ( 0.0, 3.0, 0.0 ), ( 0.0, 2.0, 1.0 ), + ( 1.0, 2.0, 1.0 ), ( 1.0, 3.0, 1.0 ), ( 0.0, 3.0, 1.0 ), + ( 2.0, 2.0, 0.0 ), ( 3.0, 2.0, 0.0 ), ( 3.0, 3.0, 0.0 ), + ( 2.0, 3.0, 0.0 ), ( 2.0, 2.0, 1.0 ), ( 3.0, 2.0, 1.0 ), + ( 3.0, 3.0, 1.0 ), ( 2.0, 3.0, 1.0 ) ] +for point_voxel in points_voxels_coords: + points_voxels.InsertNextPoint( point_voxel ) + +voxel1: vtkVoxel = vtkVoxel() +voxel1.GetPointIds().SetId( 0, 0 ) +voxel1.GetPointIds().SetId( 1, 1 ) +voxel1.GetPointIds().SetId( 2, 2 ) +voxel1.GetPointIds().SetId( 3, 3 ) +voxel1.GetPointIds().SetId( 4, 4 ) +voxel1.GetPointIds().SetId( 5, 5 ) +voxel1.GetPointIds().SetId( 6, 6 ) +voxel1.GetPointIds().SetId( 7, 7 ) + +voxel2: vtkVoxel = vtkVoxel() +voxel2.GetPointIds().SetId( 0, 8 ) +voxel2.GetPointIds().SetId( 1, 9 ) +voxel2.GetPointIds().SetId( 2, 10 ) +voxel2.GetPointIds().SetId( 3, 11 ) +voxel2.GetPointIds().SetId( 4, 12 ) +voxel2.GetPointIds().SetId( 5, 13 ) +voxel2.GetPointIds().SetId( 6, 14 ) +voxel2.GetPointIds().SetId( 7, 15 ) + +voxel3: vtkVoxel = vtkVoxel() +voxel3.GetPointIds().SetId( 0, 16 ) +voxel3.GetPointIds().SetId( 1, 17 ) +voxel3.GetPointIds().SetId( 2, 18 ) +voxel3.GetPointIds().SetId( 3, 19 ) +voxel3.GetPointIds().SetId( 4, 20 ) +voxel3.GetPointIds().SetId( 5, 21 ) +voxel3.GetPointIds().SetId( 6, 22 ) +voxel3.GetPointIds().SetId( 7, 23 ) + +voxel4: vtkVoxel = vtkVoxel() +voxel4.GetPointIds().SetId( 0, 24 ) +voxel4.GetPointIds().SetId( 1, 25 ) +voxel4.GetPointIds().SetId( 2, 26 ) +voxel4.GetPointIds().SetId( 3, 27 ) +voxel4.GetPointIds().SetId( 4, 28 ) +voxel4.GetPointIds().SetId( 5, 29 ) +voxel4.GetPointIds().SetId( 6, 30 ) +voxel4.GetPointIds().SetId( 7, 31 ) + +voxels_cells: vtkCellArray = vtkCellArray() +voxels_cells.InsertNextCell( voxel1 ) +voxels_cells.InsertNextCell( voxel2 ) +voxels_cells.InsertNextCell( voxel3 ) +voxels_cells.InsertNextCell( voxel4 ) + +voxels_grid: vtkUnstructuredGrid = vtkUnstructuredGrid() +voxels_grid.SetPoints( points_voxels ) +voxels_grid.SetCells( VTK_VOXEL, voxels_cells ) +""" +4 wedges +""" +points_wedges: vtkPoints = vtkPoints() +points_wedges_coords: list[ tuple[ float ] ] = [ ( 0.5, 0.0, 0.0 ), ( 1.5, 0.0, 0.0 ), ( 2.5, 0.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 2.0, 1.0, 0.0 ), + ( 0.5, 0.0, 1.0 ), ( 1.5, 0.0, 1.0 ), ( 2.5, 0.0, 1.0 ), + ( 0.0, 1.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( 2.0, 1.0, 1.0 ) ] +for point_wedge in points_wedges_coords: + points_wedges.InsertNextPoint( point_wedge ) + +wedge1: vtkWedge = vtkWedge() +wedge1.GetPointIds().SetId( 0, 9 ) +wedge1.GetPointIds().SetId( 1, 6 ) +wedge1.GetPointIds().SetId( 2, 10 ) +wedge1.GetPointIds().SetId( 3, 3 ) +wedge1.GetPointIds().SetId( 4, 0 ) +wedge1.GetPointIds().SetId( 5, 4 ) + +wedge2: vtkWedge = vtkWedge() +wedge2.GetPointIds().SetId( 0, 7 ) +wedge2.GetPointIds().SetId( 1, 10 ) +wedge2.GetPointIds().SetId( 2, 6 ) +wedge2.GetPointIds().SetId( 3, 1 ) +wedge2.GetPointIds().SetId( 4, 4 ) +wedge2.GetPointIds().SetId( 5, 0 ) + +wedge3: vtkWedge = vtkWedge() +wedge3.GetPointIds().SetId( 0, 10 ) +wedge3.GetPointIds().SetId( 1, 7 ) +wedge3.GetPointIds().SetId( 2, 11 ) +wedge3.GetPointIds().SetId( 3, 4 ) +wedge3.GetPointIds().SetId( 4, 1 ) +wedge3.GetPointIds().SetId( 5, 5 ) + +wedge4: vtkWedge = vtkWedge() +wedge4.GetPointIds().SetId( 0, 8 ) +wedge4.GetPointIds().SetId( 1, 11 ) +wedge4.GetPointIds().SetId( 2, 7 ) +wedge4.GetPointIds().SetId( 3, 2 ) +wedge4.GetPointIds().SetId( 4, 5 ) +wedge4.GetPointIds().SetId( 5, 1 ) + +wedges_cells: vtkCellArray = vtkCellArray() +wedges_cells.InsertNextCell( wedge1 ) +wedges_cells.InsertNextCell( wedge2 ) +wedges_cells.InsertNextCell( wedge3 ) +wedges_cells.InsertNextCell( wedge4 ) + +wedges_grid = vtkUnstructuredGrid() +wedges_grid.SetPoints( points_wedges ) +wedges_grid.SetCells( VTK_WEDGE, wedges_cells ) + +# one of every other wedge has invalid ordering +wedges_grid_invalid = vtkUnstructuredGrid() +wedges_grid_invalid.DeepCopy( wedges_grid ) +for i in range( 2 ): + reorder_cell_nodes( wedges_grid_invalid, i * 2 + 1, to_change_order[ "Wedge" ] ) +""" +4 pentagonal prisms +""" +points_penta_prisms: vtkPoints = vtkPoints() +points_penta_prisms_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.5, 0.5, 0.0 ), + ( 0.5, 1.0, 0.0 ), ( -0.5, 0.5, 0.0 ), ( 0.0, 0.0, 1.0 ), + ( 1.0, 0.0, 1.0 ), ( 1.5, 0.5, 1.0 ), ( 0.5, 1.0, 1.0 ), + ( -0.5, 0.5, 1.0 ), ( 2.0, 0.0, 0.0 ), ( 3.0, 0.0, 0.0 ), + ( 3.5, 0.5, 0.0 ), ( 2.5, 1.0, 0.0 ), ( 1.5, 0.5, 0.0 ), + ( 2.0, 0.0, 1.0 ), ( 3.0, 0.0, 1.0 ), ( 3.5, 0.5, 1.0 ), + ( 2.5, 1.0, 1.0 ), ( 1.5, 0.5, 1.0 ), ( 0.0, 2.0, 0.0 ), + ( 1.0, 2.0, 0.0 ), ( 1.5, 2.5, 0.0 ), ( 0.5, 3.0, 0.0 ), + ( -0.5, 2.5, 0.0 ), ( 0.0, 2.0, 1.0 ), ( 1.0, 2.0, 1.0 ), + ( 1.5, 2.5, 1.0 ), ( 0.5, 3.0, 1.0 ), ( -0.5, 2.5, 1.0 ), + ( 2.0, 2.0, 0.0 ), ( 3.0, 2.0, 0.0 ), ( 3.5, 2.5, 0.0 ), + ( 2.5, 3.0, 0.0 ), ( 1.5, 2.5, 0.0 ), ( 2.0, 2.0, 1.0 ), + ( 3.0, 2.0, 1.0 ), ( 3.5, 2.5, 1.0 ), ( 2.5, 3.0, 1.0 ), + ( 1.5, 2.5, 1.0 ) ] +for point_penta_prism in points_penta_prisms_coords: + points_penta_prisms.InsertNextPoint( point_penta_prism ) + +penta_prism1: vtkPentagonalPrism = vtkPentagonalPrism() +penta_prism1.GetPointIds().SetId( 0, 0 ) +penta_prism1.GetPointIds().SetId( 1, 1 ) +penta_prism1.GetPointIds().SetId( 2, 2 ) +penta_prism1.GetPointIds().SetId( 3, 3 ) +penta_prism1.GetPointIds().SetId( 4, 4 ) +penta_prism1.GetPointIds().SetId( 5, 5 ) +penta_prism1.GetPointIds().SetId( 6, 6 ) +penta_prism1.GetPointIds().SetId( 7, 7 ) +penta_prism1.GetPointIds().SetId( 8, 8 ) +penta_prism1.GetPointIds().SetId( 9, 9 ) + +penta_prism2: vtkPentagonalPrism = vtkPentagonalPrism() +penta_prism2.GetPointIds().SetId( 0, 10 ) +penta_prism2.GetPointIds().SetId( 1, 11 ) +penta_prism2.GetPointIds().SetId( 2, 12 ) +penta_prism2.GetPointIds().SetId( 3, 13 ) +penta_prism2.GetPointIds().SetId( 4, 14 ) +penta_prism2.GetPointIds().SetId( 5, 15 ) +penta_prism2.GetPointIds().SetId( 6, 16 ) +penta_prism2.GetPointIds().SetId( 7, 17 ) +penta_prism2.GetPointIds().SetId( 8, 18 ) +penta_prism2.GetPointIds().SetId( 9, 19 ) + +penta_prism3: vtkPentagonalPrism = vtkPentagonalPrism() +penta_prism3.GetPointIds().SetId( 0, 20 ) +penta_prism3.GetPointIds().SetId( 1, 21 ) +penta_prism3.GetPointIds().SetId( 2, 22 ) +penta_prism3.GetPointIds().SetId( 3, 23 ) +penta_prism3.GetPointIds().SetId( 4, 24 ) +penta_prism3.GetPointIds().SetId( 5, 25 ) +penta_prism3.GetPointIds().SetId( 6, 26 ) +penta_prism3.GetPointIds().SetId( 7, 27 ) +penta_prism3.GetPointIds().SetId( 8, 28 ) +penta_prism3.GetPointIds().SetId( 9, 29 ) + +penta_prism4: vtkPentagonalPrism = vtkPentagonalPrism() +penta_prism4.GetPointIds().SetId( 0, 30 ) +penta_prism4.GetPointIds().SetId( 1, 31 ) +penta_prism4.GetPointIds().SetId( 2, 32 ) +penta_prism4.GetPointIds().SetId( 3, 33 ) +penta_prism4.GetPointIds().SetId( 4, 34 ) +penta_prism4.GetPointIds().SetId( 5, 35 ) +penta_prism4.GetPointIds().SetId( 6, 36 ) +penta_prism4.GetPointIds().SetId( 7, 37 ) +penta_prism4.GetPointIds().SetId( 8, 38 ) +penta_prism4.GetPointIds().SetId( 9, 39 ) + +penta_prism_cells = vtkCellArray() +penta_prism_cells.InsertNextCell( penta_prism1 ) +penta_prism_cells.InsertNextCell( penta_prism2 ) +penta_prism_cells.InsertNextCell( penta_prism3 ) +penta_prism_cells.InsertNextCell( penta_prism4 ) + +penta_prism_grid = vtkUnstructuredGrid() +penta_prism_grid.SetPoints( points_penta_prisms ) +penta_prism_grid.SetCells( VTK_PENTAGONAL_PRISM, penta_prism_cells ) + +# one of every other pentagonal prism has invalid ordering +penta_prism_grid_invalid = vtkUnstructuredGrid() +penta_prism_grid_invalid.DeepCopy( penta_prism_grid ) +for i in range( 2 ): + reorder_cell_nodes( penta_prism_grid_invalid, i * 2 + 1, to_change_order[ "Prism5" ] ) +""" +4 hexagonal prisms +""" +points_hexa_prisms: vtkPoints = vtkPoints() +points_hexa_prisms_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.5, 0.5, 0.0 ), + ( 1.0, 1.0, 0.0 ), ( 0.0, 1.0, 0.0 ), ( -0.5, 0.5, 0.0 ), + ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), ( 1.5, 0.5, 1.0 ), + ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ), ( -0.5, 0.5, 1.0 ), + ( 2.0, 0.0, 0.0 ), ( 3.0, 0.0, 0.0 ), ( 3.5, 0.5, 0.0 ), + ( 3.0, 1.0, 0.0 ), ( 2.0, 1.0, 0.0 ), ( 1.5, 0.5, 0.0 ), + ( 2.0, 0.0, 1.0 ), ( 3.0, 0.0, 1.0 ), ( 3.5, 0.5, 1.0 ), + ( 3.0, 1.0, 1.0 ), ( 2.0, 1.0, 1.0 ), ( 1.5, 0.5, 1.0 ), + ( 0.0, 2.0, 0.0 ), ( 1.0, 2.0, 0.0 ), ( 1.5, 2.5, 0.0 ), + ( 1.0, 3.0, 0.0 ), ( 0.0, 3.0, 0.0 ), ( -0.5, 2.5, 0.0 ), + ( 0.0, 2.0, 1.0 ), ( 1.0, 2.0, 1.0 ), ( 1.5, 2.5, 1.0 ), + ( 1.0, 3.0, 1.0 ), ( 0.0, 3.0, 1.0 ), ( -0.5, 2.5, 1.0 ), + ( 2.0, 2.0, 0.0 ), ( 3.0, 2.0, 0.0 ), ( 3.5, 2.5, 0.0 ), + ( 3.0, 3.0, 0.0 ), ( 2.0, 3.0, 0.0 ), ( 1.5, 2.5, 0.0 ), + ( 2.0, 2.0, 1.0 ), ( 3.0, 2.0, 1.0 ), ( 3.5, 2.5, 1.0 ), + ( 3.0, 3.0, 1.0 ), ( 2.0, 3.0, 1.0 ), ( 1.5, 2.5, 1.0 ) ] +for point_hexa_prism in points_hexa_prisms_coords: + points_hexa_prisms.InsertNextPoint( point_hexa_prism ) + +hexa_prism1: vtkHexagonalPrism = vtkHexagonalPrism() +for i in range( 12 ): + hexa_prism1.GetPointIds().SetId( i, i ) + +hexa_prism2: vtkHexagonalPrism = vtkHexagonalPrism() +for i in range( 12 ): + hexa_prism2.GetPointIds().SetId( i, i + 12 ) + +hexa_prism3: vtkHexagonalPrism = vtkHexagonalPrism() +for i in range( 12 ): + hexa_prism3.GetPointIds().SetId( i, i + 24 ) + +hexa_prism4: vtkHexagonalPrism = vtkHexagonalPrism() +for i in range( 12 ): + hexa_prism4.GetPointIds().SetId( i, i + 36 ) + +hexa_prism_cells = vtkCellArray() +hexa_prism_cells.InsertNextCell( hexa_prism1 ) +hexa_prism_cells.InsertNextCell( hexa_prism2 ) +hexa_prism_cells.InsertNextCell( hexa_prism3 ) +hexa_prism_cells.InsertNextCell( hexa_prism4 ) + +hexa_prism_grid = vtkUnstructuredGrid() +hexa_prism_grid.SetPoints( points_hexa_prisms ) +hexa_prism_grid.SetCells( VTK_HEXAGONAL_PRISM, hexa_prism_cells ) + +# one of every other hexagonal prism has invalid ordering +hexa_prism_grid_invalid = vtkUnstructuredGrid() +hexa_prism_grid_invalid.DeepCopy( hexa_prism_grid ) +for i in range( 2 ): + reorder_cell_nodes( hexa_prism_grid_invalid, i * 2 + 1, to_change_order[ "Prism6" ] ) +""" +2 hexahedrons, 2 tetrahedrons, 2 wedges, 2 pyramids, 2 voxels, 2 pentagonal prisms and 2 hexagonal prisms +""" +points_mix: vtkPoints = vtkPoints() +points_mix_coords: list[ tuple[ float ] ] = [ + ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 2.0, 0.0, 0.0 ), ( 2.5, -0.5, 0.0 ), ( 3.0, 0.0, 0.0 ), ( 3.5, -0.5, 0.0 ), + ( 4.0, 0.0, 0.0 ), ( 4.5, -0.5, 0.0 ), ( 5.0, 0.0, 0.0 ), ( 5.5, -0.5, 0.0 ), ( 6.0, 0.5, 0.0 ), ( 0.0, 1.0, 0.0 ), + ( 1.0, 1.0, 0.0 ), ( 2.0, 1.0, 0.0 ), ( 2.5, 1.5, 0.0 ), ( 3.0, 1.0, 0.0 ), ( 4.0, 1.0, 0.0 ), ( 5.0, 1.0, 0.0 ), + ( 5.5, 1.5, 0.0 ), ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), ( 2.0, 0.0, 1.0 ), ( 2.5, -0.5, 1.0 ), ( 3.0, 0.0, 1.0 ), + ( 3.5, -0.5, 1.0 ), ( 4.0, 0.0, 1.0 ), ( 4.5, -0.5, 1.0 ), ( 5.0, 0.0, 1.0 ), ( 5.5, -0.5, 1.0 ), ( 6.0, 0.5, 1.0 ), + ( 0.0, 1.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( 2.0, 1.0, 1.0 ), ( 2.5, 1.5, 1.0 ), ( 3.0, 1.0, 1.0 ), ( 4.0, 1.0, 1.0 ), + ( 5.0, 1.0, 1.0 ), ( 5.5, 1.5, 1.0 ), ( 0.5, 0.5, 2.0 ), ( 0.5, 1.5, 2.0 ), ( 1.5, 0.5, 2.0 ), ( 1.5, 1.5, 2.0 ), + ( 2.0, 0.0, 2.0 ), ( 2.5, -0.5, 2.0 ), ( 3.0, 0.0, 2.0 ), ( 3.0, 1.0, 2.0 ), ( 2.5, 1.5, 2.0 ), ( 2.0, 1.0, 2.0 ), + ( 5.0, 0.0, 2.0 ), ( 5.5, -0.5, 2.0 ), ( 6.0, 0.5, 2.0 ), ( 5.5, 1.5, 2.0 ), ( 5.0, 1.0, 2.0 ) +] +for point_mix in points_mix_coords: + points_mix.InsertNextPoint( point_mix ) + +mix_hex1: vtkHexahedron = vtkHexahedron() +mix_hex1.GetPointIds().SetId( 0, 0 ) +mix_hex1.GetPointIds().SetId( 1, 1 ) +mix_hex1.GetPointIds().SetId( 2, 12 ) +mix_hex1.GetPointIds().SetId( 3, 11 ) +mix_hex1.GetPointIds().SetId( 4, 19 ) +mix_hex1.GetPointIds().SetId( 5, 20 ) +mix_hex1.GetPointIds().SetId( 6, 31 ) +mix_hex1.GetPointIds().SetId( 7, 30 ) + +mix_hex2: vtkHexahedron = vtkHexahedron() +mix_hex2.GetPointIds().SetId( 0, 1 ) +mix_hex2.GetPointIds().SetId( 1, 2 ) +mix_hex2.GetPointIds().SetId( 2, 13 ) +mix_hex2.GetPointIds().SetId( 3, 12 ) +mix_hex2.GetPointIds().SetId( 4, 20 ) +mix_hex2.GetPointIds().SetId( 5, 21 ) +mix_hex2.GetPointIds().SetId( 6, 32 ) +mix_hex2.GetPointIds().SetId( 7, 31 ) + +mix_hex3: vtkHexahedron = vtkHexahedron() +mix_hex3.GetPointIds().SetId( 0, 4 ) +mix_hex3.GetPointIds().SetId( 1, 6 ) +mix_hex3.GetPointIds().SetId( 2, 16 ) +mix_hex3.GetPointIds().SetId( 3, 15 ) +mix_hex3.GetPointIds().SetId( 4, 23 ) +mix_hex3.GetPointIds().SetId( 5, 25 ) +mix_hex3.GetPointIds().SetId( 6, 35 ) +mix_hex3.GetPointIds().SetId( 7, 34 ) + +mix_hex4: vtkHexahedron = vtkHexahedron() +mix_hex4.GetPointIds().SetId( 0, 6 ) +mix_hex4.GetPointIds().SetId( 1, 8 ) +mix_hex4.GetPointIds().SetId( 2, 17 ) +mix_hex4.GetPointIds().SetId( 3, 16 ) +mix_hex4.GetPointIds().SetId( 4, 25 ) +mix_hex4.GetPointIds().SetId( 5, 27 ) +mix_hex4.GetPointIds().SetId( 6, 36 ) +mix_hex4.GetPointIds().SetId( 7, 35 ) + +mix_pyram1: vtkPyramid = vtkPyramid() +mix_pyram1.GetPointIds().SetId( 0, 19 ) +mix_pyram1.GetPointIds().SetId( 1, 20 ) +mix_pyram1.GetPointIds().SetId( 2, 31 ) +mix_pyram1.GetPointIds().SetId( 3, 30 ) +mix_pyram1.GetPointIds().SetId( 4, 38 ) + +mix_pyram2: vtkPyramid = vtkPyramid() +mix_pyram2.GetPointIds().SetId( 0, 20 ) +mix_pyram2.GetPointIds().SetId( 1, 21 ) +mix_pyram2.GetPointIds().SetId( 2, 32 ) +mix_pyram2.GetPointIds().SetId( 3, 31 ) +mix_pyram2.GetPointIds().SetId( 4, 40 ) + +mix_tetra1: vtkTetra = vtkTetra() +mix_tetra1.GetPointIds().SetId( 0, 31 ) +mix_tetra1.GetPointIds().SetId( 1, 39 ) +mix_tetra1.GetPointIds().SetId( 2, 30 ) +mix_tetra1.GetPointIds().SetId( 3, 38 ) + +mix_tetra2: vtkTetra = vtkTetra() +mix_tetra2.GetPointIds().SetId( 0, 32 ) +mix_tetra2.GetPointIds().SetId( 1, 41 ) +mix_tetra2.GetPointIds().SetId( 2, 31 ) +mix_tetra2.GetPointIds().SetId( 3, 40 ) + +mix_hex_prism1: vtkHexagonalPrism = vtkHexagonalPrism() +mix_hex_prism1.GetPointIds().SetId( 0, 2 ) +mix_hex_prism1.GetPointIds().SetId( 1, 3 ) +mix_hex_prism1.GetPointIds().SetId( 2, 4 ) +mix_hex_prism1.GetPointIds().SetId( 3, 15 ) +mix_hex_prism1.GetPointIds().SetId( 4, 14 ) +mix_hex_prism1.GetPointIds().SetId( 5, 13 ) +mix_hex_prism1.GetPointIds().SetId( 6, 21 ) +mix_hex_prism1.GetPointIds().SetId( 7, 22 ) +mix_hex_prism1.GetPointIds().SetId( 8, 23 ) +mix_hex_prism1.GetPointIds().SetId( 9, 34 ) +mix_hex_prism1.GetPointIds().SetId( 10, 33 ) +mix_hex_prism1.GetPointIds().SetId( 11, 32 ) + +mix_hex_prism2: vtkHexagonalPrism = vtkHexagonalPrism() +mix_hex_prism2.GetPointIds().SetId( 0, 21 ) +mix_hex_prism2.GetPointIds().SetId( 1, 22 ) +mix_hex_prism2.GetPointIds().SetId( 2, 23 ) +mix_hex_prism2.GetPointIds().SetId( 3, 34 ) +mix_hex_prism2.GetPointIds().SetId( 4, 33 ) +mix_hex_prism2.GetPointIds().SetId( 5, 32 ) +mix_hex_prism2.GetPointIds().SetId( 6, 42 ) +mix_hex_prism2.GetPointIds().SetId( 7, 43 ) +mix_hex_prism2.GetPointIds().SetId( 8, 44 ) +mix_hex_prism2.GetPointIds().SetId( 9, 45 ) +mix_hex_prism2.GetPointIds().SetId( 10, 46 ) +mix_hex_prism2.GetPointIds().SetId( 11, 47 ) + +mix_wedge1: vtkWedge = vtkWedge() +mix_wedge1.GetPointIds().SetId( 0, 23 ) +mix_wedge1.GetPointIds().SetId( 1, 24 ) +mix_wedge1.GetPointIds().SetId( 2, 25 ) +mix_wedge1.GetPointIds().SetId( 3, 4 ) +mix_wedge1.GetPointIds().SetId( 4, 5 ) +mix_wedge1.GetPointIds().SetId( 5, 6 ) + +mix_wedge2: vtkWedge = vtkWedge() +mix_wedge2.GetPointIds().SetId( 0, 25 ) +mix_wedge2.GetPointIds().SetId( 1, 26 ) +mix_wedge2.GetPointIds().SetId( 2, 27 ) +mix_wedge2.GetPointIds().SetId( 3, 6 ) +mix_wedge2.GetPointIds().SetId( 4, 7 ) +mix_wedge2.GetPointIds().SetId( 5, 8 ) + +mix_penta_prism1: vtkPentagonalPrism = vtkPentagonalPrism() +mix_penta_prism1.GetPointIds().SetId( 0, 8 ) +mix_penta_prism1.GetPointIds().SetId( 1, 9 ) +mix_penta_prism1.GetPointIds().SetId( 2, 10 ) +mix_penta_prism1.GetPointIds().SetId( 3, 18 ) +mix_penta_prism1.GetPointIds().SetId( 4, 17 ) +mix_penta_prism1.GetPointIds().SetId( 5, 27 ) +mix_penta_prism1.GetPointIds().SetId( 6, 28 ) +mix_penta_prism1.GetPointIds().SetId( 7, 29 ) +mix_penta_prism1.GetPointIds().SetId( 8, 37 ) +mix_penta_prism1.GetPointIds().SetId( 9, 36 ) + +mix_penta_prism2: vtkPentagonalPrism = vtkPentagonalPrism() +mix_penta_prism2.GetPointIds().SetId( 0, 27 ) +mix_penta_prism2.GetPointIds().SetId( 1, 28 ) +mix_penta_prism2.GetPointIds().SetId( 2, 29 ) +mix_penta_prism2.GetPointIds().SetId( 3, 37 ) +mix_penta_prism2.GetPointIds().SetId( 4, 36 ) +mix_penta_prism2.GetPointIds().SetId( 5, 48 ) +mix_penta_prism2.GetPointIds().SetId( 6, 49 ) +mix_penta_prism2.GetPointIds().SetId( 7, 50 ) +mix_penta_prism2.GetPointIds().SetId( 8, 51 ) +mix_penta_prism2.GetPointIds().SetId( 9, 52 ) + +# this mix grid has only valid cell volumes +mix_grid = vtkUnstructuredGrid() +mix_grid.SetPoints( points_mix ) +all_cell_types_mix_grid = [ + VTK_HEXAHEDRON, VTK_HEXAHEDRON, VTK_PYRAMID, VTK_PYRAMID, VTK_TETRA, VTK_TETRA, VTK_HEXAGONAL_PRISM, + VTK_HEXAGONAL_PRISM, VTK_HEXAHEDRON, VTK_HEXAHEDRON, VTK_WEDGE, VTK_WEDGE, VTK_PENTAGONAL_PRISM, + VTK_PENTAGONAL_PRISM +] +all_cell_names_mix_grid = [ + "Hexahedron", "Hexahedron", "Pyramid", "Pyramid", "Tetrahedron", "Tetrahedron", "Prism6", "Prism6", "Hexahedron", + "Hexahedron", "Wedge", "Wedge", "Prism5", "Prism5" +] +all_cells_mix_grid = [ + mix_hex1, mix_hex2, mix_pyram1, mix_pyram2, mix_tetra1, mix_tetra2, mix_hex_prism1, mix_hex_prism2, mix_hex3, + mix_hex4, mix_wedge1, mix_wedge2, mix_penta_prism1, mix_penta_prism2 +] +for cell_type, cell in zip( all_cell_types_mix_grid, all_cells_mix_grid ): + mix_grid.InsertNextCell( cell_type, cell.GetPointIds() ) + +# this mix grid has one invalid cell for each different element type +mix_grid_invalid = vtkUnstructuredGrid() +mix_grid_invalid.DeepCopy( mix_grid ) +for i in range( len( all_cell_types_mix_grid ) // 2 ): + reorder_cell_nodes( mix_grid_invalid, i * 2 + 1, to_change_order[ all_cell_names_mix_grid[ i * 2 + 1 ] ] ) + + +class TestClass: + + def test_reorder_cell_nodes( self ): + expected_nodes_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), + ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ) ] + for i in range( one_hex.GetCell( 0 ).GetNumberOfPoints() ): + assert one_hex.GetCell( 0 ).GetPoints().GetPoint( i ) == expected_nodes_coords[ i ] + + # reorder the cell to make it invalid + reorder_cell_nodes( one_hex, 0, to_change_order[ "Hexahedron" ] ) + expected_nodes_coords_modified = [ ( 0.0, 0.0, 0.0 ), ( 0.0, 1.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 1.0, 0.0, 0.0 ), + ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ) ] + for i in range( one_hex.GetCell( 0 ).GetNumberOfPoints() ): + assert one_hex.GetCell( 0 ).GetPoints().GetPoint( i ) == expected_nodes_coords_modified[ i ] + + # reorder the cell again to make it valid again + reorder_cell_nodes( one_hex, 0, to_change_order[ "Hexahedron" ] ) + expected_nodes_coords_modified2 = [ ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 0.0, 1.0, 0.0 ), + ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ) ] + for i in range( one_hex.GetCell( 0 ).GetNumberOfPoints() ): + assert one_hex.GetCell( 0 ).GetPoints().GetPoint( i ) == expected_nodes_coords_modified2[ i ] + + def test_compute_mesh_cells_volume( self ): + grid_volumes = { + hexahedrons_grid: [ 1.0, 1.0, 1.0, 1.0 ], + hexahedrons_grid_invalid: [ 1.0, -0.333, 1.0, -0.333 ], + tetras_grid: [ 0.167, 0.167, 0.167, 0.167 ], + tetras_grid_invalid: [ 0.167, -0.167, 0.167, -0.167 ], + pyramids_grid: [ 0.333, 0.333, 0.333, 0.333 ], + pyramids_grid_invalid: [ 0.333, -0.333, 0.333, -0.333 ], + voxels_grid: [ 1.0, 1.0, 1.0, 1.0 ], + wedges_grid: [ 0.5, 0.5, 0.5, 0.5 ], + wedges_grid_invalid: [ 0.5, -0.167, 0.5, -0.167 ], + penta_prism_grid: [ 1.25, 1.25, 1.25, 1.25 ], + penta_prism_grid_invalid: [ 1.25, -0.083, 1.25, -0.083 ], + hexa_prism_grid: [ 1.5, 1.5, 1.5, 1.5 ], + hexa_prism_grid_invalid: [ 1.5, -0.333, 1.5, -0.333 ], + mix_grid: [ 1.0, 1.0, 0.333, 0.333, 0.167, 0.167, 1.5, 1.5, 1.0, 1.0, 0.25, 0.25, 1.25, 1.25 ], + mix_grid_invalid: + [ 1.0, -0.333, 0.333, -0.333, 0.167, -0.167, 1.5, -0.333, 1.0, -0.333, 0.25, -0.083, 1.25, -0.083 ] + } + for grid, volumes_expected in grid_volumes.items(): + volumes_computed = feo.compute_mesh_cells_volume( grid ) + for i in range( len( volumes_computed ) ): + assert round( float( volumes_computed[ i ] ), 3 ) == volumes_expected[ i ] + + def test_is_cell_to_reorder( self ): + grid_needs_ordering = { + hexahedrons_grid: [ False ] * 4, + hexahedrons_grid_invalid: [ i % 2 != 0 for i in range( 4 ) ], + tetras_grid: [ False ] * 4, + tetras_grid_invalid: [ i % 2 != 0 for i in range( 4 ) ], + pyramids_grid: [ False ] * 4, + pyramids_grid_invalid: [ i % 2 != 0 for i in range( 4 ) ], + wedges_grid: [ False ] * 4, + wedges_grid_invalid: [ i % 2 != 0 for i in range( 4 ) ], + penta_prism_grid: [ False ] * 4, + penta_prism_grid_invalid: [ i % 2 != 0 for i in range( 4 ) ], + hexa_prism_grid: [ False ] * 4, + hexa_prism_grid_invalid: [ i % 2 != 0 for i in range( 4 ) ], + mix_grid: [ False ] * 14, + mix_grid_invalid: [ i % 2 != 0 for i in range( 14 ) ] + } + for grid, needs_ordering in grid_needs_ordering.items(): + volumes = feo.compute_mesh_cells_volume( grid ) + for i in range( len( volumes ) ): + options = opt( out, to_change_order, "negative" ) + assert feo.is_cell_to_reorder( volumes[ i ], options ) == needs_ordering[ i ] + options = opt( out, to_change_order, "positive" ) + assert feo.is_cell_to_reorder( volumes[ i ], options ) != needs_ordering[ i ] + options = opt( out, to_change_order, "all" ) + assert feo.is_cell_to_reorder( volumes[ i ], options ) == True + + def test_get_cell_types_and_number( self ): + assert feo.get_cell_types_and_number( hexahedrons_grid ) == ( [ VTK_HEXAHEDRON ], [ 4 ] ) + assert feo.get_cell_types_and_number( hexahedrons_grid_invalid ) == ( [ VTK_HEXAHEDRON ], [ 4 ] ) + assert feo.get_cell_types_and_number( tetras_grid ) == ( [ VTK_TETRA ], [ 4 ] ) + assert feo.get_cell_types_and_number( tetras_grid_invalid ) == ( [ VTK_TETRA ], [ 4 ] ) + assert feo.get_cell_types_and_number( pyramids_grid ) == ( [ VTK_PYRAMID ], [ 4 ] ) + assert feo.get_cell_types_and_number( pyramids_grid_invalid ) == ( [ VTK_PYRAMID ], [ 4 ] ) + assert feo.get_cell_types_and_number( wedges_grid ) == ( [ VTK_WEDGE ], [ 4 ] ) + assert feo.get_cell_types_and_number( wedges_grid_invalid ) == ( [ VTK_WEDGE ], [ 4 ] ) + assert feo.get_cell_types_and_number( penta_prism_grid ) == ( [ VTK_PENTAGONAL_PRISM ], [ 4 ] ) + assert feo.get_cell_types_and_number( penta_prism_grid_invalid ) == ( [ VTK_PENTAGONAL_PRISM ], [ 4 ] ) + assert feo.get_cell_types_and_number( hexa_prism_grid ) == ( [ VTK_HEXAGONAL_PRISM ], [ 4 ] ) + assert feo.get_cell_types_and_number( hexa_prism_grid_invalid ) == ( [ VTK_HEXAGONAL_PRISM ], [ 4 ] ) + assert feo.get_cell_types_and_number( mix_grid ) == ( [ + VTK_TETRA, VTK_HEXAHEDRON, VTK_WEDGE, VTK_PYRAMID, VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM + ], [ 2, 4, 2, 2, 2, 2 ] ) + assert feo.get_cell_types_and_number( mix_grid_invalid ) == ( [ + VTK_TETRA, VTK_HEXAHEDRON, VTK_WEDGE, VTK_PYRAMID, VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM + ], [ 2, 4, 2, 2, 2, 2 ] ) + expected_error: str = f"Invalid type '11' for GEOS is in the mesh. Dying ..." + with pytest.raises( ValueError, match=expected_error ): + feo.get_cell_types_and_number( voxels_grid ) + + def test_reorder_nodes_to_new_mesh( self ): + options = opt( out, to_change_order, "negative" ) + # single element grids except voxels because it is an invalid cell type for GEOS + grid_cell_type = { + hexahedrons_grid_invalid: VTK_HEXAHEDRON, + tetras_grid_invalid: VTK_TETRA, + pyramids_grid_invalid: VTK_PYRAMID, + wedges_grid_invalid: VTK_WEDGE, + penta_prism_grid_invalid: VTK_PENTAGONAL_PRISM, + hexa_prism_grid_invalid: VTK_HEXAGONAL_PRISM + } + for grid, cell_type in grid_cell_type.items(): + new_invalid = vtkUnstructuredGrid() + new_invalid.DeepCopy( grid ) + not_use_invalid, reorder_stats = feo.reorder_nodes_to_new_mesh( new_invalid, options ) + expected = { + "Types reordered": [ VTK_TYPE_TO_NAME[ cell_type ] ], + "Number of cells reordered": [ 2 ], + "Types non reordered": [ VTK_TYPE_TO_NAME[ cell_type ] ], + "Number of cells non reordered": [ 2 ] + } + for prop in expected.keys(): + assert reorder_stats[ prop ] == expected[ prop ] + + # mix elements grid + mix_invalid = vtkUnstructuredGrid() + mix_invalid.DeepCopy( mix_grid_invalid ) + not_use_invalid, mix_stats = feo.reorder_nodes_to_new_mesh( mix_invalid, options ) + expected = { + "Types reordered": [ VTK_TYPE_TO_NAME[ cell_type ] for cell_type in [ 12, 15, 16, 14, 10, 13 ] ], + "Number of cells reordered": [ 2, 1, 1, 1, 1, 1 ], + "Types non reordered": [ VTK_TYPE_TO_NAME[ cell_type ] for cell_type in [ 12, 15, 16, 14, 10, 13 ] ], + "Number of cells non reordered": [ 2, 1, 1, 1, 1, 1 ] + } + for prop in expected.keys(): + assert mix_stats[ prop ] == expected[ prop ] + + def test_fix_elements_orderings_execution( self ): + # for mix_grid_invalid mesh, checks that reordered mesh was created and that reoredring_stats are valid + write_mesh( mix_grid_invalid, test_file ) + invalidTest = False + command = [ + "python", MESH_DOCTOR_FILEPATH, "-v", "-i", test_file.output, "fix_elements_orderings", "--Hexahedron", + str( to_change_order[ "Hexahedron" ] ).replace( "[", "" ).replace( "]", "" ), "--Tetrahedron", + str( to_change_order[ "Tetrahedron" ] ).replace( "[", "" ).replace( "]", "" ), "--Pyramid", + str( to_change_order[ "Pyramid" ] ).replace( "[", "" ).replace( "]", "" ), "--Wedge", + str( to_change_order[ "Wedge" ] ).replace( "[", "" ).replace( "]", "" ), "--Prism5", + str( to_change_order[ "Prism5" ] ).replace( "[", "" ).replace( "]", "" ), "--Prism6", + str( to_change_order[ "Prism6" ] ).replace( "[", "" ).replace( "]", "" ), "--volume_to_reorder", "negative", + "--data-mode", "binary", "--output", filepath_reordered_mesh + ] + try: + result = subprocess.run( command, shell=True, stderr=subprocess.PIPE, universal_newlines=True ) + os.remove( filepath_reordered_mesh ) + stderr = result.stderr + assert result.returncode == 0 + raw_stderr = r"{}".format( stderr ) + pattern = r"\[.*?\]\[.*?\] (.*)" + matches = re.findall( pattern, raw_stderr ) + no_log = "\n".join( matches ) + reordering_stats: str = no_log[ no_log.index( "Number of cells reordered" ): ] + expected_stats: str = ( "Number of cells reordered:\n" + "\tCellType\tNumber\n" + f"\tHexahedron\t\t2\n" + + f"\tPrism5\t\t1\n" + f"\tPrism6\t\t1\n" + f"\tPyramid\t\t1\n" + + f"\tTetrahedron\t\t1\n" + f"\tWedge\t\t1\n" + "Number of cells non reordered:\n" + + "\tCellType\tNumber\n" + f"\tHexahedron\t\t2\n" + f"\tPrism5\t\t1\n" + + f"\tPrism6\t\t1\n" + f"\tPyramid\t\t1\n" + f"\tTetrahedron\t\t1\n" + + f"\tWedge\t\t1" ) + assert reordering_stats == expected_stats + except Exception as e: + logging.error( "Invalid command input. Test has failed." ) + logging.error( e ) + invalidTest = True + os.remove( test_file.output ) + if invalidTest: + raise ValueError( "test_fix_elements_orderings_execution has failed." ) \ No newline at end of file From 8284615a08087aef78c45ed1135e04b9e0cadded Mon Sep 17 00:00:00 2001 From: alexbenedicto Date: Mon, 16 Sep 2024 17:58:31 -0700 Subject: [PATCH 2/5] Remove sys exit to replace by ValueError --- .../src/geos/mesh/doctor/checks/vtk_utils.py | 68 ++++++++----------- 1 file changed, 29 insertions(+), 39 deletions(-) diff --git a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py index c27872ac..8f52cbb8 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py @@ -1,25 +1,11 @@ -from dataclasses import dataclass -from os import path, access, R_OK, W_OK +import os import logging -from sys import exit -from typing import ( - Any, - Iterator, - Optional, -) - -from vtkmodules.vtkCommonCore import ( - vtkIdList, ) -from vtkmodules.vtkCommonDataModel import ( - vtkUnstructuredGrid, ) -from vtkmodules.vtkIOLegacy import ( - vtkUnstructuredGridWriter, - vtkUnstructuredGridReader, -) -from vtkmodules.vtkIOXML import ( - vtkXMLUnstructuredGridReader, - vtkXMLUnstructuredGridWriter, -) +from dataclasses import dataclass +from typing import Iterator, Optional +from vtkmodules.vtkCommonCore import vtkIdList +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from vtkmodules.vtkIOLegacy import vtkUnstructuredGridWriter, vtkUnstructuredGridReader +from vtkmodules.vtkIOXML import vtkXMLUnstructuredGridReader, vtkXMLUnstructuredGridWriter @dataclass( frozen=True ) @@ -36,7 +22,7 @@ def to_vtk_id_list( data ) -> vtkIdList: return result -def vtk_iter( l ) -> Iterator[ Any ]: +def vtk_iter( l ) -> Iterator[ any ]: """ Utility function transforming a vtk "container" (e.g. vtkIdList) into an iterable to be used for building built-ins python containers. :param l: A vtk container. @@ -62,22 +48,22 @@ def is_filepath_valid( filepath: str, is_input=True ) -> bool: Returns: bool: False if invalid, True instead. """ - dirname = path.dirname( filepath ) - if not path.isdir( dirname ): + dirname = os.path.dirname( filepath ) + if not os.path.isdir( dirname ): logging.error( f"The directory '{dirname}' specified does not exist." ) return False if is_input: - if not access( dirname, R_OK ): + if not os.access( dirname, os.R_OK ): logging.error( f"You do not have rights to read in directory '{dirname}' specified for the file." ) return False - if not path.exists( filepath ): + if not os.path.exists( filepath ): logging.error( f"The file specified does not exist." ) return False else: - if not access( dirname, W_OK ): + if not os.access( dirname, os.W_OK ): logging.error( f"You do not have rights to write in directory '{dirname}' specified for the file." ) return False - if path.exists( filepath ): + if os.path.exists( filepath ): logging.error( f"A file with the same name already exists. No over-writing possible." ) return False return True @@ -116,10 +102,11 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid: If first guess does not work, eventually all the others reader available will be tested. :return: A unstructured grid. """ - if not is_filepath_valid( vtk_input_file, True ): - logging.error( f"Invalid filepath to read. Dying ..." ) - exit( 1 ) - file_extension = path.splitext( vtk_input_file )[ -1 ] + if not is_filepath_valid( vtk_input_file ): + err_msg: str = f"Invalid file path. Could not read \"{vtk_input_file}\". Dying..." + logging.error( err_msg ) + raise ValueError( err_msg ) + file_extension = os.path.splitext( vtk_input_file )[ -1 ] extension_to_reader = { ".vtk": __read_vtk, ".vtu": __read_vtu } # Testing first the reader that should match if file_extension in extension_to_reader: @@ -132,8 +119,9 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid: if output_mesh: return output_mesh # No reader did work. Dying. - logging.critical( f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\". Dying..." ) - exit( 1 ) + err_msg = f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\". Dying..." + logging.error( err_msg ) + raise ValueError( err_msg ) def __write_vtk( mesh: vtkUnstructuredGrid, output: str ) -> int: @@ -162,15 +150,17 @@ def write_mesh( mesh: vtkUnstructuredGrid, vtk_output: VtkOutput ) -> int: :return: 0 in case of success. """ if not is_filepath_valid( vtk_output.output, False ): - logging.error( f"Invalid filepath to write. Dying ..." ) - exit( 1 ) - file_extension = path.splitext( vtk_output.output )[ -1 ] + err_msg = "Invalid filepath to write. Dying ..." + logging.error( err_msg ) + raise ValueError( err_msg ) + file_extension = os.path.splitext( vtk_output.output )[ -1 ] if file_extension == ".vtk": success_code = __write_vtk( mesh, vtk_output.output ) elif file_extension == ".vtu": success_code = __write_vtu( mesh, vtk_output.output, vtk_output.is_data_mode_binary ) else: # No writer found did work. Dying. - logging.critical( f"Could not find the appropriate VTK writer for extension \"{file_extension}\". Dying..." ) - exit( 1 ) + err_msg = f"Could not find the appropriate VTK writer for extension \"{file_extension}\". Dying..." + logging.error( err_msg ) + raise ValueError( err_msg ) return 0 if success_code else 2 # the Write member function return 1 in case of success, 0 otherwise. From 5a593a4cbf92df19977c88ce6b347f7193a70993 Mon Sep 17 00:00:00 2001 From: alexbenedicto Date: Thu, 19 Sep 2024 19:18:53 -0700 Subject: [PATCH 3/5] Node ordering correction is automatic, user does not need to enter the node ordering himself. Just needs to specify the name of elements to check. --- docs/geos-mesh.rst | 38 +++---- .../doctor/checks/fix_elements_orderings.py | 107 +++++++++++++++--- .../src/geos/mesh/doctor/checks/vtk_utils.py | 6 +- .../parsing/fix_elements_orderings_parsing.py | 56 ++++----- .../tests/test_fix_elements_orderings.py | 32 +++--- 5 files changed, 148 insertions(+), 91 deletions(-) diff --git a/docs/geos-mesh.rst b/docs/geos-mesh.rst index 739d4936..243592ed 100644 --- a/docs/geos-mesh.rst +++ b/docs/geos-mesh.rst @@ -117,41 +117,35 @@ Cells with negative volumes will typically be an issue for ``geos`` and should b """""""""""""""""""""""""" It sometimes happens that an exported mesh does not abide by the ``vtk`` orderings. -The ``fix_elements_orderings`` module can rearrange the nodes of given types of elements. +The ``fix_elements_orderings`` module can rearrange the nodes of given types of elements to respect VTK convention. This can be convenient if you cannot regenerate the mesh. .. code-block:: $ python src/geos/mesh/doctor/mesh_doctor.py fix_elements_orderings --help - usage: mesh_doctor.py fix_elements_orderings [-h] [--Tetrahedron 2,0,3,1] [--Pyramid 3,4,0,2,1] [--Wedge 3,5,4,0,2,1] - [--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] - --volume_to_reorder all,positive,negative --output OUTPUT [--data-mode binary, ascii] + usage: mesh_doctor.py fix_elements_orderings [-h] [--cell_names Hexahedron, Tetrahedron, Pyramid, Wedge, Prism5, Prism6] + [--volume_to_reorder all, positive, negative] + --output OUTPUT [--data-mode binary, ascii] options: - -h, --help show this help message and exit - --Tetrahedron 2,0,3,1 [list of integers]: node permutation for "Tetrahedron". - --Pyramid 3,4,0,2,1 [list of integers]: node permutation for "Pyramid". - --Wedge 3,5,4,0,2,1 [list of integers]: node permutation for "Wedge". - --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". - --volume_to_reorder all,positive,negative - [str]: Select which element volume is invalid and needs reordering. 'all' will allow reordering of nodes for every element, regarding of their volume. - 'positive' or 'negative' will only reorder the element with the corresponding volume. - --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. + -h, --help show this help message and exit + --cell_names Hexahedron, Tetrahedron, Pyramid, Wedge, Prism5, Prism6 + [list of str]: Cell names that can be reordered in your grid. You can use multiple names.Defaults to all cell names being used. + --volume_to_reorder all, positive, negative + [str]: Select which element volume is invalid and needs reordering. 'all' will allow reordering of nodes for every element, regarding of their volume. + 'positive' or 'negative' will only reorder the element with the corresponding volume. Defaults to 'negative'. + --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. For example, assume that you have a mesh that contains 3 different cell types like Hexahedrons, Pyramids and Tetrahedrons. After checking ``element_volumes``, you found that all your Hexahedrons and half of your Tetrahedrons have a negative volume. -To correct that, assuming you know the correct node ordering to correct these volumes, you can use this command: +To correct that, you can use this command: .. code-block:: - $ python src/geos/mesh/doctor/mesh_doctor.py -i /path/to/your/mesh/file fix_elements_orderings --Hexahedron 1,6,5,4,7,0,2,3 - --Tetrahedron 2,0,3,1 --volume_to_reorder negative --output /new/path/for/your/new/mesh/reordered/file --data-mode binary + $ python src/geos/mesh/doctor/mesh_doctor.py -i /path/to/your/mesh/file fix_elements_orderings --cell_names Hexahedron,Tetrahedron + --volume_to_reorder negative --output /new/path/for/your/new/mesh/reordered/file --data-mode binary ``generate_cube`` """"""""""""""""" diff --git a/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py b/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py index bc693a99..19e428b7 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py @@ -1,18 +1,18 @@ -from dataclasses import dataclass import numpy as np import logging +from dataclasses import dataclass +from itertools import permutations from vtk import vtkCellSizeFilter -from vtkmodules.vtkCommonCore import vtkIdList from vtkmodules.util.numpy_support import vtk_to_numpy -from vtkmodules.vtkCommonDataModel import ( vtkDataSet, VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, +from vtkmodules.vtkCommonDataModel import ( vtkDataSet, vtkCell, VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM ) -from .vtk_utils import VtkOutput, to_vtk_id_list, write_mesh, read_mesh +from geos.mesh.doctor.checks.vtk_utils import VtkOutput, vtk_iter, to_vtk_id_list, write_mesh, read_mesh @dataclass( frozen=True ) class Options: vtk_output: VtkOutput - cell_name_to_ordering: dict[ str, list[ int ] ] + cell_names_to_reorder: tuple[ str ] volume_to_reorder: str @@ -103,6 +103,90 @@ def is_cell_to_reorder( cell_volume: str, options: Options ) -> bool: return False +def are_face_nodes_counterclockwise( face: any ) -> bool: + """Checks if the nodes of a face are ordered counterclockwise when looking at the plan created by the face. + + Args: + face (any): Face of a vtkCell. + + Returns: + bool: True if counterclockwise, False instead. + """ + face_points = face.GetPoints() + number_points = face_points.GetNumberOfPoints() + if number_points < 3: + err_msg = f"The face has less than 3 nodes which is invalid." + logging.error( err_msg ) + raise ValueError( err_msg ) + # first calculate the normal vector of the face + a = np.array( face_points.GetPoint( 0 ) ) + b = np.array( face_points.GetPoint( 1 ) ) + c = np.array( face_points.GetPoint( 2 ) ) + AB = b - a + AC = c - a + normal = np.cross(AB, AC) + + # then calculate the vector cross products sum of all points of the face with a random point P + # from discussion https://math.stackexchange.com/questions/2152623/determine-the-order-of-a-3d-polygon + P = np.array( [ 0.0, 0.0, 0.0 ] ) # P position does not change the value of sum + all_points = [ np.array( face_points.GetPoint( v ) ) for v in range( number_points ) ] + vector_sum = np.array( [ 0.0, 0.0, 0.0 ] ) + for i in range( number_points ): + PAi = all_points[ i % number_points ] - P + PAiplus1 = all_points[ ( i + 1 ) % number_points ] - P + vector_sum += np.cross( PAi, PAiplus1 ) + vector_sum = vector_sum / 2 # needs to be half + + # if dot product is positive, the nodes are ordered counterclockwise + if np.dot( vector_sum, normal ) > 0.0: + return True + return False + + +def valid_cell_point_ids_ordering( cell: vtkCell ) -> list[ int ]: + """Returns the valid order of point ids of a cell that respect counterclockwise convention of each face. + + Args: + cell (vtkCell): A cell of vtk grid. + + Raises: + ValueError: "No node ids permutation made the face nodes to be ordered counterclockwise." + + Returns: + list[ int ]: [ pt_id0, pt_id2, pt_id1, ..., pt_idN ] + """ + initial_ids_order: list[ int ] = list( vtk_iter( cell.GetPointIds() ) ) + valid_points_id: list[ int ] = initial_ids_order.copy() + + for face_id in range( cell.GetNumberOfFaces() ): + face = cell.GetFace( face_id ) + if are_face_nodes_counterclockwise( face ): + continue # the face nodes respect the convention so continue to next face + + reordered = False + initial_ids: list[ int ] = [ face.GetPointId( i ) for i in range( face.GetNumberOfPoints() ) ] + initial_coords: list[ float ] = [ face.GetPoints().GetPoint( v ) for v in range( face.GetNumberOfPoints() ) ] + for permutation_ids in permutations( initial_ids ): + for i, node_id in enumerate( permutation_ids ): + face.GetPointIds().SetId( i, node_id ) + face.GetPoints().SetPoint( i, initial_coords[ initial_ids.index( node_id ) ] ) + + if are_face_nodes_counterclockwise( face ): # the correct permutation was found + for j, initial_id in enumerate(initial_ids): + if initial_id != permutation_ids[ j ]: + valid_points_id[ initial_ids_order.index( initial_id ) ] = permutation_ids[ j ] + reordered = True + break + + if not reordered: + err_msg = ( f"For face with nodes id '{initial_ids}', that corresponds to coordinates '{initial_coords}'" + + ", no node ids permutation made the face nodes to be ordered counterclockwise." ) + logging.error( err_msg ) + raise ValueError( err_msg ) + + return valid_points_id + + def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple: """From an old mesh, creates a new one where certain cell nodes are reordered to obtain a correct volume. @@ -118,7 +202,7 @@ def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple names_with_totals: dict[ str, int ] = { n: v for n, v in zip( unique_cell_names, total_per_cell_types ) } # sorted dict allow for sorted output of reordering_stats names_with_totals_sorted: dict[ str, int ] = dict( sorted( names_with_totals.items() ) ) - useful_VTK_TYPEs: list[ int ] = [ NAME_TO_VTK_TYPE[ name ] for name in options.cell_name_to_ordering.keys() ] + useful_VTK_TYPEs: list[ int ] = [ NAME_TO_VTK_TYPE[ name ] for name in options.cell_names_to_reorder ] all_cells_volume: np.array = compute_mesh_cells_volume( old_mesh ) # a new mesh with the same data is created from the old mesh new_mesh: vtkDataSet = old_mesh.NewInstance() @@ -132,7 +216,7 @@ def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple "Types non reordered": list( names_with_totals_sorted.keys() ), "Number of cells non reordered": list( names_with_totals_sorted.values() ) } - counter_cells_reordered: dict[ str, int ] = { name: 0 for name in options.cell_name_to_ordering.keys() } + counter_cells_reordered: dict[ str, int ] = { name: 0 for name in options.cell_names_to_reorder } # sorted dict allow for sorted output of reordering_stats ounter_cells_reordered_sorted: dict[ str, int ] = dict( sorted( counter_cells_reordered.items() ) ) # Reordering operations @@ -141,13 +225,8 @@ def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple if vtk_type in useful_VTK_TYPEs: if is_cell_to_reorder( float( all_cells_volume[ cell_id ] ), options ): cell_name: str = VTK_TYPE_TO_NAME[ vtk_type ] - support_point_ids = vtkIdList() - cells.GetCellAtId( cell_id, support_point_ids ) - new_support_point_ids: list[ int ] = list() - node_ordering: list[ int ] = options.cell_name_to_ordering[ cell_name ] - for i in range( len( node_ordering ) ): - new_support_point_ids.append( support_point_ids.GetId( node_ordering[ i ] ) ) - cells.ReplaceCellAtId( cell_id, to_vtk_id_list( new_support_point_ids ) ) + point_ids_ordering: list[ int ] = valid_cell_point_ids_ordering( new_mesh.GetCell( cell_id ) ) + cells.ReplaceCellAtId( cell_id, to_vtk_id_list( point_ids_ordering ) ) ounter_cells_reordered_sorted[ cell_name ] += 1 # Calculation of stats for cell_name_reordered, amount in ounter_cells_reordered_sorted.items(): diff --git a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py index 8f52cbb8..5d384049 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py @@ -103,7 +103,7 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid: :return: A unstructured grid. """ if not is_filepath_valid( vtk_input_file ): - err_msg: str = f"Invalid file path. Could not read \"{vtk_input_file}\". Dying..." + err_msg: str = f"Invalid file path. Could not read \"{vtk_input_file}\"." logging.error( err_msg ) raise ValueError( err_msg ) file_extension = os.path.splitext( vtk_input_file )[ -1 ] @@ -119,7 +119,7 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid: if output_mesh: return output_mesh # No reader did work. Dying. - err_msg = f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\". Dying..." + err_msg = f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\"." logging.error( err_msg ) raise ValueError( err_msg ) @@ -160,7 +160,7 @@ def write_mesh( mesh: vtkUnstructuredGrid, vtk_output: VtkOutput ) -> int: success_code = __write_vtu( mesh, vtk_output.output, vtk_output.is_data_mode_binary ) else: # No writer found did work. Dying. - err_msg = f"Could not find the appropriate VTK writer for extension \"{file_extension}\". Dying..." + err_msg = f"Could not find the appropriate VTK writer for extension \"{file_extension}\"." logging.error( err_msg ) raise ValueError( err_msg ) return 0 if success_code else 2 # the Write member function return 1 in case of success, 0 otherwise. diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py index 641a9505..04a3bf5b 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py +++ b/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py @@ -1,41 +1,32 @@ import logging import random -from geos.mesh.doctor.checks.fix_elements_orderings import Options, Result +from geos.mesh.doctor.checks.fix_elements_orderings import Options, Result, NAME_TO_VTK_TYPE from . import vtk_output_parsing, FIX_ELEMENTS_ORDERINGS -__CELL_NAME_WITH_NUMBER_NODES = { - "Tetrahedron": 4, - "Pyramid": 5, - "Wedge": 6, - "Hexahedron": 8, - "Prism5": 10, - "Prism6": 12 -} +__CELL_NAMES = "cell_names" +__CELL_NAMES_CHOICES = list( NAME_TO_VTK_TYPE.keys() ) __VOLUME_TO_REORDER = "volume_to_reorder" -__VOLUME_TO_REORDER_DEFAULT = "all" +__VOLUME_TO_REORDER_DEFAULT = "negative" __VOLUME_TO_REORDER_CHOICES = [ "all", "positive", "negative" ] def fill_subparser( subparsers ) -> None: p = subparsers.add_parser( FIX_ELEMENTS_ORDERINGS, help="Reorders the support nodes for the given cell types." ) - for element_name, size in __CELL_NAME_WITH_NUMBER_NODES.items(): - tmp = list( range( size ) ) - random.Random( 4 ).shuffle( tmp ) - p.add_argument( '--' + element_name, - type=str, - metavar=",".join( map( str, tmp ) ), - default=None, - required=False, - help=f"[list of integers]: node permutation for \"{element_name}\"." ) + p.add_argument( '--' + __CELL_NAMES, + type=str, + metavar=", ".join( map( str, __CELL_NAMES_CHOICES ) ), + default=", ".join( map( str, __CELL_NAMES_CHOICES ) ), + help=f"[list of str]: Cell names that can be reordered in your grid. You can use multiple names." + + "Defaults to all cell names being used." ) p.add_argument( '--' + __VOLUME_TO_REORDER, type=str, default=__VOLUME_TO_REORDER_DEFAULT, - metavar=",".join( map( str, __VOLUME_TO_REORDER_CHOICES ) ), - required=True, + metavar=", ".join( __VOLUME_TO_REORDER_CHOICES ), help="[str]: Select which element volume is invalid and needs reordering." + " 'all' will allow reordering of nodes for every element, regarding of their volume." + - " 'positive' or 'negative' will only reorder the element with the corresponding volume." ) + " 'positive' or 'negative' will only reorder the element with the corresponding volume." + + " Defaults to 'negative'." ) vtk_output_parsing.fill_vtk_output_subparser( p ) @@ -45,22 +36,19 @@ def convert( parsed_options ) -> Options: :param options_str: Parsed cli options. :return: Options instance. """ - cell_name_to_ordering: dict[ str, list[ int ] ] = {} - for element_name, size in __CELL_NAME_WITH_NUMBER_NODES.items(): - raw_mapping = parsed_options[ element_name ] - if raw_mapping: - nodes_ordering = tuple( map( int, raw_mapping.split( "," ) ) ) - if not set( nodes_ordering ) == set( range( size ) ): - err_msg: str = f"Permutation {raw_mapping} for type {element_name} is not valid." - logging.error( err_msg ) - raise ValueError( err_msg ) - cell_name_to_ordering[ element_name ] = nodes_ordering + raw_mapping = parsed_options[ __CELL_NAMES ] + cell_names_to_reorder = tuple( raw_mapping.split( "," ) ) + for cell_name in cell_names_to_reorder: + if cell_name not in __CELL_NAMES_CHOICES: + raise ValueError( f"Please choose names between these options for --{__CELL_NAMES_CHOICES}:" + + f" {__CELL_NAMES_CHOICES}." ) vtk_output = vtk_output_parsing.convert( parsed_options ) volume_to_reorder: str = parsed_options[ __VOLUME_TO_REORDER ] if volume_to_reorder.lower() not in __VOLUME_TO_REORDER_CHOICES: - raise ValueError( f"Please use one of these options for --volume_to_reorder: {__VOLUME_TO_REORDER_CHOICES}." ) + raise ValueError( f"Please use one of these options for --{__VOLUME_TO_REORDER}:" + + f" {__VOLUME_TO_REORDER_CHOICES}." ) return Options( vtk_output=vtk_output, - cell_name_to_ordering=cell_name_to_ordering, + cell_names_to_reorder=cell_names_to_reorder, volume_to_reorder=volume_to_reorder ) diff --git a/geos-mesh/tests/test_fix_elements_orderings.py b/geos-mesh/tests/test_fix_elements_orderings.py index 0f1ec40d..060c06f1 100644 --- a/geos-mesh/tests/test_fix_elements_orderings.py +++ b/geos-mesh/tests/test_fix_elements_orderings.py @@ -49,14 +49,15 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int Dict used to apply false nodes orderings for test purposes """ to_change_order: dict[ int, list[ int ] ] = { - "Hexahedron": [ 0, 3, 2, 1, 4, 5, 6, 7 ], - "Tetrahedron": [ 0, 2, 1, 3 ], - "Pyramid": [ 0, 3, 2, 1, 4 ], - "Wedge": [ 0, 2, 1, 3, 4, 5 ], - "Prism5": [ 0, 4, 3, 2, 1, 5, 6, 7, 8, 9 ], - "Prism6": [ 0, 1, 4, 2, 3, 5, 11, 10, 9, 8, 7, 6 ] + "Hexahedron": ( 0, 3, 2, 1, 4, 5, 6, 7 ), + "Tetrahedron": ( 0, 2, 1, 3 ), + "Pyramid": ( 0, 3, 2, 1, 4 ), + "Wedge": ( 0, 2, 1, 3, 4, 5 ), + "Prism5": ( 0, 4, 3, 2, 1, 5, 6, 7, 8, 9 ), + "Prism6": ( 0, 1, 4, 2, 3, 5, 11, 10, 9, 8, 7, 6 ) } to_change_order = dict( sorted( to_change_order.items() ) ) +cell_names = list( VTK_TYPE_TO_NAME.values() ) """ 1 Hexahedron: no invalid ordering """ @@ -707,11 +708,11 @@ def test_is_cell_to_reorder( self ): for grid, needs_ordering in grid_needs_ordering.items(): volumes = feo.compute_mesh_cells_volume( grid ) for i in range( len( volumes ) ): - options = opt( out, to_change_order, "negative" ) + options = opt( out, cell_names, "negative" ) assert feo.is_cell_to_reorder( volumes[ i ], options ) == needs_ordering[ i ] - options = opt( out, to_change_order, "positive" ) + options = opt( out, cell_names, "positive" ) assert feo.is_cell_to_reorder( volumes[ i ], options ) != needs_ordering[ i ] - options = opt( out, to_change_order, "all" ) + options = opt( out, cell_names, "all" ) assert feo.is_cell_to_reorder( volumes[ i ], options ) == True def test_get_cell_types_and_number( self ): @@ -738,7 +739,7 @@ def test_get_cell_types_and_number( self ): feo.get_cell_types_and_number( voxels_grid ) def test_reorder_nodes_to_new_mesh( self ): - options = opt( out, to_change_order, "negative" ) + options = opt( out, cell_names, "negative" ) # single element grids except voxels because it is an invalid cell type for GEOS grid_cell_type = { hexahedrons_grid_invalid: VTK_HEXAHEDRON, @@ -779,14 +780,9 @@ def test_fix_elements_orderings_execution( self ): write_mesh( mix_grid_invalid, test_file ) invalidTest = False command = [ - "python", MESH_DOCTOR_FILEPATH, "-v", "-i", test_file.output, "fix_elements_orderings", "--Hexahedron", - str( to_change_order[ "Hexahedron" ] ).replace( "[", "" ).replace( "]", "" ), "--Tetrahedron", - str( to_change_order[ "Tetrahedron" ] ).replace( "[", "" ).replace( "]", "" ), "--Pyramid", - str( to_change_order[ "Pyramid" ] ).replace( "[", "" ).replace( "]", "" ), "--Wedge", - str( to_change_order[ "Wedge" ] ).replace( "[", "" ).replace( "]", "" ), "--Prism5", - str( to_change_order[ "Prism5" ] ).replace( "[", "" ).replace( "]", "" ), "--Prism6", - str( to_change_order[ "Prism6" ] ).replace( "[", "" ).replace( "]", "" ), "--volume_to_reorder", "negative", - "--data-mode", "binary", "--output", filepath_reordered_mesh + "python", MESH_DOCTOR_FILEPATH, "-v", "-i", test_file.output, "fix_elements_orderings", "--cell_names", + ",".join( map( str, cell_names ) ), "--volume_to_reorder", "negative", "--data-mode", "binary", "--output", + filepath_reordered_mesh ] try: result = subprocess.run( command, shell=True, stderr=subprocess.PIPE, universal_newlines=True ) From d35509df878cdfed75cd646c2b244e142293d3b5 Mon Sep 17 00:00:00 2001 From: alexbenedicto Date: Thu, 19 Sep 2024 19:24:24 -0700 Subject: [PATCH 4/5] Node ordering correction is automatic, user does not need to enter the node ordering himself. Just needs to specify the name of elements to check. --- docs/geos-mesh.rst | 38 +++---- .../doctor/checks/fix_elements_orderings.py | 107 +++++++++++++++--- .../src/geos/mesh/doctor/checks/vtk_utils.py | 6 +- .../parsing/fix_elements_orderings_parsing.py | 56 ++++----- .../tests/test_fix_elements_orderings.py | 32 +++--- 5 files changed, 148 insertions(+), 91 deletions(-) diff --git a/docs/geos-mesh.rst b/docs/geos-mesh.rst index 739d4936..243592ed 100644 --- a/docs/geos-mesh.rst +++ b/docs/geos-mesh.rst @@ -117,41 +117,35 @@ Cells with negative volumes will typically be an issue for ``geos`` and should b """""""""""""""""""""""""" It sometimes happens that an exported mesh does not abide by the ``vtk`` orderings. -The ``fix_elements_orderings`` module can rearrange the nodes of given types of elements. +The ``fix_elements_orderings`` module can rearrange the nodes of given types of elements to respect VTK convention. This can be convenient if you cannot regenerate the mesh. .. code-block:: $ python src/geos/mesh/doctor/mesh_doctor.py fix_elements_orderings --help - usage: mesh_doctor.py fix_elements_orderings [-h] [--Tetrahedron 2,0,3,1] [--Pyramid 3,4,0,2,1] [--Wedge 3,5,4,0,2,1] - [--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] - --volume_to_reorder all,positive,negative --output OUTPUT [--data-mode binary, ascii] + usage: mesh_doctor.py fix_elements_orderings [-h] [--cell_names Hexahedron, Tetrahedron, Pyramid, Wedge, Prism5, Prism6] + [--volume_to_reorder all, positive, negative] + --output OUTPUT [--data-mode binary, ascii] options: - -h, --help show this help message and exit - --Tetrahedron 2,0,3,1 [list of integers]: node permutation for "Tetrahedron". - --Pyramid 3,4,0,2,1 [list of integers]: node permutation for "Pyramid". - --Wedge 3,5,4,0,2,1 [list of integers]: node permutation for "Wedge". - --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". - --volume_to_reorder all,positive,negative - [str]: Select which element volume is invalid and needs reordering. 'all' will allow reordering of nodes for every element, regarding of their volume. - 'positive' or 'negative' will only reorder the element with the corresponding volume. - --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. + -h, --help show this help message and exit + --cell_names Hexahedron, Tetrahedron, Pyramid, Wedge, Prism5, Prism6 + [list of str]: Cell names that can be reordered in your grid. You can use multiple names.Defaults to all cell names being used. + --volume_to_reorder all, positive, negative + [str]: Select which element volume is invalid and needs reordering. 'all' will allow reordering of nodes for every element, regarding of their volume. + 'positive' or 'negative' will only reorder the element with the corresponding volume. Defaults to 'negative'. + --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. For example, assume that you have a mesh that contains 3 different cell types like Hexahedrons, Pyramids and Tetrahedrons. After checking ``element_volumes``, you found that all your Hexahedrons and half of your Tetrahedrons have a negative volume. -To correct that, assuming you know the correct node ordering to correct these volumes, you can use this command: +To correct that, you can use this command: .. code-block:: - $ python src/geos/mesh/doctor/mesh_doctor.py -i /path/to/your/mesh/file fix_elements_orderings --Hexahedron 1,6,5,4,7,0,2,3 - --Tetrahedron 2,0,3,1 --volume_to_reorder negative --output /new/path/for/your/new/mesh/reordered/file --data-mode binary + $ python src/geos/mesh/doctor/mesh_doctor.py -i /path/to/your/mesh/file fix_elements_orderings --cell_names Hexahedron,Tetrahedron + --volume_to_reorder negative --output /new/path/for/your/new/mesh/reordered/file --data-mode binary ``generate_cube`` """"""""""""""""" diff --git a/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py b/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py index bc693a99..988a9197 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py @@ -1,18 +1,18 @@ -from dataclasses import dataclass import numpy as np import logging +from dataclasses import dataclass +from itertools import permutations from vtk import vtkCellSizeFilter -from vtkmodules.vtkCommonCore import vtkIdList from vtkmodules.util.numpy_support import vtk_to_numpy -from vtkmodules.vtkCommonDataModel import ( vtkDataSet, VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, +from vtkmodules.vtkCommonDataModel import ( vtkDataSet, vtkCell, VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM ) -from .vtk_utils import VtkOutput, to_vtk_id_list, write_mesh, read_mesh +from geos.mesh.doctor.checks.vtk_utils import VtkOutput, vtk_iter, to_vtk_id_list, write_mesh, read_mesh @dataclass( frozen=True ) class Options: vtk_output: VtkOutput - cell_name_to_ordering: dict[ str, list[ int ] ] + cell_names_to_reorder: tuple[ str ] volume_to_reorder: str @@ -103,6 +103,90 @@ def is_cell_to_reorder( cell_volume: str, options: Options ) -> bool: return False +def are_face_nodes_counterclockwise( face: any ) -> bool: + """Checks if the nodes of a face are ordered counterclockwise when looking at the plan created by the face. + + Args: + face (any): Face of a vtkCell. + + Returns: + bool: True if counterclockwise, False instead. + """ + face_points = face.GetPoints() + number_points = face_points.GetNumberOfPoints() + if number_points < 3: + err_msg = f"The face has less than 3 nodes which is invalid." + logging.error( err_msg ) + raise ValueError( err_msg ) + # first calculate the normal vector of the face + a = np.array( face_points.GetPoint( 0 ) ) + b = np.array( face_points.GetPoint( 1 ) ) + c = np.array( face_points.GetPoint( 2 ) ) + AB = b - a + AC = c - a + normal = np.cross( AB, AC ) + + # then calculate the vector cross products sum of all points of the face with a random point P + # from discussion https://math.stackexchange.com/questions/2152623/determine-the-order-of-a-3d-polygon + P = np.array( [ 0.0, 0.0, 0.0 ] ) # P position does not change the value of sum + all_points = [ np.array( face_points.GetPoint( v ) ) for v in range( number_points ) ] + vector_sum = np.array( [ 0.0, 0.0, 0.0 ] ) + for i in range( number_points ): + PAi = all_points[ i % number_points ] - P + PAiplus1 = all_points[ ( i + 1 ) % number_points ] - P + vector_sum += np.cross( PAi, PAiplus1 ) + vector_sum = vector_sum / 2 # needs to be half + + # if dot product is positive, the nodes are ordered counterclockwise + if np.dot( vector_sum, normal ) > 0.0: + return True + return False + + +def valid_cell_point_ids_ordering( cell: vtkCell ) -> list[ int ]: + """Returns the valid order of point ids of a cell that respect counterclockwise convention of each face. + + Args: + cell (vtkCell): A cell of vtk grid. + + Raises: + ValueError: "No node ids permutation made the face nodes to be ordered counterclockwise." + + Returns: + list[ int ]: [ pt_id0, pt_id2, pt_id1, ..., pt_idN ] + """ + initial_ids_order: list[ int ] = list( vtk_iter( cell.GetPointIds() ) ) + valid_points_id: list[ int ] = initial_ids_order.copy() + + for face_id in range( cell.GetNumberOfFaces() ): + face = cell.GetFace( face_id ) + if are_face_nodes_counterclockwise( face ): + continue # the face nodes respect the convention so continue to next face + + reordered = False + initial_ids: list[ int ] = [ face.GetPointId( i ) for i in range( face.GetNumberOfPoints() ) ] + initial_coords: list[ float ] = [ face.GetPoints().GetPoint( v ) for v in range( face.GetNumberOfPoints() ) ] + for permutation_ids in permutations( initial_ids ): + for i, node_id in enumerate( permutation_ids ): + face.GetPointIds().SetId( i, node_id ) + face.GetPoints().SetPoint( i, initial_coords[ initial_ids.index( node_id ) ] ) + + if are_face_nodes_counterclockwise( face ): # the correct permutation was found + for j, initial_id in enumerate( initial_ids ): + if initial_id != permutation_ids[ j ]: + valid_points_id[ initial_ids_order.index( initial_id ) ] = permutation_ids[ j ] + reordered = True + break + + if not reordered: + err_msg = ( f"For face with nodes id '{initial_ids}', that corresponds to coordinates '{initial_coords}'" + + ", no node ids permutation made the face nodes to be ordered counterclockwise." ) + logging.error( err_msg ) + raise ValueError( err_msg ) + + return valid_points_id + + def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple: """From an old mesh, creates a new one where certain cell nodes are reordered to obtain a correct volume. @@ -118,7 +202,7 @@ def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple names_with_totals: dict[ str, int ] = { n: v for n, v in zip( unique_cell_names, total_per_cell_types ) } # sorted dict allow for sorted output of reordering_stats names_with_totals_sorted: dict[ str, int ] = dict( sorted( names_with_totals.items() ) ) - useful_VTK_TYPEs: list[ int ] = [ NAME_TO_VTK_TYPE[ name ] for name in options.cell_name_to_ordering.keys() ] + useful_VTK_TYPEs: list[ int ] = [ NAME_TO_VTK_TYPE[ name ] for name in options.cell_names_to_reorder ] all_cells_volume: np.array = compute_mesh_cells_volume( old_mesh ) # a new mesh with the same data is created from the old mesh new_mesh: vtkDataSet = old_mesh.NewInstance() @@ -132,7 +216,7 @@ def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple "Types non reordered": list( names_with_totals_sorted.keys() ), "Number of cells non reordered": list( names_with_totals_sorted.values() ) } - counter_cells_reordered: dict[ str, int ] = { name: 0 for name in options.cell_name_to_ordering.keys() } + counter_cells_reordered: dict[ str, int ] = { name: 0 for name in options.cell_names_to_reorder } # sorted dict allow for sorted output of reordering_stats ounter_cells_reordered_sorted: dict[ str, int ] = dict( sorted( counter_cells_reordered.items() ) ) # Reordering operations @@ -141,13 +225,8 @@ def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple if vtk_type in useful_VTK_TYPEs: if is_cell_to_reorder( float( all_cells_volume[ cell_id ] ), options ): cell_name: str = VTK_TYPE_TO_NAME[ vtk_type ] - support_point_ids = vtkIdList() - cells.GetCellAtId( cell_id, support_point_ids ) - new_support_point_ids: list[ int ] = list() - node_ordering: list[ int ] = options.cell_name_to_ordering[ cell_name ] - for i in range( len( node_ordering ) ): - new_support_point_ids.append( support_point_ids.GetId( node_ordering[ i ] ) ) - cells.ReplaceCellAtId( cell_id, to_vtk_id_list( new_support_point_ids ) ) + point_ids_ordering: list[ int ] = valid_cell_point_ids_ordering( new_mesh.GetCell( cell_id ) ) + cells.ReplaceCellAtId( cell_id, to_vtk_id_list( point_ids_ordering ) ) ounter_cells_reordered_sorted[ cell_name ] += 1 # Calculation of stats for cell_name_reordered, amount in ounter_cells_reordered_sorted.items(): diff --git a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py index 8f52cbb8..5d384049 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py @@ -103,7 +103,7 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid: :return: A unstructured grid. """ if not is_filepath_valid( vtk_input_file ): - err_msg: str = f"Invalid file path. Could not read \"{vtk_input_file}\". Dying..." + err_msg: str = f"Invalid file path. Could not read \"{vtk_input_file}\"." logging.error( err_msg ) raise ValueError( err_msg ) file_extension = os.path.splitext( vtk_input_file )[ -1 ] @@ -119,7 +119,7 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid: if output_mesh: return output_mesh # No reader did work. Dying. - err_msg = f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\". Dying..." + err_msg = f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\"." logging.error( err_msg ) raise ValueError( err_msg ) @@ -160,7 +160,7 @@ def write_mesh( mesh: vtkUnstructuredGrid, vtk_output: VtkOutput ) -> int: success_code = __write_vtu( mesh, vtk_output.output, vtk_output.is_data_mode_binary ) else: # No writer found did work. Dying. - err_msg = f"Could not find the appropriate VTK writer for extension \"{file_extension}\". Dying..." + err_msg = f"Could not find the appropriate VTK writer for extension \"{file_extension}\"." logging.error( err_msg ) raise ValueError( err_msg ) return 0 if success_code else 2 # the Write member function return 1 in case of success, 0 otherwise. diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py index 641a9505..8fe86906 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py +++ b/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py @@ -1,41 +1,32 @@ import logging import random -from geos.mesh.doctor.checks.fix_elements_orderings import Options, Result +from geos.mesh.doctor.checks.fix_elements_orderings import Options, Result, NAME_TO_VTK_TYPE from . import vtk_output_parsing, FIX_ELEMENTS_ORDERINGS -__CELL_NAME_WITH_NUMBER_NODES = { - "Tetrahedron": 4, - "Pyramid": 5, - "Wedge": 6, - "Hexahedron": 8, - "Prism5": 10, - "Prism6": 12 -} +__CELL_NAMES = "cell_names" +__CELL_NAMES_CHOICES = list( NAME_TO_VTK_TYPE.keys() ) __VOLUME_TO_REORDER = "volume_to_reorder" -__VOLUME_TO_REORDER_DEFAULT = "all" +__VOLUME_TO_REORDER_DEFAULT = "negative" __VOLUME_TO_REORDER_CHOICES = [ "all", "positive", "negative" ] def fill_subparser( subparsers ) -> None: p = subparsers.add_parser( FIX_ELEMENTS_ORDERINGS, help="Reorders the support nodes for the given cell types." ) - for element_name, size in __CELL_NAME_WITH_NUMBER_NODES.items(): - tmp = list( range( size ) ) - random.Random( 4 ).shuffle( tmp ) - p.add_argument( '--' + element_name, - type=str, - metavar=",".join( map( str, tmp ) ), - default=None, - required=False, - help=f"[list of integers]: node permutation for \"{element_name}\"." ) + p.add_argument( '--' + __CELL_NAMES, + type=str, + metavar=", ".join( map( str, __CELL_NAMES_CHOICES ) ), + default=", ".join( map( str, __CELL_NAMES_CHOICES ) ), + help=f"[list of str]: Cell names that can be reordered in your grid. You can use multiple names." + + "Defaults to all cell names being used." ) p.add_argument( '--' + __VOLUME_TO_REORDER, type=str, default=__VOLUME_TO_REORDER_DEFAULT, - metavar=",".join( map( str, __VOLUME_TO_REORDER_CHOICES ) ), - required=True, + metavar=", ".join( __VOLUME_TO_REORDER_CHOICES ), help="[str]: Select which element volume is invalid and needs reordering." + " 'all' will allow reordering of nodes for every element, regarding of their volume." + - " 'positive' or 'negative' will only reorder the element with the corresponding volume." ) + " 'positive' or 'negative' will only reorder the element with the corresponding volume." + + " Defaults to 'negative'." ) vtk_output_parsing.fill_vtk_output_subparser( p ) @@ -45,22 +36,19 @@ def convert( parsed_options ) -> Options: :param options_str: Parsed cli options. :return: Options instance. """ - cell_name_to_ordering: dict[ str, list[ int ] ] = {} - for element_name, size in __CELL_NAME_WITH_NUMBER_NODES.items(): - raw_mapping = parsed_options[ element_name ] - if raw_mapping: - nodes_ordering = tuple( map( int, raw_mapping.split( "," ) ) ) - if not set( nodes_ordering ) == set( range( size ) ): - err_msg: str = f"Permutation {raw_mapping} for type {element_name} is not valid." - logging.error( err_msg ) - raise ValueError( err_msg ) - cell_name_to_ordering[ element_name ] = nodes_ordering + raw_mapping = parsed_options[ __CELL_NAMES ] + cell_names_to_reorder = tuple( raw_mapping.split( "," ) ) + for cell_name in cell_names_to_reorder: + if cell_name not in __CELL_NAMES_CHOICES: + raise ValueError( f"Please choose names between these options for --{__CELL_NAMES_CHOICES}:" + + f" {__CELL_NAMES_CHOICES}." ) vtk_output = vtk_output_parsing.convert( parsed_options ) volume_to_reorder: str = parsed_options[ __VOLUME_TO_REORDER ] if volume_to_reorder.lower() not in __VOLUME_TO_REORDER_CHOICES: - raise ValueError( f"Please use one of these options for --volume_to_reorder: {__VOLUME_TO_REORDER_CHOICES}." ) + raise ValueError( f"Please use one of these options for --{__VOLUME_TO_REORDER}:" + + f" {__VOLUME_TO_REORDER_CHOICES}." ) return Options( vtk_output=vtk_output, - cell_name_to_ordering=cell_name_to_ordering, + cell_names_to_reorder=cell_names_to_reorder, volume_to_reorder=volume_to_reorder ) diff --git a/geos-mesh/tests/test_fix_elements_orderings.py b/geos-mesh/tests/test_fix_elements_orderings.py index 0f1ec40d..060c06f1 100644 --- a/geos-mesh/tests/test_fix_elements_orderings.py +++ b/geos-mesh/tests/test_fix_elements_orderings.py @@ -49,14 +49,15 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int Dict used to apply false nodes orderings for test purposes """ to_change_order: dict[ int, list[ int ] ] = { - "Hexahedron": [ 0, 3, 2, 1, 4, 5, 6, 7 ], - "Tetrahedron": [ 0, 2, 1, 3 ], - "Pyramid": [ 0, 3, 2, 1, 4 ], - "Wedge": [ 0, 2, 1, 3, 4, 5 ], - "Prism5": [ 0, 4, 3, 2, 1, 5, 6, 7, 8, 9 ], - "Prism6": [ 0, 1, 4, 2, 3, 5, 11, 10, 9, 8, 7, 6 ] + "Hexahedron": ( 0, 3, 2, 1, 4, 5, 6, 7 ), + "Tetrahedron": ( 0, 2, 1, 3 ), + "Pyramid": ( 0, 3, 2, 1, 4 ), + "Wedge": ( 0, 2, 1, 3, 4, 5 ), + "Prism5": ( 0, 4, 3, 2, 1, 5, 6, 7, 8, 9 ), + "Prism6": ( 0, 1, 4, 2, 3, 5, 11, 10, 9, 8, 7, 6 ) } to_change_order = dict( sorted( to_change_order.items() ) ) +cell_names = list( VTK_TYPE_TO_NAME.values() ) """ 1 Hexahedron: no invalid ordering """ @@ -707,11 +708,11 @@ def test_is_cell_to_reorder( self ): for grid, needs_ordering in grid_needs_ordering.items(): volumes = feo.compute_mesh_cells_volume( grid ) for i in range( len( volumes ) ): - options = opt( out, to_change_order, "negative" ) + options = opt( out, cell_names, "negative" ) assert feo.is_cell_to_reorder( volumes[ i ], options ) == needs_ordering[ i ] - options = opt( out, to_change_order, "positive" ) + options = opt( out, cell_names, "positive" ) assert feo.is_cell_to_reorder( volumes[ i ], options ) != needs_ordering[ i ] - options = opt( out, to_change_order, "all" ) + options = opt( out, cell_names, "all" ) assert feo.is_cell_to_reorder( volumes[ i ], options ) == True def test_get_cell_types_and_number( self ): @@ -738,7 +739,7 @@ def test_get_cell_types_and_number( self ): feo.get_cell_types_and_number( voxels_grid ) def test_reorder_nodes_to_new_mesh( self ): - options = opt( out, to_change_order, "negative" ) + options = opt( out, cell_names, "negative" ) # single element grids except voxels because it is an invalid cell type for GEOS grid_cell_type = { hexahedrons_grid_invalid: VTK_HEXAHEDRON, @@ -779,14 +780,9 @@ def test_fix_elements_orderings_execution( self ): write_mesh( mix_grid_invalid, test_file ) invalidTest = False command = [ - "python", MESH_DOCTOR_FILEPATH, "-v", "-i", test_file.output, "fix_elements_orderings", "--Hexahedron", - str( to_change_order[ "Hexahedron" ] ).replace( "[", "" ).replace( "]", "" ), "--Tetrahedron", - str( to_change_order[ "Tetrahedron" ] ).replace( "[", "" ).replace( "]", "" ), "--Pyramid", - str( to_change_order[ "Pyramid" ] ).replace( "[", "" ).replace( "]", "" ), "--Wedge", - str( to_change_order[ "Wedge" ] ).replace( "[", "" ).replace( "]", "" ), "--Prism5", - str( to_change_order[ "Prism5" ] ).replace( "[", "" ).replace( "]", "" ), "--Prism6", - str( to_change_order[ "Prism6" ] ).replace( "[", "" ).replace( "]", "" ), "--volume_to_reorder", "negative", - "--data-mode", "binary", "--output", filepath_reordered_mesh + "python", MESH_DOCTOR_FILEPATH, "-v", "-i", test_file.output, "fix_elements_orderings", "--cell_names", + ",".join( map( str, cell_names ) ), "--volume_to_reorder", "negative", "--data-mode", "binary", "--output", + filepath_reordered_mesh ] try: result = subprocess.run( command, shell=True, stderr=subprocess.PIPE, universal_newlines=True ) From ec38eb821bbf4ed8e42bd9105b4f3ad5c791bfea Mon Sep 17 00:00:00 2001 From: alexbenedicto Date: Fri, 11 Oct 2024 11:18:45 -0700 Subject: [PATCH 5/5] Automation of reordering is redefined with limits to its scope. Degenerated cells cannot be treated because the problem lies more in the generation of the mesh itself than with a convention error. --- .../doctor/checks/fix_elements_orderings.py | 678 ++++++++++++---- .../parsing/fix_elements_orderings_parsing.py | 48 +- .../tests/test_fix_elements_orderings.py | 750 ++++++++++++++---- 3 files changed, 1186 insertions(+), 290 deletions(-) diff --git a/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py b/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py index be26164a..5550768c 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py @@ -1,12 +1,14 @@ import numpy as np import logging +from math import sqrt from dataclasses import dataclass from itertools import permutations -from vtk import vtkCellSizeFilter +from typing import TypeAlias +from vtk import vtkCellSizeFilter, vtkIdList from vtkmodules.util.numpy_support import vtk_to_numpy from vtkmodules.vtkCommonDataModel import ( vtkDataSet, vtkCell, VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM ) -from geos.mesh.doctor.checks.vtk_utils import VtkOutput, vtk_iter, to_vtk_id_list, write_mesh, read_mesh +from geos.mesh.doctor.checks.vtk_utils import VtkOutput, to_vtk_id_list, write_mesh, read_mesh @dataclass( frozen=True ) @@ -22,7 +24,19 @@ class Result: reordering_stats: dict[ str, list[ int ] ] +Coordinates3D: TypeAlias = tuple[ float, float, float ] +Points3D: TypeAlias = tuple[ Coordinates3D ] +NodeOrdering: TypeAlias = tuple[ int ] + GEOS_ACCEPTED_TYPES = [ VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM ] +VTK_TYPES_NUMBER_POINTS: dict[ int, tuple[ int, str ] ] = { + VTK_TETRA: ( 4, "Tetrahedron" ), + VTK_HEXAHEDRON: ( 8, "Hexahedron" ), + VTK_PYRAMID: ( 5, "Pyramid" ), + VTK_WEDGE: ( 6, "Wedge" ), + VTK_PENTAGONAL_PRISM: ( 10, "Pentagonal prism" ), + VTK_HEXAGONAL_PRISM: ( 12, "Hexagonal prism" ) +} # the number of different nodes that needs to be entered in parsing when dealing with a specific vtk element NAME_TO_VTK_TYPE = { "Hexahedron": VTK_HEXAHEDRON, @@ -55,32 +69,6 @@ def compute_mesh_cells_volume( mesh: vtkDataSet ) -> np.array: return vtk_to_numpy( cell_size_filter.GetOutput().GetCellData().GetArray( "Volume" ) ) -def get_cell_types_and_number( mesh: vtkDataSet ) -> tuple[ list[ int ] ]: - """Gets the cell type for every cell of a mesh and the amount for each cell type. - - Args: - mesh (vtkDataSet): A vtk grid. - - Raises: - ValueError: "Invalid type '{cell_type}' for GEOS is in the mesh. Dying ..." - - Returns: - tuple[ list[ int ] ]: ( unique_cell_types, total_per_cell_types ) - """ - number_cells: int = mesh.GetNumberOfCells() - all_cells_type: np.array = np.ones( number_cells, dtype=int ) - for cell_id in range( number_cells ): - all_cells_type[ cell_id ] = mesh.GetCellType( cell_id ) - unique_cell_types, total_per_cell_types = np.unique( all_cells_type, return_counts=True ) - unique_cell_types, total_per_cell_types = unique_cell_types.tolist(), total_per_cell_types.tolist() - for cell_type in unique_cell_types: - if cell_type not in GEOS_ACCEPTED_TYPES: - err_msg: str = f"Invalid type '{cell_type}' for GEOS is in the mesh. Dying ..." - logging.error( err_msg ) - raise ValueError( err_msg ) - return ( unique_cell_types, total_per_cell_types ) - - def is_cell_to_reorder( cell_volume: str, options: Options ) -> bool: """Check if the volume of vtkCell qualifies the cell to be reordered. @@ -103,146 +91,562 @@ def is_cell_to_reorder( cell_volume: str, options: Options ) -> bool: return False -def are_face_nodes_counterclockwise( face: any ) -> bool: - """Checks if the nodes of a face are ordered counterclockwise when looking at the plan created by the face. +def cell_ids_to_reorder_by_cell_type( mesh: vtkDataSet, options: Options ) -> dict[ int, np.array ]: + """Create an array of cell_ids for each vtk_type chosen by the user when the cell pointed by the cell_id + needs to be reorder. Args: - face (any): Face of a vtkCell. + mesh (vtkDataSet): A vtk grid. + options (Options): Options defined by the user. Returns: - bool: True if counterclockwise, False instead. + dict[ int, np.array ]: { vtk_type1: np.array, ..., vtk_typeN: np.array } """ - face_points = face.GetPoints() - number_points = face_points.GetNumberOfPoints() - if number_points < 3: - err_msg = f"The face has less than 3 nodes which is invalid." - logging.error( err_msg ) - raise ValueError( err_msg ) - # first calculate the normal vector of the face - a = np.array( face_points.GetPoint( 0 ) ) - b = np.array( face_points.GetPoint( 1 ) ) - c = np.array( face_points.GetPoint( 2 ) ) - AB = b - a - AC = c - a - normal = np.cross( AB, AC ) - - # then calculate the vector cross products sum of all points of the face with a random point P - # from discussion https://math.stackexchange.com/questions/2152623/determine-the-order-of-a-3d-polygon - P = np.array( [ 0.0, 0.0, 0.0 ] ) # P position does not change the value of sum - all_points = [ np.array( face_points.GetPoint( v ) ) for v in range( number_points ) ] - vector_sum = np.array( [ 0.0, 0.0, 0.0 ] ) - for i in range( number_points ): - PAi = all_points[ i % number_points ] - P - PAiplus1 = all_points[ ( i + 1 ) % number_points ] - P - vector_sum += np.cross( PAi, PAiplus1 ) - vector_sum = vector_sum / 2 # needs to be half - - # if dot product is positive, the nodes are ordered counterclockwise - if np.dot( vector_sum, normal ) > 0.0: - return True - return False + all_cells_volume: np.array = compute_mesh_cells_volume( mesh ) + number_cells: int = mesh.GetNumberOfCells() + useful_VTK_TYPEs: list[ int ] = [ NAME_TO_VTK_TYPE[ name ] for name in options.cell_names_to_reorder ] + all_cells: np.array = np.zeros( ( number_cells, 3 ), dtype=int ) # col0: cell_ids, col1: vtk_type, col2: 0 or 1 + for cell_id in range( number_cells ): + vtk_type: int = mesh.GetCellType( cell_id ) + all_cells[ cell_id ][ 0 ] = cell_id + all_cells[ cell_id ][ 1 ] = vtk_type + if vtk_type in useful_VTK_TYPEs: + if is_cell_to_reorder( float( all_cells_volume[ cell_id ] ), options ): + all_cells[ cell_id ][ 2 ] = 1 + continue + all_cells[ cell_id ][ 2 ] = 0 + to_reorder: np.array = all_cells[ :, 2 ] == 1 # We need to remove rows where col2 == 0 + cells_to_reorder: np.array = all_cells[ to_reorder ][ :, [ 0, 1 ] ] # col0: cell_ids, col1: vtk_type + unique_values: np.array = np.unique( cells_to_reorder[ :, 1 ] ) + # Create a dictionary to store arrays of cells_id to reorder for each unique value + cells_to_reorder_by_cell_type: dict[ int, np.array ] = { + int( value ): cells_to_reorder[ cells_to_reorder[ :, 1 ] == value ][ :, 0 ] + for value in unique_values + } + return cells_to_reorder_by_cell_type -def valid_cell_point_ids_ordering( cell: vtkCell ) -> list[ int ]: - """Returns the valid order of point ids of a cell that respect counterclockwise convention of each face. +def is_polygon_counterclockwise( points: Points3D, ref_pt: np.array ) -> bool: + """Determines if a set of points in 3D are being ordered counterclockwise or clockwise + with respect to a reference point of observation. Args: - cell (vtkCell): A cell of vtk grid. + points (Points3D): A set of points in 3D coordinates. + ref_pt (np.array): A point in 3D coordinates. Raises: - ValueError: "No node ids permutation made the face nodes to be ordered counterclockwise." + ValueError: "Polygon checked is concave with points: {points}" Returns: - list[ int ]: [ pt_id0, pt_id2, pt_id1, ..., pt_idN ] + bool: True if counterclockwise, False if clockwise. """ - initial_ids_order: list[ int ] = list( vtk_iter( cell.GetPointIds() ) ) - valid_points_id: list[ int ] = initial_ids_order.copy() - - for face_id in range( cell.GetNumberOfFaces() ): - face = cell.GetFace( face_id ) - if are_face_nodes_counterclockwise( face ): - continue # the face nodes respect the convention so continue to next face - - reordered = False - initial_ids: list[ int ] = [ face.GetPointId( i ) for i in range( face.GetNumberOfPoints() ) ] - initial_coords: list[ float ] = [ face.GetPoints().GetPoint( v ) for v in range( face.GetNumberOfPoints() ) ] - for permutation_ids in permutations( initial_ids ): - for i, node_id in enumerate( permutation_ids ): - face.GetPointIds().SetId( i, node_id ) - face.GetPoints().SetPoint( i, initial_coords[ initial_ids.index( node_id ) ] ) - - if are_face_nodes_counterclockwise( face ): # the correct permutation was found - for j, initial_id in enumerate( initial_ids ): - if initial_id != permutation_ids[ j ]: - valid_points_id[ initial_ids_order.index( initial_id ) ] = permutation_ids[ j ] - reordered = True - break + nbr_pts: int = len( points ) + assert nbr_pts > 2 + points_array = np.array( points ) + polygon_centroid = np.mean( points_array, axis=0 ) + towards_ref_vect = ref_pt - polygon_centroid + # Determine the projection plane + abs_vect = np.abs( towards_ref_vect ) + if abs_vect[ 0 ] > abs_vect[ 1 ] and abs_vect[ 0 ] > abs_vect[ 2 ]: + # Project onto the YZ-plane + projected_points: list[ Coordinates3D ] = [ list( points[ i ] )[ 1: ] for i in range( nbr_pts ) ] + is_sign_positive: bool = True if towards_ref_vect[ 0 ] < 0.0 else False + elif abs_vect[ 1 ] > abs_vect[ 0 ] and abs_vect[ 1 ] > abs_vect[ 2 ]: + # Project onto the XZ-plane + projected_points = [ [ points[ i ][ 0 ], points[ i ][ 2 ] ] for i in range( nbr_pts ) ] + is_sign_positive = True if towards_ref_vect[ 1 ] > 0.0 else False + elif abs_vect[ 2 ] > abs_vect[ 0 ] and abs_vect[ 2 ] > abs_vect[ 1 ]: + # Project onto the XY-plane + projected_points = [ list( points[ i ] )[ :2 ] for i in range( nbr_pts ) ] + is_sign_positive = True if towards_ref_vect[ 2 ] < 0.0 else False + elif abs_vect[ 0 ] > abs_vect[ 2 ] and abs_vect[ 1 ] > abs_vect[ 2 ] and abs_vect[ 0 ] == abs_vect[ 1 ]: + # Project onto the XZ-plane + projected_points = [ [ points[ i ][ 0 ], points[ i ][ 2 ] ] for i in range( nbr_pts ) ] + is_sign_positive = True if towards_ref_vect[ 1 ] > 0.0 else False + elif abs_vect[ 0 ] > abs_vect[ 1 ] and abs_vect[ 2 ] > abs_vect[ 1 ] and abs_vect[ 0 ] == abs_vect[ 2 ]: + # Project onto the YZ-plane + projected_points = [ list( points[ i ] )[ 1: ] for i in range( nbr_pts ) ] + is_sign_positive = True if towards_ref_vect[ 0 ] < 0.0 else False + elif abs_vect[ 1 ] > abs_vect[ 0 ] and abs_vect[ 2 ] > abs_vect[ 0 ] and abs_vect[ 1 ] == abs_vect[ 2 ]: + # Project onto the XY-plane + projected_points = [ list( points[ i ] )[ :2 ] for i in range( nbr_pts ) ] + is_sign_positive = True if towards_ref_vect[ 2 ] < 0.0 else False + + # translate the projected points to be with positive coordinates for det calculation + min_x = min( [ x for x, y in projected_points ] ) + min_y = min( [ y for x, y in projected_points ] ) + if min_x < 0.0: + for i in range( nbr_pts ): + projected_points[ i ][ 0 ] += abs( min_x ) + if min_y < 0.0: + for i in range( nbr_pts ): + projected_points[ i ][ 1 ] += abs( min_y ) + + turning_direction = None + for v in range( nbr_pts ): + A = projected_points[ v ] + B = projected_points[ ( v + 1 ) % nbr_pts ] + C = projected_points[ ( v + 2 ) % nbr_pts ] + det = ( B[ 0 ] - A[ 0 ] ) * ( C[ 1 ] - A[ 1 ] ) - ( C[ 0 ] - A[ 0 ] ) * ( B[ 1 ] - A[ 1 ] ) + if det != 0: + current_turn = 'ccw' if det > 0 else 'cw' + if turning_direction is None: + turning_direction = current_turn + elif turning_direction != current_turn: + raise ValueError( f"Polygon checked is concave with points: {points}." ) + if turning_direction == 'cw' and is_sign_positive: + return True + elif turning_direction == 'ccw' and not is_sign_positive: + return True + return False - if not reordered: - err_msg = ( f"For face with nodes id '{initial_ids}', that corresponds to coordinates '{initial_coords}'" + - ", no node ids permutation made the face nodes to be ordered counterclockwise." ) - logging.error( err_msg ) - raise ValueError( err_msg ) - return valid_points_id +def choose_correct_polygon( polygons: list[ Points3D ] ) -> Points3D: + """When choosing between polygons, whose points are ordered counterclockwise (ccw) or clockwise (cw) + AND that have more than 4 points, most of them will have an inappropriate shape that can generate concavity while + still being ordered correctly ccw / cw (like with a set of 5 points supposed to look like a pentagon but end being + the shape of a star). + Therefore, to choose the correct polygon, we calculate the distance between each pair of "adjacent" points and try + to minimize it to only choose the appropriate shape. + Args: + polygons (list[ Points3D ]): All possible permutations of a set of 5 points forming a ccw / cw polygon. -def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple: - """From an old mesh, creates a new one where certain cell nodes are reordered to obtain a correct volume. + Returns: + Points3D: The correct polygon. + """ + min_sum_distance: float = 1e10 + min_index: int = 0 + for i, polygon_pts in enumerate( polygons ): + sum_distance: float = 0.0 + number_points: int = len( polygon_pts ) + for j in range( number_points ): + A, B = polygon_pts[ j ], polygon_pts[ ( j + 1 ) % number_points ] + sum_distance += sqrt( ( B[ 0 ] - A[ 0 ] )**2 + ( B[ 1 ] - A[ 1 ] )**2 + ( B[ 2 ] - A[ 2 ] )**2 ) + if sum_distance < min_sum_distance: + min_sum_distance = sum_distance + min_index = i + return polygons[ min_index ] + + +def check_points_to_reorder( vtk_type: int ): + """Check that the input tuple of points for any reorder function are of the right type, + the correct number of points wrt to a VTK_TYPE and do not contain duplicates. Args: - old_mesh (vtkDataSet): A vtk grid needing nodes to be reordered. - options (Options): Options defined by the user. + vtk_type (int): VTK + """ + + def decorator( reorder_func ): + + def wrapper( arg ): + # Perform a check of the points to be reordered + if not isinstance( arg, tuple ): + raise ValueError( f"Points {arg} are not a tuple." ) + elif not isinstance( arg[ 0 ], tuple ): + raise ValueError( f"Points {arg} are not a tuple[ tuple ]." ) + elif not isinstance( arg[ 0 ][ 0 ], float ): + raise ValueError( f"Points {arg} are not a tuple[ tuple[ float ] ]." ) + elif len( arg[ 0 ] ) != 3: + raise ValueError( f"Points {arg} coordinates are not in 3D." ) + + try: + number_points, element_name = VTK_TYPES_NUMBER_POINTS[ vtk_type ] + except KeyError: + raise ValueError( f"This VTK_TYPE '{vtk_type}' is not available for fix_elements_reorderings." ) + if len( set( arg ) ) != number_points: + raise ValueError( f"Duplicated points were found in the cell {element_name} with points '{arg}'." + + " Or you meant to use a reordering function for another type of VTK cell." + + " Cannot perform reordering in this condition." ) + + # Call the original function + return reorder_func( arg ) + + return wrapper + + return decorator + + +@check_points_to_reorder( VTK_TETRA ) +def reorder_tetrahedron( points: Points3D ) -> Points3D: + points_array: np.array = np.array( points ) + centroid: np.array = np.mean( points_array, axis=0 ) + face0_pts: Points3D = ( points[ 0 ], points[ 1 ], points[ 3 ] ) # face 0 in convention + face1_pts: Points3D = ( points[ 1 ], points[ 2 ], points[ 3 ] ) # face 1 in convention + is_face0_ccw: bool = is_polygon_counterclockwise( face0_pts, centroid ) + is_face1_ccw: bool = is_polygon_counterclockwise( face1_pts, centroid ) + if not is_face0_ccw and not is_face1_ccw: + return points + else: + perms = list( permutations( points ) )[ 1: ] # first permutation is not different from input, no need to use it + for perm in perms: + face0_pts = ( perm[ 0 ], perm[ 1 ], perm[ 3 ] ) + face1_pts = ( perm[ 1 ], perm[ 2 ], perm[ 3 ] ) + is_face0_ccw = is_polygon_counterclockwise( face0_pts, centroid ) + is_face1_ccw = is_polygon_counterclockwise( face1_pts, centroid ) + if is_face0_ccw or is_face1_ccw: + continue + else: + correct_ordering: Points3D = ( face0_pts[ 0 ], face0_pts[ 1 ], face1_pts[ 1 ], face1_pts[ 2 ] ) + if len( set( correct_ordering ) ) == 4: + return correct_ordering + raise ValueError( "Error when reordering the tetrahedron." ) + + +# BIG HYPOTHESIS: The first 4 points of the pyramid represent its basis, while the 5th represent its apex. +@check_points_to_reorder( VTK_PYRAMID ) +def reorder_pyramid( points: Points3D ) -> Points3D: + points_array: np.array = np.array( points ) + centroid: np.array = np.mean( points_array, axis=0 ) + quad_pts: Points3D = ( points[ 0 ], points[ 1 ], points[ 2 ], points[ 3 ] ) # face 0 in convention + # first we verify the verify the hypothesis that the 4 first points of the pyramid represent its basis + try: + is_quad_ccw: bool = is_polygon_counterclockwise( quad_pts, centroid ) + except ValueError: + raise ValueError( "The first 4 points of your pyramid do not represent its base. No ordering possible." ) + + if is_quad_ccw: + return points + else: + perms = list( permutations( points[ :4 ] ) )[ 1: ] + for perm in perms: + try: + is_quad_ccw = is_polygon_counterclockwise( perm, centroid ) + except ValueError: # polygon created was concave + continue + if not is_quad_ccw: + continue + else: + correct_ordering: Points3D = perm + ( points[ 4 ], ) + return correct_ordering + raise ValueError( "Error when reordering the pyramid." ) + + +# BIG HYPOTHESIS: The first 3 points define the first triangle face and the next 3 define the other triangle face. +# If it is not the case, the volume of the element created will be negative +@check_points_to_reorder( VTK_WEDGE ) +def reorder_wedge( points: Points3D ) -> Points3D: + points_array: np.array = np.array( points ) + centroid: np.array = np.mean( points_array, axis=0 ) + # we check that the two supposed triangle basis are oriented the same wrt to the centroid of the wedge + initial_triangle0_pts: Points3D = ( points[ 0 ], points[ 1 ], points[ 2 ] ) # face 0 in convention + initial_triangle1_pts: Points3D = ( points[ 3 ], points[ 4 ], points[ 5 ] ) # face 1 in convention + is_triangle0_ccw: bool = is_polygon_counterclockwise( initial_triangle0_pts, centroid ) + is_triangle1_ccw: bool = is_polygon_counterclockwise( initial_triangle1_pts, centroid ) + if is_triangle0_ccw == is_triangle1_ccw: # when correct, one should be true and the other false + raise ValueError( "When looking at a wedge cell for reordering, we need to construct the two triangle faces" + + " that represent the basis. With respect to the centroid of the wedge, the faces are both" + + f" oriented in the same direction with points0 '{initial_triangle0_pts}' and with points1" + + f" '{initial_triangle1_pts}'. When respecting VTK convention, they should be oriented in" + + " opposite direction. This create a degenerated wedge that cannot be reordered." ) + # We check that the 3 quad faces are not concave to validate our hypothesis for triangle basis definition + try: + quad0 = ( points[ 0 ], points[ 3 ], points[ 4 ], points[ 1 ] ) + quad1 = ( points[ 1 ], points[ 4 ], points[ 5 ], points[ 2 ] ) + quad2 = ( points[ 2 ], points[ 5 ], points[ 3 ], points[ 0 ] ) + counter = 0 + for quad in ( quad0, quad1, quad2 ): + if not is_polygon_counterclockwise( quad, centroid ): + counter += 1 + if counter == 3: + return points # the wedge is already ordered correctly + except ValueError: # quad is concave + raise ValueError( "When looking at a wedge cell for reordering, we need to construct the two triangle faces" + + " that represent the basis. When checking its geometry, the first 3 points" + + f"'{initial_triangle0_pts}' and/or last 3 points '{initial_triangle1_pts}' cannot" + + " represent the wedge basis because they created quad faces that are concave." + + " This create a degenerated wedge that cannot be reordered." ) + # We can now reorder one triangle face base to be counterclockwise. + # Once we find the right ordering for the first one, we can deduce it for the second + # We first need to find just one correct ordering of the first triangle face of the wedge + if is_triangle0_ccw: + perms = list( permutations( points[ :3 ] ) )[ 1: ] + for perm in perms: + is_triangle_ccw = is_polygon_counterclockwise( perm, centroid ) + if is_triangle_ccw: + continue + else: + correct_triangle0_pts = perm + break + correct_triangle1_pts = list() + for pt in correct_triangle0_pts: + correct_index = initial_triangle0_pts.index( pt ) + correct_triangle1_pts.append( initial_triangle1_pts[ correct_index ] ) + correct_triangle1_pts = tuple( correct_triangle1_pts ) + # we verify that the 2nd base has the correct ordering + if is_polygon_counterclockwise( correct_triangle1_pts, centroid ): + # you just need to add the two triangle points together to form the wedge + correct_ordering = correct_triangle0_pts + correct_triangle1_pts + return correct_ordering + else: + raise ValueError( "Error when reordering the wedge." ) + + +# BIG HYPOTHESIS: The first 4 points define a quad and the next 4 define another quad. +@check_points_to_reorder( VTK_HEXAHEDRON ) +def reorder_hexahedron( points: Points3D ) -> Points3D: + points_array: np.array = np.array( points ) + centroid: np.array = np.mean( points_array, axis=0 ) + initial_quad0_pts: Points3D = ( points[ 0 ], points[ 3 ], points[ 2 ], points[ 1 ] ) # face 4 in convention + initial_quad1_pts: Points3D = ( points[ 4 ], points[ 5 ], points[ 6 ], points[ 7 ] ) # face 5 in convention + try: + is_quad0_ccw: bool = is_polygon_counterclockwise( initial_quad0_pts, centroid ) + is_quad1_ccw: bool = is_polygon_counterclockwise( initial_quad1_pts, centroid ) + except ValueError: # quad is concave + raise ValueError( "When looking at a hexahedron cell for reordering, we need to construct two quad faces" + + " that represent two faces that do not have a point common. When checking its geometry," + + f" the first 4 points '{initial_quad0_pts}' and/or last 4 points '{initial_quad1_pts}'" + + " cannot represent two hexahedron quad faces because they are concave." + + " This create a degenerated hexahedron that cannot be reordered." ) + if not is_quad0_ccw and not is_quad1_ccw: + return points # the hexahedron is already correctly ordered + if is_quad0_ccw != is_quad1_ccw: # when correct, both should be false or true + raise ValueError( "When looking at a hexahedron cell for reordering, we need to construct two quad faces" + + " that represent two faces that do not have a point common. With respect to the centroid" + + " of the hexahedron, the faces are not both oriented in the same direction with points0" + + f" '{initial_quad0_pts}' and with points1 '{initial_quad1_pts}'. When respecting VTK" + + " convention, they both should be oriented in the same direction. This create a degenerated" + + " hexahedron that cannot be reordered." ) + # We can now reorder the first quad face base to be counterclockwise. + # Once we find the right ordering for the first one, we can deduce it for the second + # We first need to find just one correct ordering of the first quad face of the hexahedron + if is_quad0_ccw: + perms = list( permutations( ( points[ 0 ], points[ 3 ], points[ 2 ], points[ 1 ] ) ) )[ 1: ] + for perm in perms: + try: + is_quad0_ccw = is_polygon_counterclockwise( ( perm[ 0 ], perm[ 3 ], perm[ 2 ], perm[ 1 ] ), centroid ) + except ValueError: + continue + if is_quad0_ccw: + continue + else: + correct_quad0_pts = ( perm[ 3 ], perm[ 0 ], perm[ 1 ], perm[ 2 ] ) + break + correct_quad1_pts = list() + correspondance_table = { 0: 0, 3: 1, 2: 2, 1: 3 } + for pt in correct_quad0_pts: + correct_index = initial_quad0_pts.index( pt ) + correct_quad1_pts.append( initial_quad1_pts[ correspondance_table[ correct_index ] ] ) + correct_quad1_pts = tuple( correct_quad1_pts ) + # we verify that the 2nd quad has the correct ordering + if not is_polygon_counterclockwise( correct_quad1_pts, centroid ): + # you just need to add the two quad points together to form the hexahedron + correct_ordering = correct_quad0_pts + correct_quad1_pts + return correct_ordering + else: + raise ValueError( "Error when reordering the hexahedron." ) + + +# BIG HYPOTHESIS: The first 5 points define a pentagon face and the next 5 define another pentagon face. +@check_points_to_reorder( VTK_PENTAGONAL_PRISM ) +def reorder_pentagonal_prism( points: Points3D ) -> Points3D: + points_array: np.array = np.array( points ) + centroid: np.array = np.mean( points_array, axis=0 ) + initial_penta0_pts: Points3D = ( points[ 0 ], points[ 4 ], points[ 3 ], points[ 2 ], points[ 1 ] ) # face 0 + initial_penta1_pts: Points3D = ( points[ 5 ], points[ 6 ], points[ 7 ], points[ 8 ], points[ 9 ] ) # face 1 + try: + is_penta0_ccw: bool = is_polygon_counterclockwise( initial_penta0_pts, centroid ) + is_penta1_ccw: bool = is_polygon_counterclockwise( initial_penta1_pts, centroid ) + except ValueError: + raise ValueError( "When looking at a pentagonal prism cell for reordering, we need to construct the two" + + " pentagon faces that represent the basis. When checking its geometry, the first 5 points" + + f"'{initial_penta0_pts}' and/or last 5 points '{initial_penta1_pts}' cannot" + + " represent the pentagonal prism basis because they created pentagon faces that are" + " concave. This create a degenerated pentagonal prism that cannot be reordered." ) + if not is_penta0_ccw and not is_penta1_ccw: + return points # the pentagonal prism is already correctly ordered + if is_penta0_ccw != is_penta1_ccw: # when correct, both should be false or true + raise ValueError( "When looking at a pentagonal prism cell for reordering, we need to construct the two" + + " pentagon faces that represent the basis. With respect to the centroid of the wedge, the" + + f" faces are oriented in opposite direction with points0 '{initial_penta0_pts}' and" + + f" with points1 '{initial_penta1_pts}'. When respecting VTK convention, they should be" + + " oriented in the same direction. This create a degenerated pentagonal prism that cannot be" + + " reordered." ) + # We can now reorder the first penta face base to be counterclockwise. + # Once we find the right ordering for the first one, we can deduce it for the second + # We first need to find just one correct ordering of the first penta face of the pentagonal prism + possible_penta0_polygons: list[ Points3D ] = list() + if is_penta0_ccw: + perms = list( permutations( ( points[ 0 ], points[ 4 ], points[ 3 ], points[ 2 ], points[ 1 ] ) ) )[ 1: ] + for perm in perms: + try: + is_penta0_ccw = is_polygon_counterclockwise( ( perm[ 0 ], perm[ 4 ], perm[ 3 ], perm[ 2 ], perm[ 1 ] ), + centroid ) + except ValueError: + continue + if is_penta0_ccw: + continue + else: + possible_penta0_polygons.append( ( perm[ 0 ], perm[ 3 ], perm[ 1 ], perm[ 4 ], perm[ 2 ] ) ) + correct_penta0_pts: Points3D = choose_correct_polygon( possible_penta0_polygons ) + correct_penta1_pts = list() + correspondance_table = { 0: 0, 4: 1, 3: 2, 2: 3, 1: 4 } + for pt in correct_penta0_pts: + correct_index = initial_penta0_pts.index( pt ) + correct_penta1_pts.append( initial_penta1_pts[ correspondance_table[ correct_index ] ] ) + correct_penta1_pts = tuple( correct_penta1_pts ) + # we verify that the 2nd penta has the correct ordering + if not is_polygon_counterclockwise( correct_penta1_pts, centroid ): + # you just need to add the two penta points together to form the pentagonal prism + correct_ordering = correct_penta0_pts + correct_penta1_pts + return correct_ordering + else: + raise ValueError( "Error when reordering the pentagonal prism." ) + + +# BIG HYPOTHESIS: The first 6 points define a hexagon face and the next 6 define another hexagon face. +@check_points_to_reorder( VTK_HEXAGONAL_PRISM ) +def reorder_hexagonal_prism( points: Points3D ) -> Points3D: + points_array: np.array = np.array( points ) + centroid: np.array = np.mean( points_array, axis=0 ) + initial_hexa0_pts: Points3D = ( points[ 0 ], points[ 5 ], points[ 4 ], points[ 3 ], points[ 2 ], points[ 1 ] ) + initial_hexa1_pts: Points3D = ( points[ 6 ], points[ 7 ], points[ 8 ], points[ 9 ], points[ 10 ], points[ 11 ] ) + try: + is_hexa0_ccw: bool = is_polygon_counterclockwise( initial_hexa0_pts, centroid ) + is_hexa1_ccw: bool = is_polygon_counterclockwise( initial_hexa1_pts, centroid ) + except ValueError: + raise ValueError( "When looking at a hexagonal prism cell for reordering, we need to construct the two" + + " hexagon faces that represent the basis. When checking its geometry, the first 6 points" + + f"'{initial_hexa0_pts}' and/or last 6 points '{initial_hexa1_pts}' cannot" + + " represent the hexagonal prism basis because they created hexagon faces that are" + + " concave. This create a degenerated hexagonal prism that cannot be reordered." ) + if not is_hexa0_ccw and not is_hexa1_ccw: + return points # the hexagonal prism is already correctly ordered + if is_hexa0_ccw != is_hexa1_ccw: # when correct, both should be false or true + raise ValueError( "When looking at a hexagonal prism cell for reordering, we need to construct the two" + + " hexagon faces that represent the basis. With respect to the centroid of the wedge, the" + + f" faces are oriented in opposite direction with points0 '{initial_hexa0_pts}' and" + + f" with points1 '{initial_hexa1_pts}'. When respecting VTK convention, they should be" + + " oriented in the same direction. This create a degenerated hexagonal prism that cannot be" + + " reordered." ) + # We can now reorder the first hexagon face base to be counterclockwise. + # Once we find the right ordering for the first one, we can deduce it for the second + # We first need to find just one correct ordering of the first hexagon face of the hexagonal prism + possible_hexa0_polygons: list[ Points3D ] = list() + if is_hexa0_ccw: + perms = list( permutations( + ( points[ 0 ], points[ 5 ], points[ 4 ], points[ 3 ], points[ 2 ], points[ 1 ] ) ) )[ 1: ] + for perm in perms: + try: + is_hexa0_ccw = is_polygon_counterclockwise( + ( perm[ 0 ], perm[ 5 ], perm[ 4 ], perm[ 3 ], perm[ 2 ], perm[ 1 ] ), centroid ) + except ValueError: + continue + if is_hexa0_ccw: + continue + else: + possible_hexa0_polygons.append( ( perm[ 0 ], perm[ 3 ], perm[ 1 ], perm[ 4 ], perm[ 2 ], perm[ 5 ] ) ) + correct_hexa0_pts: Points3D = choose_correct_polygon( possible_hexa0_polygons ) + correct_hexa1_pts = list() + correspondance_table = { 0: 0, 5: 1, 4: 2, 3: 3, 2: 4, 1: 5 } + for pt in correct_hexa0_pts: + correct_index = initial_hexa0_pts.index( pt ) + correct_hexa1_pts.append( initial_hexa1_pts[ correspondance_table[ correct_index ] ] ) + correct_hexa1_pts = tuple( correct_hexa1_pts ) + # we verify that the 2nd hexagon has the correct ordering + if not is_polygon_counterclockwise( correct_hexa1_pts, centroid ): + # you just need to add the two hexagon points together to form the hexagonal prism + correct_ordering = correct_hexa0_pts + correct_hexa1_pts + return correct_ordering + else: + raise ValueError( "Error when reordering the hexagonal prism." ) + + +REORDERING_FUNCTION: dict[ int, any ] = { + VTK_TETRA: reorder_tetrahedron, + VTK_PYRAMID: reorder_pyramid, + VTK_WEDGE: reorder_wedge, + VTK_HEXAHEDRON: reorder_hexahedron, + VTK_PENTAGONAL_PRISM: reorder_pentagonal_prism, + VTK_HEXAGONAL_PRISM: reorder_hexagonal_prism +} + + +def cell_point_ids_ordering_method( cell: vtkCell, vtk_type: int ) -> tuple[ tuple[ int ], bool ]: + """For a vtkCell, gives back the correct ordering of point ids respecting the VTK convention. + If the ordering is the same as the one given as input, specifies it by returning False as a second value ; True + if the ordering is different and therefore needs to be applied later on in other algorithms. + + Args: + cell (vtkCell): A vtk cell. + vtk_type (int): Int value describing the type of vtk cell to reorder. Returns: - tuple: ( vtkDataSet, reordering_stats ) + tuple[ tuple[ int ], bool ]: ( The ordering method, if_ordering_is_different ) """ - unique_cell_types, total_per_cell_types = get_cell_types_and_number( old_mesh ) - unique_cell_names: list[ str ] = [ VTK_TYPE_TO_NAME[ vtk_type ] for vtk_type in unique_cell_types ] - names_with_totals: dict[ str, int ] = { n: v for n, v in zip( unique_cell_names, total_per_cell_types ) } - # sorted dict allow for sorted output of reordering_stats - names_with_totals_sorted: dict[ str, int ] = dict( sorted( names_with_totals.items() ) ) - useful_VTK_TYPEs: list[ int ] = [ NAME_TO_VTK_TYPE[ name ] for name in options.cell_names_to_reorder ] - all_cells_volume: np.array = compute_mesh_cells_volume( old_mesh ) - # a new mesh with the same data is created from the old mesh - new_mesh: vtkDataSet = old_mesh.NewInstance() - new_mesh.CopyStructure( old_mesh ) - new_mesh.CopyAttributes( old_mesh ) - cells = new_mesh.GetCells() - # Statistics on how many cells have or have not been reordered + points: list[ Coordinates3D ] = list() + for v in range( cell.GetNumberOfPoints() ): + points.append( cell.GetPoints().GetPoint( v ) ) + initial_cell_points: Points3D = tuple( points ) + reordered_points: Points3D = REORDERING_FUNCTION[ vtk_type ]( initial_cell_points ) + reordering_method: list[ int ] = list() + for reorder_point in reordered_points: + matching_id: int = initial_cell_points.index( reorder_point ) + reordering_method.append( matching_id ) + no_reordering = [ i for i in range( len( reordering_method ) ) ] + change_order = not no_reordering == reordering_method + return ( tuple( reordering_method ), change_order ) + + +def reorder_points_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple: reordering_stats: dict[ str, list[ any ] ] = { "Types reordered": list(), "Number of cells reordered": list(), - "Types non reordered": list( names_with_totals_sorted.keys() ), - "Number of cells non reordered": list( names_with_totals_sorted.values() ) + "Types non reordered because ordering is already correct": list(), + "Number of cells non reordered because ordering is already correct": list(), + "Types non reordered because of errors": list(), + "Number of cells non reordered because of errors": list(), + "Error message given": list() } - counter_cells_reordered: dict[ str, int ] = { name: 0 for name in options.cell_names_to_reorder } - # sorted dict allow for sorted output of reordering_stats - ounter_cells_reordered_sorted: dict[ str, int ] = dict( sorted( counter_cells_reordered.items() ) ) - # Reordering operations - for cell_id in range( new_mesh.GetNumberOfCells() ): - vtk_type: int = new_mesh.GetCellType( cell_id ) - if vtk_type in useful_VTK_TYPEs: # only cell types specified are checked - if is_cell_to_reorder( float( all_cells_volume[ cell_id ] ), options ): - cell_name: str = VTK_TYPE_TO_NAME[ vtk_type ] - point_ids_ordering: list[ int ] = valid_cell_point_ids_ordering( new_mesh.GetCell( cell_id ) ) - cells.ReplaceCellAtId( cell_id, to_vtk_id_list( point_ids_ordering ) ) - ounter_cells_reordered_sorted[ cell_name ] += 1 - # Calculation of stats - for cell_name_reordered, amount in ounter_cells_reordered_sorted.items(): - if amount > 0: - reordering_stats[ "Types reordered" ].append( cell_name_reordered ) - reordering_stats[ "Number of cells reordered" ].append( amount ) - index_non_reordered: int = reordering_stats[ "Types non reordered" ].index( cell_name_reordered ) - reordering_stats[ "Number of cells non reordered" ][ index_non_reordered ] -= amount - if reordering_stats[ "Number of cells non reordered" ][ index_non_reordered ] == 0: - reordering_stats[ "Types non reordered" ].pop( index_non_reordered ) - reordering_stats[ "Number of cells non reordered" ].pop( index_non_reordered ) + # 1st step: find the correct ordering method for each vtk type. This type should be unique because a mesh cannot + # be built with cells that use different points ordering, only one per cell type + cell_ids_to_reorder: dict[ int, np.array ] = cell_ids_to_reorder_by_cell_type( old_mesh, options ) + ordering_method_per_vtk_type: dict[ int, tuple[ int ] ] = dict() + for vtk_type, cell_ids in cell_ids_to_reorder.items(): + cell_to_check: vtkCell = old_mesh.GetCell( cell_ids[ 0 ] ) + cell_counter = 1 + try: + ordering_method, change_order = cell_point_ids_ordering_method( cell_to_check, vtk_type ) + while change_order == False and cell_counter < cell_ids.size: + cell_to_check = old_mesh.GetCell( cell_ids[ cell_counter ] ) + ordering_method, change_order = cell_point_ids_ordering_method( cell_to_check, vtk_type ) + cell_counter += 1 + if change_order: + ordering_method_per_vtk_type[ vtk_type ] = ordering_method + reordering_stats[ "Types reordered" ].append( VTK_TYPE_TO_NAME[ vtk_type ] ) + reordering_stats[ "Number of cells reordered" ].append( cell_ids.size ) + else: + ordering_method_per_vtk_type[ vtk_type ] = ordering_method + # yapf: disable + reordering_stats[ "Types non reordered because ordering is already correct" ].append( VTK_TYPE_TO_NAME[ vtk_type ] ) + reordering_stats[ "Number of cells non reordered because ordering is already correct" ].append( cell_ids.size ) + # yapf: enable + except ValueError as err_msg: + reordering_stats[ "Types non reordered because of errors" ].append( VTK_TYPE_TO_NAME[ vtk_type ] ) + reordering_stats[ "Number of cells non reordered because of errors" ].append( cell_ids.size ) + reordering_stats[ "Error message given" ].append( err_msg ) + + # 2nd step: once the ordering has been found for each vtk type, we can modify the ordering of the cells + # of this type if they have to be reordered + new_mesh = old_mesh.NewInstance() + new_mesh.CopyStructure( old_mesh ) + new_mesh.CopyAttributes( old_mesh ) + cells = new_mesh.GetCells() + for vtk_type, new_ordering in ordering_method_per_vtk_type.items(): + cell_ids: np.array = cell_ids_to_reorder[ vtk_type ] + for cell_id in cell_ids: + support_point_ids = vtkIdList() + cells.GetCellAtId( cell_id, support_point_ids ) + new_support_point_ids: list[ int ] = list() + for matching_id in new_ordering: + new_support_point_ids.append( support_point_ids.GetId( matching_id ) ) + cells.ReplaceCellAtId( cell_id, to_vtk_id_list( new_support_point_ids ) ) + return ( new_mesh, reordering_stats ) def __check( mesh, options: Options ) -> Result: - output_mesh, reordering_stats = reorder_nodes_to_new_mesh( mesh, options ) + output_mesh, reordering_stats = reorder_points_to_new_mesh( mesh, options ) write_mesh( output_mesh, options.vtk_output ) return Result( output=options.vtk_output.output, reordering_stats=reordering_stats ) diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py index 8fe86906..e5d5a628 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py +++ b/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py @@ -1,5 +1,4 @@ import logging -import random from geos.mesh.doctor.checks.fix_elements_orderings import Options, Result, NAME_TO_VTK_TYPE from . import vtk_output_parsing, FIX_ELEMENTS_ORDERINGS @@ -17,16 +16,16 @@ def fill_subparser( subparsers ) -> None: type=str, metavar=", ".join( map( str, __CELL_NAMES_CHOICES ) ), default=", ".join( map( str, __CELL_NAMES_CHOICES ) ), - help=f"[list of str]: Cell names that can be reordered in your grid. You can use multiple names." + - "Defaults to all cell names being used." ) + help=( "[list of str]: Cell names that can be reordered in your grid. You can use multiple names." + + " Defaults to all cell names being used." ) ) p.add_argument( '--' + __VOLUME_TO_REORDER, type=str, default=__VOLUME_TO_REORDER_DEFAULT, metavar=", ".join( __VOLUME_TO_REORDER_CHOICES ), - help="[str]: Select which element volume is invalid and needs reordering." + - " 'all' will allow reordering of nodes for every element, regarding of their volume." + - " 'positive' or 'negative' will only reorder the element with the corresponding volume." + - " Defaults to 'negative'." ) + help=( "[str]: Select which element volume is invalid and needs reordering." + + " 'all' will allow reordering of nodes for every element, regarding of their volume." + + " 'positive' or 'negative' will only reorder the element with the corresponding volume." + + " Defaults to 'negative'." ) ) vtk_output_parsing.fill_vtk_output_subparser( p ) @@ -40,7 +39,7 @@ def convert( parsed_options ) -> Options: cell_names_to_reorder = tuple( raw_mapping.split( "," ) ) for cell_name in cell_names_to_reorder: if cell_name not in __CELL_NAMES_CHOICES: - raise ValueError( f"Please choose names between these options for --{__CELL_NAMES_CHOICES}:" + + raise ValueError( f"Please choose names between these options for --{__CELL_NAMES}:" + f" {__CELL_NAMES_CHOICES}." ) vtk_output = vtk_output_parsing.convert( parsed_options ) volume_to_reorder: str = parsed_options[ __VOLUME_TO_REORDER ] @@ -57,13 +56,26 @@ def display_results( options: Options, result: Result ): logging.info( f"New mesh was written to file '{result.output}'" ) else: logging.info( "No output file was written." ) - logging.info( f"Number of cells reordered:" ) - logging.info( f"\tCellType\tNumber" ) - for i in range( len( result.reordering_stats[ "Types reordered" ] ) ): - logging.info( f"\t{result.reordering_stats[ 'Types reordered' ][ i ]}" + - f"\t\t{result.reordering_stats[ 'Number of cells reordered' ][ i ]}" ) - logging.info( f"Number of cells non reordered:" ) - logging.info( f"\tCellType\tNumber" ) - for i in range( len( result.reordering_stats[ "Types non reordered" ] ) ): - logging.info( f"\t{result.reordering_stats[ 'Types non reordered' ][ i ]}" + - f"\t\t{result.reordering_stats[ 'Number of cells non reordered' ][ i ]}" ) + if len( result.reordering_stats[ "Types reordered" ] ) > 0: + logging.info( "Number of cells reordered:" ) + logging.info( "\tCellType\tNumber" ) + for i in range( len( result.reordering_stats[ "Types reordered" ] ) ): + type_r = result.reordering_stats[ "Types reordered" ][ i ] + number = result.reordering_stats[ "Number of cells reordered" ][ i ] + logging.info( f"\t{type_r}\t\t{number}" ) + if len( result.reordering_stats[ "Types non reordered because ordering is already correct" ] ) > 0: + logging.info( "Number of cells non reordered because ordering is already correct:" ) + logging.info( "\tCellType\tNumber" ) + for i in range( len( result.reordering_stats[ "Types non reordered because ordering is already correct" ] ) ): + type_nr = result.reordering_stats[ "Types non reordered because ordering is already correct" ][ i ] + number = result.reordering_stats[ "Number of cells non reordered because ordering is already correct" ][ i ] + logging.info( f"\t{type_nr}\t\t{number}" ) + if len( result.reordering_stats[ "Types non reordered because of errors" ] ) > 0: + logging.info( "Number of cells non reordered because of errors:" ) + logging.info( "\tCellType\tNumber" ) + for i in range( len( result.reordering_stats[ "Types non reordered because of errors" ] ) ): + type_nr = result.reordering_stats[ "Types non reordered because of errors" ][ i ] + number = result.reordering_stats[ "Number of cells non reordered because of errors" ][ i ] + err_msg = result.reordering_stats[ "Error message given" ][ i ] + logging.info( f"\t{type_nr}\t\t{number}" ) + logging.info( f"\tError message: {err_msg}" ) \ No newline at end of file diff --git a/geos-mesh/tests/test_fix_elements_orderings.py b/geos-mesh/tests/test_fix_elements_orderings.py index 060c06f1..c53c4e8f 100644 --- a/geos-mesh/tests/test_fix_elements_orderings.py +++ b/geos-mesh/tests/test_fix_elements_orderings.py @@ -4,6 +4,7 @@ import logging import subprocess import numpy as np +from math import pi from geos.mesh.doctor.mesh_doctor import MESH_DOCTOR_FILEPATH from geos.mesh.doctor.checks import fix_elements_orderings as feo from geos.mesh.doctor.checks.generate_cube import Options, __build @@ -16,7 +17,27 @@ VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM ) -def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int ] ): +def degrees_to_radians( degrees ): + return degrees * ( pi / 180 ) + + +def radians_to_degrees( radians ): + return radians * ( 180 / pi ) + + +# yapf: disable +def rotate_polygon_around_point( polygon: np.array, rotation_point: np.array, a: float, b: float, g: float ) -> np.array: + rotation_matrix = np.array([ [ np.cos(b)*np.cos(g), np.sin(a)*np.sin(b)*np.cos(g) - np.cos(a)*np.sin(g), np.cos(a)*np.sin(b)*np.cos(g) + np.sin(a)*np.sin(g) ], + [ np.cos(b)*np.sin(g), np.sin(a)*np.sin(b)*np.sin(g) + np.cos(a)*np.cos(g), np.cos(a)*np.sin(b)*np.sin(g) - np.sin(a)*np.cos(g) ], + [ -np.sin(b), np.sin(a)*np.cos(b), np.cos(a)*np.cos(b) ] ]) + translated_points = polygon - rotation_point + rotated_points = np.dot( translated_points, rotation_matrix.T ) + rotated_polygon = rotated_points + rotation_point + return rotated_polygon +# yapf: enable + + +def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: tuple[ int ] ): """Utility function to reorder the nodes of one cell for test purposes. Args: @@ -26,7 +47,7 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int """ if mesh.GetCell( cell_id ).GetNumberOfPoints() != len( node_ordering ): raise ValueError( f"The cell to reorder needs to have '{mesh.GetCell( cell_id ).GetNumberOfPoints()}'" + - "nodes in reordering." ) + " nodes in reordering." ) cells = mesh.GetCells() support_point_ids = vtkIdList() cells.GetCellAtId( cell_id, support_point_ids ) @@ -45,16 +66,26 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int filepath_non_ordered_mesh: str = os.path.join( dir_name, "to_reorder_mesh.vtu" ) filepath_reordered_mesh: str = os.path.join( dir_name, "reordered_mesh.vtu" ) test_file: VtkOutput = VtkOutput( filepath_non_ordered_mesh, True ) +filepath_non_ordered_mesh2: str = os.path.join( dir_name, "to_reorder_mesh2.vtu" ) +filepath_reordered_mesh2: str = os.path.join( dir_name, "reordered_mesh2.vtu" ) +test_file2: VtkOutput = VtkOutput( filepath_non_ordered_mesh2, True ) """ Dict used to apply false nodes orderings for test purposes """ -to_change_order: dict[ int, list[ int ] ] = { - "Hexahedron": ( 0, 3, 2, 1, 4, 5, 6, 7 ), +to_change_order: dict[ int, tuple[ int ] ] = { + "Hexahedron": ( 0, 3, 2, 1, 4, 7, 6, 5 ), "Tetrahedron": ( 0, 2, 1, 3 ), "Pyramid": ( 0, 3, 2, 1, 4 ), - "Wedge": ( 0, 2, 1, 3, 4, 5 ), - "Prism5": ( 0, 4, 3, 2, 1, 5, 6, 7, 8, 9 ), - "Prism6": ( 0, 1, 4, 2, 3, 5, 11, 10, 9, 8, 7, 6 ) + "Wedge": ( 0, 2, 1, 3, 5, 4 ), + "Prism5": ( 0, 4, 3, 2, 1, 5, 9, 8, 7, 6 ), + "Prism6": ( 0, 5, 4, 3, 2, 1, 6, 11, 10, 9, 8, 7 ) +} +to_change_order = dict( sorted( to_change_order.items() ) ) +to_change_order_for_degenerated_cells: dict[ int, tuple[ int ] ] = { + "Hexahedron": ( 0, 2, 1, 3, 4, 6, 5, 7 ), + "Wedge": ( 0, 1, 2, 5, 4, 3 ), + "Prism5": ( 0, 1, 2, 3, 4, 7, 6, 5, 9, 8 ), + "Prism6": ( 0, 1, 2, 3, 4, 5, 11, 10, 9, 8, 7, 6 ) } to_change_order = dict( sorted( to_change_order.items() ) ) cell_names = list( VTK_TYPE_TO_NAME.values() ) @@ -94,13 +125,23 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int hexahedrons_grid_invalid: vtkDataSet = __build( options_hexahedrons_grid ) for i in range( 2 ): reorder_cell_nodes( hexahedrons_grid_invalid, i * 2 + 1, to_change_order[ "Hexahedron" ] ) + +opt_hexahedrons_grid = opt( test_file, [ "Hexahedron" ], "negative" ) +opt_hexahedrons_grid_invalid = opt( test_file, [ "Hexahedron" ], "negative" ) """ 4 tetrahedrons """ points_tetras: vtkPoints = vtkPoints() -points_tetras_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ), - ( 0.0, 1.0, 0.0 ), ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), - ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ) ] +# yapf: disable +points_tetras_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), # point0 + ( 1.0, 0.0, 0.0 ), + ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), + ( 0.0, 0.0, 1.0 ), + ( 1.0, 0.0, 1.0 ), # point5 + ( 1.0, 1.0, 1.0 ), + ( 0.0, 1.0, 1.0 ) ] +# yapf: enable for point_tetra in points_tetras_coords: points_tetras.InsertNextPoint( point_tetra ) @@ -143,17 +184,35 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int tetras_grid_invalid.DeepCopy( tetras_grid ) for i in range( 2 ): reorder_cell_nodes( tetras_grid_invalid, i * 2 + 1, to_change_order[ "Tetrahedron" ] ) + +opt_tetras_grid = opt( test_file, [ "Tetrahedron" ], "negative" ) +opt_tetras_grid_invalid = opt( test_file, [ "Tetrahedron" ], "negative" ) """ 4 pyramids """ points_pyramids: vtkPoints = vtkPoints() -points_pyramids_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ), - ( 0.0, 1.0, 0.0 ), ( 0.5, 0.5, 1.0 ), ( 2.0, 0.0, 0.0 ), - ( 3.0, 0.0, 0.0 ), ( 3.0, 1.0, 0.0 ), ( 2.0, 1.0, 0.0 ), - ( 2.5, 0.5, 1.0 ), ( 0.0, 2.0, 0.0 ), ( 1.0, 2.0, 0.0 ), - ( 1.0, 3.0, 0.0 ), ( 0.0, 3.0, 0.0 ), ( 0.5, 2.5, 1.0 ), - ( 2.0, 2.0, 0.0 ), ( 3.0, 2.0, 0.0 ), ( 3.0, 3.0, 0.0 ), - ( 2.0, 3.0, 0.0 ), ( 2.5, 2.5, 1.0 ) ] +# yapf: disable +points_pyramids_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), # point0 + ( 1.0, 0.0, 0.0 ), + ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), + ( 0.5, 0.5, 1.0 ), + ( 2.0, 0.0, 0.0 ), # point5 + ( 3.0, 0.0, 0.0 ), + ( 3.0, 1.0, 0.0 ), + ( 2.0, 1.0, 0.0 ), + ( 2.5, 0.5, 1.0 ), + ( 0.0, 2.0, 0.0 ), # point10 + ( 1.0, 2.0, 0.0 ), + ( 1.0, 3.0, 0.0 ), + ( 0.0, 3.0, 0.0 ), + ( 0.5, 2.5, 1.0 ), + ( 2.0, 2.0, 0.0 ), # point15 + ( 3.0, 2.0, 0.0 ), + ( 3.0, 3.0, 0.0 ), + ( 2.0, 3.0, 0.0 ), + ( 2.5, 2.5, 1.0 ) ] +# yapf: enable for point_pyramid in points_pyramids_coords: points_pyramids.InsertNextPoint( point_pyramid ) @@ -200,21 +259,47 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int pyramids_grid_invalid.DeepCopy( pyramids_grid ) for i in range( 2 ): reorder_cell_nodes( pyramids_grid_invalid, i * 2 + 1, to_change_order[ "Pyramid" ] ) + +opt_pyramids_grid = opt( test_file, [ "Pyramid" ], "negative" ) +opt_pyramids_grid_invalid = opt( test_file, [ "Pyramid" ], "negative" ) """ 4 voxels: this type of element cannot be used in GEOS, we just test that the feature rejects them """ points_voxels: vtkPoints = vtkPoints() -points_voxels_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ), - ( 0.0, 1.0, 0.0 ), ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), - ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ), ( 2.0, 0.0, 0.0 ), - ( 3.0, 0.0, 0.0 ), ( 3.0, 1.0, 0.0 ), ( 2.0, 1.0, 0.0 ), - ( 2.0, 0.0, 1.0 ), ( 3.0, 0.0, 1.0 ), ( 3.0, 1.0, 1.0 ), - ( 2.0, 1.0, 1.0 ), ( 0.0, 2.0, 0.0 ), ( 1.0, 2.0, 0.0 ), - ( 1.0, 3.0, 0.0 ), ( 0.0, 3.0, 0.0 ), ( 0.0, 2.0, 1.0 ), - ( 1.0, 2.0, 1.0 ), ( 1.0, 3.0, 1.0 ), ( 0.0, 3.0, 1.0 ), - ( 2.0, 2.0, 0.0 ), ( 3.0, 2.0, 0.0 ), ( 3.0, 3.0, 0.0 ), - ( 2.0, 3.0, 0.0 ), ( 2.0, 2.0, 1.0 ), ( 3.0, 2.0, 1.0 ), - ( 3.0, 3.0, 1.0 ), ( 2.0, 3.0, 1.0 ) ] +# yapf: disable +points_voxels_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), # point0 + ( 1.0, 0.0, 0.0 ), + ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), + ( 0.0, 0.0, 1.0 ), + ( 1.0, 0.0, 1.0 ), # point5 + ( 1.0, 1.0, 1.0 ), + ( 0.0, 1.0, 1.0 ), + ( 2.0, 0.0, 0.0 ), + ( 3.0, 0.0, 0.0 ), + ( 3.0, 1.0, 0.0 ), # point10 + ( 2.0, 1.0, 0.0 ), + ( 2.0, 0.0, 1.0 ), + ( 3.0, 0.0, 1.0 ), + ( 3.0, 1.0, 1.0 ), + ( 2.0, 1.0, 1.0 ), # point15 + ( 0.0, 2.0, 0.0 ), + ( 1.0, 2.0, 0.0 ), + ( 1.0, 3.0, 0.0 ), + ( 0.0, 3.0, 0.0 ), + ( 0.0, 2.0, 1.0 ), # point20 + ( 1.0, 2.0, 1.0 ), + ( 1.0, 3.0, 1.0 ), + ( 0.0, 3.0, 1.0 ), + ( 2.0, 2.0, 0.0 ), + ( 3.0, 2.0, 0.0 ), # point25 + ( 3.0, 3.0, 0.0 ), + ( 2.0, 3.0, 0.0 ), + ( 2.0, 2.0, 1.0 ), + ( 3.0, 2.0, 1.0 ), + ( 3.0, 3.0, 1.0 ), # point30 + ( 2.0, 3.0, 1.0 ) ] +# yapf: enable for point_voxel in points_voxels_coords: points_voxels.InsertNextPoint( point_voxel ) @@ -267,14 +352,27 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int voxels_grid: vtkUnstructuredGrid = vtkUnstructuredGrid() voxels_grid.SetPoints( points_voxels ) voxels_grid.SetCells( VTK_VOXEL, voxels_cells ) + +opt_voxels_grid = opt( test_file, [ "Voxel" ], "negative" ) +opt_voxels_grid_invalid = opt( test_file, [ "Voxel" ], "negative" ) """ 4 wedges """ points_wedges: vtkPoints = vtkPoints() -points_wedges_coords: list[ tuple[ float ] ] = [ ( 0.5, 0.0, 0.0 ), ( 1.5, 0.0, 0.0 ), ( 2.5, 0.0, 0.0 ), - ( 0.0, 1.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 2.0, 1.0, 0.0 ), - ( 0.5, 0.0, 1.0 ), ( 1.5, 0.0, 1.0 ), ( 2.5, 0.0, 1.0 ), - ( 0.0, 1.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( 2.0, 1.0, 1.0 ) ] +# yapf: disable +points_wedges_coords: list[ tuple[ float ] ] = [ ( 0.5, 0.0, 0.0 ), # point0 + ( 1.5, 0.0, 0.0 ), + ( 2.5, 0.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), + ( 1.0, 1.0, 0.0 ), + ( 2.0, 1.0, 0.0 ), # point5 + ( 0.5, 0.0, 1.0 ), + ( 1.5, 0.0, 1.0 ), + ( 2.5, 0.0, 1.0 ), + ( 0.0, 1.0, 1.0 ), + ( 1.0, 1.0, 1.0 ), # point10 + ( 2.0, 1.0, 1.0 ) ] +# yapf: enable for point_wedge in points_wedges_coords: points_wedges.InsertNextPoint( point_wedge ) @@ -325,24 +423,55 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int wedges_grid_invalid.DeepCopy( wedges_grid ) for i in range( 2 ): reorder_cell_nodes( wedges_grid_invalid, i * 2 + 1, to_change_order[ "Wedge" ] ) + +opt_wedges_grid = opt( test_file, [ "Wedge" ], "negative" ) +opt_wedges_grid_invalid = opt( test_file, [ "Wedge" ], "negative" ) """ 4 pentagonal prisms """ points_penta_prisms: vtkPoints = vtkPoints() -points_penta_prisms_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.5, 0.5, 0.0 ), - ( 0.5, 1.0, 0.0 ), ( -0.5, 0.5, 0.0 ), ( 0.0, 0.0, 1.0 ), - ( 1.0, 0.0, 1.0 ), ( 1.5, 0.5, 1.0 ), ( 0.5, 1.0, 1.0 ), - ( -0.5, 0.5, 1.0 ), ( 2.0, 0.0, 0.0 ), ( 3.0, 0.0, 0.0 ), - ( 3.5, 0.5, 0.0 ), ( 2.5, 1.0, 0.0 ), ( 1.5, 0.5, 0.0 ), - ( 2.0, 0.0, 1.0 ), ( 3.0, 0.0, 1.0 ), ( 3.5, 0.5, 1.0 ), - ( 2.5, 1.0, 1.0 ), ( 1.5, 0.5, 1.0 ), ( 0.0, 2.0, 0.0 ), - ( 1.0, 2.0, 0.0 ), ( 1.5, 2.5, 0.0 ), ( 0.5, 3.0, 0.0 ), - ( -0.5, 2.5, 0.0 ), ( 0.0, 2.0, 1.0 ), ( 1.0, 2.0, 1.0 ), - ( 1.5, 2.5, 1.0 ), ( 0.5, 3.0, 1.0 ), ( -0.5, 2.5, 1.0 ), - ( 2.0, 2.0, 0.0 ), ( 3.0, 2.0, 0.0 ), ( 3.5, 2.5, 0.0 ), - ( 2.5, 3.0, 0.0 ), ( 1.5, 2.5, 0.0 ), ( 2.0, 2.0, 1.0 ), - ( 3.0, 2.0, 1.0 ), ( 3.5, 2.5, 1.0 ), ( 2.5, 3.0, 1.0 ), +# yapf: disable +points_penta_prisms_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), # point0 + ( 1.0, 0.0, 0.0 ), + ( 1.5, 0.5, 0.0 ), + ( 0.5, 1.0, 0.0 ), + ( -0.5, 0.5, 0.0 ), + ( 0.0, 0.0, 1.0 ), # point5 + ( 1.0, 0.0, 1.0 ), + ( 1.5, 0.5, 1.0 ), + ( 0.5, 1.0, 1.0 ), + ( -0.5, 0.5, 1.0 ), + ( 2.0, 0.0, 0.0 ), # point10 + ( 3.0, 0.0, 0.0 ), + ( 3.5, 0.5, 0.0 ), + ( 2.5, 1.0, 0.0 ), + ( 1.5, 0.5, 0.0 ), + ( 2.0, 0.0, 1.0 ), # point15 + ( 3.0, 0.0, 1.0 ), + ( 3.5, 0.5, 1.0 ), + ( 2.5, 1.0, 1.0 ), + ( 1.5, 0.5, 1.0 ), + ( 0.0, 2.0, 0.0 ), # point20 + ( 1.0, 2.0, 0.0 ), + ( 1.5, 2.5, 0.0 ), + ( 0.5, 3.0, 0.0 ), + ( -0.5, 2.5, 0.0 ), + ( 0.0, 2.0, 1.0 ), # point25 + ( 1.0, 2.0, 1.0 ), + ( 1.5, 2.5, 1.0 ), + ( 0.5, 3.0, 1.0 ), + ( -0.5, 2.5, 1.0 ), + ( 2.0, 2.0, 0.0 ), # point30 + ( 3.0, 2.0, 0.0 ), + ( 3.5, 2.5, 0.0 ), + ( 2.5, 3.0, 0.0 ), + ( 1.5, 2.5, 0.0 ), + ( 2.0, 2.0, 1.0 ), # point35 + ( 3.0, 2.0, 1.0 ), + ( 3.5, 2.5, 1.0 ), + ( 2.5, 3.0, 1.0 ), ( 1.5, 2.5, 1.0 ) ] +# yapf: enable for point_penta_prism in points_penta_prisms_coords: points_penta_prisms.InsertNextPoint( point_penta_prism ) @@ -409,26 +538,63 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int penta_prism_grid_invalid.DeepCopy( penta_prism_grid ) for i in range( 2 ): reorder_cell_nodes( penta_prism_grid_invalid, i * 2 + 1, to_change_order[ "Prism5" ] ) + +opt_penta_prism_grid = opt( test_file, [ "Prism5" ], "negative" ) +opt_penta_prism_grid_invalid = opt( test_file, [ "Prism5" ], "negative" ) """ 4 hexagonal prisms """ points_hexa_prisms: vtkPoints = vtkPoints() -points_hexa_prisms_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.5, 0.5, 0.0 ), - ( 1.0, 1.0, 0.0 ), ( 0.0, 1.0, 0.0 ), ( -0.5, 0.5, 0.0 ), - ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), ( 1.5, 0.5, 1.0 ), - ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ), ( -0.5, 0.5, 1.0 ), - ( 2.0, 0.0, 0.0 ), ( 3.0, 0.0, 0.0 ), ( 3.5, 0.5, 0.0 ), - ( 3.0, 1.0, 0.0 ), ( 2.0, 1.0, 0.0 ), ( 1.5, 0.5, 0.0 ), - ( 2.0, 0.0, 1.0 ), ( 3.0, 0.0, 1.0 ), ( 3.5, 0.5, 1.0 ), - ( 3.0, 1.0, 1.0 ), ( 2.0, 1.0, 1.0 ), ( 1.5, 0.5, 1.0 ), - ( 0.0, 2.0, 0.0 ), ( 1.0, 2.0, 0.0 ), ( 1.5, 2.5, 0.0 ), - ( 1.0, 3.0, 0.0 ), ( 0.0, 3.0, 0.0 ), ( -0.5, 2.5, 0.0 ), - ( 0.0, 2.0, 1.0 ), ( 1.0, 2.0, 1.0 ), ( 1.5, 2.5, 1.0 ), - ( 1.0, 3.0, 1.0 ), ( 0.0, 3.0, 1.0 ), ( -0.5, 2.5, 1.0 ), - ( 2.0, 2.0, 0.0 ), ( 3.0, 2.0, 0.0 ), ( 3.5, 2.5, 0.0 ), - ( 3.0, 3.0, 0.0 ), ( 2.0, 3.0, 0.0 ), ( 1.5, 2.5, 0.0 ), - ( 2.0, 2.0, 1.0 ), ( 3.0, 2.0, 1.0 ), ( 3.5, 2.5, 1.0 ), - ( 3.0, 3.0, 1.0 ), ( 2.0, 3.0, 1.0 ), ( 1.5, 2.5, 1.0 ) ] +# yapf: disable +points_hexa_prisms_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), # point0 + ( 1.0, 0.0, 0.0 ), + ( 1.5, 0.5, 0.0 ), + ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), + ( -0.5, 0.5, 0.0 ), # point5 + ( 0.0, 0.0, 1.0 ), + ( 1.0, 0.0, 1.0 ), + ( 1.5, 0.5, 1.0 ), + ( 1.0, 1.0, 1.0 ), + ( 0.0, 1.0, 1.0 ), # point10 + ( -0.5, 0.5, 1.0 ), + ( 2.0, 0.0, 0.0 ), + ( 3.0, 0.0, 0.0 ), + ( 3.5, 0.5, 0.0 ), + ( 3.0, 1.0, 0.0 ), # point15 + ( 2.0, 1.0, 0.0 ), + ( 1.5, 0.5, 0.0 ), + ( 2.0, 0.0, 1.0 ), + ( 3.0, 0.0, 1.0 ), + ( 3.5, 0.5, 1.0 ), # point20 + ( 3.0, 1.0, 1.0 ), + ( 2.0, 1.0, 1.0 ), + ( 1.5, 0.5, 1.0 ), + ( 0.0, 2.0, 0.0 ), + ( 1.0, 2.0, 0.0 ), # point25 + ( 1.5, 2.5, 0.0 ), + ( 1.0, 3.0, 0.0 ), + ( 0.0, 3.0, 0.0 ), + ( -0.5, 2.5, 0.0 ), + ( 0.0, 2.0, 1.0 ), # point30 + ( 1.0, 2.0, 1.0 ), + ( 1.5, 2.5, 1.0 ), + ( 1.0, 3.0, 1.0 ), + ( 0.0, 3.0, 1.0 ), + ( -0.5, 2.5, 1.0 ), # point35 + ( 2.0, 2.0, 0.0 ), + ( 3.0, 2.0, 0.0 ), + ( 3.5, 2.5, 0.0 ), + ( 3.0, 3.0, 0.0 ), + ( 2.0, 3.0, 0.0 ), # point40 + ( 1.5, 2.5, 0.0 ), + ( 2.0, 2.0, 1.0 ), + ( 3.0, 2.0, 1.0 ), + ( 3.5, 2.5, 1.0 ), + ( 3.0, 3.0, 1.0 ), # point45 + ( 2.0, 3.0, 1.0 ), + ( 1.5, 2.5, 1.0 ) ] +# yapf: enable for point_hexa_prism in points_hexa_prisms_coords: points_hexa_prisms.InsertNextPoint( point_hexa_prism ) @@ -463,21 +629,68 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int hexa_prism_grid_invalid.DeepCopy( hexa_prism_grid ) for i in range( 2 ): reorder_cell_nodes( hexa_prism_grid_invalid, i * 2 + 1, to_change_order[ "Prism6" ] ) + +opt_hexa_prism_grid = opt( test_file, [ "Prism6" ], "negative" ) +opt_hexa_prism_grid_invalid = opt( test_file, [ "Prism6" ], "negative" ) """ 2 hexahedrons, 2 tetrahedrons, 2 wedges, 2 pyramids, 2 voxels, 2 pentagonal prisms and 2 hexagonal prisms """ points_mix: vtkPoints = vtkPoints() -points_mix_coords: list[ tuple[ float ] ] = [ - ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 2.0, 0.0, 0.0 ), ( 2.5, -0.5, 0.0 ), ( 3.0, 0.0, 0.0 ), ( 3.5, -0.5, 0.0 ), - ( 4.0, 0.0, 0.0 ), ( 4.5, -0.5, 0.0 ), ( 5.0, 0.0, 0.0 ), ( 5.5, -0.5, 0.0 ), ( 6.0, 0.5, 0.0 ), ( 0.0, 1.0, 0.0 ), - ( 1.0, 1.0, 0.0 ), ( 2.0, 1.0, 0.0 ), ( 2.5, 1.5, 0.0 ), ( 3.0, 1.0, 0.0 ), ( 4.0, 1.0, 0.0 ), ( 5.0, 1.0, 0.0 ), - ( 5.5, 1.5, 0.0 ), ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), ( 2.0, 0.0, 1.0 ), ( 2.5, -0.5, 1.0 ), ( 3.0, 0.0, 1.0 ), - ( 3.5, -0.5, 1.0 ), ( 4.0, 0.0, 1.0 ), ( 4.5, -0.5, 1.0 ), ( 5.0, 0.0, 1.0 ), ( 5.5, -0.5, 1.0 ), ( 6.0, 0.5, 1.0 ), - ( 0.0, 1.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( 2.0, 1.0, 1.0 ), ( 2.5, 1.5, 1.0 ), ( 3.0, 1.0, 1.0 ), ( 4.0, 1.0, 1.0 ), - ( 5.0, 1.0, 1.0 ), ( 5.5, 1.5, 1.0 ), ( 0.5, 0.5, 2.0 ), ( 0.5, 1.5, 2.0 ), ( 1.5, 0.5, 2.0 ), ( 1.5, 1.5, 2.0 ), - ( 2.0, 0.0, 2.0 ), ( 2.5, -0.5, 2.0 ), ( 3.0, 0.0, 2.0 ), ( 3.0, 1.0, 2.0 ), ( 2.5, 1.5, 2.0 ), ( 2.0, 1.0, 2.0 ), - ( 5.0, 0.0, 2.0 ), ( 5.5, -0.5, 2.0 ), ( 6.0, 0.5, 2.0 ), ( 5.5, 1.5, 2.0 ), ( 5.0, 1.0, 2.0 ) -] +# yapf: disable +points_mix_coords: list[ tuple[ float ] ] = [ ( 0.0, 0.0, 0.0 ), # point0 + ( 1.0, 0.0, 0.0 ), + ( 2.0, 0.0, 0.0 ), + ( 2.5, -0.5, 0.0 ), + ( 3.0, 0.0, 0.0 ), + ( 3.5, -0.5, 0.0 ), # point5 + ( 4.0, 0.0, 0.0 ), + ( 4.5, -0.5, 0.0 ), + ( 5.0, 0.0, 0.0 ), + ( 5.5, -0.5, 0.0 ), + ( 6.0, 0.5, 0.0 ), # point10 + ( 0.0, 1.0, 0.0 ), + ( 1.0, 1.0, 0.0 ), + ( 2.0, 1.0, 0.0 ), + ( 2.5, 1.5, 0.0 ), + ( 3.0, 1.0, 0.0 ), # point15 + ( 4.0, 1.0, 0.0 ), + ( 5.0, 1.0, 0.0 ), + ( 5.5, 1.5, 0.0 ), + ( 0.0, 0.0, 1.0 ), + ( 1.0, 0.0, 1.0 ), # point20 + ( 2.0, 0.0, 1.0 ), + ( 2.5, -0.5, 1.0 ), + ( 3.0, 0.0, 1.0 ), + ( 3.5, -0.5, 1.0 ), + ( 4.0, 0.0, 1.0 ), # point25 + ( 4.5, -0.5, 1.0 ), + ( 5.0, 0.0, 1.0 ), + ( 5.5, -0.5, 1.0 ), + ( 6.0, 0.5, 1.0 ), + ( 0.0, 1.0, 1.0 ), # point30 + ( 1.0, 1.0, 1.0 ), + ( 2.0, 1.0, 1.0 ), + ( 2.5, 1.5, 1.0 ), + ( 3.0, 1.0, 1.0 ), + ( 4.0, 1.0, 1.0 ), # point35 + ( 5.0, 1.0, 1.0 ), + ( 5.5, 1.5, 1.0 ), + ( 0.5, 0.5, 2.0 ), + ( 0.5, 1.5, 2.0 ), + ( 1.5, 0.5, 2.0 ), # point40 + ( 1.5, 1.5, 2.0 ), + ( 2.0, 0.0, 2.0 ), + ( 2.5, -0.5, 2.0 ), + ( 3.0, 0.0, 2.0 ), + ( 3.0, 1.0, 2.0 ), # point45 + ( 2.5, 1.5, 2.0 ), + ( 2.0, 1.0, 2.0 ), + ( 5.0, 0.0, 2.0 ), + ( 5.5, -0.5, 2.0 ), + ( 6.0, 0.5, 2.0 ), # point50 + ( 5.5, 1.5, 2.0 ), + ( 5.0, 1.0, 2.0 ) ] +# yapf: enable for point_mix in points_mix_coords: points_mix.InsertNextPoint( point_mix ) @@ -640,6 +853,16 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int for i in range( len( all_cell_types_mix_grid ) // 2 ): reorder_cell_nodes( mix_grid_invalid, i * 2 + 1, to_change_order[ all_cell_names_mix_grid[ i * 2 + 1 ] ] ) +opt_mix_grid = opt( test_file, tuple( ( "Hexahedron", "Tetrahedron", "Pyramid", "Wedge", "Prism5", "Prism6" ) ), + "negative" ) +opt_mix_grid_inv = opt( test_file, tuple( ( "Hexahedron", "Tetrahedron", "Pyramid", "Wedge", "Prism5", "Prism6" ) ), + "negative" ) + +mix_grid_degenerated = vtkUnstructuredGrid() +mix_grid_degenerated.DeepCopy( mix_grid ) +for i in [ 1, 7, 11, 13 ]: + reorder_cell_nodes( mix_grid_degenerated, i, to_change_order_for_degenerated_cells[ all_cell_names_mix_grid[ i ] ] ) + class TestClass: @@ -653,7 +876,7 @@ def test_reorder_cell_nodes( self ): # reorder the cell to make it invalid reorder_cell_nodes( one_hex, 0, to_change_order[ "Hexahedron" ] ) expected_nodes_coords_modified = [ ( 0.0, 0.0, 0.0 ), ( 0.0, 1.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 1.0, 0.0, 0.0 ), - ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ) ] + ( 0.0, 0.0, 1.0 ), ( 0.0, 1.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( 1.0, 0.0, 1.0 ) ] for i in range( one_hex.GetCell( 0 ).GetNumberOfPoints() ): assert one_hex.GetCell( 0 ).GetPoints().GetPoint( i ) == expected_nodes_coords_modified[ i ] @@ -665,24 +888,25 @@ def test_reorder_cell_nodes( self ): assert one_hex.GetCell( 0 ).GetPoints().GetPoint( i ) == expected_nodes_coords_modified2[ i ] def test_compute_mesh_cells_volume( self ): + # yapf: disable grid_volumes = { hexahedrons_grid: [ 1.0, 1.0, 1.0, 1.0 ], - hexahedrons_grid_invalid: [ 1.0, -0.333, 1.0, -0.333 ], + hexahedrons_grid_invalid: [ 1.0, -1.0, 1.0, -1.0 ], tetras_grid: [ 0.167, 0.167, 0.167, 0.167 ], tetras_grid_invalid: [ 0.167, -0.167, 0.167, -0.167 ], pyramids_grid: [ 0.333, 0.333, 0.333, 0.333 ], pyramids_grid_invalid: [ 0.333, -0.333, 0.333, -0.333 ], voxels_grid: [ 1.0, 1.0, 1.0, 1.0 ], wedges_grid: [ 0.5, 0.5, 0.5, 0.5 ], - wedges_grid_invalid: [ 0.5, -0.167, 0.5, -0.167 ], + wedges_grid_invalid: [ 0.5, -0.5, 0.5, -0.5 ], penta_prism_grid: [ 1.25, 1.25, 1.25, 1.25 ], - penta_prism_grid_invalid: [ 1.25, -0.083, 1.25, -0.083 ], + penta_prism_grid_invalid: [ 1.25, -1.25, 1.25, -1.25 ], hexa_prism_grid: [ 1.5, 1.5, 1.5, 1.5 ], - hexa_prism_grid_invalid: [ 1.5, -0.333, 1.5, -0.333 ], + hexa_prism_grid_invalid: [ 1.5, -1.5, 1.5, -1.5 ], mix_grid: [ 1.0, 1.0, 0.333, 0.333, 0.167, 0.167, 1.5, 1.5, 1.0, 1.0, 0.25, 0.25, 1.25, 1.25 ], - mix_grid_invalid: - [ 1.0, -0.333, 0.333, -0.333, 0.167, -0.167, 1.5, -0.333, 1.0, -0.333, 0.25, -0.083, 1.25, -0.083 ] + mix_grid_invalid: [ 1.0, -1.0, 0.333, -0.333, 0.167, -0.167, 1.5, -1.5, 1.0, -1.0, 0.25, -0.25, 1.25, -1.25 ] } + # yapf: enable for grid, volumes_expected in grid_volumes.items(): volumes_computed = feo.compute_mesh_cells_volume( grid ) for i in range( len( volumes_computed ) ): @@ -705,75 +929,291 @@ def test_is_cell_to_reorder( self ): mix_grid: [ False ] * 14, mix_grid_invalid: [ i % 2 != 0 for i in range( 14 ) ] } + grid_options = { + hexahedrons_grid: opt_hexahedrons_grid, + hexahedrons_grid_invalid: opt_hexahedrons_grid_invalid, + tetras_grid: opt_tetras_grid, + tetras_grid_invalid: opt_tetras_grid_invalid, + pyramids_grid: opt_pyramids_grid, + pyramids_grid_invalid: opt_pyramids_grid_invalid, + wedges_grid: opt_wedges_grid, + wedges_grid_invalid: opt_wedges_grid_invalid, + penta_prism_grid: opt_penta_prism_grid, + penta_prism_grid_invalid: opt_penta_prism_grid, + hexa_prism_grid: opt_hexa_prism_grid, + hexa_prism_grid_invalid: opt_hexa_prism_grid_invalid, + mix_grid: opt_mix_grid, + mix_grid_invalid: opt_mix_grid_inv + } for grid, needs_ordering in grid_needs_ordering.items(): volumes = feo.compute_mesh_cells_volume( grid ) + opt_to_use = grid_options[ grid ] for i in range( len( volumes ) ): - options = opt( out, cell_names, "negative" ) - assert feo.is_cell_to_reorder( volumes[ i ], options ) == needs_ordering[ i ] - options = opt( out, cell_names, "positive" ) - assert feo.is_cell_to_reorder( volumes[ i ], options ) != needs_ordering[ i ] - options = opt( out, cell_names, "all" ) - assert feo.is_cell_to_reorder( volumes[ i ], options ) == True - - def test_get_cell_types_and_number( self ): - assert feo.get_cell_types_and_number( hexahedrons_grid ) == ( [ VTK_HEXAHEDRON ], [ 4 ] ) - assert feo.get_cell_types_and_number( hexahedrons_grid_invalid ) == ( [ VTK_HEXAHEDRON ], [ 4 ] ) - assert feo.get_cell_types_and_number( tetras_grid ) == ( [ VTK_TETRA ], [ 4 ] ) - assert feo.get_cell_types_and_number( tetras_grid_invalid ) == ( [ VTK_TETRA ], [ 4 ] ) - assert feo.get_cell_types_and_number( pyramids_grid ) == ( [ VTK_PYRAMID ], [ 4 ] ) - assert feo.get_cell_types_and_number( pyramids_grid_invalid ) == ( [ VTK_PYRAMID ], [ 4 ] ) - assert feo.get_cell_types_and_number( wedges_grid ) == ( [ VTK_WEDGE ], [ 4 ] ) - assert feo.get_cell_types_and_number( wedges_grid_invalid ) == ( [ VTK_WEDGE ], [ 4 ] ) - assert feo.get_cell_types_and_number( penta_prism_grid ) == ( [ VTK_PENTAGONAL_PRISM ], [ 4 ] ) - assert feo.get_cell_types_and_number( penta_prism_grid_invalid ) == ( [ VTK_PENTAGONAL_PRISM ], [ 4 ] ) - assert feo.get_cell_types_and_number( hexa_prism_grid ) == ( [ VTK_HEXAGONAL_PRISM ], [ 4 ] ) - assert feo.get_cell_types_and_number( hexa_prism_grid_invalid ) == ( [ VTK_HEXAGONAL_PRISM ], [ 4 ] ) - assert feo.get_cell_types_and_number( mix_grid ) == ( [ - VTK_TETRA, VTK_HEXAHEDRON, VTK_WEDGE, VTK_PYRAMID, VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM - ], [ 2, 4, 2, 2, 2, 2 ] ) - assert feo.get_cell_types_and_number( mix_grid_invalid ) == ( [ - VTK_TETRA, VTK_HEXAHEDRON, VTK_WEDGE, VTK_PYRAMID, VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM - ], [ 2, 4, 2, 2, 2, 2 ] ) - expected_error: str = f"Invalid type '11' for GEOS is in the mesh. Dying ..." - with pytest.raises( ValueError, match=expected_error ): - feo.get_cell_types_and_number( voxels_grid ) - - def test_reorder_nodes_to_new_mesh( self ): + assert feo.is_cell_to_reorder( volumes[ i ], opt_to_use ) == needs_ordering[ i ] + + def test_cell_ids_to_reorder_by_cell_type( self ): + options_per_grid = { + hexahedrons_grid: opt_hexahedrons_grid, + hexahedrons_grid_invalid: opt_hexahedrons_grid_invalid, + tetras_grid: opt_tetras_grid, + tetras_grid_invalid: opt_tetras_grid_invalid, + pyramids_grid: opt_pyramids_grid, + pyramids_grid_invalid: opt_pyramids_grid_invalid, + wedges_grid: opt_wedges_grid, + wedges_grid_invalid: opt_wedges_grid_invalid, + penta_prism_grid: opt_penta_prism_grid, + penta_prism_grid_invalid: opt_penta_prism_grid_invalid, + hexa_prism_grid: opt_hexa_prism_grid, + hexa_prism_grid_invalid: opt_hexa_prism_grid_invalid, + mix_grid: opt_mix_grid, + mix_grid_invalid: opt_mix_grid_inv + } + # yapf: disable + expected_per_grid = { + hexahedrons_grid: {}, + hexahedrons_grid_invalid: { 12: np.array( [ 1, 3 ] ) }, + tetras_grid: {}, + tetras_grid_invalid: { 10: np.array( [ 1, 3 ] ) }, + pyramids_grid: {}, + pyramids_grid_invalid: { 14: np.array( [ 1, 3 ] ) }, + wedges_grid: {}, + wedges_grid_invalid: { 13: np.array( [ 1, 3 ] ) }, + penta_prism_grid: {}, + penta_prism_grid_invalid: { 15: np.array( [ 1, 3 ] ) }, + hexa_prism_grid: {}, + hexa_prism_grid_invalid: { 16: np.array( [ 1, 3 ] ) }, + mix_grid: {}, + mix_grid_invalid: { 12: np.array( [ 1, 9 ] ), 14: np.array( [ 3 ] ), 10: np.array( [ 5 ] ), + 16: np.array( [ 7 ] ), 13: np.array( [ 11 ] ), 15: np.array( [ 13 ] ) } + } + # yapf: enable + for grid, options in options_per_grid.items(): + result = feo.cell_ids_to_reorder_by_cell_type( grid, options ) + for vtk_type, array in result.items(): + assert np.array_equal( array, expected_per_grid[ grid ][ vtk_type ] ) + + def test_is_polygon_counterclockwise( self ): + # yapf: disable + polygons_available = { + "triangle_ccw": np.array([ [0.0, 0.5, 0.33], [0.0, 1.5, 0.33], [0.0, 1.0, 0.67] ]), + "triangle_cw": np.array([ [0.0, 0.5, 0.33], [0.0, 1.0, 0.67], [0.0, 1.5, 0.33] ]), + "quad_ccw": np.array([ [0.0, 0.5, 0.33], [0.0, 1.5, 0.33], [0.0, 1.5, 0.67], [0.0, 0.5, 0.67] ]), + "quad_cw": np.array([ [0.0, 0.5, 0.33], [0.0, 0.5, 0.67], [0.0, 1.5, 0.67], [0.0, 1.5, 0.33] ]), + "hexagon_ccw": np.array([ [0.0, 0.5, 0.33], [0.0, 1.5, 0.33], [0.0, 2.0, 0.5], [0, 1.5, 0.67], [0, 0.5, 0.67], [0.0, 0.0, 0.5] ]), + "hexagon_cw": np.array([ [0.0, 0.5, 0.33], [0.0, 0.0, 0.5], [0, 0.5, 0.67], [0, 1.5, 0.67], [0.0, 2.0, 0.5], [0.0, 1.5, 0.33] ]), + } + rotation_point = np.array([1, 1, 1]) + + for polygon_name, polygon in polygons_available.items(): + subdivisions = 24 + all_alphas = [ degrees_to_radians( ( 360/subdivisions )*i ) for i in range( subdivisions ) ] + all_betas = all_alphas.copy() + all_gammas = all_alphas.copy() + + polygons_are_clockwise = list() + for alpha_id in range( len( all_alphas ) ): + for beta_id in range( len( all_betas ) ): + for gamma_id in range( len( all_betas ) ): + rotated_polygon = rotate_polygon_around_point( polygon, rotation_point, all_alphas[ alpha_id ], + all_betas[ beta_id ], all_gammas[ gamma_id ] ) + is_clockwise = feo.is_polygon_counterclockwise( rotated_polygon.tolist(), rotation_point.tolist() ) + polygons_are_clockwise.append( is_clockwise ) + + if "ccw" in polygon_name: + assert False not in polygons_are_clockwise + else: + assert True not in polygons_are_clockwise + + concave_polygon = ( ( 0.0, 0.5, 0.33 ), ( 0.0, 1.5, 0.67 ), ( 0.0, 1.5, 0.33 ), ( 0.0, 0.5, 0.67 ) ) + with pytest.raises( ValueError ) as err_msg: + feo.is_polygon_counterclockwise( concave_polygon, rotation_point.tolist() ) + assert str( err_msg.value ) == f"Polygon checked is concave with points: {concave_polygon}." + # yapf: enable + + def test_reordering_functions( self ): + # yapf: disable + points_tetras_coords = [ + ( ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 0.0, 0.0, 1.0 ) ), # valid + ( ( 0.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 0.0, 0.0, 1.0 ) ), # invalid ordering, can be reordered + ] + assert feo.reorder_tetrahedron( points_tetras_coords[ 0 ] ) == ((0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 1.0, 0.0), (0.0, 0.0, 1.0)) + assert feo.reorder_tetrahedron( points_tetras_coords[ 1 ] ) == ((0.0, 0.0, 0.0), (1.0, 1.0, 0.0), (0.0, 0.0, 1.0), (1.0, 0.0, 0.0)) + + points_pyrams_coords = [ + ( ( 1.0, 1.0, 0.0 ), ( 0.0, 1.0, 0.0 ), ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 0.5, 0.5, 1.0 ) ), # valid + ( ( 1.0, 1.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 0.0, 0.0, 0.0 ), ( 0.0, 1.0, 0.0 ), ( 0.5, 0.5, 1.0 ) ), # invalid quad base ordering, can be reordered + ( ( 1.0, 1.0, 0.0 ), ( 0.0, 0.0, 0.0 ), ( 0.0, 1.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 0.5, 0.5, 1.0 ) ), # invalid base definition, cannot reorder + ] + assert feo.reorder_pyramid( points_pyrams_coords[ 0 ] ) == ((1.0, 1.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (0.5, 0.5, 1.0)) + assert feo.reorder_pyramid( points_pyrams_coords[ 1 ] ) == ((1.0, 1.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (0.5, 0.5, 1.0)) + with pytest.raises( ValueError ) as err_msg: + feo.reorder_pyramid( points_pyrams_coords[ 2 ] ) + assert str( err_msg.value ) == "The first 4 points of your pyramid do not represent its base. No ordering possible." + + points_wedges_coords = [ + ( ( 0.0, 0.5, 1.0 ), ( 0.0, 0.0, 0.0 ), ( 0.0, 1.0, 0.0 ), ( 1.0, 0.5, 1.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ) ), # valid + ( ( 0.0, 0.5, 1.0 ), ( 0.0, 1.0, 0.0 ), ( 0.0, 0.0, 0.0 ), ( 1.0, 0.5, 1.0 ), ( 1.0, 1.0, 0.0 ), ( 1.0, 0.0, 0.0 ) ), # two invalid bases ordering in the same way, can be reordered + ( ( 0.0, 0.5, 1.0 ), ( 0.0, 1.0, 0.0 ), ( 0.0, 0.0, 0.0 ), ( 1.0, 0.5, 1.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ) ), # one invalid base ordering creating degenerated wedge + ( ( 0.0, 0.5, 1.0 ), ( 0.0, 0.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 0.0, 1.0, 0.0 ), ( 1.0, 0.5, 1.0 ), ( 1.0, 0.0, 0.0 ) ), # the first / last 3 points do not represent the triangle base + ] + assert feo.reorder_wedge( points_wedges_coords[ 0 ] ) == ((0.0, 0.5, 1.0), (0.0, 1.0, 0.0), (0.0, 0.0, 0.0), (1.0, 0.5, 1.0), (1.0, 1.0, 0.0), (1.0, 0.0, 0.0)) + assert feo.reorder_wedge( points_wedges_coords[ 1 ] ) == ((0.0, 0.5, 1.0), (0.0, 1.0, 0.0), (0.0, 0.0, 0.0), (1.0, 0.5, 1.0), (1.0, 1.0, 0.0), (1.0, 0.0, 0.0)) + with pytest.raises( ValueError ) as err_msg: + feo.reorder_wedge( points_wedges_coords[ 2 ] ) + assert str( err_msg.value ) == ( "When looking at a wedge cell for reordering, we need to construct the two triangle faces that represent the basis. With respect to the centroid of the wedge, the faces are both oriented in the same direction with points0 '((0.0, 0.5, 1.0), (0.0, 1.0, 0.0), (0.0, 0.0, 0.0))' and with points1 '((1.0, 0.5, 1.0), (1.0, 0.0, 0.0), (1.0, 1.0, 0.0))'. When respecting VTK convention, they should be oriented in opposite direction. This create a degenerated wedge that cannot be reordered." ) + with pytest.raises( ValueError ) as err_msg: + feo.reorder_wedge( points_wedges_coords[ 3 ] ) + assert str( err_msg.value ) == ( "When looking at a wedge cell for reordering, we need to construct the two triangle faces that represent the basis. When checking its geometry, the first 3 points'((0.0, 0.5, 1.0), (0.0, 0.0, 0.0), (1.0, 1.0, 0.0))' and/or last 3 points '((0.0, 1.0, 0.0), (1.0, 0.5, 1.0), (1.0, 0.0, 0.0))' cannot represent the wedge basis because they created quad faces that are concave. This create a degenerated wedge that cannot be reordered." ) + + points_hexas_coords = [ + ( ( 0.0, 1.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ), ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 0.0, 1.0 ), ( 0.0, 0.0, 1.0 ) ), # valid + ( ( 0.0, 1.0, 0.0 ), ( 0.0, 1.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( 1.0, 1.0, 0.0 ), ( 0.0, 0.0, 0.0 ), ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), ( 1.0, 0.0, 0.0 ) ), # two invalid bases ordering in the same way, can be reordered + ( ( 0.0, 1.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( 1.0, 1.0, 0.0 ), ( 0.0, 1.0, 0.0 ), ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 0.0, 1.0 ), ( 0.0, 0.0, 1.0 ) ), # one invalid base ordering creating degenerated hexahedron + ( ( 0.0, 1.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 0.0, 1.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.0, 0.0, 1.0 ), ( 0.0, 0.0, 1.0 ) ), # the first / last 4 points do not represent the quad base + ] + assert feo.reorder_hexahedron( points_hexas_coords[ 0 ] ) == ((0.0, 1.0, 0.0), (1.0, 1.0, 0.0), (1.0, 1.0, 1.0), (0.0, 1.0, 1.0), (0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (0.0, 0.0, 1.0)) + assert feo.reorder_hexahedron( points_hexas_coords[ 1 ] ) == ((0.0, 1.0, 0.0), (1.0, 1.0, 0.0), (1.0, 1.0, 1.0), (0.0, 1.0, 1.0), (0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (0.0, 0.0, 1.0)) + with pytest.raises( ValueError ) as err_msg: + feo.reorder_hexahedron( points_hexas_coords[ 2 ] ) + assert str( err_msg.value ) == ( "When looking at a hexahedron cell for reordering, we need to construct two quad faces that represent two faces that do not have a point common. With respect to the centroid of the hexahedron, the faces are not both oriented in the same direction with points0 '((0.0, 1.0, 1.0), (0.0, 1.0, 0.0), (1.0, 1.0, 0.0), (1.0, 1.0, 1.0))' and with points1 '((0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (0.0, 0.0, 1.0))'. When respecting VTK convention, they both should be oriented in the same direction. This create a degenerated hexahedron that cannot be reordered." ) + with pytest.raises( ValueError ) as err_msg: + feo.reorder_hexahedron( points_hexas_coords[ 3 ] ) + assert str( err_msg.value ) == ( "When looking at a hexahedron cell for reordering, we need to construct two quad faces that represent two faces that do not have a point common. When checking its geometry, the first 4 points '((0.0, 1.0, 0.0), (1.0, 1.0, 1.0), (0.0, 1.0, 1.0), (1.0, 1.0, 0.0))' and/or last 4 points '((0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.0, 1.0), (0.0, 0.0, 1.0))' cannot represent two hexahedron quad faces because they are concave. This create a degenerated hexahedron that cannot be reordered." ) + + points_penta_prisms_coords = [ + ( ( 0.0, 1.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 1.5, 1.0, 0.5 ), ( 0.5, 1.0, 1.0 ), ( -0.5, 1.0, 0.5 ), ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.5, 0.0, 0.5 ), ( 0.5, 0.0, 1.0 ), ( -0.5, 0.0, 0.5 ) ), # valid + ( ( 0.0, 1.0, 0.0 ), ( -0.5, 1.0, 0.5 ), ( 0.5, 1.0, 1.0 ), ( 1.5, 1.0, 0.5 ), ( 1.0, 1.0, 0.0 ), ( 0.0, 0.0, 0.0 ), ( -0.5, 0.0, 0.5 ), ( 0.5, 0.0, 1.0 ), ( 1.5, 0.0, 0.5 ), ( 1.0, 0.0, 0.0 ) ), # two invalid bases ordering in the same way, can be reordered + ( ( 0.0, 1.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 1.5, 1.0, 0.5 ), ( 0.5, 1.0, 1.0 ), ( -0.5, 1.0, 0.5 ), ( 0.0, 0.0, 0.0 ), ( -0.5, 0.0, 0.5 ), ( 0.5, 0.0, 1.0 ), ( 1.5, 0.0, 0.5 ), ( 1.0, 0.0, 0.0 ) ), # one invalid base ordering creating degenerated pentagonal prism + ( ( 0.0, 1.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 1.5, 1.0, 0.5 ), ( -0.5, 1.0, 0.5 ), ( 0.5, 1.0, 1.0 ), ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.5, 0.0, 0.5 ), ( 0.5, 0.0, 1.0 ), ( -0.5, 0.0, 0.5 ) ), # the first / last 5 points do not represent the pentagon base + ] + assert feo.reorder_pentagonal_prism( points_penta_prisms_coords[ 0 ] ) == ((0.0, 1.0, 0.0), (1.0, 1.0, 0.0), (1.5, 1.0, 0.5), (0.5, 1.0, 1.0), (-0.5, 1.0, 0.5), (0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.5, 0.0, 0.5), (0.5, 0.0, 1.0), (-0.5, 0.0, 0.5)) + assert feo.reorder_pentagonal_prism( points_penta_prisms_coords[ 1 ] ) == ((0.0, 1.0, 0.0), (1.0, 1.0, 0.0), (1.5, 1.0, 0.5), (0.5, 1.0, 1.0), (-0.5, 1.0, 0.5), (0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.5, 0.0, 0.5), (0.5, 0.0, 1.0), (-0.5, 0.0, 0.5)) + with pytest.raises( ValueError ) as err_msg: + feo.reorder_pentagonal_prism( points_penta_prisms_coords[ 2 ] ) + assert str( err_msg.value ) == ( "When looking at a pentagonal prism cell for reordering, we need to construct the two pentagon faces that represent the basis. With respect to the centroid of the wedge, the faces are oriented in opposite direction with points0 '((0.0, 1.0, 0.0), (-0.5, 1.0, 0.5), (0.5, 1.0, 1.0), (1.5, 1.0, 0.5), (1.0, 1.0, 0.0))' and with points1 '((0.0, 0.0, 0.0), (-0.5, 0.0, 0.5), (0.5, 0.0, 1.0), (1.5, 0.0, 0.5), (1.0, 0.0, 0.0))'. When respecting VTK convention, they should be oriented in the same direction. This create a degenerated pentagonal prism that cannot be reordered." ) + with pytest.raises( ValueError ) as err_msg: + feo.reorder_pentagonal_prism( points_penta_prisms_coords[ 3 ] ) + assert str( err_msg.value ) == ( "When looking at a pentagonal prism cell for reordering, we need to construct the two pentagon faces that represent the basis. When checking its geometry, the first 5 points'((0.0, 1.0, 0.0), (0.5, 1.0, 1.0), (-0.5, 1.0, 0.5), (1.5, 1.0, 0.5), (1.0, 1.0, 0.0))' and/or last 5 points '((0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.5, 0.0, 0.5), (0.5, 0.0, 1.0), (-0.5, 0.0, 0.5))' cannot represent the pentagonal prism basis because they created pentagon faces that are concave. This create a degenerated pentagonal prism that cannot be reordered." ) + + points_hexa_prisms_coords = [ + ( ( 0.0, 1.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 1.5, 1.0, 0.5 ), ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ), ( -0.5, 1.0, 0.5 ), ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.5, 0.0, 0.5 ), ( 1.0, 0.0, 1.0 ), ( 0.0, 0.0, 1.0 ), ( -0.5, 0.0, 0.5 ) ), # valid + ( ( 0.0, 1.0, 0.0 ), ( -0.5, 1.0, 0.5 ), ( 0.0, 1.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( 1.5, 1.0, 0.5 ), ( 1.0, 1.0, 0.0 ), ( 0.0, 0.0, 0.0 ), ( -0.5, 0.0, 0.5 ), ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), ( 1.5, 0.0, 0.5 ), ( 1.0, 0.0, 0.0 ) ), # two invalid bases ordering in the same way, can be reordered + ( ( 0.0, 1.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 1.5, 1.0, 0.5 ), ( 1.0, 1.0, 1.0 ), ( 0.0, 1.0, 1.0 ), ( -0.5, 1.0, 0.5 ), ( 0.0, 0.0, 0.0 ), ( -0.5, 0.0, 0.5 ), ( 0.0, 0.0, 1.0 ), ( 1.0, 0.0, 1.0 ), ( 1.5, 0.0, 0.5 ), ( 1.0, 0.0, 0.0 ) ), # one invalid base ordering creating degenerated hexagonal prism + ( ( 0.0, 1.0, 0.0 ), ( 1.0, 1.0, 0.0 ), ( 1.5, 1.0, 0.5 ), ( 0.0, 1.0, 1.0 ), ( 1.0, 1.0, 1.0 ), ( -0.5, 1.0, 0.5 ), ( 0.0, 0.0, 0.0 ), ( 1.0, 0.0, 0.0 ), ( 1.5, 0.0, 0.5 ), ( 1.0, 0.0, 1.0 ), ( 0.0, 0.0, 1.0 ), ( -0.5, 0.0, 0.5 ) ), # the first / last 5 points do not represent the hexagon base + ] + assert feo.reorder_hexagonal_prism( points_hexa_prisms_coords[ 0 ] ) == ((0.0, 1.0, 0.0), (1.0, 1.0, 0.0), (1.5, 1.0, 0.5), (1.0, 1.0, 1.0), (0.0, 1.0, 1.0), (-0.5, 1.0, 0.5), (0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.5, 0.0, 0.5), (1.0, 0.0, 1.0), (0.0, 0.0, 1.0), (-0.5, 0.0, 0.5)) + assert feo.reorder_hexagonal_prism( points_hexa_prisms_coords[ 1 ] ) == ((1.0, 1.0, 0.0), (1.5, 1.0, 0.5), (1.0, 1.0, 1.0), (0.0, 1.0, 1.0), (-0.5, 1.0, 0.5), (0.0, 1.0, 0.0), (1.0, 0.0, 0.0), (1.5, 0.0, 0.5), (1.0, 0.0, 1.0), (0.0, 0.0, 1.0), (-0.5, 0.0, 0.5), (0.0, 0.0, 0.0)) + with pytest.raises( ValueError ) as err_msg: + feo.reorder_hexagonal_prism( points_hexa_prisms_coords[ 2 ] ) + assert str( err_msg.value ) == ( "When looking at a hexagonal prism cell for reordering, we need to construct the two hexagon faces that represent the basis. With respect to the centroid of the wedge, the faces are oriented in opposite direction with points0 '((0.0, 1.0, 0.0), (-0.5, 1.0, 0.5), (0.0, 1.0, 1.0), (1.0, 1.0, 1.0), (1.5, 1.0, 0.5), (1.0, 1.0, 0.0))' and with points1 '((0.0, 0.0, 0.0), (-0.5, 0.0, 0.5), (0.0, 0.0, 1.0), (1.0, 0.0, 1.0), (1.5, 0.0, 0.5), (1.0, 0.0, 0.0))'. When respecting VTK convention, they should be oriented in the same direction. This create a degenerated hexagonal prism that cannot be reordered." ) + with pytest.raises( ValueError ) as err_msg: + feo.reorder_hexagonal_prism( points_hexa_prisms_coords[ 3 ] ) + assert str( err_msg.value ) == ( "When looking at a hexagonal prism cell for reordering, we need to construct the two hexagon faces that represent the basis. When checking its geometry, the first 6 points'((0.0, 1.0, 0.0), (-0.5, 1.0, 0.5), (1.0, 1.0, 1.0), (0.0, 1.0, 1.0), (1.5, 1.0, 0.5), (1.0, 1.0, 0.0))' and/or last 6 points '((0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.5, 0.0, 0.5), (1.0, 0.0, 1.0), (0.0, 0.0, 1.0), (-0.5, 0.0, 0.5))' cannot represent the hexagonal prism basis because they created hexagon faces that are concave. This create a degenerated hexagonal prism that cannot be reordered." ) + # yapf: enable + + def test_cell_point_ids_ordering_method( self ): + # yapf: disable + assert feo.cell_point_ids_ordering_method( tetras_grid.GetCell( 0 ), 10 ) == ( ( 0, 1, 2, 3 ), False ) + assert feo.cell_point_ids_ordering_method( tetras_grid_invalid.GetCell( 0 ), 10 ) == ( ( 0, 1, 2, 3 ), False ) + assert feo.cell_point_ids_ordering_method( tetras_grid.GetCell( 1 ), 10 ) == ( ( 0, 1, 2, 3 ), False ) + assert feo.cell_point_ids_ordering_method( tetras_grid_invalid.GetCell( 1 ), 10 ) == ( ( 0, 1, 3, 2 ), True ) + + assert feo.cell_point_ids_ordering_method( pyramids_grid.GetCell( 0 ), 14 ) == ( ( 0, 1, 2, 3, 4 ), False ) + assert feo.cell_point_ids_ordering_method( pyramids_grid_invalid.GetCell( 0 ), 14 ) == ( ( 0, 1, 2, 3, 4 ), False ) + assert feo.cell_point_ids_ordering_method( pyramids_grid.GetCell( 1 ), 14 ) == ( ( 0, 1, 2, 3, 4 ), False ) + assert feo.cell_point_ids_ordering_method( pyramids_grid_invalid.GetCell( 1 ), 14 ) == ( ( 0, 3, 2, 1, 4 ), True ) + + assert feo.cell_point_ids_ordering_method( wedges_grid.GetCell( 0 ), 13 ) == ( ( 0, 1, 2, 3, 4, 5 ), False ) + assert feo.cell_point_ids_ordering_method( wedges_grid_invalid.GetCell( 0 ), 13 ) == ( ( 0, 1, 2, 3, 4, 5 ), False ) + assert feo.cell_point_ids_ordering_method( wedges_grid.GetCell( 1 ), 13 ) == ( ( 0, 1, 2, 3, 4, 5 ), False ) + assert feo.cell_point_ids_ordering_method( wedges_grid_invalid.GetCell( 1 ), 13 ) == ( ( 0, 2, 1, 3, 5, 4 ), True ) + + assert feo.cell_point_ids_ordering_method( hexahedrons_grid.GetCell( 0 ), 12 ) == ( ( 0, 1, 2, 3, 4, 5, 6, 7 ), False ) + assert feo.cell_point_ids_ordering_method( hexahedrons_grid_invalid.GetCell( 0 ), 12 ) == ( ( 0, 1, 2, 3, 4, 5, 6, 7 ), False ) + assert feo.cell_point_ids_ordering_method( hexahedrons_grid.GetCell( 1 ), 12 ) == ( ( 0, 1, 2, 3, 4, 5, 6, 7 ), False ) + assert feo.cell_point_ids_ordering_method( hexahedrons_grid_invalid.GetCell( 1 ), 12 ) == ( ( 0, 3, 2, 1, 4, 7, 6, 5 ), True ) + + assert feo.cell_point_ids_ordering_method( penta_prism_grid.GetCell( 0 ), 15 ) == ( ( 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ), False ) + assert feo.cell_point_ids_ordering_method( penta_prism_grid_invalid.GetCell( 0 ), 15 ) == ( ( 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ), False ) + assert feo.cell_point_ids_ordering_method( penta_prism_grid.GetCell( 1 ), 15 ) == ( ( 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 ), False ) + assert feo.cell_point_ids_ordering_method( penta_prism_grid_invalid.GetCell( 1 ), 15 ) == ( ( 0, 4, 3, 2, 1, 5, 9, 8, 7, 6 ), True ) + + assert feo.cell_point_ids_ordering_method( hexa_prism_grid.GetCell( 0 ), 16 ) == ( ( 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 ), False ) + assert feo.cell_point_ids_ordering_method( hexa_prism_grid_invalid.GetCell( 0 ), 16 ) == ( ( 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 ), False ) + assert feo.cell_point_ids_ordering_method( hexa_prism_grid.GetCell( 1 ), 16 ) == ( ( 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 ), False ) + assert feo.cell_point_ids_ordering_method( hexa_prism_grid_invalid.GetCell( 1 ), 16 ) == ( ( 5, 4, 3, 2, 1, 0, 11, 10, 9, 8, 7, 6 ), True ) + # yapf: enable + + def test_reorder_points_to_new_mesh( self ): options = opt( out, cell_names, "negative" ) # single element grids except voxels because it is an invalid cell type for GEOS - grid_cell_type = { - hexahedrons_grid_invalid: VTK_HEXAHEDRON, - tetras_grid_invalid: VTK_TETRA, - pyramids_grid_invalid: VTK_PYRAMID, - wedges_grid_invalid: VTK_WEDGE, - penta_prism_grid_invalid: VTK_PENTAGONAL_PRISM, - hexa_prism_grid_invalid: VTK_HEXAGONAL_PRISM + grid_cell_type_and_expected_volumes = { + hexahedrons_grid_invalid: ( VTK_HEXAHEDRON, [ 1.0, 1.0, 1.0, 1.0 ] ), + tetras_grid_invalid: ( VTK_TETRA, [ 0.167, 0.167, 0.167, 0.167 ] ), + pyramids_grid_invalid: ( VTK_PYRAMID, [ 0.333, 0.333, 0.333, 0.333 ] ), + wedges_grid_invalid: ( VTK_WEDGE, [ 0.5, 0.5, 0.5, 0.5 ] ), + penta_prism_grid_invalid: ( VTK_PENTAGONAL_PRISM, [ 1.25, 1.25, 1.25, 1.25 ] ), + hexa_prism_grid_invalid: ( VTK_HEXAGONAL_PRISM, [ 1.5, 1.5, 1.5, 1.5 ] ) } - for grid, cell_type in grid_cell_type.items(): + for grid, cell_type_and_volumes in grid_cell_type_and_expected_volumes.items(): new_invalid = vtkUnstructuredGrid() new_invalid.DeepCopy( grid ) - not_use_invalid, reorder_stats = feo.reorder_nodes_to_new_mesh( new_invalid, options ) + new_corrected, reorder_stats = feo.reorder_points_to_new_mesh( new_invalid, options ) expected = { - "Types reordered": [ VTK_TYPE_TO_NAME[ cell_type ] ], + "Types reordered": [ VTK_TYPE_TO_NAME[ cell_type_and_volumes[ 0 ] ] ], "Number of cells reordered": [ 2 ], - "Types non reordered": [ VTK_TYPE_TO_NAME[ cell_type ] ], - "Number of cells non reordered": [ 2 ] + 'Types non reordered because of errors': list(), + 'Number of cells non reordered because of errors': list(), + 'Error message given': list() } for prop in expected.keys(): assert reorder_stats[ prop ] == expected[ prop ] + volumes = feo.compute_mesh_cells_volume( new_corrected ) + expected_volumes = cell_type_and_volumes[ 1 ] + for i in range( len( volumes ) ): + assert round( float( volumes[ i ] ), 3 ) == expected_volumes[ i ] # mix elements grid mix_invalid = vtkUnstructuredGrid() mix_invalid.DeepCopy( mix_grid_invalid ) - not_use_invalid, mix_stats = feo.reorder_nodes_to_new_mesh( mix_invalid, options ) + new_mix_corrected, mix_stats = feo.reorder_points_to_new_mesh( mix_invalid, options ) expected = { - "Types reordered": [ VTK_TYPE_TO_NAME[ cell_type ] for cell_type in [ 12, 15, 16, 14, 10, 13 ] ], - "Number of cells reordered": [ 2, 1, 1, 1, 1, 1 ], - "Types non reordered": [ VTK_TYPE_TO_NAME[ cell_type ] for cell_type in [ 12, 15, 16, 14, 10, 13 ] ], - "Number of cells non reordered": [ 2, 1, 1, 1, 1, 1 ] + "Types reordered": [ VTK_TYPE_TO_NAME[ cell_type ] for cell_type in [ 10, 12, 13, 14, 15, 16 ] ], + "Number of cells reordered": [ 1, 2, 1, 1, 1, 1 ], + "Types non reordered because ordering is already correct": list(), + "Number of cells non reordered because ordering is already correct": list(), + "Types non reordered because of errors": list(), + "Number of cells non reordered because of errors": list(), + "Error message given": list() } for prop in expected.keys(): assert mix_stats[ prop ] == expected[ prop ] + volumes = feo.compute_mesh_cells_volume( new_mix_corrected ) + expected_volumes = [ 1.0, 1.0, 0.333, 0.333, 0.167, 0.167, 1.5, 1.5, 1.0, 1.0, 0.25, 0.25, 1.25, 1.25 ] + for i in range( len( volumes ) ): + assert round( float( volumes[ i ] ), 3 ) == expected_volumes[ i ] + + # mix elements grid + options = opt( out, cell_names, "all" ) + mix_degenerated = vtkUnstructuredGrid() + mix_degenerated.DeepCopy( mix_grid_degenerated ) + not_used, mix_degen_stats = feo.reorder_points_to_new_mesh( mix_degenerated, options ) + # yapf: disable + expected = { + "Types reordered": list(), + "Number of cells reordered": list(), + "Types non reordered because ordering is already correct": [ VTK_TYPE_TO_NAME[ cell_type ] for cell_type in [ 10, 14 ] ], + "Number of cells non reordered because ordering is already correct": [ 2, 2 ], + "Types non reordered because of errors": [ VTK_TYPE_TO_NAME[ cell_type ] for cell_type in [ 12, 13, 15, 16 ] ], + "Number of cells non reordered because of errors": [ 4, 2, 2, 2 ], + "Error message given": list() + } + # yapf: enable + for prop in expected.keys(): + if prop == "Error message given": + assert len( mix_degen_stats[ prop ] ) == 4 + else: + assert mix_degen_stats[ prop ] == expected[ prop ] def test_fix_elements_orderings_execution( self ): # for mix_grid_invalid mesh, checks that reordered mesh was created and that reoredring_stats are valid @@ -781,7 +1221,7 @@ def test_fix_elements_orderings_execution( self ): invalidTest = False command = [ "python", MESH_DOCTOR_FILEPATH, "-v", "-i", test_file.output, "fix_elements_orderings", "--cell_names", - ",".join( map( str, cell_names ) ), "--volume_to_reorder", "negative", "--data-mode", "binary", "--output", + ",".join( map( str, cell_names ) ), "--volume_to_reorder", "all", "--data-mode", "binary", "--output", filepath_reordered_mesh ] try: @@ -794,17 +1234,57 @@ def test_fix_elements_orderings_execution( self ): matches = re.findall( pattern, raw_stderr ) no_log = "\n".join( matches ) reordering_stats: str = no_log[ no_log.index( "Number of cells reordered" ): ] - expected_stats: str = ( "Number of cells reordered:\n" + "\tCellType\tNumber\n" + f"\tHexahedron\t\t2\n" + - f"\tPrism5\t\t1\n" + f"\tPrism6\t\t1\n" + f"\tPyramid\t\t1\n" + - f"\tTetrahedron\t\t1\n" + f"\tWedge\t\t1\n" + "Number of cells non reordered:\n" + - "\tCellType\tNumber\n" + f"\tHexahedron\t\t2\n" + f"\tPrism5\t\t1\n" + - f"\tPrism6\t\t1\n" + f"\tPyramid\t\t1\n" + f"\tTetrahedron\t\t1\n" + - f"\tWedge\t\t1" ) + # yapf: disable + expected_stats: str = ( "Number of cells reordered:\n" + + "\tCellType\tNumber\n" + + "\tTetrahedron\t\t2\n" + + "\tHexahedron\t\t4\n" + + "\tWedge\t\t2\n" + + "\tPyramid\t\t2\n" + + "\tPrism5\t\t2\n" + + "\tPrism6\t\t2" ) + # yapf: enable assert reordering_stats == expected_stats except Exception as e: logging.error( "Invalid command input. Test has failed." ) logging.error( e ) invalidTest = True os.remove( test_file.output ) + + write_mesh( mix_grid_degenerated, test_file2 ) + command = [ + "python", MESH_DOCTOR_FILEPATH, "-v", "-i", test_file2.output, "fix_elements_orderings", "--cell_names", + ",".join( map( str, cell_names ) ), "--volume_to_reorder", "all", "--data-mode", "binary", "--output", + filepath_reordered_mesh2 + ] + try: + result = subprocess.run( command, shell=True, stderr=subprocess.PIPE, universal_newlines=True ) + os.remove( filepath_reordered_mesh2 ) + stderr = result.stderr + assert result.returncode == 0 + raw_stderr = r"{}".format( stderr ) + pattern = r"\[.*?\]\[.*?\] (.*)" + matches = re.findall( pattern, raw_stderr ) + no_log = "\n".join( matches ) + reordering_stats: str = no_log[ no_log.index( "Number of cells non reordered because of errors" ): ] + # yapf: disable + expected_stats: str = ( "Number of cells non reordered because of errors:\n" + + "\tCellType\tNumber\n" + + "\tHexahedron\t\t4\n" + + "\tError message: When looking at a hexahedron cell for reordering, we need to construct two quad faces that represent two faces that do not have a point common. When checking its geometry, the first 4 points '((1.0, 0.0, 0.0), (1.0, 1.0, 0.0), (2.0, 0.0, 0.0), (2.0, 1.0, 0.0))' and/or last 4 points '((1.0, 0.0, 1.0), (2.0, 1.0, 1.0), (2.0, 0.0, 1.0), (1.0, 1.0, 1.0))' cannot represent two hexahedron quad faces because they are concave. This create a degenerated hexahedron that cannot be reordered.\n" + + "\tWedge\t\t2\n" + + "\tError message: When looking at a wedge cell for reordering, we need to construct the two triangle faces that represent the basis. With respect to the centroid of the wedge, the faces are both oriented in the same direction with points0 '((4.0, 0.0, 1.0), (4.5, -0.5, 1.0), (5.0, 0.0, 1.0))' and with points1 '((5.0, 0.0, 0.0), (4.5, -0.5, 0.0), (4.0, 0.0, 0.0))'. When respecting VTK convention, they should be oriented in opposite direction. This create a degenerated wedge that cannot be reordered.\n" + + "\tPrism5\t\t2\n" + + "\tError message: When looking at a pentagonal prism cell for reordering, we need to construct the two pentagon faces that represent the basis. With respect to the centroid of the wedge, the faces are oriented in opposite direction with points0 '((5.0, 0.0, 1.0), (5.0, 1.0, 1.0), (5.5, 1.5, 1.0), (6.0, 0.5, 1.0), (5.5, -0.5, 1.0))' and with points1 '((6.0, 0.5, 2.0), (5.5, -0.5, 2.0), (5.0, 0.0, 2.0), (5.0, 1.0, 2.0), (5.5, 1.5, 2.0))'. When respecting VTK convention, they should be oriented in the same direction. This create a degenerated pentagonal prism that cannot be reordered.\n" + + "\tPrism6\t\t2\n" + + "\tError message: When looking at a hexagonal prism cell for reordering, we need to construct the two hexagon faces that represent the basis. With respect to the centroid of the wedge, the faces are oriented in opposite direction with points0 '((2.0, 0.0, 1.0), (2.0, 1.0, 1.0), (2.5, 1.5, 1.0), (3.0, 1.0, 1.0), (3.0, 0.0, 1.0), (2.5, -0.5, 1.0))' and with points1 '((2.0, 1.0, 2.0), (2.5, 1.5, 2.0), (3.0, 1.0, 2.0), (3.0, 0.0, 2.0), (2.5, -0.5, 2.0), (2.0, 0.0, 2.0))'. When respecting VTK convention, they should be oriented in the same direction. This create a degenerated hexagonal prism that cannot be reordered." ) + # yapf: enable + assert reordering_stats == expected_stats + except Exception as e: + logging.error( "Invalid command input. Test has failed." ) + logging.error( e ) + invalidTest = True + os.remove( test_file2.output ) + if invalidTest: - raise ValueError( "test_fix_elements_orderings_execution has failed." ) \ No newline at end of file + raise ValueError( "test_fix_elements_orderings_execution has failed." )