diff --git a/docs/geos-mesh.rst b/docs/geos-mesh.rst index cbdea0e9..243592ed 100644 --- a/docs/geos-mesh.rst +++ b/docs/geos-mesh.rst @@ -117,33 +117,36 @@ 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] [--Hexahedron 1,6,5,4,7,0,2,3] [--Prism5 8,2,0,7,6,9,5,1,4,3] - [--Prism6 11,2,8,10,5,0,9,7,6,1,4,3] [--Pyramid 3,4,0,2,1] - [--Tetrahedron 2,0,3,1] [--Voxel 1,6,5,4,7,0,2,3] - [--Wedge 3,5,4,0,2,1] --output OUTPUT [--data-mode binary, ascii] + 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 - --Hexahedron 1,6,5,4,7,0,2,3 - [list of integers]: node permutation for "Hexahedron". - --Prism5 8,2,0,7,6,9,5,1,4,3 - [list of integers]: node permutation for "Prism5". - --Prism6 11,2,8,10,5,0,9,7,6,1,4,3 - [list of integers]: node permutation for "Prism6". - --Pyramid 3,4,0,2,1 [list of integers]: node permutation for "Pyramid". - --Tetrahedron 2,0,3,1 [list of integers]: node permutation for "Tetrahedron". - --Voxel 1,6,5,4,7,0,2,3 [list of integers]: node permutation for "Voxel". - --Wedge 3,5,4,0,2,1 [list of integers]: node permutation for "Wedge". + --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, you can use this command: + +.. code-block:: + + $ 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`` """"""""""""""""" @@ -340,8 +343,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..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,65 +1,656 @@ -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 math import sqrt +from dataclasses import dataclass +from itertools import permutations +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, 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_names_to_reorder: tuple[ str ] + volume_to_reorder: str @dataclass( frozen=True ) class Result: output: str - unchanged_cell_types: FrozenSet[ int ] + reordering_stats: dict[ str, list[ int ] ] -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: +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, + "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 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 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: + mesh (vtkDataSet): A vtk grid. + options (Options): Options defined by the user. + + Returns: + dict[ int, np.array ]: { vtk_type1: np.array, ..., vtk_typeN: np.array } + """ + 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 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: + points (Points3D): A set of points in 3D coordinates. + ref_pt (np.array): A point in 3D coordinates. + + Raises: + ValueError: "Polygon checked is concave with points: {points}" + + Returns: + bool: True if counterclockwise, False if clockwise. + """ + 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 + + +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. + + 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: + 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[ tuple[ int ], bool ]: ( The ordering method, if_ordering_is_different ) + """ + 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 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() + } + # 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_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 ) ) + 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_points_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 3b581b75..da1ceacc 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py @@ -1,4 +1,4 @@ -import os.path +import os import logging from dataclasses import dataclass from typing import Iterator, Optional @@ -36,6 +36,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 = 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 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 os.path.exists( filepath ): + logging.error( f"The file specified does not exist." ) + return False + else: + 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 os.path.exists( filepath ): + logging.error( f"A file with the same name already exists. No over-writing possible." ) + return False + return True + + def has_invalid_field( mesh: vtkUnstructuredGrid, invalid_fields: list[ str ] ) -> bool: """Checks if a mesh contains at least a data arrays within its cell, field or point data having a certain name. If so, returns True, else False. @@ -101,8 +134,8 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid: If first guess does not work, eventually all the others reader available will be tested. :return: A unstructured grid. """ - if not os.path.exists( vtk_input_file ): - err_msg: str = f"Invalid file path. Could not read \"{vtk_input_file}\". Dying..." + if not is_filepath_valid( vtk_input_file ): + 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 ] @@ -118,7 +151,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 ) @@ -148,9 +181,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 + if not is_filepath_valid( vtk_output.output, False ): + 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 ) @@ -158,7 +192,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/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..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,52 +1,31 @@ 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 geos.mesh.doctor.checks.fix_elements_orderings import Options, Result, NAME_TO_VTK_TYPE 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_NAMES = "cell_names" +__CELL_NAMES_CHOICES = list( NAME_TO_VTK_TYPE.keys() ) -__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 = "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 key, vtk_key in __CELL_TYPE_MAPPING.items(): - tmp = list( range( __CELL_TYPE_SUPPORT_SIZE[ vtk_key ] ) ) - random.Random( 4 ).shuffle( tmp ) - p.add_argument( '--' + key, - type=str, - metavar=",".join( map( str, tmp ) ), - default=None, - required=False, - help=f"[list of integers]: node permutation for \"{key}\"." ) + p.add_argument( '--' + __CELL_NAMES, + type=str, + metavar=", ".join( map( str, __CELL_NAMES_CHOICES ) ), + default=", ".join( map( str, __CELL_NAMES_CHOICES ) ), + 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'." ) ) vtk_output_parsing.fill_vtk_output_subparser( p ) @@ -56,26 +35,47 @@ 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 ] - 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." - logging.error( err_msg ) - raise ValueError( err_msg ) - cell_type_to_ordering[ vtk_key ] = tmp + 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}:" + + f" {__CELL_NAMES_CHOICES}." ) 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}:" + + f" {__VOLUME_TO_REORDER_CHOICES}." ) + return Options( vtk_output=vtk_output, + cell_names_to_reorder=cell_names_to_reorder, + 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." ) + 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 new file mode 100644 index 00000000..c53c4e8f --- /dev/null +++ b/geos-mesh/tests/test_fix_elements_orderings.py @@ -0,0 +1,1290 @@ +import os +import re +import pytest +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 +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 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: + 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 ) +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, 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, 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() ) +""" +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" ] ) + +opt_hexahedrons_grid = opt( test_file, [ "Hexahedron" ], "negative" ) +opt_hexahedrons_grid_invalid = opt( test_file, [ "Hexahedron" ], "negative" ) +""" +4 tetrahedrons +""" +points_tetras: vtkPoints = vtkPoints() +# 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 ) + +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" ] ) + +opt_tetras_grid = opt( test_file, [ "Tetrahedron" ], "negative" ) +opt_tetras_grid_invalid = opt( test_file, [ "Tetrahedron" ], "negative" ) +""" +4 pyramids +""" +points_pyramids: vtkPoints = vtkPoints() +# 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 ) + +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" ] ) + +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() +# 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 ) + +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 ) + +opt_voxels_grid = opt( test_file, [ "Voxel" ], "negative" ) +opt_voxels_grid_invalid = opt( test_file, [ "Voxel" ], "negative" ) +""" +4 wedges +""" +points_wedges: vtkPoints = vtkPoints() +# 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 ) + +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" ] ) + +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() +# 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 ) + +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" ] ) + +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() +# 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 ) + +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" ] ) + +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() +# 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 ) + +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 ] ] ) + +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: + + 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 ), ( 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 ] + + # 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 ): + # yapf: disable + grid_volumes = { + hexahedrons_grid: [ 1.0, 1.0, 1.0, 1.0 ], + 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.5, 0.5, -0.5 ], + penta_prism_grid: [ 1.25, 1.25, 1.25, 1.25 ], + 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, -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, -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 ) ): + 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 ) ] + } + 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 ) ): + 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_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_and_volumes in grid_cell_type_and_expected_volumes.items(): + new_invalid = vtkUnstructuredGrid() + new_invalid.DeepCopy( grid ) + new_corrected, reorder_stats = feo.reorder_points_to_new_mesh( new_invalid, options ) + expected = { + "Types reordered": [ VTK_TYPE_TO_NAME[ cell_type_and_volumes[ 0 ] ] ], + "Number of cells 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 ) + 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 [ 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 + write_mesh( mix_grid_invalid, test_file ) + 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", "all", "--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" ): ] + # 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." )