diff --git a/docs/geos-mesh.rst b/docs/geos-mesh.rst index 739d493..243592e 100644 --- a/docs/geos-mesh.rst +++ b/docs/geos-mesh.rst @@ -117,41 +117,35 @@ Cells with negative volumes will typically be an issue for ``geos`` and should b """""""""""""""""""""""""" It sometimes happens that an exported mesh does not abide by the ``vtk`` orderings. -The ``fix_elements_orderings`` module can rearrange the nodes of given types of elements. +The ``fix_elements_orderings`` module can rearrange the nodes of given types of elements to respect VTK convention. This can be convenient if you cannot regenerate the mesh. .. code-block:: $ python src/geos/mesh/doctor/mesh_doctor.py fix_elements_orderings --help - usage: mesh_doctor.py fix_elements_orderings [-h] [--Tetrahedron 2,0,3,1] [--Pyramid 3,4,0,2,1] [--Wedge 3,5,4,0,2,1] - [--Hexahedron 1,6,5,4,7,0,2,3] [--Prism5 8,2,0,7,6,9,5,1,4,3] [--Prism6 11,2,8,10,5,0,9,7,6,1,4,3] - --volume_to_reorder all,positive,negative --output OUTPUT [--data-mode binary, ascii] + usage: mesh_doctor.py fix_elements_orderings [-h] [--cell_names Hexahedron, Tetrahedron, Pyramid, Wedge, Prism5, Prism6] + [--volume_to_reorder all, positive, negative] + --output OUTPUT [--data-mode binary, ascii] options: - -h, --help show this help message and exit - --Tetrahedron 2,0,3,1 [list of integers]: node permutation for "Tetrahedron". - --Pyramid 3,4,0,2,1 [list of integers]: node permutation for "Pyramid". - --Wedge 3,5,4,0,2,1 [list of integers]: node permutation for "Wedge". - --Hexahedron 1,6,5,4,7,0,2,3 - [list of integers]: node permutation for "Hexahedron". - --Prism5 8,2,0,7,6,9,5,1,4,3 - [list of integers]: node permutation for "Prism5". - --Prism6 11,2,8,10,5,0,9,7,6,1,4,3 - [list of integers]: node permutation for "Prism6". - --volume_to_reorder all,positive,negative - [str]: Select which element volume is invalid and needs reordering. 'all' will allow reordering of nodes for every element, regarding of their volume. - 'positive' or 'negative' will only reorder the element with the corresponding volume. - --output OUTPUT [string]: The vtk output file destination. - --data-mode binary, ascii [string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary. + -h, --help show this help message and exit + --cell_names Hexahedron, Tetrahedron, Pyramid, Wedge, Prism5, Prism6 + [list of str]: Cell names that can be reordered in your grid. You can use multiple names.Defaults to all cell names being used. + --volume_to_reorder all, positive, negative + [str]: Select which element volume is invalid and needs reordering. 'all' will allow reordering of nodes for every element, regarding of their volume. + 'positive' or 'negative' will only reorder the element with the corresponding volume. Defaults to 'negative'. + --output OUTPUT [string]: The vtk output file destination. + --data-mode binary, ascii + [string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary. For example, assume that you have a mesh that contains 3 different cell types like Hexahedrons, Pyramids and Tetrahedrons. After checking ``element_volumes``, you found that all your Hexahedrons and half of your Tetrahedrons have a negative volume. -To correct that, assuming you know the correct node ordering to correct these volumes, you can use this command: +To correct that, you can use this command: .. code-block:: - $ python src/geos/mesh/doctor/mesh_doctor.py -i /path/to/your/mesh/file fix_elements_orderings --Hexahedron 1,6,5,4,7,0,2,3 - --Tetrahedron 2,0,3,1 --volume_to_reorder negative --output /new/path/for/your/new/mesh/reordered/file --data-mode binary + $ python src/geos/mesh/doctor/mesh_doctor.py -i /path/to/your/mesh/file fix_elements_orderings --cell_names Hexahedron,Tetrahedron + --volume_to_reorder negative --output /new/path/for/your/new/mesh/reordered/file --data-mode binary ``generate_cube`` """"""""""""""""" diff --git a/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py b/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py index bc693a9..19e428b 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py @@ -1,18 +1,18 @@ -from dataclasses import dataclass import numpy as np import logging +from dataclasses import dataclass +from itertools import permutations from vtk import vtkCellSizeFilter -from vtkmodules.vtkCommonCore import vtkIdList from vtkmodules.util.numpy_support import vtk_to_numpy -from vtkmodules.vtkCommonDataModel import ( vtkDataSet, VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, +from vtkmodules.vtkCommonDataModel import ( vtkDataSet, vtkCell, VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM ) -from .vtk_utils import VtkOutput, to_vtk_id_list, write_mesh, read_mesh +from geos.mesh.doctor.checks.vtk_utils import VtkOutput, vtk_iter, to_vtk_id_list, write_mesh, read_mesh @dataclass( frozen=True ) class Options: vtk_output: VtkOutput - cell_name_to_ordering: dict[ str, list[ int ] ] + cell_names_to_reorder: tuple[ str ] volume_to_reorder: str @@ -103,6 +103,90 @@ def is_cell_to_reorder( cell_volume: str, options: Options ) -> bool: return False +def are_face_nodes_counterclockwise( face: any ) -> bool: + """Checks if the nodes of a face are ordered counterclockwise when looking at the plan created by the face. + + Args: + face (any): Face of a vtkCell. + + Returns: + bool: True if counterclockwise, False instead. + """ + face_points = face.GetPoints() + number_points = face_points.GetNumberOfPoints() + if number_points < 3: + err_msg = f"The face has less than 3 nodes which is invalid." + logging.error( err_msg ) + raise ValueError( err_msg ) + # first calculate the normal vector of the face + a = np.array( face_points.GetPoint( 0 ) ) + b = np.array( face_points.GetPoint( 1 ) ) + c = np.array( face_points.GetPoint( 2 ) ) + AB = b - a + AC = c - a + normal = np.cross(AB, AC) + + # then calculate the vector cross products sum of all points of the face with a random point P + # from discussion https://math.stackexchange.com/questions/2152623/determine-the-order-of-a-3d-polygon + P = np.array( [ 0.0, 0.0, 0.0 ] ) # P position does not change the value of sum + all_points = [ np.array( face_points.GetPoint( v ) ) for v in range( number_points ) ] + vector_sum = np.array( [ 0.0, 0.0, 0.0 ] ) + for i in range( number_points ): + PAi = all_points[ i % number_points ] - P + PAiplus1 = all_points[ ( i + 1 ) % number_points ] - P + vector_sum += np.cross( PAi, PAiplus1 ) + vector_sum = vector_sum / 2 # needs to be half + + # if dot product is positive, the nodes are ordered counterclockwise + if np.dot( vector_sum, normal ) > 0.0: + return True + return False + + +def valid_cell_point_ids_ordering( cell: vtkCell ) -> list[ int ]: + """Returns the valid order of point ids of a cell that respect counterclockwise convention of each face. + + Args: + cell (vtkCell): A cell of vtk grid. + + Raises: + ValueError: "No node ids permutation made the face nodes to be ordered counterclockwise." + + Returns: + list[ int ]: [ pt_id0, pt_id2, pt_id1, ..., pt_idN ] + """ + initial_ids_order: list[ int ] = list( vtk_iter( cell.GetPointIds() ) ) + valid_points_id: list[ int ] = initial_ids_order.copy() + + for face_id in range( cell.GetNumberOfFaces() ): + face = cell.GetFace( face_id ) + if are_face_nodes_counterclockwise( face ): + continue # the face nodes respect the convention so continue to next face + + reordered = False + initial_ids: list[ int ] = [ face.GetPointId( i ) for i in range( face.GetNumberOfPoints() ) ] + initial_coords: list[ float ] = [ face.GetPoints().GetPoint( v ) for v in range( face.GetNumberOfPoints() ) ] + for permutation_ids in permutations( initial_ids ): + for i, node_id in enumerate( permutation_ids ): + face.GetPointIds().SetId( i, node_id ) + face.GetPoints().SetPoint( i, initial_coords[ initial_ids.index( node_id ) ] ) + + if are_face_nodes_counterclockwise( face ): # the correct permutation was found + for j, initial_id in enumerate(initial_ids): + if initial_id != permutation_ids[ j ]: + valid_points_id[ initial_ids_order.index( initial_id ) ] = permutation_ids[ j ] + reordered = True + break + + if not reordered: + err_msg = ( f"For face with nodes id '{initial_ids}', that corresponds to coordinates '{initial_coords}'" + + ", no node ids permutation made the face nodes to be ordered counterclockwise." ) + logging.error( err_msg ) + raise ValueError( err_msg ) + + return valid_points_id + + def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple: """From an old mesh, creates a new one where certain cell nodes are reordered to obtain a correct volume. @@ -118,7 +202,7 @@ def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple names_with_totals: dict[ str, int ] = { n: v for n, v in zip( unique_cell_names, total_per_cell_types ) } # sorted dict allow for sorted output of reordering_stats names_with_totals_sorted: dict[ str, int ] = dict( sorted( names_with_totals.items() ) ) - useful_VTK_TYPEs: list[ int ] = [ NAME_TO_VTK_TYPE[ name ] for name in options.cell_name_to_ordering.keys() ] + useful_VTK_TYPEs: list[ int ] = [ NAME_TO_VTK_TYPE[ name ] for name in options.cell_names_to_reorder ] all_cells_volume: np.array = compute_mesh_cells_volume( old_mesh ) # a new mesh with the same data is created from the old mesh new_mesh: vtkDataSet = old_mesh.NewInstance() @@ -132,7 +216,7 @@ def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple "Types non reordered": list( names_with_totals_sorted.keys() ), "Number of cells non reordered": list( names_with_totals_sorted.values() ) } - counter_cells_reordered: dict[ str, int ] = { name: 0 for name in options.cell_name_to_ordering.keys() } + counter_cells_reordered: dict[ str, int ] = { name: 0 for name in options.cell_names_to_reorder } # sorted dict allow for sorted output of reordering_stats ounter_cells_reordered_sorted: dict[ str, int ] = dict( sorted( counter_cells_reordered.items() ) ) # Reordering operations @@ -141,13 +225,8 @@ def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple if vtk_type in useful_VTK_TYPEs: if is_cell_to_reorder( float( all_cells_volume[ cell_id ] ), options ): cell_name: str = VTK_TYPE_TO_NAME[ vtk_type ] - support_point_ids = vtkIdList() - cells.GetCellAtId( cell_id, support_point_ids ) - new_support_point_ids: list[ int ] = list() - node_ordering: list[ int ] = options.cell_name_to_ordering[ cell_name ] - for i in range( len( node_ordering ) ): - new_support_point_ids.append( support_point_ids.GetId( node_ordering[ i ] ) ) - cells.ReplaceCellAtId( cell_id, to_vtk_id_list( new_support_point_ids ) ) + point_ids_ordering: list[ int ] = valid_cell_point_ids_ordering( new_mesh.GetCell( cell_id ) ) + cells.ReplaceCellAtId( cell_id, to_vtk_id_list( point_ids_ordering ) ) ounter_cells_reordered_sorted[ cell_name ] += 1 # Calculation of stats for cell_name_reordered, amount in ounter_cells_reordered_sorted.items(): diff --git a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py index 8f52cbb..5d38404 100644 --- a/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py +++ b/geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py @@ -103,7 +103,7 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid: :return: A unstructured grid. """ if not is_filepath_valid( vtk_input_file ): - err_msg: str = f"Invalid file path. Could not read \"{vtk_input_file}\". Dying..." + err_msg: str = f"Invalid file path. Could not read \"{vtk_input_file}\"." logging.error( err_msg ) raise ValueError( err_msg ) file_extension = os.path.splitext( vtk_input_file )[ -1 ] @@ -119,7 +119,7 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid: if output_mesh: return output_mesh # No reader did work. Dying. - err_msg = f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\". Dying..." + err_msg = f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\"." logging.error( err_msg ) raise ValueError( err_msg ) @@ -160,7 +160,7 @@ def write_mesh( mesh: vtkUnstructuredGrid, vtk_output: VtkOutput ) -> int: success_code = __write_vtu( mesh, vtk_output.output, vtk_output.is_data_mode_binary ) else: # No writer found did work. Dying. - err_msg = f"Could not find the appropriate VTK writer for extension \"{file_extension}\". Dying..." + err_msg = f"Could not find the appropriate VTK writer for extension \"{file_extension}\"." logging.error( err_msg ) raise ValueError( err_msg ) return 0 if success_code else 2 # the Write member function return 1 in case of success, 0 otherwise. diff --git a/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py b/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py index 641a950..04a3bf5 100644 --- a/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py +++ b/geos-mesh/src/geos/mesh/doctor/parsing/fix_elements_orderings_parsing.py @@ -1,41 +1,32 @@ import logging import random -from geos.mesh.doctor.checks.fix_elements_orderings import Options, Result +from geos.mesh.doctor.checks.fix_elements_orderings import Options, Result, NAME_TO_VTK_TYPE from . import vtk_output_parsing, FIX_ELEMENTS_ORDERINGS -__CELL_NAME_WITH_NUMBER_NODES = { - "Tetrahedron": 4, - "Pyramid": 5, - "Wedge": 6, - "Hexahedron": 8, - "Prism5": 10, - "Prism6": 12 -} +__CELL_NAMES = "cell_names" +__CELL_NAMES_CHOICES = list( NAME_TO_VTK_TYPE.keys() ) __VOLUME_TO_REORDER = "volume_to_reorder" -__VOLUME_TO_REORDER_DEFAULT = "all" +__VOLUME_TO_REORDER_DEFAULT = "negative" __VOLUME_TO_REORDER_CHOICES = [ "all", "positive", "negative" ] def fill_subparser( subparsers ) -> None: p = subparsers.add_parser( FIX_ELEMENTS_ORDERINGS, help="Reorders the support nodes for the given cell types." ) - for element_name, size in __CELL_NAME_WITH_NUMBER_NODES.items(): - tmp = list( range( size ) ) - random.Random( 4 ).shuffle( tmp ) - p.add_argument( '--' + element_name, - type=str, - metavar=",".join( map( str, tmp ) ), - default=None, - required=False, - help=f"[list of integers]: node permutation for \"{element_name}\"." ) + p.add_argument( '--' + __CELL_NAMES, + type=str, + metavar=", ".join( map( str, __CELL_NAMES_CHOICES ) ), + default=", ".join( map( str, __CELL_NAMES_CHOICES ) ), + help=f"[list of str]: Cell names that can be reordered in your grid. You can use multiple names." + + "Defaults to all cell names being used." ) p.add_argument( '--' + __VOLUME_TO_REORDER, type=str, default=__VOLUME_TO_REORDER_DEFAULT, - metavar=",".join( map( str, __VOLUME_TO_REORDER_CHOICES ) ), - required=True, + metavar=", ".join( __VOLUME_TO_REORDER_CHOICES ), help="[str]: Select which element volume is invalid and needs reordering." + " 'all' will allow reordering of nodes for every element, regarding of their volume." + - " 'positive' or 'negative' will only reorder the element with the corresponding volume." ) + " 'positive' or 'negative' will only reorder the element with the corresponding volume." + + " Defaults to 'negative'." ) vtk_output_parsing.fill_vtk_output_subparser( p ) @@ -45,22 +36,19 @@ def convert( parsed_options ) -> Options: :param options_str: Parsed cli options. :return: Options instance. """ - cell_name_to_ordering: dict[ str, list[ int ] ] = {} - for element_name, size in __CELL_NAME_WITH_NUMBER_NODES.items(): - raw_mapping = parsed_options[ element_name ] - if raw_mapping: - nodes_ordering = tuple( map( int, raw_mapping.split( "," ) ) ) - if not set( nodes_ordering ) == set( range( size ) ): - err_msg: str = f"Permutation {raw_mapping} for type {element_name} is not valid." - logging.error( err_msg ) - raise ValueError( err_msg ) - cell_name_to_ordering[ element_name ] = nodes_ordering + raw_mapping = parsed_options[ __CELL_NAMES ] + cell_names_to_reorder = tuple( raw_mapping.split( "," ) ) + for cell_name in cell_names_to_reorder: + if cell_name not in __CELL_NAMES_CHOICES: + raise ValueError( f"Please choose names between these options for --{__CELL_NAMES_CHOICES}:" + + f" {__CELL_NAMES_CHOICES}." ) vtk_output = vtk_output_parsing.convert( parsed_options ) volume_to_reorder: str = parsed_options[ __VOLUME_TO_REORDER ] if volume_to_reorder.lower() not in __VOLUME_TO_REORDER_CHOICES: - raise ValueError( f"Please use one of these options for --volume_to_reorder: {__VOLUME_TO_REORDER_CHOICES}." ) + raise ValueError( f"Please use one of these options for --{__VOLUME_TO_REORDER}:" + + f" {__VOLUME_TO_REORDER_CHOICES}." ) return Options( vtk_output=vtk_output, - cell_name_to_ordering=cell_name_to_ordering, + cell_names_to_reorder=cell_names_to_reorder, volume_to_reorder=volume_to_reorder ) diff --git a/geos-mesh/tests/test_fix_elements_orderings.py b/geos-mesh/tests/test_fix_elements_orderings.py index 0f1ec40..060c06f 100644 --- a/geos-mesh/tests/test_fix_elements_orderings.py +++ b/geos-mesh/tests/test_fix_elements_orderings.py @@ -49,14 +49,15 @@ def reorder_cell_nodes( mesh: vtkDataSet, cell_id: int, node_ordering: list[ int Dict used to apply false nodes orderings for test purposes """ to_change_order: dict[ int, list[ int ] ] = { - "Hexahedron": [ 0, 3, 2, 1, 4, 5, 6, 7 ], - "Tetrahedron": [ 0, 2, 1, 3 ], - "Pyramid": [ 0, 3, 2, 1, 4 ], - "Wedge": [ 0, 2, 1, 3, 4, 5 ], - "Prism5": [ 0, 4, 3, 2, 1, 5, 6, 7, 8, 9 ], - "Prism6": [ 0, 1, 4, 2, 3, 5, 11, 10, 9, 8, 7, 6 ] + "Hexahedron": ( 0, 3, 2, 1, 4, 5, 6, 7 ), + "Tetrahedron": ( 0, 2, 1, 3 ), + "Pyramid": ( 0, 3, 2, 1, 4 ), + "Wedge": ( 0, 2, 1, 3, 4, 5 ), + "Prism5": ( 0, 4, 3, 2, 1, 5, 6, 7, 8, 9 ), + "Prism6": ( 0, 1, 4, 2, 3, 5, 11, 10, 9, 8, 7, 6 ) } to_change_order = dict( sorted( to_change_order.items() ) ) +cell_names = list( VTK_TYPE_TO_NAME.values() ) """ 1 Hexahedron: no invalid ordering """ @@ -707,11 +708,11 @@ def test_is_cell_to_reorder( self ): for grid, needs_ordering in grid_needs_ordering.items(): volumes = feo.compute_mesh_cells_volume( grid ) for i in range( len( volumes ) ): - options = opt( out, to_change_order, "negative" ) + options = opt( out, cell_names, "negative" ) assert feo.is_cell_to_reorder( volumes[ i ], options ) == needs_ordering[ i ] - options = opt( out, to_change_order, "positive" ) + options = opt( out, cell_names, "positive" ) assert feo.is_cell_to_reorder( volumes[ i ], options ) != needs_ordering[ i ] - options = opt( out, to_change_order, "all" ) + options = opt( out, cell_names, "all" ) assert feo.is_cell_to_reorder( volumes[ i ], options ) == True def test_get_cell_types_and_number( self ): @@ -738,7 +739,7 @@ def test_get_cell_types_and_number( self ): feo.get_cell_types_and_number( voxels_grid ) def test_reorder_nodes_to_new_mesh( self ): - options = opt( out, to_change_order, "negative" ) + options = opt( out, cell_names, "negative" ) # single element grids except voxels because it is an invalid cell type for GEOS grid_cell_type = { hexahedrons_grid_invalid: VTK_HEXAHEDRON, @@ -779,14 +780,9 @@ def test_fix_elements_orderings_execution( self ): write_mesh( mix_grid_invalid, test_file ) invalidTest = False command = [ - "python", MESH_DOCTOR_FILEPATH, "-v", "-i", test_file.output, "fix_elements_orderings", "--Hexahedron", - str( to_change_order[ "Hexahedron" ] ).replace( "[", "" ).replace( "]", "" ), "--Tetrahedron", - str( to_change_order[ "Tetrahedron" ] ).replace( "[", "" ).replace( "]", "" ), "--Pyramid", - str( to_change_order[ "Pyramid" ] ).replace( "[", "" ).replace( "]", "" ), "--Wedge", - str( to_change_order[ "Wedge" ] ).replace( "[", "" ).replace( "]", "" ), "--Prism5", - str( to_change_order[ "Prism5" ] ).replace( "[", "" ).replace( "]", "" ), "--Prism6", - str( to_change_order[ "Prism6" ] ).replace( "[", "" ).replace( "]", "" ), "--volume_to_reorder", "negative", - "--data-mode", "binary", "--output", filepath_reordered_mesh + "python", MESH_DOCTOR_FILEPATH, "-v", "-i", test_file.output, "fix_elements_orderings", "--cell_names", + ",".join( map( str, cell_names ) ), "--volume_to_reorder", "negative", "--data-mode", "binary", "--output", + filepath_reordered_mesh ] try: result = subprocess.run( command, shell=True, stderr=subprocess.PIPE, universal_newlines=True )