Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix elements orderings to use cells volume #39

Closed
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 25 additions & 20 deletions docs/geos-mesh.rst
Original file line number Diff line number Diff line change
Expand Up @@ -123,26 +123,35 @@ 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] [--Tetrahedron 2,0,3,1] [--Pyramid 3,4,0,2,1] [--Wedge 3,5,4,0,2,1]
[--Hexahedron 1,6,5,4,7,0,2,3] [--Prism5 8,2,0,7,6,9,5,1,4,3] [--Prism6 11,2,8,10,5,0,9,7,6,1,4,3]
--volume_to_reorder all,positive,negative --output OUTPUT [--data-mode binary, ascii]

options:
-h, --help show this help message and exit
-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".
[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".
[list of integers]: node permutation for "Prism5".
--Prism6 11,2,8,10,5,0,9,7,6,1,4,3
[list of integers]: node permutation for "Prism6".
--Pyramid 3,4,0,2,1 [list of integers]: node permutation for "Pyramid".
--Tetrahedron 2,0,3,1 [list of integers]: node permutation for "Tetrahedron".
--Voxel 1,6,5,4,7,0,2,3 [list of integers]: node permutation for "Voxel".
--Wedge 3,5,4,0,2,1 [list of integers]: node permutation for "Wedge".
--output OUTPUT [string]: The vtk output file destination.
--data-mode binary, ascii
[string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary.
[list of integers]: node permutation for "Prism6".
--volume_to_reorder all,positive,negative
[str]: Select which element volume is invalid and needs reordering. 'all' will allow reordering of nodes for every element, regarding of their volume.
'positive' or 'negative' will only reorder the element with the corresponding volume.
--output OUTPUT [string]: The vtk output file destination.
--data-mode binary, ascii [string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary.

For example, assume that you have a mesh that contains 3 different cell types like Hexahedrons, Pyramids and Tetrahedrons.
After checking ``element_volumes``, you found that all your Hexahedrons and half of your Tetrahedrons have a negative volume.
To correct that, assuming you know the correct node ordering to correct these volumes, you can use this command:

.. code-block::

$ python src/geos/mesh/doctor/mesh_doctor.py -i /path/to/your/mesh/file fix_elements_orderings --Hexahedron 1,6,5,4,7,0,2,3
--Tetrahedron 2,0,3,1 --volume_to_reorder negative --output /new/path/for/your/new/mesh/reordered/file --data-mode binary

``generate_cube``
"""""""""""""""""
Expand Down Expand Up @@ -340,8 +349,4 @@ API
^^^

.. automodule:: geos.mesh.conversion.abaqus_converter
:members:




:members:
200 changes: 154 additions & 46 deletions geos-mesh/src/geos/mesh/doctor/checks/fix_elements_orderings.py
Original file line number Diff line number Diff line change
@@ -1,65 +1,173 @@
from dataclasses import dataclass
import numpy as np
import logging
from typing import (
List,
Dict,
Set,
FrozenSet,
)

from vtkmodules.vtkCommonCore import (
vtkIdList, )

from . import vtk_utils
from .vtk_utils import (
to_vtk_id_list,
VtkOutput,
)
from vtk import vtkCellSizeFilter
from vtkmodules.vtkCommonCore import vtkIdList
from vtkmodules.util.numpy_support import vtk_to_numpy
from vtkmodules.vtkCommonDataModel import ( vtkDataSet, VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE,
VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM )
from .vtk_utils import VtkOutput, to_vtk_id_list, write_mesh, read_mesh


@dataclass( frozen=True )
class Options:
vtk_output: VtkOutput
cell_type_to_ordering: Dict[ int, List[ int ] ]
cell_name_to_ordering: dict[ str, list[ int ] ]
volume_to_reorder: str


@dataclass( frozen=True )
class Result:
output: str
unchanged_cell_types: FrozenSet[ int ]
reordering_stats: dict[ str, list[ int ] ]


GEOS_ACCEPTED_TYPES = [ VTK_HEXAHEDRON, VTK_TETRA, VTK_PYRAMID, VTK_WEDGE, VTK_PENTAGONAL_PRISM, VTK_HEXAGONAL_PRISM ]
# the number of different nodes that needs to be entered in parsing when dealing with a specific vtk element
NAME_TO_VTK_TYPE = {
"Hexahedron": VTK_HEXAHEDRON,
"Tetrahedron": VTK_TETRA,
"Pyramid": VTK_PYRAMID,
"Wedge": VTK_WEDGE,
"Prism5": VTK_PENTAGONAL_PRISM,
"Prism6": VTK_HEXAGONAL_PRISM
}
VTK_TYPE_TO_NAME = { vtk_type: name for name, vtk_type in NAME_TO_VTK_TYPE.items() }


# Knowing the calculation of cell volumes in vtk was discussed there: https://github.com/GEOS-DEV/GEOS/issues/2253
# Here, we do not need to have the exact volumes matching the simulation softwares results
# because we are mostly interested in knowing the sign of the volume for the rest of the workflow.
# Therefore, there is no need to use vtkMeshQuality for specific cell types when vtkCellSizeFilter works with all types.
def compute_mesh_cells_volume( mesh: vtkDataSet ) -> np.array:
"""Generates a volume array that was calculated on all cells of a mesh.

Args:
mesh (vtkDataSet): A vtk grid.

Returns:
np.array: Volume for every cell of a mesh.
"""
cell_size_filter = vtkCellSizeFilter()
cell_size_filter.SetInputData( mesh )
cell_size_filter.SetComputeVolume( True )
cell_size_filter.Update()
return vtk_to_numpy( cell_size_filter.GetOutput().GetCellData().GetArray( "Volume" ) )


def get_cell_types_and_number( mesh: vtkDataSet ) -> tuple[ list[ int ] ]:
"""Gets the cell type for every cell of a mesh and the amount for each cell type.

Args:
mesh (vtkDataSet): A vtk grid.

Raises:
ValueError: "Invalid type '{cell_type}' for GEOS is in the mesh. Dying ..."

Returns:
tuple[ list[ int ] ]: ( unique_cell_types, total_per_cell_types )
"""
number_cells: int = mesh.GetNumberOfCells()
all_cells_type: np.array = np.ones( number_cells, dtype=int )
for cell_id in range( number_cells ):
all_cells_type[ cell_id ] = mesh.GetCellType( cell_id )
unique_cell_types, total_per_cell_types = np.unique( all_cells_type, return_counts=True )
unique_cell_types, total_per_cell_types = unique_cell_types.tolist(), total_per_cell_types.tolist()
for cell_type in unique_cell_types:
if cell_type not in GEOS_ACCEPTED_TYPES:
err_msg: str = f"Invalid type '{cell_type}' for GEOS is in the mesh. Dying ..."
logging.error( err_msg )
raise ValueError( err_msg )
return ( unique_cell_types, total_per_cell_types )


def is_cell_to_reorder( cell_volume: str, options: Options ) -> bool:
"""Check if the volume of vtkCell qualifies the cell to be reordered.

Args:
cell_volume (float): The volume of a vtkCell.
options (Options): Options defined by the user.

Returns:
bool: True if cell needs to be reordered
"""
if options.volume_to_reorder == "all":
return True
if cell_volume == 0.0:
return True
sign_of_cell_volume: int = int( cell_volume / abs( cell_volume ) )
if options.volume_to_reorder == "positive" and sign_of_cell_volume == 1:
return True
elif options.volume_to_reorder == "negative" and sign_of_cell_volume == -1:
return True
return False


def reorder_nodes_to_new_mesh( old_mesh: vtkDataSet, options: Options ) -> tuple:
"""From an old mesh, creates a new one where certain cell nodes are reordered to obtain a correct volume.

Args:
old_mesh (vtkDataSet): A vtk grid needing nodes to be reordered.
options (Options): Options defined by the user.

Returns:
tuple: ( vtkDataSet, reordering_stats )
"""
unique_cell_types, total_per_cell_types = get_cell_types_and_number( old_mesh )
unique_cell_names: list[ str ] = [ VTK_TYPE_TO_NAME[ vtk_type ] for vtk_type in unique_cell_types ]
names_with_totals: dict[ str, int ] = { n: v for n, v in zip( unique_cell_names, total_per_cell_types ) }
# sorted dict allow for sorted output of reordering_stats
names_with_totals_sorted: dict[ str, int ] = dict( sorted( names_with_totals.items() ) )
useful_VTK_TYPEs: list[ int ] = [ NAME_TO_VTK_TYPE[ name ] for name in options.cell_name_to_ordering.keys() ]
all_cells_volume: np.array = compute_mesh_cells_volume( old_mesh )
# a new mesh with the same data is created from the old mesh
new_mesh: vtkDataSet = old_mesh.NewInstance()
new_mesh.CopyStructure( old_mesh )
new_mesh.CopyAttributes( old_mesh )
cells = new_mesh.GetCells()
# Statistics on how many cells have or have not been reordered
reordering_stats: dict[ str, list[ any ] ] = {
"Types reordered": list(),
"Number of cells reordered": list(),
"Types non reordered": list( names_with_totals_sorted.keys() ),
"Number of cells non reordered": list( names_with_totals_sorted.values() )
}
counter_cells_reordered: dict[ str, int ] = { name: 0 for name in options.cell_name_to_ordering.keys() }
# sorted dict allow for sorted output of reordering_stats
ounter_cells_reordered_sorted: dict[ str, int ] = dict( sorted( counter_cells_reordered.items() ) )
# Reordering operations
for cell_id in range( new_mesh.GetNumberOfCells() ):
vtk_type: int = new_mesh.GetCellType( cell_id )
if vtk_type in useful_VTK_TYPEs:
if is_cell_to_reorder( float( all_cells_volume[ cell_id ] ), options ):
cell_name: str = VTK_TYPE_TO_NAME[ vtk_type ]
support_point_ids = vtkIdList()
cells.GetCellAtId( cell_id, support_point_ids )
new_support_point_ids: list[ int ] = list()
node_ordering: list[ int ] = options.cell_name_to_ordering[ cell_name ]
for i in range( len( node_ordering ) ):
new_support_point_ids.append( support_point_ids.GetId( node_ordering[ i ] ) )
cells.ReplaceCellAtId( cell_id, to_vtk_id_list( new_support_point_ids ) )
ounter_cells_reordered_sorted[ cell_name ] += 1
# Calculation of stats
for cell_name_reordered, amount in ounter_cells_reordered_sorted.items():
if amount > 0:
reordering_stats[ "Types reordered" ].append( cell_name_reordered )
reordering_stats[ "Number of cells reordered" ].append( amount )
index_non_reordered: int = reordering_stats[ "Types non reordered" ].index( cell_name_reordered )
reordering_stats[ "Number of cells non reordered" ][ index_non_reordered ] -= amount
if reordering_stats[ "Number of cells non reordered" ][ index_non_reordered ] == 0:
reordering_stats[ "Types non reordered" ].pop( index_non_reordered )
reordering_stats[ "Number of cells non reordered" ].pop( index_non_reordered )
return ( new_mesh, reordering_stats )


def __check( mesh, options: Options ) -> Result:
# The vtk cell type is an int and will be the key of the following mapping,
# that will point to the relevant permutation.
cell_type_to_ordering: Dict[ int, List[ int ] ] = options.cell_type_to_ordering
unchanged_cell_types: Set[ int ] = set() # For logging purpose

# Preparing the output mesh by first keeping the same instance type.
output_mesh = mesh.NewInstance()
output_mesh.CopyStructure( mesh )
output_mesh.CopyAttributes( mesh )

# `output_mesh` now contains a full copy of the input mesh.
# We'll now modify the support nodes orderings in place if needed.
cells = output_mesh.GetCells()
for cell_idx in range( output_mesh.GetNumberOfCells() ):
cell_type: int = output_mesh.GetCell( cell_idx ).GetCellType()
new_ordering = cell_type_to_ordering.get( cell_type )
if new_ordering:
support_point_ids = vtkIdList()
cells.GetCellAtId( cell_idx, support_point_ids )
new_support_point_ids = []
for i, v in enumerate( new_ordering ):
new_support_point_ids.append( support_point_ids.GetId( new_ordering[ i ] ) )
cells.ReplaceCellAtId( cell_idx, to_vtk_id_list( new_support_point_ids ) )
else:
unchanged_cell_types.add( cell_type )
is_written_error = vtk_utils.write_mesh( output_mesh, options.vtk_output )
return Result( output=options.vtk_output.output if not is_written_error else "",
unchanged_cell_types=frozenset( unchanged_cell_types ) )
output_mesh, reordering_stats = reorder_nodes_to_new_mesh( mesh, options )
write_mesh( output_mesh, options.vtk_output )
return Result( output=options.vtk_output.output, reordering_stats=reordering_stats )


def check( vtk_input_file: str, options: Options ) -> Result:
mesh = vtk_utils.read_mesh( vtk_input_file )
mesh = read_mesh( vtk_input_file )
return __check( mesh, options )
54 changes: 45 additions & 9 deletions geos-mesh/src/geos/mesh/doctor/checks/vtk_utils.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
from dataclasses import dataclass
import os.path
from os import path, access, R_OK, W_OK
import logging
import sys
from sys import exit
from typing import (
Any,
Iterator,
Expand Down Expand Up @@ -50,6 +50,39 @@ def vtk_iter( l ) -> Iterator[ Any ]:
yield l.GetCellType( i )


# Inspired from https://stackoverflow.com/a/78870363
def is_filepath_valid( filepath: str, is_input=True ) -> bool:
"""Checks if a filepath can be used to read or to create a new file.

Args:
filepath (str): A filepath.
is_input (bool, optional): If the filepath is used to read a file, use True.
If the filepath is used to create a new file, use False. Defaults to True.

Returns:
bool: False if invalid, True instead.
"""
dirname = path.dirname( filepath )
if not path.isdir( dirname ):
logging.error( f"The directory '{dirname}' specified does not exist." )
return False
if is_input:
if not access( dirname, R_OK ):
logging.error( f"You do not have rights to read in directory '{dirname}' specified for the file." )
return False
if not path.exists( filepath ):
logging.error( f"The file specified does not exist." )
return False
else:
if not access( dirname, W_OK ):
logging.error( f"You do not have rights to write in directory '{dirname}' specified for the file." )
return False
if path.exists( filepath ):
logging.error( f"A file with the same name already exists. No over-writing possible." )
return False
return True


def __read_vtk( vtk_input_file: str ) -> Optional[ vtkUnstructuredGrid ]:
reader = vtkUnstructuredGridReader()
logging.info( f"Testing file format \"{vtk_input_file}\" using legacy format reader..." )
Expand Down Expand Up @@ -83,7 +116,10 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid:
If first guess does not work, eventually all the others reader available will be tested.
:return: A unstructured grid.
"""
file_extension = os.path.splitext( vtk_input_file )[ -1 ]
if not is_filepath_valid( vtk_input_file, True ):
logging.error( f"Invalid filepath to read. Dying ..." )
exit( 1 )
file_extension = path.splitext( vtk_input_file )[ -1 ]
extension_to_reader = { ".vtk": __read_vtk, ".vtu": __read_vtu }
# Testing first the reader that should match
if file_extension in extension_to_reader:
Expand All @@ -97,7 +133,7 @@ def read_mesh( vtk_input_file: str ) -> vtkUnstructuredGrid:
return output_mesh
# No reader did work. Dying.
logging.critical( f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\". Dying..." )
sys.exit( 1 )
exit( 1 )


def __write_vtk( mesh: vtkUnstructuredGrid, output: str ) -> int:
Expand Down Expand Up @@ -125,16 +161,16 @@ def write_mesh( mesh: vtkUnstructuredGrid, vtk_output: VtkOutput ) -> int:
:param vtk_output: Where to write. The file extension will be used to select the VTK file format.
:return: 0 in case of success.
"""
if os.path.exists( vtk_output.output ):
logging.error( f"File \"{vtk_output.output}\" already exists, nothing done." )
return 1
file_extension = os.path.splitext( vtk_output.output )[ -1 ]
if not is_filepath_valid( vtk_output.output, False ):
logging.error( f"Invalid filepath to write. Dying ..." )
exit( 1 )
file_extension = path.splitext( vtk_output.output )[ -1 ]
if file_extension == ".vtk":
success_code = __write_vtk( mesh, vtk_output.output )
elif file_extension == ".vtu":
success_code = __write_vtu( mesh, vtk_output.output, vtk_output.is_data_mode_binary )
else:
# No writer found did work. Dying.
logging.critical( f"Could not find the appropriate VTK writer for extension \"{file_extension}\". Dying..." )
sys.exit( 1 )
exit( 1 )
return 0 if success_code else 2 # the Write member function return 1 in case of success, 0 otherwise.
12 changes: 6 additions & 6 deletions geos-mesh/src/geos/mesh/doctor/mesh_doctor.py
Original file line number Diff line number Diff line change
@@ -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():
Expand Down
Loading
Loading