diff --git a/docs/geos-mesh.rst b/docs/geos-mesh.rst index 061f596d..2221cc94 100644 --- a/docs/geos-mesh.rst +++ b/docs/geos-mesh.rst @@ -9,6 +9,8 @@ GEOS Mesh tools ./geos_mesh_docs/doctor + ./geos_mesh_docs/doctor_filters + ./geos_mesh_docs/converter ./geos_mesh_docs/io diff --git a/docs/geos_mesh_docs/doctor.rst b/docs/geos_mesh_docs/doctor.rst index 0e66d84f..6e07353a 100644 --- a/docs/geos_mesh_docs/doctor.rst +++ b/docs/geos_mesh_docs/doctor.rst @@ -1,9 +1,9 @@ Mesh Doctor ---------------- +----------- -``mesh-doctor`` is a ``python`` executable that can be used through the command line to perform various checks, validations, and tiny fixes to the ``vtk`` mesh that are meant to be used in ``geos``. -``mesh-doctor`` is organized as a collection of modules with their dedicated sets of options. -The current page will introduce those modules, but the details and all the arguments can be retrieved by using the ``--help`` option for each module. +| ``mesh-doctor`` is a ``python`` executable that can be used through the command line to perform various checks, validations, and tiny fixes to the ``vtkUnstructuredGrid`` mesh that are meant to be used in ``geos``. + ``mesh-doctor`` is organized as a collection of modules with their dedicated sets of options. +| The current page will introduce those modules, but the details and all the arguments can be retrieved by using the ``--help`` option for each module. Prerequisites ^^^^^^^^^^^^^ @@ -310,8 +310,55 @@ It will also verify that the ``VTK_POLYHEDRON`` cells can effectively get conver .. code-block:: $ mesh-doctor supported_elements --help - usage: mesh_doctor.py supported_elements [-h] [--chunck_size 1] [--nproc 8] + usage: mesh_doctor.py supported_elements [-h] [--chunk_size 1] [--nproc 8] options: -h, --help show this help message and exit - --chunck_size 1 [int]: Defaults chunk size for parallel processing to 1 - --nproc 8 [int]: Number of threads used for parallel processing. Defaults to your CPU count 8. \ No newline at end of file + --chunk_size 1 [int]: Defaults chunk size for parallel processing to 1 + --nproc 8 [int]: Number of threads used for parallel processing. Defaults to your CPU count 8. + + +Why only use vtkUnstructuredGrid? +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The mesh doctor is designed specifically for unstructured meshes used in GEOS. +| All input files are expected to be ``.vtu`` (VTK Unstructured Grid) format. +| What about other formats? + +VTK Hierarchy +""""""""""""" + +Supposedly, other grid types that are part of the following VTK hierarchy could be used:: + + vtkDataObject + └── vtkDataSet + └── vtkCartesianGrid + └── vtkRectilinearGrid + └── vtkImageData + └── vtkStructuredPoints + └── vtkUniformGrid + └── vtkPointSet + └── vtkExplicitStructuredGrid + └── vtkPolyData + └── vtkStructuredGrid + └── vtkUnstructuredGrid + +And when looking at specific methods used in mesh-doctor, it could suggest that other formats could be used: + +* Points access: mesh.GetPoints() - Available in all vtkPointSet subclasses ✓ +* Cell iteration: mesh.GetNumberOfCells(), mesh.GetCell() - Available in all vtkDataSet subclasses ✓ +* Cell types: mesh.GetCellType() - Available in all vtkDataSet subclasses ✓ +* Cell/Point data: mesh.GetCellData(), mesh.GetPointData() - Available in all vtkDataSet subclasses ✓ + +VTK Filter Compatibility +"""""""""""""""""""""""" + +| vtkCellSizeFilter, vtkMeshQuality, and other VTK filters used in the actions expect vtkDataSet or its subclasses + vtkUnstructuredGrid is compatible with all VTK filters used. +| vtkPolyData has a different data structure, not suitable for 3D volumetric meshes. + +Specific Operations Require vtkUnstructuredGrid +""""""""""""""""""""""""""""""""""""""""""""""" + +* GetCellNeighbors() - Only available in vtkUnstructuredGrid +* GetFaceStream() - Only available in vtkUnstructuredGrid (for polyhedron support) +* GetDistinctCellTypesArray() - Only available in vtkUnstructuredGrid \ No newline at end of file diff --git a/docs/geos_mesh_docs/doctor_filters.rst b/docs/geos_mesh_docs/doctor_filters.rst new file mode 100644 index 00000000..119ee508 --- /dev/null +++ b/docs/geos_mesh_docs/doctor_filters.rst @@ -0,0 +1,19 @@ +VTK Filters +=========== + +In addition to the command-line interface, mesh-doctor functionality is also available as VTK filters for programmatic use in Python. These filters provide the same mesh validation and processing capabilities but can be integrated directly into Python workflows and visualization pipelines. + +The VTK filters offer several advantages: + +* **Integration**: Easy integration with existing VTK/ParaView workflows +* **Scripting**: Programmatic access for batch processing and automation +* **Visualization**: Direct output suitable for ParaView visualization +* **Flexibility**: Configurable parameters and output options +* **Chaining**: Ability to combine multiple filters in processing pipelines + +For detailed documentation on available filters, their parameters, and usage examples, see: + +.. toctree:: + :maxdepth: 1 + + filters/index \ No newline at end of file diff --git a/docs/geos_mesh_docs/filters/AllChecks.rst b/docs/geos_mesh_docs/filters/AllChecks.rst new file mode 100644 index 00000000..d7dd6fd5 --- /dev/null +++ b/docs/geos_mesh_docs/filters/AllChecks.rst @@ -0,0 +1,88 @@ +AllChecks Filter +================ + +.. autoclass:: geos.mesh.doctor.filters.Checks.AllChecks + :members: + :undoc-members: + :show-inheritance: + +Overview +-------- + +The AllChecks filter performs comprehensive mesh validation by running all available quality checks on a vtkUnstructuredGrid. This filter is part of the mesh doctor toolkit and provides detailed analysis of mesh quality, topology, and geometric integrity. + +Features +-------- + +* Comprehensive mesh validation with all available quality checks +* Configurable check parameters for customized validation +* Detailed reporting of found issues +* Integration with mesh doctor parsing system +* Support for both individual check parameter customization and global parameter setting + +Usage Example +------------- + +.. code-block:: python + + from geos.mesh.doctor.filters.Checks import AllChecks + + # Instantiate the filter for all available checks + allChecksFilter = AllChecks() + + # Set input mesh + allChecksFilter.SetInputData(mesh) + + # Optionally customize check parameters + allChecksFilter.setCheckParameter("collocated_nodes", "tolerance", 1e-6) + allChecksFilter.setGlobalParameter("tolerance", 1e-6) # applies to all checks with tolerance parameter + + # Execute the checks and get output + output_mesh = allChecksFilter.getGrid() + + # Get check results + check_results = allChecksFilter.getCheckResults() + + # Write the output mesh + allChecksFilter.writeGrid("output/mesh_with_check_results.vtu") + +Parameters +---------- + +The AllChecks filter supports parameter customization for individual checks: + +* **setCheckParameter(check_name, parameter_name, value)**: Set specific parameter for a named check +* **setGlobalParameter(parameter_name, value)**: Apply parameter to all checks that support it + +Common parameters include: + +* **tolerance**: Distance tolerance for geometric checks (e.g., collocated nodes, non-conformal interfaces) +* **angle_tolerance**: Angular tolerance for orientation checks +* **min_volume**: Minimum acceptable element volume + +Available Checks +---------------- + +The AllChecks filter includes all checks available in the mesh doctor system: + +* Collocated nodes detection +* Element volume validation +* Self-intersecting elements detection +* Non-conformal interface detection +* Supported element type validation +* And many more quality checks + +Output +------ + +* **Input**: vtkUnstructuredGrid +* **Output**: vtkUnstructuredGrid (copy of input with potential additional arrays marking issues) +* **Check Results**: Detailed dictionary with results from all performed checks + +See Also +-------- + +* :doc:`MainChecks ` - Subset of most important checks +* :doc:`CollocatedNodes ` - Individual collocated nodes check +* :doc:`ElementVolumes ` - Individual element volume check +* :doc:`SelfIntersectingElements ` - Individual self-intersection check diff --git a/docs/geos_mesh_docs/filters/CollocatedNodes.rst b/docs/geos_mesh_docs/filters/CollocatedNodes.rst new file mode 100644 index 00000000..96edae1a --- /dev/null +++ b/docs/geos_mesh_docs/filters/CollocatedNodes.rst @@ -0,0 +1,131 @@ +CollocatedNodes Filter +====================== + +.. automodule:: geos.mesh.doctor.filters.CollocatedNodes + :members: + :undoc-members: + :show-inheritance: + +Overview +-------- + +The CollocatedNodes filter identifies and handles duplicated/collocated nodes in a vtkUnstructuredGrid. Collocated nodes are nodes that are positioned very close to each other (within a specified tolerance), which can indicate mesh quality issues or modeling problems. + +Features +-------- + +* Detection of collocated/duplicated nodes within specified tolerance +* Identification of elements with wrong support (nodes appearing multiple times) +* Optional marking of problematic elements in output mesh +* Configurable tolerance for distance-based node comparison +* Detailed reporting of found collocated node groups + +Usage Example +------------- + +.. code-block:: python + + from geos.mesh.doctor.filters.CollocatedNodes import CollocatedNodes + + # Instantiate the filter + collocatedNodesFilter = CollocatedNodes() + + # Set the tolerance for detecting collocated nodes + collocatedNodesFilter.setTolerance(1e-6) + + # Optionally enable painting of wrong support elements + collocatedNodesFilter.setPaintWrongSupportElements(1) # 1 to enable, 0 to disable + + # Set input mesh + collocatedNodesFilter.SetInputData(mesh) + + # Execute the filter and get output + output_mesh = collocatedNodesFilter.getGrid() + + # Get results + collocated_buckets = collocatedNodesFilter.getCollocatedNodeBuckets() # list of tuples with collocated node indices + wrong_support_elements = collocatedNodesFilter.getWrongSupportElements() # list of problematic element indices + + # Write the output mesh + collocatedNodesFilter.writeGrid("output/mesh_with_collocated_info.vtu") + +Parameters +---------- + +setTolerance(tolerance) + Set the distance tolerance for determining if two nodes are collocated. + + * **tolerance** (float): Distance threshold below which nodes are considered collocated + * **Default**: 0.0 + +setPaintWrongSupportElements(choice) + Enable/disable creation of array marking elements with wrong support nodes. + + * **choice** (int): 1 to enable marking, 0 to disable + * **Default**: 0 + +Results Access +-------------- + +getCollocatedNodeBuckets() + Returns groups of collocated node indices. + + * **Returns**: list[tuple[int]] - Each tuple contains indices of nodes that are collocated + +getWrongSupportElements() + Returns element indices that have support nodes appearing more than once. + + * **Returns**: list[int] - Element indices with problematic support nodes + +Understanding the Results +------------------------- + +**Collocated Node Buckets** + +Each bucket is a tuple containing node indices that are within the specified tolerance of each other: + +.. code-block:: python + + # Example result: [(0, 15, 23), (7, 42), (100, 101, 102, 103)] + # This means: + # - Nodes 0, 15, and 23 are collocated + # - Nodes 7 and 42 are collocated + # - Nodes 100, 101, 102, and 103 are collocated + +**Wrong Support Elements** + +Elements where the same node appears multiple times in the element's connectivity. This usually indicates: + +* Degenerate elements +* Mesh generation errors +* Topology problems + +Common Use Cases +---------------- + +* **Mesh Quality Assessment**: Identify potential mesh issues before simulation +* **Mesh Preprocessing**: Clean up meshes by detecting node duplicates +* **Debugging**: Understand why meshes might have connectivity problems +* **Validation**: Ensure mesh meets quality standards for specific applications + +Output +------ + +* **Input**: vtkUnstructuredGrid +* **Output**: vtkUnstructuredGrid with optional arrays marking problematic elements +* **Additional Data**: When painting is enabled, adds "WrongSupportElements" array to cell data + +Best Practices +-------------- + +* Set tolerance based on mesh scale and precision requirements +* Use smaller tolerances for high-precision meshes +* Enable painting to visualize problematic areas in ParaView +* Check both collocated buckets and wrong support elements for comprehensive analysis + +See Also +-------- + +* :doc:`AllChecks ` - Includes collocated nodes check among others +* :doc:`MainChecks ` - Includes collocated nodes check in main set +* :doc:`SelfIntersectingElements ` - Related geometric validation diff --git a/docs/geos_mesh_docs/filters/ElementVolumes.rst b/docs/geos_mesh_docs/filters/ElementVolumes.rst new file mode 100644 index 00000000..2a9910a6 --- /dev/null +++ b/docs/geos_mesh_docs/filters/ElementVolumes.rst @@ -0,0 +1,160 @@ +ElementVolumes Filter +===================== + +.. automodule:: geos.mesh.doctor.filters.ElementVolumes + :members: + :undoc-members: + :show-inheritance: + +Overview +-------- + +The ElementVolumes filter calculates the volumes of all elements in a vtkUnstructuredGrid and can identify elements with negative or zero volumes. Such elements typically indicate serious mesh quality issues including inverted elements, degenerate cells, or incorrect node ordering. + +Features +-------- + +* Volume calculation for all element types +* Detection of negative and zero volume elements +* Detailed reporting of problematic elements with their volumes +* Integration with VTK's cell size computation +* Optional filtering to return only problematic elements + +Usage Example +------------- + +.. code-block:: python + + from geos.mesh.doctor.filters.ElementVolumes import ElementVolumes + + # Instantiate the filter + elementVolumesFilter = ElementVolumes() + + # Optionally enable detection of negative/zero volume elements + elementVolumesFilter.setReturnNegativeZeroVolumes(True) + + # Set input mesh + elementVolumesFilter.SetInputData(mesh) + + # Execute the filter and get output + output_mesh = elementVolumesFilter.getGrid() + + # Get problematic elements (if enabled) + if elementVolumesFilter.m_returnNegativeZeroVolumes: + negative_zero_volumes = elementVolumesFilter.getNegativeZeroVolumes() + # Returns numpy array with shape (n, 2) where first column is element index, second is volume + + # Write the output mesh with volume information + elementVolumesFilter.writeGrid("output/mesh_with_volumes.vtu") + +Parameters +---------- + +setReturnNegativeZeroVolumes(returnNegativeZeroVolumes) + Enable/disable specific detection and return of elements with negative or zero volumes. + + * **returnNegativeZeroVolumes** (bool): True to enable detection, False to disable + * **Default**: False + +Results Access +-------------- + +getNegativeZeroVolumes() + Returns detailed information about elements with negative or zero volumes. + + * **Returns**: numpy.ndarray with shape (n, 2) + + * Column 0: Element indices with problematic volumes + * Column 1: Corresponding volume values (≤ 0) + + * **Note**: Only available when returnNegativeZeroVolumes is enabled + +Understanding Volume Issues +--------------------------- + +**Negative Volumes** + +Indicate elements with inverted geometry: + +* **Tetrahedra**: Nodes ordered incorrectly (clockwise instead of counter-clockwise) +* **Hexahedra**: Faces oriented inward instead of outward +* **Other elements**: Similar orientation/ordering issues + +**Zero Volumes** + +Indicate degenerate elements: + +* **Collapsed elements**: All nodes lie in the same plane (3D) or line (2D) +* **Duplicate nodes**: Multiple nodes at the same location within the element +* **Extreme aspect ratios**: Elements stretched to near-zero thickness + +**Impact on Simulations** + +Elements with non-positive volumes can cause: + +* Numerical instabilities +* Convergence problems +* Incorrect physical results +* Solver failures + +Common Fixes +------------ + +**For Negative Volumes:** + +1. **Reorder nodes**: Correct the connectivity order +2. **Flip elements**: Use mesh repair tools +3. **Regenerate mesh**: If issues are widespread + +**For Zero Volumes:** + +1. **Remove degenerate elements**: Delete problematic cells +2. **Merge collocated nodes**: Fix duplicate node issues +3. **Improve mesh quality**: Regenerate with better settings + +Output +------ + +* **Input**: vtkUnstructuredGrid +* **Output**: vtkUnstructuredGrid with volume information added as cell data +* **Volume Array**: "Volume" array added to cell data containing computed volumes +* **Additional Data**: When returnNegativeZeroVolumes is enabled, provides detailed problematic element information + +Integration with Other Filters +------------------------------ + +The ElementVolumes filter works well in combination with: + +* **CollocatedNodes**: Fix node duplication that can cause zero volumes +* **SelfIntersectingElements**: Detect other geometric problems +* **AllChecks/MainChecks**: Comprehensive mesh validation including volume checks + +Example Workflow +---------------- + +.. code-block:: python + + # Complete volume analysis workflow + volumeFilter = ElementVolumes() + volumeFilter.setReturnNegativeZeroVolumes(True) + volumeFilter.SetInputData(mesh) + + output_mesh = volumeFilter.getGrid() + + # Analyze results + problematic = volumeFilter.getNegativeZeroVolumes() + + if len(problematic) > 0: + print(f"Found {len(problematic)} elements with non-positive volumes:") + for idx, volume in problematic: + print(f" Element {idx}: volume = {volume}") + else: + print("All elements have positive volumes - mesh is good!") + +See Also +-------- + +* :doc:`AllChecks ` - Includes element volume check among others +* :doc:`MainChecks ` - Includes element volume check in main set +* :doc:`CollocatedNodes ` - Can help fix zero volume issues +* :doc:`SelfIntersectingElements ` - Related geometric validation diff --git a/docs/geos_mesh_docs/filters/GenerateFractures.rst b/docs/geos_mesh_docs/filters/GenerateFractures.rst new file mode 100644 index 00000000..78c6feab --- /dev/null +++ b/docs/geos_mesh_docs/filters/GenerateFractures.rst @@ -0,0 +1,220 @@ +GenerateFractures Filter +======================== + +.. automodule:: geos.mesh.doctor.filters.GenerateFractures + :members: + :undoc-members: + :show-inheritance: + +Overview +-------- + +The GenerateFractures filter splits a vtkUnstructuredGrid along non-embedded fractures. When a fracture plane is defined between two cells, the nodes of the shared face are duplicated to create a discontinuity. The filter generates both the split main mesh and separate fracture meshes. + +Features +-------- + +* Mesh splitting along fracture planes with node duplication +* Multiple fracture policy support (internal vs boundary fractures) +* Configurable fracture identification via field values +* Generation of separate fracture mesh outputs +* Flexible output data modes (ASCII/binary) +* Automatic fracture mesh file management + +Usage Example +------------- + +.. code-block:: python + + from geos.mesh.doctor.filters.GenerateFractures import GenerateFractures + + # Instantiate the filter + generateFracturesFilter = GenerateFractures() + + # Set the field name that defines fracture regions + generateFracturesFilter.setFieldName("fracture_field") + + # Set the field values that identify fracture boundaries + generateFracturesFilter.setFieldValues("1,2") # comma-separated values + + # Set fracture policy (0 for internal fractures, 1 for boundary fractures) + generateFracturesFilter.setPolicy(1) + + # Set output directory for fracture meshes + generateFracturesFilter.setFracturesOutputDirectory("./fractures/") + + # Optionally set data mode (0 for ASCII, 1 for binary) + generateFracturesFilter.setOutputDataMode(1) + generateFracturesFilter.setFracturesDataMode(1) + + # Set input mesh + generateFracturesFilter.SetInputData(mesh) + + # Execute the filter + generateFracturesFilter.Update() + + # Get the split mesh and fracture meshes + split_mesh, fracture_meshes = generateFracturesFilter.getAllGrids() + + # Write all meshes + generateFracturesFilter.writeMeshes("output/split_mesh.vtu", is_data_mode_binary=True) + +Parameters +---------- + +setFieldName(field_name) + Set the name of the cell data field that defines fracture regions. + + * **field_name** (str): Name of the field in cell data + +setFieldValues(field_values) + Set the field values that identify fracture boundaries. + + * **field_values** (str): Comma-separated list of values (e.g., "1,2,3") + +setPolicy(choice) + Set the fracture generation policy. + + * **choice** (int): + + * 0: Internal fractures policy + * 1: Boundary fractures policy + +setFracturesOutputDirectory(directory) + Set the output directory for individual fracture mesh files. + + * **directory** (str): Path to output directory + +setOutputDataMode(choice) + Set the data format for the main output mesh. + + * **choice** (int): + + * 0: ASCII format + * 1: Binary format + +setFracturesDataMode(choice) + Set the data format for fracture mesh files. + + * **choice** (int): + + * 0: ASCII format + * 1: Binary format + +Fracture Policies +----------------- + +**Internal Fractures Policy (0)** + +* Creates fractures within the mesh interior +* Duplicates nodes at internal interfaces +* Suitable for modeling embedded fracture networks + +**Boundary Fractures Policy (1)** + +* Creates fractures at mesh boundaries +* Handles fractures that extend to domain edges +* Suitable for modeling fault systems extending beyond the domain + +Results Access +-------------- + +getAllGrids() + Returns both the split mesh and individual fracture meshes. + + * **Returns**: tuple (split_mesh, fracture_meshes) + + * **split_mesh**: vtkUnstructuredGrid - Main mesh with fractures applied + * **fracture_meshes**: list[vtkUnstructuredGrid] - Individual fracture surfaces + +writeMeshes(filepath, is_data_mode_binary, canOverwrite) + Write all generated meshes to files. + + * **filepath** (str): Path for main split mesh + * **is_data_mode_binary** (bool): Use binary format + * **canOverwrite** (bool): Allow overwriting existing files + +Understanding Fracture Generation +--------------------------------- + +**Input Requirements** + +1. **Fracture Field**: Cell data array identifying regions separated by fractures +2. **Field Values**: Specific values indicating fracture boundaries +3. **Policy**: How to handle fracture creation + +**Process** + +1. **Identification**: Find interfaces between cells with different field values +2. **Node Duplication**: Create separate nodes for each side of the fracture +3. **Mesh Splitting**: Update connectivity to use duplicated nodes +4. **Fracture Extraction**: Generate separate meshes for fracture surfaces + +**Output Structure** + +* **Split Mesh**: Original mesh with fractures as discontinuities +* **Fracture Meshes**: Individual surface meshes for each fracture + +Common Use Cases +---------------- + +* **Geomechanics**: Modeling fault systems in geological domains +* **Fluid Flow**: Creating discrete fracture networks +* **Contact Mechanics**: Preparing meshes for contact simulations +* **Multi-physics**: Coupling different physics across fracture interfaces + +Example Workflow +---------------- + +.. code-block:: python + + # Complete fracture generation workflow + fracture_filter = GenerateFractures() + + # Configure fracture detection + fracture_filter.setFieldName("material_id") + fracture_filter.setFieldValues("1,2") # Fracture between materials 1 and 2 + fracture_filter.setPolicy(1) # Boundary fractures + + # Configure output + fracture_filter.setFracturesOutputDirectory("./fractures/") + fracture_filter.setOutputDataMode(1) # Binary for efficiency + fracture_filter.setFracturesDataMode(1) + + # Process mesh + fracture_filter.SetInputData(original_mesh) + fracture_filter.Update() + + # Get results + split_mesh, fracture_surfaces = fracture_filter.getAllGrids() + + print(f"Generated {len(fracture_surfaces)} fracture surfaces") + + # Write all outputs + fracture_filter.writeMeshes("output/domain_with_fractures.vtu") + +Output +------ + +* **Input**: vtkUnstructuredGrid with fracture identification field +* **Outputs**: + + * Split mesh with fractures as discontinuities + * Individual fracture surface meshes +* **File Output**: Automatic writing of fracture meshes to specified directory + +Best Practices +-------------- + +* Ensure fracture field values are properly assigned to cells +* Use appropriate policy based on fracture geometry +* Check that fractures form continuous surfaces +* Verify mesh quality after fracture generation +* Use binary format for large meshes to improve I/O performance + +See Also +-------- + +* :doc:`GenerateRectilinearGrid ` - Basic mesh generation +* :doc:`CollocatedNodes ` - May be needed after fracture generation +* :doc:`ElementVolumes ` - Quality check after splitting diff --git a/docs/geos_mesh_docs/filters/GenerateRectilinearGrid.rst b/docs/geos_mesh_docs/filters/GenerateRectilinearGrid.rst new file mode 100644 index 00000000..a29540ee --- /dev/null +++ b/docs/geos_mesh_docs/filters/GenerateRectilinearGrid.rst @@ -0,0 +1,218 @@ +GenerateRectilinearGrid Filter +============================== + +.. automodule:: geos.mesh.doctor.filters.GenerateRectilinearGrid + :members: + :undoc-members: + :show-inheritance: + +Overview +-------- + +The GenerateRectilinearGrid filter creates simple rectilinear (structured) grids as vtkUnstructuredGrid objects. This filter is useful for generating regular meshes for testing, validation, or as starting points for more complex mesh generation workflows. + +Features +-------- + +* Generation of 3D rectilinear grids with customizable dimensions +* Flexible block-based coordinate specification +* Variable element density per block +* Optional global ID generation for points and cells +* Customizable field generation with arbitrary dimensions +* Direct mesh generation without input requirements + +Usage Example +------------- + +.. code-block:: python + + from geos.mesh.doctor.filters.GenerateRectilinearGrid import GenerateRectilinearGrid + from geos.mesh.doctor.actions.generate_cube import FieldInfo + + # Instantiate the filter + generateRectilinearGridFilter = GenerateRectilinearGrid() + + # Set the coordinates of each block border for the X, Y and Z axis + generateRectilinearGridFilter.setCoordinates([0.0, 5.0, 10.0], [0.0, 5.0, 10.0], [0.0, 10.0]) + + # For each block defined, specify the number of cells that they should contain in the X, Y, Z axis + generateRectilinearGridFilter.setNumberElements([5, 5], [5, 5], [10]) + + # To add the GlobalIds for cells and points, set to True the generate global ids options + generateRectilinearGridFilter.setGenerateCellsGlobalIds(True) + generateRectilinearGridFilter.setGeneratePointsGlobalIds(True) + + # To create new arrays with a specific dimension, you can use the following commands + cells_dim1 = FieldInfo("cell1", 1, "CELLS") # array "cell1" of shape (number of cells, 1) + cells_dim3 = FieldInfo("cell3", 3, "CELLS") # array "cell3" of shape (number of cells, 3) + points_dim1 = FieldInfo("point1", 1, "POINTS") # array "point1" of shape (number of points, 1) + points_dim3 = FieldInfo("point3", 3, "POINTS") # array "point3" of shape (number of points, 3) + generateRectilinearGridFilter.setFields([cells_dim1, cells_dim3, points_dim1, points_dim3]) + + # Then, to obtain the constructed mesh out of all these operations, 2 solutions are available + + # Solution1 (recommended) + mesh = generateRectilinearGridFilter.getGrid() + + # Solution2 (manual) + generateRectilinearGridFilter.Update() + mesh = generateRectilinearGridFilter.GetOutputDataObject(0) + + # Finally, you can write the mesh at a specific destination with: + generateRectilinearGridFilter.writeGrid("output/filepath/of/your/grid.vtu") + +Parameters +---------- + +setCoordinates(coordsX, coordsY, coordsZ) + Set the coordinates that define block boundaries along each axis. + + * **coordsX** (Sequence[float]): X-coordinates of block boundaries + * **coordsY** (Sequence[float]): Y-coordinates of block boundaries + * **coordsZ** (Sequence[float]): Z-coordinates of block boundaries + +setNumberElements(numberElementsX, numberElementsY, numberElementsZ) + Set the number of elements in each block along each axis. + + * **numberElementsX** (Sequence[int]): Number of elements in each X-block + * **numberElementsY** (Sequence[int]): Number of elements in each Y-block + * **numberElementsZ** (Sequence[int]): Number of elements in each Z-block + +setGenerateCellsGlobalIds(generate) + Enable/disable generation of global cell IDs. + + * **generate** (bool): True to generate global cell IDs + +setGeneratePointsGlobalIds(generate) + Enable/disable generation of global point IDs. + + * **generate** (bool): True to generate global point IDs + +setFields(fields) + Specify additional cell or point arrays to be added to the grid. + + * **fields** (Iterable[FieldInfo]): List of field specifications + +Field Specification +------------------- + +Fields are specified using FieldInfo objects: + +.. code-block:: python + + from geos.mesh.doctor.actions.generate_cube import FieldInfo + + # Create a field specification + field = FieldInfo(name, dimension, location) + +**Parameters:** + +* **name** (str): Name of the array +* **dimension** (int): Number of components (1 for scalars, 3 for vectors, etc.) +* **location** (str): "CELLS" for cell data, "POINTS" for point data + +**Examples:** + +.. code-block:: python + + # Scalar cell field + pressure = FieldInfo("pressure", 1, "CELLS") + + # Vector point field + velocity = FieldInfo("velocity", 3, "POINTS") + + # Tensor cell field + stress = FieldInfo("stress", 6, "CELLS") # 6 components for symmetric tensor + +Grid Construction Logic +----------------------- + +**Coordinate Specification** + +Coordinates define block boundaries. For example: + +.. code-block:: python + + coordsX = [0.0, 5.0, 10.0] # Creates 2 blocks in X: [0,5] and [5,10] + numberElementsX = [5, 10] # 5 elements in first block, 10 in second + +**Element Distribution** + +Each block can have different element densities: + +* Block 1: 5 elements distributed uniformly in [0.0, 5.0] +* Block 2: 10 elements distributed uniformly in [5.0, 10.0] + +**Total Grid Size** + +* Total elements: numberElementsX[0] × numberElementsY[0] × numberElementsZ[0] + ... (for each block combination) +* Total points: (sum(numberElementsX) + len(numberElementsX)) × (sum(numberElementsY) + len(numberElementsY)) × (sum(numberElementsZ) + len(numberElementsZ)) + +Example: Multi-Block Grid +------------------------- + +.. code-block:: python + + # Create a grid with 3 blocks in X, 2 in Y, 1 in Z + filter = GenerateRectilinearGrid() + + # Define block boundaries + filter.setCoordinates( + [0.0, 1.0, 3.0, 5.0], # 3 blocks: [0,1], [1,3], [3,5] + [0.0, 2.0, 4.0], # 2 blocks: [0,2], [2,4] + [0.0, 1.0] # 1 block: [0,1] + ) + + # Define element counts per block + filter.setNumberElements( + [10, 20, 10], # 10, 20, 10 elements in X blocks + [15, 15], # 15, 15 elements in Y blocks + [5] # 5 elements in Z block + ) + + # Add global IDs and custom fields + filter.setGenerateCellsGlobalIds(True) + filter.setGeneratePointsGlobalIds(True) + + material_id = FieldInfo("material", 1, "CELLS") + coordinates = FieldInfo("coords", 3, "POINTS") + filter.setFields([material_id, coordinates]) + + # Generate the grid + mesh = filter.getGrid() + +Output +------ + +* **Input**: None (generator filter) +* **Output**: vtkUnstructuredGrid with hexahedral elements +* **Additional Arrays**: + + * Global IDs (if enabled) + * Custom fields (if specified) + * All arrays filled with constant value 1.0 + +Use Cases +--------- + +* **Testing**: Create simple grids for algorithm testing +* **Validation**: Generate known geometries for code validation +* **Prototyping**: Quick mesh generation for initial simulations +* **Benchmarking**: Standard grids for performance testing +* **Education**: Simple examples for learning mesh-based methods + +Best Practices +-------------- + +* Start with simple single-block grids before using multi-block configurations +* Ensure coordinate sequences are monotonically increasing +* Match the length of numberElements arrays with coordinate block counts +* Use meaningful field names for better mesh organization +* Enable global IDs when mesh will be used in parallel computations + +See Also +-------- + +* :doc:`GenerateFractures ` - Advanced mesh modification +* :doc:`ElementVolumes ` - Quality assessment for generated meshes +* :doc:`AllChecks ` - Comprehensive validation of generated meshes diff --git a/docs/geos_mesh_docs/filters/MainChecks.rst b/docs/geos_mesh_docs/filters/MainChecks.rst new file mode 100644 index 00000000..722e9818 --- /dev/null +++ b/docs/geos_mesh_docs/filters/MainChecks.rst @@ -0,0 +1,98 @@ +MainChecks Filter +================= + +.. autoclass:: geos.mesh.doctor.filters.Checks.MainChecks + :members: + :undoc-members: + :show-inheritance: + +Overview +-------- + +The MainChecks filter performs essential mesh validation by running the most important quality checks on a vtkUnstructuredGrid. This filter provides a streamlined subset of checks that are most critical for mesh quality assessment. + +Features +-------- + +* Essential mesh validation with the most important quality checks +* Faster execution compared to AllChecks +* Configurable check parameters +* Focused on critical mesh quality issues +* Same interface as AllChecks for consistency + +Usage Example +------------- + +.. code-block:: python + + from geos.mesh.doctor.filters.Checks import MainChecks + + # Instantiate the filter for main checks only + mainChecksFilter = MainChecks() + + # Set input mesh + mainChecksFilter.SetInputData(mesh) + + # Optionally customize check parameters + mainChecksFilter.setCheckParameter("collocated_nodes", "tolerance", 1e-6) + + # Execute the checks and get output + output_mesh = mainChecksFilter.getGrid() + + # Get check results + check_results = mainChecksFilter.getCheckResults() + + # Write the output mesh + mainChecksFilter.writeGrid("output/mesh_main_checks.vtu") + +Main Checks Included +-------------------- + +The MainChecks filter includes a curated subset of the most important checks: + +* **Collocated nodes**: Detect duplicate/overlapping nodes +* **Element volumes**: Identify negative or zero volume elements +* **Self-intersecting elements**: Find geometrically invalid elements +* **Basic topology validation**: Ensure mesh connectivity is valid + +These checks cover the most common and critical mesh quality issues that can affect simulation stability and accuracy. + +When to Use MainChecks vs AllChecks +----------------------------------- + +**Use MainChecks when:** + +* You need quick mesh validation +* You're doing routine quality checks +* Performance is important +* You want to focus on critical issues only + +**Use AllChecks when:** + +* You need comprehensive mesh analysis +* You're debugging complex mesh issues +* You have time for thorough validation +* You need detailed reporting on all aspects + +Parameters +---------- + +Same parameter interface as AllChecks: + +* **setCheckParameter(check_name, parameter_name, value)**: Set specific parameter for a named check +* **setGlobalParameter(parameter_name, value)**: Apply parameter to all checks that support it + +Output +------ + +* **Input**: vtkUnstructuredGrid +* **Output**: vtkUnstructuredGrid (copy of input with potential additional arrays marking issues) +* **Check Results**: Dictionary with results from performed main checks + +See Also +-------- + +* :doc:`AllChecks ` - Comprehensive mesh validation with all checks +* :doc:`CollocatedNodes ` - Individual collocated nodes check +* :doc:`ElementVolumes ` - Individual element volume check +* :doc:`SelfIntersectingElements ` - Individual self-intersection check diff --git a/docs/geos_mesh_docs/filters/NonConformal.rst b/docs/geos_mesh_docs/filters/NonConformal.rst new file mode 100644 index 00000000..44991af8 --- /dev/null +++ b/docs/geos_mesh_docs/filters/NonConformal.rst @@ -0,0 +1,250 @@ +NonConformal Filter +=================== + +.. automodule:: geos.mesh.doctor.filters.NonConformal + :members: + :undoc-members: + :show-inheritance: + +Overview +-------- + +The NonConformal filter detects non-conformal mesh interfaces in a vtkUnstructuredGrid. Non-conformal interfaces occur when adjacent cells do not share nodes or faces properly, which can indicate mesh quality issues or intentional non-matching grid interfaces that require special handling in simulations. + +Features +-------- + +* Detection of non-conformal interfaces between mesh elements +* Configurable tolerance parameters for different geometric aspects +* Optional marking of non-conformal cells in output mesh +* Support for point, face, and angle tolerance specifications +* Detailed reporting of non-conformal cell pairs + +Usage Example +------------- + +.. code-block:: python + + from geos.mesh.doctor.filters.NonConformal import NonConformal + + # Instantiate the filter + nonConformalFilter = NonConformal() + + # Set tolerance parameters + nonConformalFilter.setPointTolerance(1e-6) # tolerance for point matching + nonConformalFilter.setFaceTolerance(1e-6) # tolerance for face matching + nonConformalFilter.setAngleTolerance(10.0) # angle tolerance in degrees + + # Optionally enable painting of non-conformal cells + nonConformalFilter.setPaintNonConformalCells(1) # 1 to enable, 0 to disable + + # Set input mesh + nonConformalFilter.SetInputData(mesh) + + # Execute the filter and get output + output_mesh = nonConformalFilter.getGrid() + + # Get non-conformal cell pairs + non_conformal_cells = nonConformalFilter.getNonConformalCells() + # Returns list of tuples with (cell1_id, cell2_id) for non-conformal interfaces + + # Write the output mesh + nonConformalFilter.writeGrid("output/mesh_with_nonconformal_info.vtu") + +Parameters +---------- + +setPointTolerance(tolerance) + Set the tolerance for point position matching. + + * **tolerance** (float): Distance below which points are considered coincident + * **Default**: 0.0 + +setFaceTolerance(tolerance) + Set the tolerance for face geometry matching. + + * **tolerance** (float): Distance tolerance for face-to-face matching + * **Default**: 0.0 + +setAngleTolerance(tolerance) + Set the tolerance for face normal angle differences. + + * **tolerance** (float): Maximum angle difference in degrees + * **Default**: 10.0 + +setPaintNonConformalCells(choice) + Enable/disable creation of array marking non-conformal cells. + + * **choice** (int): 1 to enable marking, 0 to disable + * **Default**: 0 + +Results Access +-------------- + +getNonConformalCells() + Returns pairs of cell indices that have non-conformal interfaces. + + * **Returns**: list[tuple[int, int]] - Each tuple contains (cell1_id, cell2_id) + +getAngleTolerance() + Get the current angle tolerance setting. + + * **Returns**: float - Current angle tolerance in degrees + +getFaceTolerance() + Get the current face tolerance setting. + + * **Returns**: float - Current face tolerance + +getPointTolerance() + Get the current point tolerance setting. + + * **Returns**: float - Current point tolerance + +Understanding Non-Conformal Interfaces +--------------------------------------- + +**Conformal vs Non-Conformal** + +**Conformal Interface**: Adjacent cells share exact nodes and faces + +.. code-block:: + + Cell A: nodes [1, 2, 3, 4] + Cell B: nodes [3, 4, 5, 6] # Shares nodes 3,4 with Cell A + → CONFORMAL + +**Non-Conformal Interface**: Adjacent cells do not share nodes/faces exactly + +.. code-block:: + + Cell A: nodes [1, 2, 3, 4] + Cell B: nodes [5, 6, 7, 8] # No shared nodes with Cell A but geometrically adjacent + → NON-CONFORMAL + +**Types of Non-Conformity** + +1. **Hanging Nodes**: T-junctions where one element edge intersects another element face +2. **Mismatched Boundaries**: Interfaces where element boundaries don't align +3. **Gap Interfaces**: Small gaps between elements that should be connected +4. **Overlapping Interfaces**: Elements that overlap slightly due to meshing errors + +Tolerance Parameters Explained +------------------------------ + +**Point Tolerance** + +Controls how close points must be to be considered the same: + +* **Too small**: May miss near-coincident points that should match +* **Too large**: May incorrectly group distinct points +* **Typical values**: 1e-6 to 1e-12 depending on mesh scale + +**Face Tolerance** + +Controls how closely face geometries must match: + +* **Distance-based**: Maximum separation between face centroids or boundaries +* **Affects**: Detection of faces that should be conformal but have small gaps +* **Typical values**: 1e-6 to 1e-10 + +**Angle Tolerance** + +Controls how closely face normals must align: + +* **In degrees**: Maximum angle between face normal vectors +* **Affects**: Detection of faces that should be coplanar but have slight orientation differences +* **Typical values**: 0.1 to 10.0 degrees + +Common Causes of Non-Conformity +------------------------------- + +1. **Mesh Generation Issues**: + + * Different mesh densities in adjacent regions + * Boundary misalignment during mesh merging + * Floating-point precision errors + +2. **Intentional Design**: + + * Adaptive mesh refinement interfaces + * Multi-scale coupling boundaries + * Domain decomposition interfaces + +3. **Mesh Processing Errors**: + + * Node merging tolerances too strict + * Coordinate transformation errors + * File format conversion issues + +Impact on Simulations +--------------------- + +**Potential Problems**: + +* **Gaps**: Can cause fluid/heat leakage in flow simulations +* **Overlaps**: May create artificial constraints or stress concentrations +* **Inconsistent Physics**: Different discretizations across interfaces + +**When Non-Conformity is Acceptable**: + +* **Mortar Methods**: Designed to handle non-matching grids +* **Penalty Methods**: Use constraints to enforce continuity +* **Adaptive Refinement**: Temporary non-conformity during adaptation + +Example Analysis Workflow +------------------------- + +.. code-block:: python + + # Comprehensive non-conformity analysis + nc_filter = NonConformal() + + # Configure for sensitive detection + nc_filter.setPointTolerance(1e-8) + nc_filter.setFaceTolerance(1e-8) + nc_filter.setAngleTolerance(1.0) # 1 degree tolerance + + # Enable visualization + nc_filter.setPaintNonConformalCells(1) + + # Process mesh + nc_filter.SetInputData(mesh) + output_mesh = nc_filter.getGrid() + + # Analyze results + non_conformal_pairs = nc_filter.getNonConformalCells() + + if len(non_conformal_pairs) == 0: + print("Mesh is fully conformal - all interfaces match properly") + else: + print(f"Found {len(non_conformal_pairs)} non-conformal interfaces:") + for cell1, cell2 in non_conformal_pairs[:10]: # Show first 10 + print(f" Cells {cell1} and {cell2} have non-conformal interface") + + # Write mesh with marking for visualization + nc_filter.writeGrid("output/mesh_nonconformal_marked.vtu") + +Output +------ + +* **Input**: vtkUnstructuredGrid +* **Output**: vtkUnstructuredGrid with optional marking arrays +* **Marking Array**: When painting is enabled, adds "IsNonConformal" array to cell data +* **Cell Pairs**: List of non-conformal cell index pairs + +Best Practices +-------------- + +* **Set appropriate tolerances** based on mesh precision and simulation requirements +* **Use painting** to visualize non-conformal regions in ParaView +* **Consider physics requirements** when deciding if non-conformity is acceptable +* **Combine with other checks** for comprehensive mesh validation +* **Document intentional non-conformity** for future reference + +See Also +-------- + +* :doc:`AllChecks ` - Includes non-conformal check among others +* :doc:`CollocatedNodes ` - Related to point matching issues +* :doc:`SelfIntersectingElements ` - Related geometric validation diff --git a/docs/geos_mesh_docs/filters/SelfIntersectingElements.rst b/docs/geos_mesh_docs/filters/SelfIntersectingElements.rst new file mode 100644 index 00000000..059f7e56 --- /dev/null +++ b/docs/geos_mesh_docs/filters/SelfIntersectingElements.rst @@ -0,0 +1,293 @@ +SelfIntersectingElements Filter +=============================== + +.. automodule:: geos.mesh.doctor.filters.SelfIntersectingElements + :members: + :undoc-members: + :show-inheritance: + +Overview +-------- + +The SelfIntersectingElements filter identifies various types of invalid or problematic elements in a vtkUnstructuredGrid. It performs comprehensive geometric validation to detect elements with intersecting edges, intersecting faces, non-contiguous edges, non-convex shapes, incorrectly oriented faces, and wrong number of points. + +Features +-------- + +* Detection of multiple types of geometric element problems +* Configurable minimum distance parameter for intersection detection +* Optional marking of invalid elements in output mesh +* Detailed classification of different problem types +* Comprehensive reporting of all detected issues + +Usage Example +------------- + +.. code-block:: python + + from geos.mesh.doctor.filters.SelfIntersectingElements import SelfIntersectingElements + + # Instantiate the filter + selfIntersectingElementsFilter = SelfIntersectingElements() + + # Set minimum distance parameter for intersection detection + selfIntersectingElementsFilter.setMinDistance(1e-6) + + # Optionally enable painting of invalid elements + selfIntersectingElementsFilter.setPaintInvalidElements(1) # 1 to enable, 0 to disable + + # Set input mesh + selfIntersectingElementsFilter.SetInputData(mesh) + + # Execute the filter and get output + output_mesh = selfIntersectingElementsFilter.getGrid() + + # Get different types of problematic elements + wrong_points_elements = selfIntersectingElementsFilter.getWrongNumberOfPointsElements() + intersecting_edges_elements = selfIntersectingElementsFilter.getIntersectingEdgesElements() + intersecting_faces_elements = selfIntersectingElementsFilter.getIntersectingFacesElements() + non_contiguous_edges_elements = selfIntersectingElementsFilter.getNonContiguousEdgesElements() + non_convex_elements = selfIntersectingElementsFilter.getNonConvexElements() + wrong_oriented_faces_elements = selfIntersectingElementsFilter.getFacesOrientedIncorrectlyElements() + + # Write the output mesh + selfIntersectingElementsFilter.writeGrid("output/mesh_with_invalid_elements.vtu") + +Parameters +---------- + +setMinDistance(distance) + Set the minimum distance parameter for intersection detection. + + * **distance** (float): Minimum distance threshold for geometric calculations + * **Default**: 0.0 + * **Usage**: Smaller values detect more subtle problems, larger values ignore minor issues + +setPaintInvalidElements(choice) + Enable/disable creation of arrays marking invalid elements. + + * **choice** (int): 1 to enable marking, 0 to disable + * **Default**: 0 + * **Effect**: When enabled, adds arrays to cell data identifying different problem types + +Types of Problems Detected +-------------------------- + +getWrongNumberOfPointsElements() + Elements with incorrect number of points for their cell type. + + * **Returns**: list[int] - Element indices with wrong point counts + * **Examples**: Triangle with 4 points, hexahedron with 7 points + +getIntersectingEdgesElements() + Elements where edges intersect themselves. + + * **Returns**: list[int] - Element indices with self-intersecting edges + * **Examples**: Twisted quadrilaterals, folded triangles + +getIntersectingFacesElements() + Elements where faces intersect each other. + + * **Returns**: list[int] - Element indices with self-intersecting faces + * **Examples**: Inverted tetrahedra, twisted hexahedra + +getNonContiguousEdgesElements() + Elements where edges are not properly connected. + + * **Returns**: list[int] - Element indices with connectivity issues + * **Examples**: Disconnected edge loops, gaps in element boundaries + +getNonConvexElements() + Elements that are not convex as required. + + * **Returns**: list[int] - Element indices that are non-convex + * **Examples**: Concave quadrilaterals, non-convex polygons + +getFacesOrientedIncorrectlyElements() + Elements with incorrectly oriented faces. + + * **Returns**: list[int] - Element indices with orientation problems + * **Examples**: Inward-pointing face normals, inconsistent winding + +Understanding Element Problems +------------------------------ + +**Wrong Number of Points** + +Each VTK cell type has a specific number of points: + +* Triangle: 3 points +* Quadrilateral: 4 points +* Tetrahedron: 4 points +* Hexahedron: 8 points +* etc. + +**Self-Intersecting Edges** + +Edges that cross over themselves: + +.. code-block:: + + Valid triangle: Invalid triangle (bow-tie): + A A + / \ /|\ + / \ / | \ + B-----C B--+--C + \|/ + D + +**Self-Intersecting Faces** + +3D elements where faces intersect: + +.. code-block:: + + Valid tetrahedron Invalid tetrahedron (inverted) + D D + /|\ /|\ + / | \ / | \ + A--+--C C--+--A (face ABC flipped) + \ | / \ | / + \|/ \|/ + B B + +**Non-Contiguous Edges** + +Element boundaries that don't form continuous loops: + +* Missing edges between consecutive points +* Duplicate edges +* Gaps in the boundary + +**Non-Convex Elements** + +Elements that have internal angles > 180 degrees: + +* Can cause numerical issues in finite element calculations +* May indicate mesh generation problems +* Some solvers require strictly convex elements + +**Incorrectly Oriented Faces** + +Faces with normals pointing in wrong directions: + +* Outward normals pointing inward +* Inconsistent winding order +* Can affect normal-based calculations + +Common Causes and Solutions +--------------------------- + +**Wrong Number of Points** + +* **Cause**: Mesh file corruption, wrong cell type specification +* **Solution**: Fix cell type definitions or regenerate mesh + +**Self-Intersecting Edges/Faces** + +* **Cause**: Node coordinate errors, mesh deformation, bad mesh generation +* **Solution**: Fix node coordinates, improve mesh quality settings + +**Non-Contiguous Edges** + +* **Cause**: Missing connectivity information, duplicate nodes +* **Solution**: Fix element connectivity, merge duplicate nodes + +**Non-Convex Elements** + +* **Cause**: Poor mesh quality, extreme deformation +* **Solution**: Improve mesh generation parameters, element quality checks + +**Wrong Face Orientation** + +* **Cause**: Inconsistent node ordering, mesh processing errors +* **Solution**: Fix element node ordering, use mesh repair tools + +Example Comprehensive Analysis +------------------------------ + +.. code-block:: python + + # Detailed element validation workflow + validator = SelfIntersectingElements() + validator.setMinDistance(1e-8) # Very sensitive detection + validator.setPaintInvalidElements(1) # Enable visualization + + validator.SetInputData(mesh) + output_mesh = validator.getGrid() + + # Collect all problems + problems = { + 'Wrong points': validator.getWrongNumberOfPointsElements(), + 'Intersecting edges': validator.getIntersectingEdgesElements(), + 'Intersecting faces': validator.getIntersectingFacesElements(), + 'Non-contiguous edges': validator.getNonContiguousEdgesElements(), + 'Non-convex': validator.getNonConvexElements(), + 'Wrong orientation': validator.getFacesOrientedIncorrectlyElements() + } + + # Report results + total_problems = sum(len(elements) for elements in problems.values()) + + if total_problems == 0: + print("✓ All elements are geometrically valid!") + else: + print(f"⚠ Found {total_problems} problematic elements:") + for problem_type, elements in problems.items(): + if elements: + print(f" {problem_type}: {len(elements)} elements") + print(f" Examples: {elements[:5]}") # Show first 5 + + # Save results for visualization + validator.writeGrid("output/mesh_validation_results.vtu") + +Impact on Simulations +--------------------- + +**Numerical Issues** + +* Poor convergence +* Solver instabilities +* Incorrect results +* Simulation crashes + +**Physical Accuracy** + +* Wrong material volumes +* Incorrect flow paths +* Bad stress/strain calculations +* Energy conservation violations + +**Performance Impact** + +* Slower convergence +* Need for smaller time steps +* Additional stabilization methods +* Increased computational cost + +Output +------ + +* **Input**: vtkUnstructuredGrid +* **Output**: vtkUnstructuredGrid with optional marking arrays +* **Problem Lists**: Separate lists for each type of geometric problem +* **Marking Arrays**: When painting is enabled, adds arrays identifying problem types + +Best Practices +-------------- + +* **Set appropriate minimum distance** based on mesh precision +* **Enable painting** to visualize problems in ParaView +* **Check all problem types** for comprehensive validation +* **Fix problems before simulation** to avoid numerical issues +* **Use with other validators** for complete mesh assessment +* **Document any intentionally invalid elements** if they serve a purpose + +See Also +-------- + +* :doc:`AllChecks ` - Includes self-intersection check among others +* :doc:`MainChecks ` - Includes self-intersection check in main set +* :doc:`ElementVolumes ` - Related to geometric validity +* :doc:`CollocatedNodes ` - Can help fix some geometric issues +* :doc:`NonConformal ` - Related interface validation diff --git a/docs/geos_mesh_docs/filters/SupportedElements.rst b/docs/geos_mesh_docs/filters/SupportedElements.rst new file mode 100644 index 00000000..a2fb59e6 --- /dev/null +++ b/docs/geos_mesh_docs/filters/SupportedElements.rst @@ -0,0 +1,228 @@ +SupportedElements Filter +======================== + +.. automodule:: geos.mesh.doctor.filters.SupportedElements + :members: + :undoc-members: + :show-inheritance: + +Overview +-------- + +The SupportedElements filter identifies unsupported element types and problematic polyhedron elements in a vtkUnstructuredGrid. It validates that all elements in the mesh are supported by GEOS and checks polyhedron elements for geometric correctness. + +.. note:: + This filter is currently disabled due to multiprocessing requirements that are incompatible with the VTK filter framework. The implementation exists but is commented out in the source code. + +Features (When Available) +------------------------- + +* Detection of unsupported VTK element types +* Validation of polyhedron element geometry +* Optional marking of unsupported elements in output mesh +* Integration with parallel processing for large meshes +* Detailed reporting of element type compatibility + +Intended Usage Example +---------------------- + +.. code-block:: python + + from geos.mesh.doctor.filters.SupportedElements import SupportedElements + + # Instantiate the filter + supportedElementsFilter = SupportedElements() + + # Optionally enable painting of unsupported element types + supportedElementsFilter.setPaintUnsupportedElementTypes(1) # 1 to enable, 0 to disable + + # Set input mesh + supportedElementsFilter.SetInputData(mesh) + + # Execute the filter and get output + output_mesh = supportedElementsFilter.getGrid() + + # Get unsupported elements + unsupported_elements = supportedElementsFilter.getUnsupportedElements() + + # Write the output mesh + supportedElementsFilter.writeGrid("output/mesh_with_support_info.vtu") + +GEOS Supported Element Types +---------------------------- + +GEOS supports the following VTK element types: + +**Standard Elements** + +* **VTK_VERTEX** (1): Point elements +* **VTK_LINE** (3): Line segments +* **VTK_TRIANGLE** (5): Triangular elements +* **VTK_QUAD** (9): Quadrilateral elements +* **VTK_TETRA** (10): Tetrahedral elements +* **VTK_HEXAHEDRON** (12): Hexahedral (brick) elements +* **VTK_WEDGE** (13): Wedge/prism elements +* **VTK_PYRAMID** (14): Pyramid elements + +**Higher-Order Elements** + +* **VTK_QUADRATIC_TRIANGLE** (22): 6-node triangles +* **VTK_QUADRATIC_QUAD** (23): 8-node quadrilaterals +* **VTK_QUADRATIC_TETRA** (24): 10-node tetrahedra +* **VTK_QUADRATIC_HEXAHEDRON** (25): 20-node hexahedra + +**Special Elements** + +* **VTK_POLYHEDRON** (42): General polyhedra (with validation) + +Unsupported Element Types +------------------------- + +Elements not supported by GEOS include: + +* **VTK_PIXEL** (8): Axis-aligned rectangles +* **VTK_VOXEL** (11): Axis-aligned boxes +* **VTK_PENTAGONAL_PRISM** (15): 5-sided prisms +* **VTK_HEXAGONAL_PRISM** (16): 6-sided prisms +* Various specialized or experimental VTK cell types + +Polyhedron Validation +--------------------- + +For polyhedron elements (VTK_POLYHEDRON), additional checks are performed: + +**Geometric Validation** + +* Face planarity +* Edge connectivity +* Volume calculation +* Normal consistency + +**Topological Validation** + +* Manifold surface verification +* Closed surface check +* Face orientation consistency + +**Quality Checks** + +* Aspect ratio limits +* Volume positivity +* Face area positivity + +Common Issues and Solutions +--------------------------- + +**Unsupported Standard Elements** + +* **Problem**: Mesh contains VTK_PIXEL or VTK_VOXEL elements +* **Solution**: Convert to VTK_QUAD or VTK_HEXAHEDRON respectively +* **Tools**: VTK conversion filters or mesh processing software + +**Invalid Polyhedra** + +* **Problem**: Non-manifold or self-intersecting polyhedra +* **Solution**: Use mesh repair tools or regenerate with better quality settings +* **Prevention**: Validate polyhedra during mesh generation + +**Mixed Element Types** + +* **Problem**: Mesh contains both supported and unsupported elements +* **Solution**: Selective element type conversion or mesh region separation + +Current Status and Alternatives +------------------------------- + +**Why Disabled** + +The SupportedElements filter requires multiprocessing capabilities for efficient polyhedron validation on large meshes. However, the VTK Python filter framework doesn't integrate well with multiprocessing, leading to: + +* Process spawning issues +* Memory management problems +* Inconsistent results across platforms + +**Alternative Approaches** + +1. **Command-Line Tool**: Use the ``mesh-doctor supported_elements`` command instead +2. **Direct Function Calls**: Import and use the underlying functions directly +3. **Manual Validation**: Check element types programmatically + +**Command-Line Alternative** + +.. code-block:: bash + + # Use mesh-doctor command line tool instead + mesh-doctor -i input_mesh.vtu supported_elements --help + +**Direct Function Usage** + +.. code-block:: python + + from geos.mesh.doctor.actions.supported_elements import ( + find_unsupported_std_elements_types, + find_unsupported_polyhedron_elements + ) + + # Direct function usage (not as VTK filter) + unsupported_std = find_unsupported_std_elements_types(mesh) + # Note: polyhedron validation requires multiprocessing setup + +Future Development +------------------ + +**Planned Improvements** + +* Integration with VTK's parallel processing capabilities +* Alternative implementation without multiprocessing dependency +* Better error handling and reporting +* Performance optimizations for large meshes + +**Workaround Implementation** + +Until the filter is re-enabled, users can: + +1. Use the command-line interface +2. Implement custom validation loops +3. Use external mesh validation tools +4. Perform validation in separate processes + +Example Manual Validation +------------------------- + +.. code-block:: python + + import vtk + + def check_supported_elements(mesh): + """Manual check for supported element types.""" + supported_types = { + vtk.VTK_VERTEX, vtk.VTK_LINE, vtk.VTK_TRIANGLE, vtk.VTK_QUAD, + vtk.VTK_TETRA, vtk.VTK_HEXAHEDRON, vtk.VTK_WEDGE, vtk.VTK_PYRAMID, + vtk.VTK_QUADRATIC_TRIANGLE, vtk.VTK_QUADRATIC_QUAD, + vtk.VTK_QUADRATIC_TETRA, vtk.VTK_QUADRATIC_HEXAHEDRON, + vtk.VTK_POLYHEDRON + } + + unsupported = [] + for i in range(mesh.GetNumberOfCells()): + cell_type = mesh.GetCellType(i) + if cell_type not in supported_types: + unsupported.append((i, cell_type)) + + return unsupported + + # Usage + unsupported_elements = check_supported_elements(mesh) + if unsupported_elements: + print(f"Found {len(unsupported_elements)} unsupported elements") + for cell_id, cell_type in unsupported_elements[:5]: + print(f" Cell {cell_id}: type {cell_type}") + +See Also +-------- + +* :doc:`AllChecks ` - Would include supported elements check when available +* :doc:`SelfIntersectingElements ` - Related geometric validation +* :doc:`ElementVolumes ` - Basic element validation +* GEOS documentation on supported element types +* VTK documentation on cell types diff --git a/docs/geos_mesh_docs/filters/index.rst b/docs/geos_mesh_docs/filters/index.rst new file mode 100644 index 00000000..456dd2e8 --- /dev/null +++ b/docs/geos_mesh_docs/filters/index.rst @@ -0,0 +1,192 @@ +Mesh Doctor Filters +=================== + +The mesh doctor filters provide VTK-based tools for mesh quality assessment, validation, and processing. All filters work with vtkUnstructuredGrid data and follow a consistent interface pattern. + +Quality Assessment Filters +-------------------------- + +These filters analyze existing meshes for various quality issues and geometric problems. + +.. toctree:: + :maxdepth: 1 + + AllChecks + MainChecks + CollocatedNodes + ElementVolumes + SelfIntersectingElements + NonConformal + +Mesh Generation Filters +----------------------- + +These filters create new meshes from scratch or modify existing meshes. + +.. toctree:: + :maxdepth: 1 + + GenerateRectilinearGrid + GenerateFractures + +Processing Filters +------------------ + +These filters perform specialized processing and validation tasks. + +.. toctree:: + :maxdepth: 1 + + SupportedElements + +Common Usage Pattern +==================== + +All mesh doctor filters follow a consistent usage pattern: + +.. code-block:: python + + from geos.mesh.doctor.filters.FilterName import FilterName + + # Instantiate the filter + filter = FilterName() + + # Configure filter parameters + filter.setParameter(value) + + # Set input data (for processing filters, not needed for generators) + filter.SetInputData(mesh) + + # Execute the filter and get output + output_mesh = filter.getGrid() + + # Access specific results (filter-dependent) + results = filter.getSpecificResults() + + # Write results to file + filter.writeGrid("output/result.vtu") + +Filter Categories Explained +=========================== + +Quality Assessment +------------------ + +**Purpose**: Identify mesh quality issues, geometric problems, and topology errors + +**When to use**: +- Before running simulations +- After mesh generation or modification +- During mesh debugging +- For mesh quality reporting + +**Key filters**: +- **AllChecks/MainChecks**: Comprehensive validation suites +- **CollocatedNodes**: Duplicate node detection +- **ElementVolumes**: Volume validation +- **SelfIntersectingElements**: Geometric integrity +- **NonConformal**: Interface validation + +Mesh Generation +--------------- + +**Purpose**: Create new meshes or modify existing ones + +**When to use**: +- Creating test meshes +- Generating simple geometries +- Adding fractures to existing meshes +- Prototyping mesh-based algorithms + +**Key filters**: +- **GenerateRectilinearGrid**: Simple structured grids +- **GenerateFractures**: Fracture network generation + +Processing +---------- + +**Purpose**: Specialized mesh processing and validation + +**When to use**: +- Validating element type compatibility +- Preparing meshes for specific solvers +- Advanced geometric analysis + +**Key filters**: +- **SupportedElements**: GEOS compatibility validation + +Quick Reference +=============== + +Filter Selection Guide +---------------------- + +**For routine mesh validation**: + Use :doc:`MainChecks ` for fast, essential checks + +**For comprehensive analysis**: + Use :doc:`AllChecks ` for detailed validation + +**For specific problems**: + - Duplicate nodes → :doc:`CollocatedNodes ` + - Negative volumes → :doc:`ElementVolumes ` + - Invalid geometry → :doc:`SelfIntersectingElements ` + - Interface issues → :doc:`NonConformal ` + +**For mesh generation**: + - Simple grids → :doc:`GenerateRectilinearGrid ` + - Fracture networks → :doc:`GenerateFractures ` + +**For compatibility checking**: + - GEOS support → :doc:`SupportedElements ` + +Parameter Guidelines +-------------------- + +**Tolerance Parameters**: + - High precision meshes: 1e-12 to 1e-8 + - Standard meshes: 1e-8 to 1e-6 + - Coarse meshes: 1e-6 to 1e-4 + +**Painting Options**: + - Enable (1) for visualization in ParaView + - Disable (0) for performance in batch processing + +**Output Modes**: + - Binary for large meshes and performance + - ASCII for debugging and text processing + +Best Practices +============== + +Workflow Integration +-------------------- + +1. **Start with quality assessment** using MainChecks or AllChecks +2. **Address specific issues** with targeted filters +3. **Validate fixes** by re-running quality checks +4. **Document mesh quality** for simulation reference + +Performance Considerations +-------------------------- + +- Use appropriate tolerances (not unnecessarily strict) +- Enable painting only when needed for visualization +- Use binary output for large meshes +- Run comprehensive checks during development, lighter checks in production + +Error Handling +-------------- + +- Check filter results before proceeding with simulations +- Save problematic meshes for debugging +- Document known issues and their acceptable thresholds +- Use multiple validation approaches for critical applications + +See Also +======== + +- **GEOS Documentation**: Supported element types and mesh requirements +- **VTK Documentation**: VTK data formats and cell types +- **ParaView**: Visualization of mesh quality results +- **Mesh Generation Tools**: Creating high-quality input meshes diff --git a/geos-mesh/src/geos/mesh/doctor/actions/all_checks.py b/geos-mesh/src/geos/mesh/doctor/actions/all_checks.py index 253165d9..f8f8d2c2 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/all_checks.py +++ b/geos-mesh/src/geos/mesh/doctor/actions/all_checks.py @@ -1,4 +1,6 @@ from dataclasses import dataclass +from typing import Any +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid from geos.mesh.doctor.register import __load_module_action from geos.mesh.doctor.parsing.cli_parsing import setup_logger @@ -6,21 +8,33 @@ @dataclass( frozen=True ) class Options: checks_to_perform: list[ str ] - checks_options: dict[ str, any ] - check_displays: dict[ str, any ] + checks_options: dict[ str, Any ] + check_displays: dict[ str, Any ] @dataclass( frozen=True ) class Result: - check_results: dict[ str, any ] + check_results: dict[ str, Any ] -def action( vtk_input_file: str, options: Options ) -> list[ Result ]: +def get_check_results( vtk_input: str | vtkUnstructuredGrid, options: Options ) -> dict[ str, Any ]: + isFilepath: bool = isinstance( vtk_input, str ) + isVtkUnstructuredGrid: bool = isinstance( vtk_input, vtkUnstructuredGrid ) + assert isFilepath | isVtkUnstructuredGrid, "Invalid input type, should either be a filepath to .vtu file" \ + " or a vtkUnstructuredGrid object" check_results: dict[ str, any ] = dict() for check_name in options.checks_to_perform: - check_action = __load_module_action( check_name ) + if isVtkUnstructuredGrid: # we need to call the mesh_action function that takes a vtkPointSet as input + check_action = __load_module_action( check_name, "mesh_action" ) + else: # because its a filepath, we can call the regular action function + check_action = __load_module_action( check_name ) setup_logger.info( f"Performing check '{check_name}'." ) option = options.checks_options[ check_name ] - check_result = check_action( vtk_input_file, option ) + check_result = check_action( vtk_input, option ) check_results[ check_name ] = check_result + return check_results + + +def action( vtk_input_file: str, options: Options ) -> Result: + check_results: dict[ str, Any ] = get_check_results( vtk_input_file, options ) return Result( check_results=check_results ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/collocated_nodes.py b/geos-mesh/src/geos/mesh/doctor/actions/collocated_nodes.py index 4881a1d3..aa6d2276 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/collocated_nodes.py +++ b/geos-mesh/src/geos/mesh/doctor/actions/collocated_nodes.py @@ -1,11 +1,10 @@ from collections import defaultdict from dataclasses import dataclass import numpy -from typing import Collection, Iterable from vtkmodules.vtkCommonCore import reference, vtkPoints -from vtkmodules.vtkCommonDataModel import vtkIncrementalOctreePointLocator +from vtkmodules.vtkCommonDataModel import vtkCell, vtkIncrementalOctreePointLocator, vtkUnstructuredGrid from geos.mesh.doctor.parsing.cli_parsing import setup_logger -from geos.mesh.io.vtkIO import read_mesh +from geos.mesh.io.vtkIO import read_unstructured_grid @dataclass( frozen=True ) @@ -15,23 +14,22 @@ class Options: @dataclass( frozen=True ) class Result: - nodes_buckets: Iterable[ Iterable[ int ] ] # Each bucket contains the duplicated node indices. - wrong_support_elements: Collection[ int ] # Element indices with support node indices appearing more than once. + nodes_buckets: list[ tuple[ int ] ] # Each bucket contains the duplicated node indices. + wrong_support_elements: list[ int ] # Element indices with support node indices appearing more than once. -def __action( mesh, options: Options ) -> Result: - points = mesh.GetPoints() - +def find_collocated_nodes_buckets( mesh: vtkUnstructuredGrid, tolerance: float ) -> list[ tuple[ int ] ]: + points: vtkPoints = mesh.GetPoints() locator = vtkIncrementalOctreePointLocator() - locator.SetTolerance( options.tolerance ) + locator.SetTolerance( tolerance ) output = vtkPoints() locator.InitPointInsertion( output, points.GetBounds() ) # original ids to/from filtered ids. filtered_to_original = numpy.ones( points.GetNumberOfPoints(), dtype=int ) * -1 - rejected_points = defaultdict( list ) - point_id = reference( 0 ) + rejected_points: dict[ int, list[ int ] ] = defaultdict( list ) + point_id: int = reference( 0 ) for i in range( points.GetNumberOfPoints() ): is_inserted = locator.InsertUniquePoint( points.GetPoint( i ), point_id ) if not is_inserted: @@ -48,21 +46,29 @@ def __action( mesh, options: Options ) -> Result: # original_to_filtered[i] = point_id.get() filtered_to_original[ point_id.get() ] = i - tmp = [] + collocated_nodes_buckets: list[ tuple[ int ] ] = list() for n, ns in rejected_points.items(): - tmp.append( ( n, *ns ) ) + collocated_nodes_buckets.append( ( n, *ns ) ) + return collocated_nodes_buckets + +def find_wrong_support_elements( mesh: vtkUnstructuredGrid ) -> list[ int ]: # Checking that the support node indices appear only once per element. - wrong_support_elements = [] + wrong_support_elements: list[ int ] = list() for c in range( mesh.GetNumberOfCells() ): - cell = mesh.GetCell( c ) - num_points_per_cell = cell.GetNumberOfPoints() + cell: vtkCell = mesh.GetCell( c ) + num_points_per_cell: int = cell.GetNumberOfPoints() if len( { cell.GetPointId( i ) for i in range( num_points_per_cell ) } ) != num_points_per_cell: wrong_support_elements.append( c ) + return wrong_support_elements + - return Result( nodes_buckets=tmp, wrong_support_elements=wrong_support_elements ) +def mesh_action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + collocated_nodes_buckets = find_collocated_nodes_buckets( mesh, options.tolerance ) + wrong_support_elements = find_wrong_support_elements( mesh ) + return Result( nodes_buckets=collocated_nodes_buckets, wrong_support_elements=wrong_support_elements ) def action( vtk_input_file: str, options: Options ) -> Result: - mesh = read_mesh( vtk_input_file ) - return __action( mesh, options ) + mesh: vtkUnstructuredGrid = read_unstructured_grid( vtk_input_file ) + return mesh_action( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/element_volumes.py b/geos-mesh/src/geos/mesh/doctor/actions/element_volumes.py index e5380c3c..888235e3 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/element_volumes.py +++ b/geos-mesh/src/geos/mesh/doctor/actions/element_volumes.py @@ -1,11 +1,11 @@ from dataclasses import dataclass from typing import List, Tuple import uuid -from vtkmodules.vtkCommonDataModel import VTK_HEXAHEDRON, VTK_PYRAMID, VTK_TETRA, VTK_WEDGE +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, VTK_HEXAHEDRON, VTK_PYRAMID, VTK_TETRA, VTK_WEDGE from vtkmodules.vtkFiltersVerdict import vtkCellSizeFilter, vtkMeshQuality from vtkmodules.util.numpy_support import vtk_to_numpy from geos.mesh.doctor.parsing.cli_parsing import setup_logger -from geos.mesh.io.vtkIO import read_mesh +from geos.mesh.io.vtkIO import read_unstructured_grid @dataclass( frozen=True ) @@ -18,7 +18,7 @@ class Result: element_volumes: List[ Tuple[ int, float ] ] -def __action( mesh, options: Options ) -> Result: +def mesh_action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: cs = vtkCellSizeFilter() cs.ComputeAreaOff() @@ -67,5 +67,5 @@ def __action( mesh, options: Options ) -> Result: def action( vtk_input_file: str, options: Options ) -> Result: - mesh = read_mesh( vtk_input_file ) - return __action( mesh, options ) + mesh: vtkUnstructuredGrid = read_unstructured_grid( vtk_input_file ) + return mesh_action( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/fix_elements_orderings.py b/geos-mesh/src/geos/mesh/doctor/actions/fix_elements_orderings.py index 3e00cf52..3947ce51 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/fix_elements_orderings.py +++ b/geos-mesh/src/geos/mesh/doctor/actions/fix_elements_orderings.py @@ -1,8 +1,9 @@ from dataclasses import dataclass from typing import Dict, FrozenSet, List, Set from vtkmodules.vtkCommonCore import vtkIdList +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid from geos.mesh.utils.genericHelpers import to_vtk_id_list -from geos.mesh.io.vtkIO import VtkOutput, read_mesh, write_mesh +from geos.mesh.io.vtkIO import VtkOutput, read_unstructured_grid, write_mesh @dataclass( frozen=True ) @@ -17,7 +18,7 @@ class Result: unchanged_cell_types: FrozenSet[ int ] -def __action( mesh, options: Options ) -> Result: +def mesh_action( mesh: vtkUnstructuredGrid, 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 @@ -49,5 +50,5 @@ def __action( mesh, options: Options ) -> Result: def action( vtk_input_file: str, options: Options ) -> Result: - mesh = read_mesh( vtk_input_file ) - return __action( mesh, options ) + mesh: vtkUnstructuredGrid = read_unstructured_grid( vtk_input_file ) + return mesh_action( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/generate_cube.py b/geos-mesh/src/geos/mesh/doctor/actions/generate_cube.py index f30d2089..782f1a08 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/generate_cube.py +++ b/geos-mesh/src/geos/mesh/doctor/actions/generate_cube.py @@ -1,13 +1,15 @@ from dataclasses import dataclass -import numpy +import numpy as np +import numpy.typing as npt from typing import Iterable, Sequence from vtkmodules.util.numpy_support import numpy_to_vtk from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkCommonDataModel import ( vtkCellArray, vtkHexahedron, vtkRectilinearGrid, vtkUnstructuredGrid, VTK_HEXAHEDRON ) -from geos.mesh.doctor.actions.generate_global_ids import __build_global_ids +from geos.mesh.doctor.actions.generate_global_ids import build_global_ids from geos.mesh.doctor.parsing.cli_parsing import setup_logger from geos.mesh.io.vtkIO import VtkOutput, write_mesh +from geos.mesh.utils.arrayModifiers import createConstantAttributeDataSet @dataclass( frozen=True ) @@ -38,9 +40,64 @@ class Options: @dataclass( frozen=True ) class XYZ: - x: numpy.ndarray - y: numpy.ndarray - z: numpy.ndarray + x: npt.NDArray + y: npt.NDArray + z: npt.NDArray + + +def build_coordinates( positions, num_elements ): + result = [] + it = zip( zip( positions, positions[ 1: ] ), num_elements ) + try: + coords, n = next( it ) + while True: + start, stop = coords + end_point = False + tmp = np.linspace( start=start, stop=stop, num=n + end_point, endpoint=end_point ) + coords, n = next( it ) + result.append( tmp ) + except StopIteration: + end_point = True + tmp = np.linspace( start=start, stop=stop, num=n + end_point, endpoint=end_point ) + result.append( tmp ) + return np.concatenate( result ) + + +def build_rectilinear_grid( x: npt.NDArray, y: npt.NDArray, z: npt.NDArray ) -> vtkUnstructuredGrid: + """ + Builds an unstructured vtk grid from the x,y,z coordinates. + :return: The unstructured mesh, even if it's topologically structured. + """ + rg = vtkRectilinearGrid() + rg.SetDimensions( len( x ), len( y ), len( z ) ) + rg.SetXCoordinates( numpy_to_vtk( x ) ) + rg.SetYCoordinates( numpy_to_vtk( y ) ) + rg.SetZCoordinates( numpy_to_vtk( z ) ) + + num_points = rg.GetNumberOfPoints() + num_cells = rg.GetNumberOfCells() + + points = vtkPoints() + points.Allocate( num_points ) + for i in range( rg.GetNumberOfPoints() ): + points.InsertNextPoint( rg.GetPoint( i ) ) + + cell_types = [ VTK_HEXAHEDRON ] * num_cells + cells = vtkCellArray() + cells.AllocateExact( num_cells, num_cells * 8 ) + + m = ( 0, 1, 3, 2, 4, 5, 7, 6 ) # VTK_VOXEL and VTK_HEXAHEDRON do not share the same ordering. + for i in range( rg.GetNumberOfCells() ): + c = rg.GetCell( i ) + new_cell = vtkHexahedron() + for j in range( 8 ): + new_cell.GetPointIds().SetId( j, c.GetPointId( m[ j ] ) ) + cells.InsertNextCell( new_cell ) + + mesh = vtkUnstructuredGrid() + mesh.SetPoints( points ) + mesh.SetCells( cell_types, cells ) + return mesh def build_rectilinear_blocks_mesh( xyzs: Iterable[ XYZ ] ) -> vtkUnstructuredGrid: @@ -89,50 +146,37 @@ def build_rectilinear_blocks_mesh( xyzs: Iterable[ XYZ ] ) -> vtkUnstructuredGri return mesh -def __add_fields( mesh: vtkUnstructuredGrid, fields: Iterable[ FieldInfo ] ) -> vtkUnstructuredGrid: +def add_fields( mesh: vtkUnstructuredGrid, fields: Iterable[ FieldInfo ] ) -> vtkUnstructuredGrid: + """ + Add constant fields to the mesh using arrayModifiers utilities. + Each field is filled with ones (1.0) for all components. + """ for field_info in fields: - if field_info.support == "CELLS": - data = mesh.GetCellData() - n = mesh.GetNumberOfCells() - elif field_info.support == "POINTS": - data = mesh.GetPointData() - n = mesh.GetNumberOfPoints() - array = numpy.ones( ( n, field_info.dimension ), dtype=float ) - vtk_array = numpy_to_vtk( array ) - vtk_array.SetName( field_info.name ) - data.AddArray( vtk_array ) + onPoints = field_info.support == "POINTS" + # Create list of values (all 1.0) for each component + listValues = [ 1.0 ] * field_info.dimension + # Use the robust createConstantAttributeDataSet function + success = createConstantAttributeDataSet( dataSet=mesh, + listValues=listValues, + attributeName=field_info.name, + onPoints=onPoints ) + if not success: + setup_logger.warning( f"Failed to create field {field_info.name}" ) return mesh def __build( options: Options ): - def build_coordinates( positions, num_elements ): - result = [] - it = zip( zip( positions, positions[ 1: ] ), num_elements ) - try: - coords, n = next( it ) - while True: - start, stop = coords - end_point = False - tmp = numpy.linspace( start=start, stop=stop, num=n + end_point, endpoint=end_point ) - coords, n = next( it ) - result.append( tmp ) - except StopIteration: - end_point = True - tmp = numpy.linspace( start=start, stop=stop, num=n + end_point, endpoint=end_point ) - result.append( tmp ) - return numpy.concatenate( result ) - x = build_coordinates( options.xs, options.nxs ) y = build_coordinates( options.ys, options.nys ) z = build_coordinates( options.zs, options.nzs ) cube = build_rectilinear_blocks_mesh( ( XYZ( x, y, z ), ) ) - cube = __add_fields( cube, options.fields ) - __build_global_ids( cube, options.generate_cells_global_ids, options.generate_points_global_ids ) + cube = add_fields( cube, options.fields ) + build_global_ids( cube, options.generate_cells_global_ids, options.generate_points_global_ids ) return cube -def __action( options: Options ) -> Result: +def mesh_action( options: Options ) -> Result: output_mesh = __build( options ) write_mesh( output_mesh, options.vtk_output ) return Result( info=f"Mesh was written to {options.vtk_output.output}" ) @@ -140,7 +184,7 @@ def __action( options: Options ) -> Result: def action( vtk_input_file: str, options: Options ) -> Result: try: - return __action( options ) + return mesh_action( options ) except BaseException as e: setup_logger.error( e ) return Result( info="Something went wrong." ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/generate_fractures.py b/geos-mesh/src/geos/mesh/doctor/actions/generate_fractures.py index 32f809db..953be3f1 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/generate_fractures.py +++ b/geos-mesh/src/geos/mesh/doctor/actions/generate_fractures.py @@ -15,7 +15,7 @@ from geos.mesh.doctor.parsing.cli_parsing import setup_logger from geos.mesh.utils.arrayHelpers import has_array from geos.mesh.utils.genericHelpers import to_vtk_id_list, vtk_iter -from geos.mesh.io.vtkIO import VtkOutput, read_mesh, write_mesh +from geos.mesh.io.vtkIO import VtkOutput, read_unstructured_grid, write_mesh """ TypeAliases cannot be used with Python 3.9. A simple assignment like described there will be used: https://docs.python.org/3/library/typing.html#typing.TypeAlias:~:text=through%20simple%20assignment%3A-,Vector%20%3D%20list%5Bfloat%5D,-Or%20marked%20with @@ -527,8 +527,8 @@ def __generate_fracture_mesh( old_mesh: vtkUnstructuredGrid, fracture_info: Frac return fracture_mesh -def __split_mesh_on_fractures( mesh: vtkUnstructuredGrid, - options: Options ) -> tuple[ vtkUnstructuredGrid, list[ vtkUnstructuredGrid ] ]: +def split_mesh_on_fractures( mesh: vtkUnstructuredGrid, + options: Options ) -> tuple[ vtkUnstructuredGrid, list[ vtkUnstructuredGrid ] ]: all_fracture_infos: list[ FractureInfo ] = list() for fracture_id in range( len( options.field_values_per_fracture ) ): fracture_info: FractureInfo = build_fracture_info( mesh, options, False, fracture_id ) @@ -546,8 +546,8 @@ def __split_mesh_on_fractures( mesh: vtkUnstructuredGrid, return ( output_mesh, fracture_meshes ) -def __action( mesh, options: Options ) -> Result: - output_mesh, fracture_meshes = __split_mesh_on_fractures( mesh, options ) +def mesh_action( mesh, options: Options ) -> Result: + output_mesh, fracture_meshes = split_mesh_on_fractures( mesh, options ) write_mesh( output_mesh, options.mesh_VtkOutput ) for i, fracture_mesh in enumerate( fracture_meshes ): write_mesh( fracture_mesh, options.all_fractures_VtkOutput[ i ] ) @@ -557,14 +557,14 @@ def __action( mesh, options: Options ) -> Result: def action( vtk_input_file: str, options: Options ) -> Result: try: - mesh = read_mesh( vtk_input_file ) + mesh: vtkUnstructuredGrid = read_unstructured_grid( vtk_input_file ) # Mesh cannot contain global ids before splitting. if has_array( mesh, [ "GLOBAL_IDS_POINTS", "GLOBAL_IDS_CELLS" ] ): err_msg: str = ( "The mesh cannot contain global ids for neither cells nor points. The correct procedure " + " is to split the mesh and then generate global ids for new split meshes." ) setup_logger.error( err_msg ) raise ValueError( err_msg ) - return __action( mesh, options ) + return mesh_action( mesh, options ) except BaseException as e: setup_logger.error( e ) return Result( info="Something went wrong" ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/generate_global_ids.py b/geos-mesh/src/geos/mesh/doctor/actions/generate_global_ids.py index f4df1871..e73e21ca 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/generate_global_ids.py +++ b/geos-mesh/src/geos/mesh/doctor/actions/generate_global_ids.py @@ -1,7 +1,8 @@ from dataclasses import dataclass from vtkmodules.vtkCommonCore import vtkIdTypeArray +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid from geos.mesh.doctor.parsing.cli_parsing import setup_logger -from geos.mesh.io.vtkIO import VtkOutput, read_mesh, write_mesh +from geos.mesh.io.vtkIO import VtkOutput, read_unstructured_grid, write_mesh @dataclass( frozen=True ) @@ -16,7 +17,7 @@ class Result: info: str -def __build_global_ids( mesh, generate_cells_global_ids: bool, generate_points_global_ids: bool ) -> None: +def build_global_ids( mesh, generate_cells_global_ids: bool, generate_points_global_ids: bool ) -> None: """ Adds the global ids for cells and points in place into the mesh instance. :param mesh: @@ -45,16 +46,16 @@ def __build_global_ids( mesh, generate_cells_global_ids: bool, generate_points_g mesh.GetCellData().SetGlobalIds( cells_global_ids ) -def __action( mesh, options: Options ) -> Result: - __build_global_ids( mesh, options.generate_cells_global_ids, options.generate_points_global_ids ) +def mesh_action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + build_global_ids( mesh, options.generate_cells_global_ids, options.generate_points_global_ids ) write_mesh( mesh, options.vtk_output ) return Result( info=f"Mesh was written to {options.vtk_output.output}" ) def action( vtk_input_file: str, options: Options ) -> Result: try: - mesh = read_mesh( vtk_input_file ) - return __action( mesh, options ) + mesh: vtkUnstructuredGrid = read_unstructured_grid( vtk_input_file ) + return mesh_action( mesh, options ) except BaseException as e: setup_logger.error( e ) return Result( info="Something went wrong." ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/non_conformal.py b/geos-mesh/src/geos/mesh/doctor/actions/non_conformal.py index d1c83a37..728954ad 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/non_conformal.py +++ b/geos-mesh/src/geos/mesh/doctor/actions/non_conformal.py @@ -1,10 +1,11 @@ from dataclasses import dataclass import math -import numpy +import numpy as np +import numpy.typing as npt from tqdm import tqdm -from typing import List, Tuple, Any +from typing import Any, Sequence from vtk import reference as vtk_reference -from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints +from vtkmodules.vtkCommonCore import vtkDataArray, vtkIdList, vtkPoints from vtkmodules.vtkCommonDataModel import ( vtkBoundingBox, vtkCell, vtkCellArray, vtkPointSet, vtkPolyData, vtkStaticCellLocator, vtkStaticPointLocator, vtkUnstructuredGrid, VTK_POLYHEDRON ) @@ -14,7 +15,7 @@ from vtkmodules.vtkFiltersModeling import vtkCollisionDetectionFilter, vtkLinearExtrusionFilter from geos.mesh.doctor.actions import reorient_mesh, triangle_distance from geos.mesh.utils.genericHelpers import vtk_iter -from geos.mesh.io.vtkIO import read_mesh +from geos.mesh.io.vtkIO import read_unstructured_grid @dataclass( frozen=True ) @@ -26,12 +27,11 @@ class Options: @dataclass( frozen=True ) class Result: - non_conformal_cells: List[ Tuple[ int, int ] ] + non_conformal_cells: list[ tuple[ int, int ] ] class BoundaryMesh: - """ - A BoundaryMesh is the envelope of the 3d mesh on which we want to perform the simulations. + """A BoundaryMesh is the envelope of the 3d mesh on which we want to perform the simulations. It is computed by vtk. But we want to be sure that the normals of the envelope are directed outwards. The `vtkDataSetSurfaceFilter` does not have the same behavior for standard vtk cells (like tets or hexs), and for polyhedron meshes, for which the result is a bit brittle. @@ -40,9 +40,10 @@ class BoundaryMesh: """ def __init__( self, mesh: vtkUnstructuredGrid ): - """ - Builds a boundary mesh. - :param mesh: The 3d mesh. + """Builds a boundary mesh. + + Args: + mesh (vtkUnstructuredGrid): The 3d mesh. """ # Building the boundary meshes boundary_mesh, __normals, self.__original_cells = BoundaryMesh.__build_boundary_mesh( mesh ) @@ -53,13 +54,13 @@ def __init__( self, mesh: vtkUnstructuredGrid ): self.re_boundary_mesh, re_normals, _ = BoundaryMesh.__build_boundary_mesh( reoriented_mesh, consistency=False ) num_cells = boundary_mesh.GetNumberOfCells() # Precomputing the underlying cell type - self.__is_underlying_cell_type_a_polyhedron = numpy.zeros( num_cells, dtype=bool ) + self.__is_underlying_cell_type_a_polyhedron = np.zeros( num_cells, dtype=bool ) for ic in range( num_cells ): self.__is_underlying_cell_type_a_polyhedron[ ic ] = mesh.GetCell( self.__original_cells.GetValue( ic ) ).GetCellType() == VTK_POLYHEDRON # Precomputing the normals - self.__normals: numpy.ndarray = numpy.empty( ( num_cells, 3 ), dtype=numpy.double, - order='C' ) # Do not modify the storage layout + self.__normals: np.ndarray = np.empty( ( num_cells, 3 ), dtype=np.float64, + order='C' ) # Do not modify the storage layout for ic in range( num_cells ): if self.__is_underlying_cell_type_a_polyhedron[ ic ]: self.__normals[ ic, : ] = re_normals.GetTuple3( ic ) @@ -67,13 +68,16 @@ def __init__( self, mesh: vtkUnstructuredGrid ): self.__normals[ ic, : ] = __normals.GetTuple3( ic ) @staticmethod - def __build_boundary_mesh( mesh: vtkUnstructuredGrid, consistency=True ) -> Tuple[ vtkUnstructuredGrid, Any, Any ]: - """ - From a 3d mesh, build the envelope meshes. - :param mesh: The input 3d mesh. - :param consistency: The vtk option passed to the `vtkDataSetSurfaceFilter`. - :return: A tuple containing the boundary mesh, the normal vectors array, - an array that maps the id of the boundary element to the id of the 3d cell it touches. + def __build_boundary_mesh( mesh: vtkUnstructuredGrid, consistency=True ) -> tuple[ vtkUnstructuredGrid, Any, Any ]: + """From a 3d mesh, build the envelope meshes. + + Args: + mesh (vtkUnstructuredGrid): The input 3d mesh. + consistency (bool, optional): The vtk option passed to the `vtkDataSetSurfaceFilter`. Defaults to True. + + Returns: + tuple[ vtkUnstructuredGrid, Any, Any ]: A tuple containing the boundary mesh, the normal vectors array, + an array that maps the id of the boundary element to the id of the 3d cell it touches. """ f = vtkDataSetSurfaceFilter() f.PassThroughCellIdsOn() @@ -81,7 +85,7 @@ def __build_boundary_mesh( mesh: vtkUnstructuredGrid, consistency=True ) -> Tupl f.FastModeOff() # Note that we do not need the original points, but we could keep them as well if needed - original_cells_key = "ORIGINAL_CELLS" + original_cells_key: str = "ORIGINAL_CELLS" f.SetOriginalCellIdsName( original_cells_key ) boundary_mesh = vtkPolyData() @@ -94,7 +98,7 @@ def __build_boundary_mesh( mesh: vtkUnstructuredGrid, consistency=True ) -> Tupl n.ComputeCellNormalsOn() n.SetInputData( boundary_mesh ) n.Update() - normals = n.GetOutput().GetCellData().GetArray( "Normals" ) + normals: vtkDataArray = n.GetOutput().GetCellData().GetArray( "Normals" ) assert normals assert normals.GetNumberOfComponents() == 3 assert normals.GetNumberOfTuples() == boundary_mesh.GetNumberOfCells() @@ -103,74 +107,92 @@ def __build_boundary_mesh( mesh: vtkUnstructuredGrid, consistency=True ) -> Tupl return boundary_mesh, normals, original_cells def GetNumberOfCells( self ) -> int: - """ - The number of cells. - :return: An integer. + """The number of cells. + + Returns: + int """ return self.re_boundary_mesh.GetNumberOfCells() def GetNumberOfPoints( self ) -> int: - """ - The number of points. - :return: An integer. + """The number of points. + + Returns: + int """ return self.re_boundary_mesh.GetNumberOfPoints() - def bounds( self, i ) -> Tuple[ float, float, float, float, float, float ]: - """ - The boundrary box of cell `i`. - :param i: The boundary cell index. - :return: The vtk bounding box. + def bounds( self, i: int ) -> tuple[ float, float, float, float, float, float ]: + """The boundrary box of cell `i`. + + Args: + i (int): The boundary cell index. + + Returns: + tuple[ float, float, float, float, float, float ]: The vtk bounding box. """ return self.re_boundary_mesh.GetCell( i ).GetBounds() - def normals( self, i ) -> numpy.ndarray: - """ - The normal of cell `i`. This normal will be directed outwards - :param i: The boundary cell index. - :return: The normal as a length-3 numpy array. + def normals( self, i: int ) -> npt.NDArray: + """The normal of cell `i`. This normal will be directed outwards + + Args: + i (int): The boundary cell index. + + Returns: + npt.NDArray: The normal as a length-3 numpy array. """ return self.__normals[ i ] - def GetCell( self, i ) -> vtkCell: - """ - Cell i of the boundary mesh. This cell will have its normal directed outwards. - :param i: The boundary cell index. - :return: The cell instance. - :warning: This member function relies on the vtkUnstructuredGrid.GetCell member function which is not thread safe. + def GetCell( self, i: int ) -> vtkCell: + """Cell i of the boundary mesh. This cell will have its normal directed outwards. + This member function relies on the vtkUnstructuredGrid.GetCell member function which is not thread safe. + + Args: + i (int): The boundary cell index. + + Returns: + vtkCell: The cell instance. """ return self.re_boundary_mesh.GetCell( i ) - def GetPoint( self, i ) -> Tuple[ float, float, float ]: - """ - Point i of the boundary mesh. - :param i: The boundary point index. - :return: A length-3 tuple containing the coordinates of the point. - :warning: This member function relies on the vtkUnstructuredGrid.GetPoint member function which is not thread safe. + def GetPoint( self, i: int ) -> tuple[ float, float, float ]: + """Point i of the boundary mesh. + This member function relies on the vtkUnstructuredGrid.GetPoint member function which is not thread safe. + + Args: + i (int): The boundary point index. + + Returns: + tuple[ float, float, float ]: A length-3 tuple containing the coordinates of the point. """ return self.re_boundary_mesh.GetPoint( i ) @property - def original_cells( self ): - """ - Returns the 2d boundary cell to the 3d cell index of the original mesh. - :return: A 1d array. + def original_cells( self ) -> vtkDataArray: + """Returns the 2d boundary cell to the 3d cell index of the original mesh. + + Returns: + vtkDataArray: A 1d array. """ return self.__original_cells def build_poly_data_for_extrusion( i: int, boundary_mesh: BoundaryMesh ) -> vtkPolyData: - """ - Creates a vtkPolyData containing the unique cell `i` of the boundary mesh. + """Creates a vtkPolyData containing the unique cell `i` of the boundary mesh. This operation is needed to use the vtk extrusion filter. - :param i: The boundary cell index that will eventually be extruded. - :param boundary_mesh: - :return: The created vtkPolyData. + + Args: + i (int): The boundary cell index that will eventually be extruded. + boundary_mesh (BoundaryMesh) + + Returns: + vtkPolyData: The created vtkPolyData. """ cell = boundary_mesh.GetCell( i ) copied_cell = cell.NewInstance() copied_cell.DeepCopy( cell ) - points_ids_mapping = [] + points_ids_mapping: list[ int ] = list() for i in range( copied_cell.GetNumberOfPoints() ): copied_cell.GetPointIds().SetId( i, i ) points_ids_mapping.append( cell.GetPointId( i ) ) @@ -187,12 +209,15 @@ def build_poly_data_for_extrusion( i: int, boundary_mesh: BoundaryMesh ) -> vtkP def are_points_conformal( point_tolerance: float, cell_i: vtkCell, cell_j: vtkCell ) -> bool: - """ - Checks if points of cell `i` matches, one by one, the points of cell `j`. - :param point_tolerance: The point tolerance to consider that two points match. - :param cell_i: The first cell. - :param cell_j: The second cell. - :return: A boolean. + """Checks if points of cell `i` matches, one by one, the points of cell `j`. + + Args: + point_tolerance (float): The point tolerance to consider that two points match. + cell_i (vtkCell): The first cell. + cell_j (vtkCell): The second cell. + + Returns: + bool """ # In this last step, we check that the nodes are (or not) matching each other. if cell_i.GetNumberOfPoints() != cell_j.GetNumberOfPoints(): @@ -203,34 +228,34 @@ def are_points_conformal( point_tolerance: float, cell_i: vtkCell, cell_j: vtkCe points.SetPoints( cell_i.GetPoints() ) point_locator.SetDataSet( points ) point_locator.BuildLocator() - found_points = set() + found_points: set[ int ] = set() for ip in range( cell_j.GetNumberOfPoints() ): p = cell_j.GetPoints().GetPoint( ip ) squared_dist = vtk_reference( 0. ) # unused - found_point = point_locator.FindClosestPointWithinRadius( point_tolerance, p, squared_dist ) + found_point: int = point_locator.FindClosestPointWithinRadius( point_tolerance, p, squared_dist ) found_points.add( found_point ) return found_points == set( range( cell_i.GetNumberOfPoints() ) ) class Extruder: - """ - Computes and stores all the extrusions of the boundary faces. + """Computes and stores all the extrusions of the boundary faces. The main reason for this class is to be lazy and cache the extrusions. """ def __init__( self, boundary_mesh: BoundaryMesh, face_tolerance: float ): - self.__extrusions: List[ vtkPolyData ] = [ - None, - ] * boundary_mesh.GetNumberOfCells() - self.__boundary_mesh = boundary_mesh - self.__face_tolerance = face_tolerance + self.__extrusions: list[ vtkPolyData ] = [ None ] * boundary_mesh.GetNumberOfCells() + self.__boundary_mesh: BoundaryMesh = boundary_mesh + self.__face_tolerance: float = face_tolerance - def __extrude( self, polygon_poly_data, normal ) -> vtkPolyData: - """ - Extrude the polygon data to create a volume that will be used for intersection. - :param polygon_poly_data: The data to extrude - :param normal: The (uniform) direction of the extrusion. - :return: The extrusion. + def __extrude( self, polygon_poly_data: vtkPolyData, normal: Sequence[ float ] ) -> vtkPolyData: + """Extrude the polygon data to create a volume that will be used for intersection. + + Args: + polygon_poly_data (_type_): The data to extrude + normal (_type_): The (uniform) direction of the extrusion. + + Returns: + vtkPolyData: The extrusion. """ extruder = vtkLinearExtrusionFilter() extruder.SetExtrusionTypeToVectorExtrusion() @@ -240,11 +265,14 @@ def __extrude( self, polygon_poly_data, normal ) -> vtkPolyData: extruder.Update() return extruder.GetOutput() - def __getitem__( self, i ) -> vtkPolyData: - """ - Returns the vtk extrusion for boundary element i. - :param i: The cell index. - :return: The vtk instance. + def __getitem__( self, i: int ) -> vtkPolyData: + """Returns the vtk extrusion for boundary element i. + + Args: + i (int): The cell index. + + Returns: + vtkPolyData: The vtk instance. """ extrusion = self.__extrusions[ i ] if extrusion: @@ -257,14 +285,17 @@ def __getitem__( self, i ) -> vtkPolyData: def are_faces_conformal_using_extrusions( extrusions: Extruder, i: int, j: int, boundary_mesh: vtkUnstructuredGrid, point_tolerance: float ) -> bool: - """ - Tests if two boundary faces are conformal, checking for intersection between their normal extruded volumes. - :param extrusions: The extrusions cache. - :param i: The cell index of the first cell. - :param j: The cell index of the second cell. - :param boundary_mesh: The boundary mesh. - :param point_tolerance: The point tolerance to consider that two points match. - :return: A boolean. + """Tests if two boundary faces are conformal, checking for intersection between their normal extruded volumes. + + Args: + extrusions (Extruder): The extrusions cache. + i (int): The cell index of the first cell. + j (int): The cell index of the second cell. + boundary_mesh (vtkUnstructuredGrid): The boundary mesh. + point_tolerance (float): The point tolerance to consider that two points match. + + Returns: + bool """ collision = vtkCollisionDetectionFilter() collision.SetCollisionModeToFirstContact() @@ -289,21 +320,24 @@ def are_faces_conformal_using_extrusions( extrusions: Extruder, i: int, j: int, def are_faces_conformal_using_distances( i: int, j: int, boundary_mesh: vtkUnstructuredGrid, face_tolerance: float, point_tolerance: float ) -> bool: - """ - Tests if two boundary faces are conformal, checking the minimal distance between triangulated surfaces. - :param i: The cell index of the first cell. - :param j: The cell index of the second cell. - :param boundary_mesh: The boundary mesh. - :param face_tolerance: The tolerance under which we should consider the two faces "touching" each other. - :param point_tolerance: The point tolerance to consider that two points match. - :return: A boolean. + """Tests if two boundary faces are conformal, checking the minimal distance between triangulated surfaces. + + Args: + i (int): The cell index of the first cell. + j (int): The cell index of the second cell. + boundary_mesh (vtkUnstructuredGrid): The boundary mesh. + face_tolerance (float): The tolerance under which we should consider the two faces "touching" each other. + point_tolerance (float): The point tolerance to consider that two points match. + + Returns: + bool """ cp_i = boundary_mesh.GetCell( i ).NewInstance() cp_i.DeepCopy( boundary_mesh.GetCell( i ) ) cp_j = boundary_mesh.GetCell( j ).NewInstance() cp_j.DeepCopy( boundary_mesh.GetCell( j ) ) - def triangulate( cell ): + def triangulate( cell: vtkCell ): assert cell.GetCellDimension() == 2 __points_ids = vtkIdList() __points = vtkPoints() @@ -322,13 +356,13 @@ def build_numpy_triangles( points_ids ): __t = [] for __pi in points_ids[ __i:__i + 3 ]: __t.append( boundary_mesh.GetPoint( __pi ) ) - __triangles.append( numpy.array( __t, dtype=float ) ) + __triangles.append( np.array( __t, dtype=float ) ) return __triangles triangles_i = build_numpy_triangles( points_ids_i ) triangles_j = build_numpy_triangles( points_ids_j ) - min_dist = numpy.inf + min_dist = np.inf for ti, tj in [ ( ti, tj ) for ti in triangles_i for tj in triangles_j ]: # Note that here, we compute the exact distance to compare with the threshold. # We could improve by exiting the iterative distance computation as soon as @@ -344,44 +378,57 @@ def build_numpy_triangles( points_ids ): return are_points_conformal( point_tolerance, cp_i, cp_j ) -def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: - """ - Checks if the mesh is "conformal" (i.e. if some of its boundary faces may not be too close to each other without matching nodes). - :param mesh: The vtk mesh - :param options: The check options. - :return: The Result instance. - """ - boundary_mesh = BoundaryMesh( mesh ) - cos_theta = abs( math.cos( numpy.deg2rad( options.angle_tolerance ) ) ) - num_cells = boundary_mesh.GetNumberOfCells() +def compute_bounding_box( boundary_mesh: BoundaryMesh, face_tolerance: float ) -> npt.NDArray[ np.float64 ]: + # Precomputing the bounding boxes. + # The options are important to directly interact with memory in C++. + bounding_boxes = np.empty( ( boundary_mesh.GetNumberOfCells(), 6 ), dtype=np.float64, order="C" ) + for i in range( boundary_mesh.GetNumberOfCells() ): + bb = vtkBoundingBox( boundary_mesh.bounds( i ) ) + bb.Inflate( 2 * face_tolerance ) + assert bounding_boxes[ + i, : ].data.contiguous # Do not modify the storage layout since vtk deals with raw memory here. + bb.GetBounds( bounding_boxes[ i, : ] ) + return bounding_boxes + +def compute_number_cells_per_node( boundary_mesh: BoundaryMesh ) -> npt.NDArray[ np.int64 ]: # Computing the exact number of cells per node - num_cells_per_node = numpy.zeros( boundary_mesh.GetNumberOfPoints(), dtype=int ) + num_cells_per_node = np.zeros( boundary_mesh.GetNumberOfPoints(), dtype=int ) for ic in range( boundary_mesh.GetNumberOfCells() ): c = boundary_mesh.GetCell( ic ) point_ids = c.GetPointIds() for point_id in vtk_iter( point_ids ): num_cells_per_node[ point_id ] += 1 + return num_cells_per_node + +def build_cell_locator( mesh: vtkUnstructuredGrid, numberMaxCellPerNode: int ) -> vtkStaticCellLocator: cell_locator = vtkStaticCellLocator() cell_locator.Initialize() - cell_locator.SetNumberOfCellsPerNode( num_cells_per_node.max() ) - cell_locator.SetDataSet( boundary_mesh.re_boundary_mesh ) + cell_locator.SetNumberOfCellsPerNode( numberMaxCellPerNode ) + cell_locator.SetDataSet( mesh ) cell_locator.BuildLocator() + return cell_locator - # Precomputing the bounding boxes. - # The options are important to directly interact with memory in C++. - bounding_boxes = numpy.empty( ( boundary_mesh.GetNumberOfCells(), 6 ), dtype=numpy.double, order="C" ) - for i in range( boundary_mesh.GetNumberOfCells() ): - bb = vtkBoundingBox( boundary_mesh.bounds( i ) ) - bb.Inflate( 2 * options.face_tolerance ) - assert bounding_boxes[ - i, : ].data.contiguous # Do not modify the storage layout since vtk deals with raw memory here. - bb.GetBounds( bounding_boxes[ i, : ] ) - non_conformal_cells = [] +def find_non_conformal_cells( mesh: vtkUnstructuredGrid, options: Options ) -> list[ tuple[ int, int ] ]: + # Extracts the outer surface of the 3D mesh. + # Ensures that face normals are consistently oriented outward. + boundary_mesh = BoundaryMesh( mesh ) + num_cells: int = boundary_mesh.GetNumberOfCells() + + # Used to filter out face pairs that are not facing each other. + cos_theta: float = abs( math.cos( np.deg2rad( options.angle_tolerance ) ) ) + + # Prepares extruded volumes of boundary faces for intersection testing. extrusions = Extruder( boundary_mesh, options.face_tolerance ) + + num_cells_per_node = compute_number_cells_per_node( boundary_mesh ) + bounding_boxes = compute_bounding_box( boundary_mesh, options.face_tolerance ) + cell_locator = build_cell_locator( boundary_mesh.re_boundary_mesh, num_cells_per_node.max() ) + close_cells = vtkIdList() + non_conformal_cells_boundary_id: list[ tuple[ int, int ] ] = list() # Looping on all the pairs of boundary cells. We'll hopefully discard most of the pairs. for i in tqdm( range( num_cells ), desc="Non conformal elements" ): cell_locator.FindCellsWithinBounds( bounding_boxes[ i ], close_cells ) @@ -390,19 +437,34 @@ def __action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: continue # Discarding pairs that are not facing each others (with a threshold). normal_i, normal_j = boundary_mesh.normals( i ), boundary_mesh.normals( j ) - if numpy.dot( normal_i, normal_j ) > -cos_theta: # opposite directions only (can be facing or not) + if np.dot( normal_i, normal_j ) > -cos_theta: # opposite directions only (can be facing or not) continue # At this point, back-to-back and face-to-face pairs of elements are considered. if not are_faces_conformal_using_extrusions( extrusions, i, j, boundary_mesh, options.point_tolerance ): - non_conformal_cells.append( ( i, j ) ) + non_conformal_cells_boundary_id.append( ( i, j ) ) # Extracting the original 3d element index (and not the index of the boundary mesh). - tmp = [] - for i, j in non_conformal_cells: - tmp.append( ( boundary_mesh.original_cells.GetValue( i ), boundary_mesh.original_cells.GetValue( j ) ) ) + non_conformal_cells: list[ tuple[ int, int ] ] = list() + for i, j in non_conformal_cells_boundary_id: + non_conformal_cells.append( + ( boundary_mesh.original_cells.GetValue( i ), boundary_mesh.original_cells.GetValue( j ) ) ) + return non_conformal_cells - return Result( non_conformal_cells=tmp ) + +def mesh_action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + """Checks if the mesh is "conformal" + (i.e. if some of its boundary faces may not be too close to each other without matching nodes). + + Args: + mesh (vtkUnstructuredGrid): The vtk mesh + options (Options): The check options. + + Returns: + Result: The Result instance. + """ + non_conformal_cells = find_non_conformal_cells( mesh, options ) + return Result( non_conformal_cells=non_conformal_cells ) def action( vtk_input_file: str, options: Options ) -> Result: - mesh = read_mesh( vtk_input_file ) - return __action( mesh, options ) + mesh: vtkUnstructuredGrid = read_unstructured_grid( vtk_input_file ) + return mesh_action( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/self_intersecting_elements.py b/geos-mesh/src/geos/mesh/doctor/actions/self_intersecting_elements.py index 3b7d313a..fb796990 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/self_intersecting_elements.py +++ b/geos-mesh/src/geos/mesh/doctor/actions/self_intersecting_elements.py @@ -1,9 +1,10 @@ from dataclasses import dataclass -from typing import Collection, List +from typing import Collection from vtkmodules.util.numpy_support import vtk_to_numpy from vtkmodules.vtkFiltersGeneral import vtkCellValidator from vtkmodules.vtkCommonCore import vtkOutputWindow, vtkFileOutputWindow -from geos.mesh.io.vtkIO import read_mesh +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from geos.mesh.io.vtkIO import read_unstructured_grid @dataclass( frozen=True ) @@ -18,62 +19,85 @@ class Result: intersecting_faces_elements: Collection[ int ] non_contiguous_edges_elements: Collection[ int ] non_convex_elements: Collection[ int ] - faces_are_oriented_incorrectly_elements: Collection[ int ] + faces_oriented_incorrectly_elements: Collection[ int ] -def __action( mesh, options: Options ) -> Result: +def get_invalid_cell_ids( mesh: vtkUnstructuredGrid, min_distance: float ) -> dict[ str, list[ int ] ]: + """For every cell element in a vtk mesh, check if the cell is invalid regarding 6 specific criteria: + "wrong_number_of_points", "intersecting_edges", "intersecting_faces", + "non_contiguous_edges","non_convex" and "faces_oriented_incorrectly". + + If any of this criteria was met, the cell index is added to a list corresponding to this specific criteria. + The dict with the complete list of cell indices by criteria is returned. + + Args: + mesh (vtkUnstructuredGrid): A vtk grid. + min_distance (float): Minimum distance in the computation. + + Returns: + dict[ str, list[ int ] ]: + { + "wrong_number_of_points": [ 10, 34, ... ], + "intersecting_edges": [ ... ], + "intersecting_faces": [ ... ], + "non_contiguous_edges": [ ... ], + "non_convex": [ ... ], + "faces_oriented_incorrectly": [ ... ] + } + """ + # The goal of this first block is to silence the standard error output from VTK. The vtkCellValidator can be very + # verbose, printing a message for every cell it checks. We redirect the output to /dev/null to remove it. err_out = vtkFileOutputWindow() - err_out.SetFileName( "/dev/null" ) # vtkCellValidator outputs loads for each cell... + err_out.SetFileName( "/dev/null" ) vtk_std_err_out = vtkOutputWindow() vtk_std_err_out.SetInstance( err_out ) - valid = 0x0 - wrong_number_of_points = 0x01 - intersecting_edges = 0x02 - intersecting_faces = 0x04 - non_contiguous_edges = 0x08 - non_convex = 0x10 - faces_are_oriented_incorrectly = 0x20 - - wrong_number_of_points_elements: List[ int ] = [] - intersecting_edges_elements: List[ int ] = [] - intersecting_faces_elements: List[ int ] = [] - non_contiguous_edges_elements: List[ int ] = [] - non_convex_elements: List[ int ] = [] - faces_are_oriented_incorrectly_elements: List[ int ] = [] + # Different types of cell invalidity are defined as hexadecimal values, specific to vtkCellValidator + # Here NonPlanarFaces and DegenerateFaces can also be obtained. + error_masks: dict[ str, int ] = { + "wrong_number_of_points_elements": 0x01, # 0000 0001 + "intersecting_edges_elements": 0x02, # 0000 0010 + "intersecting_faces_elements": 0x04, # 0000 0100 + "non_contiguous_edges_elements": 0x08, # 0000 1000 + "non_convex_elements": 0x10, # 0001 0000 + "faces_oriented_incorrectly_elements": 0x20, # 0010 0000 + } - f = vtkCellValidator() - f.SetTolerance( options.min_distance ) + # The results can be stored in a dictionary where keys are the error names + # and values are the lists of cell indices with that error. + # We can initialize it directly from the keys of our error_masks dictionary. + invalid_cell_ids: dict[ str, list[ int ] ] = { error_name: list() for error_name in error_masks } + f = vtkCellValidator() + f.SetTolerance( min_distance ) f.SetInputData( mesh ) - f.Update() + f.Update() # executes the filter output = f.GetOutput() - validity = output.GetCellData().GetArray( "ValidityState" ) # Could not change name using the vtk interface. + validity = output.GetCellData().GetArray( "ValidityState" ) assert validity is not None + # array of np.int16 that combines the flags using a bitwise OR operation for each cell index. validity = vtk_to_numpy( validity ) - for i, v in enumerate( validity ): - if not v & valid: - if v & wrong_number_of_points: - wrong_number_of_points_elements.append( i ) - if v & intersecting_edges: - intersecting_edges_elements.append( i ) - if v & intersecting_faces: - intersecting_faces_elements.append( i ) - if v & non_contiguous_edges: - non_contiguous_edges_elements.append( i ) - if v & non_convex: - non_convex_elements.append( i ) - if v & faces_are_oriented_incorrectly: - faces_are_oriented_incorrectly_elements.append( i ) - return Result( wrong_number_of_points_elements=wrong_number_of_points_elements, - intersecting_edges_elements=intersecting_edges_elements, - intersecting_faces_elements=intersecting_faces_elements, - non_contiguous_edges_elements=non_contiguous_edges_elements, - non_convex_elements=non_convex_elements, - faces_are_oriented_incorrectly_elements=faces_are_oriented_incorrectly_elements ) + for cell_index, validity_flag in enumerate( validity ): + if not validity_flag: # Skip valid cells ( validity_flag == 0 or 0000 0000 ) + continue + for error_name, error_mask in error_masks.items(): # Check only invalid cells against all possible errors. + if validity_flag & error_mask: + invalid_cell_ids[ error_name ].append( cell_index ) + + return invalid_cell_ids + + +def mesh_action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + invalid_cell_ids: dict[ str, list[ int ] ] = get_invalid_cell_ids( mesh, options.min_distance ) + return Result( wrong_number_of_points_elements=invalid_cell_ids[ "wrong_number_of_points_elements" ], + intersecting_edges_elements=invalid_cell_ids[ "intersecting_edges_elements" ], + intersecting_faces_elements=invalid_cell_ids[ "intersecting_faces_elements" ], + non_contiguous_edges_elements=invalid_cell_ids[ "non_contiguous_edges_elements" ], + non_convex_elements=invalid_cell_ids[ "non_convex_elements" ], + faces_oriented_incorrectly_elements=invalid_cell_ids[ "faces_oriented_incorrectly_elements" ] ) def action( vtk_input_file: str, options: Options ) -> Result: - mesh = read_mesh( vtk_input_file ) - return __action( mesh, options ) + mesh: vtkUnstructuredGrid = read_unstructured_grid( vtk_input_file ) + return mesh_action( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/actions/supported_elements.py b/geos-mesh/src/geos/mesh/doctor/actions/supported_elements.py index 8d9fd46a..f19680aa 100644 --- a/geos-mesh/src/geos/mesh/doctor/actions/supported_elements.py +++ b/geos-mesh/src/geos/mesh/doctor/actions/supported_elements.py @@ -1,8 +1,9 @@ from dataclasses import dataclass import multiprocessing import networkx +from numpy import ones from tqdm import tqdm -from typing import FrozenSet, Iterable, Mapping, Optional +from typing import Iterable, Mapping, Optional from vtkmodules.util.numpy_support import vtk_to_numpy from vtkmodules.vtkCommonCore import vtkIdList from vtkmodules.vtkCommonDataModel import ( vtkCellTypes, vtkUnstructuredGrid, VTK_HEXAGONAL_PRISM, VTK_HEXAHEDRON, @@ -10,7 +11,7 @@ VTK_WEDGE ) from geos.mesh.doctor.actions.vtk_polyhedron import build_face_to_face_connectivity_through_edges, FaceStream from geos.mesh.doctor.parsing.cli_parsing import setup_logger -from geos.mesh.io.vtkIO import read_mesh +from geos.mesh.io.vtkIO import read_unstructured_grid from geos.mesh.utils.genericHelpers import vtk_iter @@ -22,32 +23,26 @@ class Options: @dataclass( frozen=True ) class Result: - unsupported_std_elements_types: FrozenSet[ int ] # list of unsupported types - unsupported_polyhedron_elements: FrozenSet[ + unsupported_std_elements_types: list[ str ] # list of unsupported types + unsupported_polyhedron_elements: frozenset[ int ] # list of polyhedron elements that could not be converted to supported std elements # for multiprocessing, vtkUnstructuredGrid cannot be pickled. Let's use a global variable instead. +# Global variable to be set in each worker process MESH: Optional[ vtkUnstructuredGrid ] = None -def init_worker_mesh( input_file_for_worker: str ): - """Initializer for multiprocessing.Pool to set the global MESH variable in each worker process. - - Args: - input_file_for_worker (str): Filepath to vtk grid - """ +def init_worker( mesh_to_init: vtkUnstructuredGrid ) -> None: + """Initializer for each worker process to set the global mesh.""" global MESH - setup_logger.debug( - f"Worker process (PID: {multiprocessing.current_process().pid}) initializing MESH from file: {input_file_for_worker}" - ) - MESH = read_mesh( input_file_for_worker ) - if MESH is None: - setup_logger.error( - f"Worker process (PID: {multiprocessing.current_process().pid}) failed to load mesh from {input_file_for_worker}" - ) - # You might want to raise an error here or ensure MESH being None is handled downstream - # For now, the assert MESH is not None in __call__ will catch this. + MESH = mesh_to_init + + +supported_cell_types: set[ int ] = { + VTK_HEXAGONAL_PRISM, VTK_HEXAHEDRON, VTK_PENTAGONAL_PRISM, VTK_POLYHEDRON, VTK_PYRAMID, VTK_TETRA, VTK_VOXEL, + VTK_WEDGE +} class IsPolyhedronConvertible: @@ -55,14 +50,11 @@ class IsPolyhedronConvertible: def __init__( self ): def build_prism_graph( n: int, name: str ) -> networkx.Graph: - """Builds the face to face connectivities (through edges) for prism graphs. - - Args: - n (int): The number of nodes of the basis (i.e. the pentagonal prims gets n = 5) - name (str): A human-readable name for logging purpose. - - Returns: - networkx.Graph: A graph instance. + """ + Builds the face to face connectivities (through edges) for prism graphs. + :param n: The number of nodes of the basis (i.e. the pentagonal prims gets n = 5) + :param name: A human-readable name for logging purpose. + :return: A graph instance. """ tmp = networkx.cycle_graph( n ) for node in range( n ): @@ -90,34 +82,26 @@ def build_prism_graph( n: int, name: str ) -> networkx.Graph: } def __is_polyhedron_supported( self, face_stream ) -> str: - """Checks if a polyhedron can be converted into a supported cell. + """ + Checks if a polyhedron can be converted into a supported cell. If so, returns the name of the type. If not, the returned name will be empty. - - Args: - face_stream (_type_): The polyhedron. - - Returns: - str: The name of the supported type or an empty string. + :param face_stream: The polyhedron. + :return: The name of the supported type or an empty string. """ cell_graph = build_face_to_face_connectivity_through_edges( face_stream, add_compatibility=True ) - if cell_graph.order() not in self.__reference_graphs: - return "" for reference_graph in self.__reference_graphs[ cell_graph.order() ]: if networkx.is_isomorphic( reference_graph, cell_graph ): return str( reference_graph.name ) return "" def __call__( self, ic: int ) -> int: - """Checks if a vtk polyhedron cell can be converted into a supported GEOSX element. - - Args: - ic (int): The index element. - - Returns: - int: -1 if the polyhedron vtk element can be converted into a supported element type. The index otherwise. + """ + Checks if a vtk polyhedron cell can be converted into a supported GEOSX element. + :param ic: The index element. + :return: -1 if the polyhedron vtk element can be converted into a supported element type. The index otherwise. """ global MESH - assert MESH is not None, f"MESH global variable not initialized in worker process (PID: {multiprocessing.current_process().pid}). This should have been set by init_worker_mesh." + assert MESH is not None if MESH.GetCellType( ic ) != VTK_POLYHEDRON: return -1 pt_ids = vtkIdList() @@ -128,50 +112,45 @@ def __call__( self, ic: int ) -> int: setup_logger.debug( f"Polyhedron cell {ic} can be converted into \"{converted_type_name}\"" ) return -1 else: - setup_logger.debug( - f"Polyhedron cell {ic} (in PID {multiprocessing.current_process().pid}) cannot be converted into any supported element." - ) + setup_logger.debug( f"Polyhedron cell {ic} cannot be converted into any supported element." ) return ic -def __action( vtk_input_file: str, options: Options ) -> Result: - # Main process loads the mesh for its own use - mesh = read_mesh( vtk_input_file ) - if mesh is None: - setup_logger.error( f"Main process failed to load mesh from {vtk_input_file}. Aborting." ) - # Return an empty/error result or raise an exception - return Result( unsupported_std_elements_types=frozenset(), unsupported_polyhedron_elements=frozenset() ) - - if hasattr( mesh, "GetDistinctCellTypesArray" ): - cell_types_numpy = vtk_to_numpy( mesh.GetDistinctCellTypesArray() ) - cell_types = set( cell_types_numpy.tolist() ) +def find_unsupported_std_elements_types( mesh: vtkUnstructuredGrid ) -> list[ str ]: + if hasattr( mesh, "GetDistinctCellTypesArray" ): # For more recent versions of vtk. + unique_cell_types = set( vtk_to_numpy( mesh.GetDistinctCellTypesArray() ) ) else: - vtk_cell_types_obj = vtkCellTypes() - mesh.GetCellTypes( vtk_cell_types_obj ) - cell_types = set( vtk_iter( vtk_cell_types_obj ) ) + vtk_cell_types = vtkCellTypes() + mesh.GetCellTypes( vtk_cell_types ) + unique_cell_types = set( vtk_iter( vtk_cell_types ) ) + result_values: set[ int ] = unique_cell_types - supported_cell_types + results = [ f"Type {i}: {vtkCellTypes.GetClassNameFromTypeId( i )}" for i in frozenset( result_values ) ] + return results - supported_cell_types = { - VTK_HEXAGONAL_PRISM, VTK_HEXAHEDRON, VTK_PENTAGONAL_PRISM, VTK_POLYHEDRON, VTK_PYRAMID, VTK_TETRA, VTK_VOXEL, - VTK_WEDGE - } - unsupported_std_elements_types = cell_types - supported_cell_types +def find_unsupported_polyhedron_elements( mesh: vtkUnstructuredGrid, options: Options ) -> list[ int ]: # Dealing with polyhedron elements. - num_cells = mesh.GetNumberOfCells() - polyhedron_converter = IsPolyhedronConvertible() + num_cells: int = mesh.GetNumberOfCells() + result = ones( num_cells, dtype=int ) * -1 + # Use the initializer to set up each worker process + # Pass the mesh to the initializer + with multiprocessing.Pool( processes=options.nproc, initializer=init_worker, initargs=( mesh, ) ) as pool: + # Pass a mesh-free instance of the class to the workers. + # The MESH global will already be set in each worker. + generator = pool.imap_unordered( IsPolyhedronConvertible(), range( num_cells ), chunksize=options.chunk_size ) + for i, val in enumerate( tqdm( generator, total=num_cells, desc="Testing support for elements" ) ): + result[ i ] = val + + return [ i for i in result if i > -1 ] - unsupported_polyhedron_indices = [] - # Pass the vtk_input_file to the initializer - with multiprocessing.Pool( processes=options.nproc, initializer=init_worker_mesh, - initargs=( vtk_input_file, ) ) as pool: # Comma makes it a tuple - generator = pool.imap_unordered( polyhedron_converter, range( num_cells ), chunksize=options.chunk_size ) - for cell_index_or_neg_one in tqdm( generator, total=num_cells, desc="Testing support for elements" ): - if cell_index_or_neg_one != -1: - unsupported_polyhedron_indices.append( cell_index_or_neg_one ) - return Result( unsupported_std_elements_types=frozenset( unsupported_std_elements_types ), - unsupported_polyhedron_elements=frozenset( unsupported_polyhedron_indices ) ) +def mesh_action( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + unsupported_std_elements_types: list[ str ] = find_unsupported_std_elements_types( mesh ) + unsupported_polyhedron_elements: list[ int ] = find_unsupported_polyhedron_elements( mesh, options ) + return Result( unsupported_std_elements_types=unsupported_std_elements_types, + unsupported_polyhedron_elements=frozenset( unsupported_polyhedron_elements ) ) def action( vtk_input_file: str, options: Options ) -> Result: - return __action( vtk_input_file, options ) + mesh: vtkUnstructuredGrid = read_unstructured_grid( vtk_input_file ) + return mesh_action( mesh, options ) diff --git a/geos-mesh/src/geos/mesh/doctor/filters/Checks.py b/geos-mesh/src/geos/mesh/doctor/filters/Checks.py new file mode 100644 index 00000000..e972aa77 --- /dev/null +++ b/geos-mesh/src/geos/mesh/doctor/filters/Checks.py @@ -0,0 +1,231 @@ +from types import SimpleNamespace +from typing_extensions import Self +from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from geos.mesh.doctor.actions.all_checks import Options, get_check_results +from geos.mesh.doctor.filters.MeshDoctorBase import MeshDoctorBase +from geos.mesh.doctor.parsing._shared_checks_parsing_logic import CheckFeature, display_results +from geos.mesh.doctor.parsing.all_checks_parsing import ( CHECK_FEATURES_CONFIG as cfc_all_checks, ORDERED_CHECK_NAMES + as ocn_all_checks ) +from geos.mesh.doctor.parsing.main_checks_parsing import ( CHECK_FEATURES_CONFIG as cfc_main_checks, ORDERED_CHECK_NAMES + as ocn_main_checks ) + +__doc__ = """ +Checks module is a vtk filter that performs comprehensive mesh validation checks on a vtkUnstructuredGrid. +This module contains AllChecks and MainChecks filters that run various quality checks including element validation, +node validation, topology checks, and geometric integrity verification. + +One filter input is vtkUnstructuredGrid, one filter output which is vtkUnstructuredGrid. + +To use the AllChecks filter: + +.. code-block:: python + + from filters.Checks import AllChecks + + # instantiate the filter for all available checks + allChecksFilter: AllChecks = AllChecks() + + # set input mesh + allChecksFilter.SetInputData(mesh) + + # optionally customize check parameters + allChecksFilter.setCheckParameter("collocated_nodes", "tolerance", 1e-6) + allChecksFilter.setGlobalParameter("tolerance", 1e-6) # applies to all checks with tolerance parameter + + # execute the checks + output_mesh: vtkUnstructuredGrid = allChecksFilter.getGrid() + + # get check results + check_results = allChecksFilter.getCheckResults() + +To use the MainChecks filter (subset of most important checks): + +.. code-block:: python + + from filters.Checks import MainChecks + + # instantiate the filter for main checks only + mainChecksFilter: MainChecks = MainChecks() + + # set input mesh and run checks + mainChecksFilter.SetInputData(mesh) + output_mesh: vtkUnstructuredGrid = mainChecksFilter.getGrid() +""" + + +class MeshDoctorChecks( MeshDoctorBase ): + + def __init__( self: Self, checks_to_perform: list[ str ], check_features_config: dict[ str, CheckFeature ], + ordered_check_names: list[ str ] ) -> None: + super().__init__() + self.m_checks_to_perform: list[ str ] = checks_to_perform + self.m_check_parameters: dict[ str, dict[ str, any ] ] = dict() # Custom parameters override + self.m_check_results: dict[ str, any ] = dict() + self.m_CHECK_FEATURES_CONFIG: dict[ str, CheckFeature ] = check_features_config + self.m_ORDERED_CHECK_NAMES: list[ str ] = ordered_check_names + + def RequestData( self: Self, request: vtkInformation, inInfoVec: list[ vtkInformationVector ], + outInfo: vtkInformationVector ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestData. + + Args: + request (vtkInformation): request + inInfoVec (list[vtkInformationVector]): input objects + outInfoVec (vtkInformationVector): output objects + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + input_mesh: vtkUnstructuredGrid = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] ) + output = vtkUnstructuredGrid.GetData( outInfo ) + + # Build the options using the parsing logic structure + options = self._buildOptions() + self.m_check_results = get_check_results( input_mesh, options ) + + results_wrapper = SimpleNamespace( check_results=self.m_check_results ) + display_results( options, results_wrapper ) + + output_mesh: vtkUnstructuredGrid = self.copyInputToOutput( input_mesh ) + output.ShallowCopy( output_mesh ) + + return 1 + + def _buildOptions( self: Self ) -> Options: + """Build Options object using the same logic as the parsing system. + + Returns: + Options: Properly configured options for all checks + """ + # Start with default parameters for all configured checks + default_params: dict[ str, dict[ str, any ] ] = { + name: feature.default_params.copy() + for name, feature in self.m_CHECK_FEATURES_CONFIG.items() + } + final_check_params: dict[ str, dict[ str, any ] ] = { + name: default_params[ name ] + for name in self.m_checks_to_perform + } + + # Apply any custom parameter overrides + for check_name in self.m_checks_to_perform: + if check_name in self.m_check_parameters: + final_check_params[ check_name ].update( self.m_check_parameters[ check_name ] ) + + # Instantiate Options objects for the selected checks + individual_check_options: dict[ str, any ] = dict() + individual_check_display: dict[ str, any ] = dict() + + for check_name in self.m_checks_to_perform: + if check_name not in self.m_CHECK_FEATURES_CONFIG: + self.m_logger.warning( f"Check '{check_name}' is not available. Skipping." ) + continue + + params = final_check_params[ check_name ] + feature_config = self.m_CHECK_FEATURES_CONFIG[ check_name ] + try: + individual_check_options[ check_name ] = feature_config.options_cls( **params ) + individual_check_display[ check_name ] = feature_config.display + except Exception as e: + self.m_logger.error( f"Failed to create options for check '{check_name}': {e}. " + f"This check will be skipped." ) + + return Options( checks_to_perform=list( individual_check_options.keys() ), + checks_options=individual_check_options, + check_displays=individual_check_display ) + + def getAvailableChecks( self: Self ) -> list[ str ]: + """Returns the list of available check names. + + Returns: + list[str]: List of available check names + """ + return self.m_ORDERED_CHECK_NAMES + + def getCheckResults( self: Self ) -> dict[ str, any ]: + """Returns the results of all performed checks. + + Args: + self (Self) + + Returns: + dict[str, any]: Dictionary mapping check names to their results + """ + return self.m_check_results + + def getDefaultParameters( self: Self, check_name: str ) -> dict[ str, any ]: + """Get the default parameters for a specific check. + + Args: + check_name (str): Name of the check + + Returns: + dict[str, any]: Dictionary of default parameters + """ + if check_name in self.m_CHECK_FEATURES_CONFIG: + return self.m_CHECK_FEATURES_CONFIG[ check_name ].default_params + return {} + + def setChecksToPerform( self: Self, checks_to_perform: list[ str ] ) -> None: + """Set which checks to perform. + + Args: + self (Self) + checks_to_perform (list[str]): List of check names to perform. + """ + self.m_checks_to_perform = checks_to_perform + self.Modified() + + def setCheckParameter( self: Self, check_name: str, parameter_name: str, value: any ) -> None: + """Set a parameter for a specific check. + + Args: + self (Self) + check_name (str): Name of the check (e.g., "collocated_nodes") + parameter_name (str): Name of the parameter (e.g., "tolerance") + value (any): Value to set for the parameter + """ + if check_name not in self.m_check_parameters: + self.m_check_parameters[ check_name ] = {} + self.m_check_parameters[ check_name ][ parameter_name ] = value + self.Modified() + + def setAllChecksParameter( self: Self, parameter_name: str, value: any ) -> None: + """Set a parameter for all checks that support it. + + Args: + self (Self) + parameter_name (str): Name of the parameter (e.g., "tolerance") + value (any): Value to set for the parameter + """ + for check_name in self.m_checks_to_perform: + if check_name in self.m_CHECK_FEATURES_CONFIG: + default_params = self.m_CHECK_FEATURES_CONFIG[ check_name ].default_params + if parameter_name in default_params: + self.setCheckParameter( check_name, parameter_name, value ) + self.Modified() + + +class AllChecks( MeshDoctorChecks ): + + def __init__( self: Self ) -> None: + """Vtk filter to ... of a vtkUnstructuredGrid. + + Output mesh is vtkUnstructuredGrid. + """ + super().__init__( checks_to_perform=ocn_all_checks, + check_features_config=cfc_all_checks, + ordered_check_names=ocn_all_checks ) + + +class MainChecks( MeshDoctorChecks ): + + def __init__( self: Self ) -> None: + """Vtk filter to ... of a vtkUnstructuredGrid. + + Output mesh is vtkUnstructuredGrid. + """ + super().__init__( checks_to_perform=ocn_main_checks, + check_features_config=cfc_main_checks, + ordered_check_names=ocn_main_checks ) diff --git a/geos-mesh/src/geos/mesh/doctor/filters/CollocatedNodes.py b/geos-mesh/src/geos/mesh/doctor/filters/CollocatedNodes.py new file mode 100644 index 00000000..a9c30dab --- /dev/null +++ b/geos-mesh/src/geos/mesh/doctor/filters/CollocatedNodes.py @@ -0,0 +1,146 @@ +import numpy as np +import numpy.typing as npt +from typing_extensions import Self +from vtkmodules.util.numpy_support import numpy_to_vtk +from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector, vtkDataArray +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from geos.mesh.doctor.actions.collocated_nodes import find_collocated_nodes_buckets, find_wrong_support_elements +from geos.mesh.doctor.filters.MeshDoctorBase import MeshDoctorBase + +__doc__ = """ +CollocatedNodes module is a vtk filter that identifies and handles duplicated/collocated nodes in a vtkUnstructuredGrid. +The filter can detect nodes that are within a specified tolerance distance and optionally identify elements +that have support nodes appearing more than once (wrong support elements). + +One filter input is vtkUnstructuredGrid, one filter output which is vtkUnstructuredGrid. + +To use the filter: + +.. code-block:: python + + from filters.CollocatedNodes import CollocatedNodes + + # instantiate the filter + collocatedNodesFilter: CollocatedNodes = CollocatedNodes() + + # set the tolerance for detecting collocated nodes + collocatedNodesFilter.setTolerance(1e-6) + + # optionally enable painting of wrong support elements + collocatedNodesFilter.setPaintWrongSupportElements(1) # 1 to enable, 0 to disable + + # set input mesh + collocatedNodesFilter.SetInputData(mesh) + + # execute the filter + output_mesh: vtkUnstructuredGrid = collocatedNodesFilter.getGrid() + + # get results + collocated_buckets = collocatedNodesFilter.getCollocatedNodeBuckets() # list of tuples with collocated node indices + wrong_support_elements = collocatedNodesFilter.getWrongSupportElements() # list of problematic element indices + + # write the output mesh + collocatedNodesFilter.writeGrid("output/mesh_with_collocated_info.vtu") +""" + + +class CollocatedNodes( MeshDoctorBase ): + + def __init__( self: Self ) -> None: + """Vtk filter to find the duplicated nodes of a vtkUnstructuredGrid. + + Output mesh is vtkUnstructuredGrid. + """ + super().__init__( nInputPorts=1, + nOutputPorts=1, + inputType='vtkUnstructuredGrid', + outputType='vtkUnstructuredGrid' ) + self.m_collocatedNodesBuckets: list[ tuple[ int ] ] = list() + self.m_paintWrongSupportElements: int = 0 + self.m_tolerance: float = 0.0 + self.m_wrongSupportElements: list[ int ] = list() + + def RequestData( self: Self, request: vtkInformation, inInfoVec: list[ vtkInformationVector ], + outInfo: vtkInformationVector ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestData. + + Args: + request (vtkInformation): request + inInfoVec (list[vtkInformationVector]): input objects + outInfoVec (vtkInformationVector): output objects + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + input_mesh: vtkUnstructuredGrid = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] ) + output = vtkUnstructuredGrid.GetData( outInfo ) + + self.m_collocatedNodesBuckets = find_collocated_nodes_buckets( input_mesh, self.m_tolerance ) + self.m_wrongSupportElements = find_wrong_support_elements( input_mesh ) + + self.m_logger.info( "The following list displays the nodes buckets that contains the duplicated node indices." ) + self.m_logger.info( self.getCollocatedNodeBuckets() ) + + self.m_logger.info( "The following list displays the indexes of the cells with support node indices " + " appearing twice or more." ) + self.m_logger.info( self.getWrongSupportElements() ) + + output_mesh: vtkUnstructuredGrid = input_mesh.NewInstance() + output_mesh.CopyStructure( input_mesh ) + output_mesh.CopyAttributes( input_mesh ) + + if self.m_paintWrongSupportElements: + arrayWSP: npt.NDArray = np.zeros( ( output_mesh.GetNumberOfCells(), 1 ), dtype=int ) + arrayWSP[ self.m_wrongSupportElements ] = 1 + vtkArrayWSP: vtkDataArray = numpy_to_vtk( arrayWSP ) + vtkArrayWSP.SetName( "HasDuplicatedNodes" ) + output_mesh.GetCellData().AddArray( vtkArrayWSP ) + + output.ShallowCopy( output_mesh ) + + return 1 + + def getCollocatedNodeBuckets( self: Self ) -> list[ tuple[ int ] ]: + """Returns the nodes buckets that contains the duplicated node indices. + + Args: + self (Self) + + Returns: + list[ tuple[ int ] ] + """ + return self.m_collocatedNodesBuckets + + def getWrongSupportElements( self: Self ) -> list[ int ]: + """Returns the element indices with support node indices appearing more than once. + + Args: + self (Self) + + Returns: + list[ int ] + """ + return self.m_wrongSupportElements + + def setPaintWrongSupportElements( self: Self, choice: int ) -> None: + """Set 0 or 1 to choose if you want to create a new "WrongSupportElements" array in your output data. + + Args: + self (Self) + choice (int): 0 or 1 + """ + if choice not in [ 0, 1 ]: + self.m_logger.error( f"setPaintWrongSupportElements: Please choose either 0 or 1 not '{choice}'." ) + else: + self.m_paintWrongSupportElements = choice + self.Modified() + + def setTolerance( self: Self, tolerance: float ) -> None: + """Set the tolerance parameter to define if two points are collocated or not. + + Args: + self (Self) + tolerance (float) + """ + self.m_tolerance = tolerance + self.Modified() diff --git a/geos-mesh/src/geos/mesh/doctor/filters/ElementVolumes.py b/geos-mesh/src/geos/mesh/doctor/filters/ElementVolumes.py new file mode 100644 index 00000000..abb6fbe0 --- /dev/null +++ b/geos-mesh/src/geos/mesh/doctor/filters/ElementVolumes.py @@ -0,0 +1,122 @@ +import numpy as np +import numpy.typing as npt +from typing_extensions import Self +from vtkmodules.util.numpy_support import vtk_to_numpy +from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector, vtkDataArray +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from vtkmodules.vtkFiltersVerdict import vtkCellSizeFilter +from geos.mesh.doctor.filters.MeshDoctorBase import MeshDoctorBase + +__doc__ = """ +ElementVolumes module is a vtk filter that calculates the volumes of all elements in a vtkUnstructuredGrid. +The filter can identify elements with negative or zero volumes, which typically indicate mesh quality issues +such as inverted elements or degenerate cells. + +One filter input is vtkUnstructuredGrid, one filter output which is vtkUnstructuredGrid. + +To use the filter: + +.. code-block:: python + + from filters.ElementVolumes import ElementVolumes + + # instantiate the filter + elementVolumesFilter: ElementVolumes = ElementVolumes() + + # optionally enable detection of negative/zero volume elements + elementVolumesFilter.setReturnNegativeZeroVolumes(True) + + # set input mesh + elementVolumesFilter.SetInputData(mesh) + + # execute the filter + output_mesh: vtkUnstructuredGrid = elementVolumesFilter.getGrid() + + # get problematic elements (if enabled) + if elementVolumesFilter.m_returnNegativeZeroVolumes: + negative_zero_volumes = elementVolumesFilter.getNegativeZeroVolumes() + # returns numpy array with shape (n, 2) where first column is element index, second is volume + + # write the output mesh with volume information + elementVolumesFilter.writeGrid("output/mesh_with_volumes.vtu") +""" + + +class ElementVolumes( MeshDoctorBase ): + + def __init__( self: Self ) -> None: + """Vtk filter to calculate the volume of every element of a vtkUnstructuredGrid. + + Output mesh is vtkUnstructuredGrid. + """ + super().__init__( nInputPorts=1, + nOutputPorts=1, + inputType='vtkUnstructuredGrid', + outputType='vtkUnstructuredGrid' ) + self.m_returnNegativeZeroVolumes: bool = False + self.m_volumes: npt.NDArray = None + + def RequestData( self: Self, request: vtkInformation, inInfoVec: list[ vtkInformationVector ], + outInfo: vtkInformationVector ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestData. + + Args: + request (vtkInformation): request + inInfoVec (list[vtkInformationVector]): input objects + outInfoVec (vtkInformationVector): output objects + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + input_mesh: vtkUnstructuredGrid = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] ) + output = vtkUnstructuredGrid.GetData( outInfo ) + + cellSize = vtkCellSizeFilter() + cellSize.ComputeAreaOff() + cellSize.ComputeLengthOff() + cellSize.ComputeSumOff() + cellSize.ComputeVertexCountOff() + cellSize.ComputeVolumeOn() + volume_array_name: str = "MESH_DOCTOR_VOLUME" + cellSize.SetVolumeArrayName( volume_array_name ) + + cellSize.SetInputData( input_mesh ) + cellSize.Update() + volumes: vtkDataArray = cellSize.GetOutput().GetCellData().GetArray( volume_array_name ) + self.m_volumes = volumes + + output_mesh: vtkUnstructuredGrid = input_mesh.NewInstance() + output_mesh.CopyStructure( input_mesh ) + output_mesh.CopyAttributes( input_mesh ) + output_mesh.GetCellData().AddArray( volumes ) + output.ShallowCopy( output_mesh ) + + if self.m_returnNegativeZeroVolumes: + self.m_logger.info( "The following table displays the indexes of the cells with a zero or negative volume" ) + self.m_logger.info( self.getNegativeZeroVolumes() ) + + return 1 + + def getNegativeZeroVolumes( self: Self ) -> npt.NDArray: + """Returns a numpy array of all the negative and zero volumes of the input vtkUnstructuredGrid. + + Args: + self (Self) + + Returns: + npt.NDArray + """ + assert self.m_volumes is not None + volumes_np: npt.NDArray = vtk_to_numpy( self.m_volumes ) + indices = np.where( volumes_np <= 0 )[ 0 ] + return np.column_stack( ( indices, volumes_np[ indices ] ) ) + + def setReturnNegativeZeroVolumes( self: Self, returnNegativeZeroVolumes: bool ) -> None: + """Set the condition to return or not the negative and Zero volumes when calculating the volumes. + + Args: + self (Self) + returnNegativeZeroVolumes (bool) + """ + self.m_returnNegativeZeroVolumes = returnNegativeZeroVolumes + self.Modified() diff --git a/geos-mesh/src/geos/mesh/doctor/filters/GenerateFractures.py b/geos-mesh/src/geos/mesh/doctor/filters/GenerateFractures.py new file mode 100644 index 00000000..be024662 --- /dev/null +++ b/geos-mesh/src/geos/mesh/doctor/filters/GenerateFractures.py @@ -0,0 +1,194 @@ +from typing_extensions import Self +from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from geos.mesh.doctor.actions.generate_fractures import Options, split_mesh_on_fractures +from geos.mesh.doctor.filters.MeshDoctorBase import MeshDoctorBase +from geos.mesh.doctor.parsing.generate_fractures_parsing import convert, convert_to_fracture_policy +from geos.mesh.doctor.parsing.generate_fractures_parsing import ( __FIELD_NAME, __FIELD_VALUES, __FRACTURES_DATA_MODE, + __FRACTURES_OUTPUT_DIR, __FRACTURES_DATA_MODE_VALUES, + __POLICIES, __POLICY ) +from geos.mesh.io.vtkIO import VtkOutput, write_mesh +from geos.mesh.utils.arrayHelpers import has_array + +__doc__ = """ +GenerateFractures module is a vtk filter that splits a vtkUnstructuredGrid along non-embedded fractures. +When a fracture plane is defined between two cells, the nodes of the shared face will be duplicated +to create a discontinuity. The filter generates both the split main mesh and separate fracture meshes. + +One filter input is vtkUnstructuredGrid, two filter outputs which are vtkUnstructuredGrid. + +To use the filter: + +.. code-block:: python + + from filters.GenerateFractures import GenerateFractures + + # instantiate the filter + generateFracturesFilter: GenerateFractures = GenerateFractures() + + # set the field name that defines fracture regions + generateFracturesFilter.setFieldName("fracture_field") + + # set the field values that identify fracture boundaries + generateFracturesFilter.setFieldValues("1,2") # comma-separated values + + # set fracture policy (0 for internal fractures, 1 for boundary fractures) + generateFracturesFilter.setPolicy(1) + + # set output directory for fracture meshes + generateFracturesFilter.setFracturesOutputDirectory("./fractures/") + + # optionally set data mode (0 for ASCII, 1 for binary) + generateFracturesFilter.setOutputDataMode(1) + generateFracturesFilter.setFracturesDataMode(1) + + # set input mesh + generateFracturesFilter.SetInputData(mesh) + + # execute the filter + generateFracturesFilter.Update() + + # get the split mesh and fracture meshes + split_mesh, fracture_meshes = generateFracturesFilter.getAllGrids() + + # write all meshes + generateFracturesFilter.writeMeshes("output/split_mesh.vtu", is_data_mode_binary=True) +""" + +FIELD_NAME = __FIELD_NAME +FIELD_VALUES = __FIELD_VALUES +FRACTURES_DATA_MODE = __FRACTURES_DATA_MODE +DATA_MODE = __FRACTURES_DATA_MODE_VALUES +FRACTURES_OUTPUT_DIR = __FRACTURES_OUTPUT_DIR +POLICIES = __POLICIES +POLICY = __POLICY + + +class GenerateFractures( MeshDoctorBase ): + + def __init__( self: Self ) -> None: + """Vtk filter to generate a simple rectilinear grid. + + Output mesh is vtkUnstructuredGrid. + """ + super().__init__( nInputPorts=1, + nOutputPorts=2, + inputType='vtkUnstructuredGrid', + outputType='vtkUnstructuredGrid' ) + self.m_policy: str = POLICIES[ 1 ] + self.m_field_name: str = None + self.m_field_values: str = None + self.m_fractures_output_dir: str = None + self.m_output_modes_binary: str = { "mesh": DATA_MODE[ 0 ], "fractures": DATA_MODE[ 1 ] } + self.m_mesh_VtkOutput: VtkOutput = None + self.m_all_fractures_VtkOutput: list[ VtkOutput ] = None + + def RequestData( self: Self, request: vtkInformation, inInfoVec: list[ vtkInformationVector ], + outInfo: list[ vtkInformationVector ] ) -> int: + input_mesh = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] ) + if has_array( input_mesh, [ "GLOBAL_IDS_POINTS", "GLOBAL_IDS_CELLS" ] ): + err_msg: str = ( "The mesh cannot contain global ids for neither cells nor points. The correct procedure " + + " is to split the mesh and then generate global ids for new split meshes." ) + self.m_logger.error( err_msg ) + return 0 + + parsed_options: dict[ str, str ] = self.getParsedOptions() + self.m_logger.critical( f"Parsed_options:\n{parsed_options}" ) + if len( parsed_options ) < 5: + self.m_logger.error( "You must set all variables before trying to create fractures." ) + return 0 + + options: Options = convert( parsed_options ) + self.m_all_fractures_VtkOutput = options.all_fractures_VtkOutput + output_mesh, fracture_meshes = split_mesh_on_fractures( input_mesh, options ) + opt = vtkUnstructuredGrid.GetData( outInfo, 0 ) + opt.ShallowCopy( output_mesh ) + + nbr_faults: int = len( fracture_meshes ) + self.SetNumberOfOutputPorts( 1 + nbr_faults ) # one output port for splitted mesh, the rest for every fault + for i in range( nbr_faults ): + opt_fault = vtkUnstructuredGrid.GetData( outInfo, i + 1 ) + opt_fault.ShallowCopy( fracture_meshes[ i ] ) + + return 1 + + def getAllGrids( self: Self ) -> tuple[ vtkUnstructuredGrid, list[ vtkUnstructuredGrid ] ]: + """Returns the vtkUnstructuredGrid with volumes. + + Args: + self (Self) + + Returns: + vtkUnstructuredGrid + """ + self.Update() # triggers RequestData + splitted_grid: vtkUnstructuredGrid = self.GetOutputDataObject( 0 ) + nbrOutputPorts: int = self.GetNumberOfOutputPorts() + fracture_meshes: list[ vtkUnstructuredGrid ] = list() + for i in range( 1, nbrOutputPorts ): + fracture_meshes.append( self.GetOutputDataObject( i ) ) + return ( splitted_grid, fracture_meshes ) + + def getParsedOptions( self: Self ) -> dict[ str, str ]: + parsed_options: dict[ str, str ] = { "output": "./mesh.vtu", "data_mode": DATA_MODE[ 0 ] } + parsed_options[ POLICY ] = self.m_policy + parsed_options[ FRACTURES_DATA_MODE ] = self.m_output_modes_binary[ "fractures" ] + if self.m_field_name: + parsed_options[ FIELD_NAME ] = self.m_field_name + else: + self.m_logger.error( "No field name provided. Please use setFieldName." ) + if self.m_field_values: + parsed_options[ FIELD_VALUES ] = self.m_field_values + else: + self.m_logger.error( "No field values provided. Please use setFieldValues." ) + if self.m_fractures_output_dir: + parsed_options[ FRACTURES_OUTPUT_DIR ] = self.m_fractures_output_dir + else: + self.m_logger.error( "No fracture output directory provided. Please use setFracturesOutputDirectory." ) + return parsed_options + + def setFieldName( self: Self, field_name: str ) -> None: + self.m_field_name = field_name + self.Modified() + + def setFieldValues( self: Self, field_values: str ) -> None: + self.m_field_values = field_values + self.Modified() + + def setFracturesDataMode( self: Self, choice: int ) -> None: + if choice not in [ 0, 1 ]: + self.m_logger.error( f"setFracturesDataMode: Please choose either 0 for {DATA_MODE[ 0 ]} or 1 for", + f" {DATA_MODE[ 1 ]}, not '{choice}'." ) + else: + self.m_output_modes_binary[ "fractures" ] = DATA_MODE[ choice ] + self.Modified() + + def setFracturesOutputDirectory( self: Self, directory: str ) -> None: + self.m_fractures_output_dir = directory + self.Modified() + + def setOutputDataMode( self: Self, choice: int ) -> None: + if choice not in [ 0, 1 ]: + self.m_logger.error( f"setOutputDataMode: Please choose either 0 for {DATA_MODE[ 0 ]} or 1 for", + f" {DATA_MODE[ 1 ]}, not '{choice}'." ) + else: + self.m_output_modes_binary[ "mesh" ] = DATA_MODE[ choice ] + self.Modified() + + def setPolicy( self: Self, choice: int ) -> None: + if choice not in [ 0, 1 ]: + self.m_logger.error( f"setPolicy: Please choose either 0 for {POLICIES[ 0 ]} or 1 for {POLICIES[ 1 ]}," + f" not '{choice}'." ) + else: + self.m_policy = convert_to_fracture_policy( POLICIES[ choice ] ) + self.Modified() + + def writeMeshes( self, filepath: str, is_data_mode_binary: bool = True, canOverwrite: bool = False ) -> None: + splitted_grid, fracture_meshes = self.getAllGrids() + if splitted_grid: + write_mesh( splitted_grid, VtkOutput( filepath, is_data_mode_binary ), canOverwrite ) + else: + self.m_logger.error( f"No output grid was built. Cannot output vtkUnstructuredGrid at {filepath}." ) + + for i, fracture_mesh in enumerate( fracture_meshes ): + write_mesh( fracture_mesh, self.m_all_fractures_VtkOutput[ i ] ) diff --git a/geos-mesh/src/geos/mesh/doctor/filters/GenerateRectilinearGrid.py b/geos-mesh/src/geos/mesh/doctor/filters/GenerateRectilinearGrid.py new file mode 100644 index 00000000..56bfe299 --- /dev/null +++ b/geos-mesh/src/geos/mesh/doctor/filters/GenerateRectilinearGrid.py @@ -0,0 +1,146 @@ +import numpy.typing as npt +from typing import Iterable, Sequence +from typing_extensions import Self +from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from geos.mesh.doctor.actions.generate_global_ids import build_global_ids +from geos.mesh.doctor.actions.generate_cube import FieldInfo, add_fields, build_coordinates, build_rectilinear_grid +from geos.mesh.doctor.filters.MeshDoctorBase import MeshDoctorGenerator + +__doc__ = """ +GenerateRectilinearGrid module is a vtk filter that allows to create a simple vtkUnstructuredGrid rectilinear grid. +GlobalIds for points and cells can be added. +You can create CellArray and PointArray of constant value = 1 and any dimension >= 1. + +No filter input and one filter output which is vtkUnstructuredGrid. + +To use the filter: + +.. code-block:: python + + from filters.GenerateRectilinearGrid import GenerateRectilinearGrid + + # instanciate the filter + generateRectilinearGridFilter: GenerateRectilinearGrid = GenerateRectilinearGrid() + + # set the coordinates of each block border for the X, Y and Z axis + generateRectilinearGridFilter.setCoordinates( [ 0.0, 5.0, 10.0 ], [ 0.0, 5.0, 10.0 ], [ 0.0, 10.0 ] ) + + # for each block defined, specify the number of cells that they should contain in the X, Y, Z axis + generateRectilinearGridFilter.setNumberElements( [ 5, 5 ], [ 5, 5 ], [ 10 ] ) + + # to add the GlobalIds for cells and points, set to True the generate global ids options + generateRectilinearGridFilter.setGenerateCellsGlobalIds( True ) + generateRectilinearGridFilter.setGeneratePointsGlobalIds( True ) + + # to create new arrays with a specific dimension, you can use the following commands + cells_dim1 = FieldInfo( "cell1", 1, "CELLS" ) # array "cell1" of shape ( number of cells, 1 ) + cells_dim3 = FieldInfo( "cell3", 3, "CELLS" ) # array "cell3" of shape ( number of cells, 3 ) + points_dim1 = FieldInfo( "point1", 1, "POINTS" ) # array "point1" of shape ( number of points, 1 ) + points_dim3 = FieldInfo( "point3", 3, "POINTS" ) # array "point3" of shape ( number of points, 3 ) + generateRectilinearGridFilter.setFields( [ cells_dim1, cells_dim3, points_dim1, points_dim3 ] ) + + # then, to obtain the constructed mesh out of all these operations, 2 solutions are available + + # solution1 + mesh: vtkUnstructuredGrid = generateRectilinearGridFilter.getGrid() + + # solution2, which calls the same method as above + generateRectilinearGridFilter.Update() + mesh: vtkUnstructuredGrid = generateRectilinearGridFilter.GetOutputDataObject( 0 ) + + # finally, you can write the mesh at a specific destination with: + generateRectilinearGridFilter.writeGrid( "output/filepath/of/your/grid.vtu" ) +""" + + +class GenerateRectilinearGrid( MeshDoctorGenerator ): + + def __init__( self: Self ) -> None: + """Vtk filter to generate a simple rectilinear grid. + + Output mesh is vtkUnstructuredGrid. + """ + super().__init__() + self.m_generateCellsGlobalIds: bool = False + self.m_generatePointsGlobalIds: bool = False + self.m_coordsX: Sequence[ float ] = None + self.m_coordsY: Sequence[ float ] = None + self.m_coordsZ: Sequence[ float ] = None + self.m_numberElementsX: Sequence[ int ] = None + self.m_numberElementsY: Sequence[ int ] = None + self.m_numberElementsZ: Sequence[ int ] = None + self.m_fields: Iterable[ FieldInfo ] = list() + + def RequestData( self: Self, request: vtkInformation, inInfo: vtkInformationVector, + outInfo: vtkInformationVector ) -> int: + opt = vtkUnstructuredGrid.GetData( outInfo ) + x: npt.NDArray = build_coordinates( self.m_coordsX, self.m_numberElementsX ) + y: npt.NDArray = build_coordinates( self.m_coordsY, self.m_numberElementsY ) + z: npt.NDArray = build_coordinates( self.m_coordsZ, self.m_numberElementsZ ) + output: vtkUnstructuredGrid = build_rectilinear_grid( x, y, z ) + output = add_fields( output, self.m_fields ) + build_global_ids( output, self.m_generateCellsGlobalIds, self.m_generatePointsGlobalIds ) + opt.ShallowCopy( output ) + return 1 + + def setCoordinates( self: Self, coordsX: Sequence[ float ], coordsY: Sequence[ float ], + coordsZ: Sequence[ float ] ) -> None: + """Set the coordinates of the block you want to have in your grid by specifying the beginning and ending + coordinates along the X, Y and Z axis. + + Args: + self (Self) + coordsX (Sequence[ float ]) + coordsY (Sequence[ float ]) + coordsZ (Sequence[ float ]) + """ + self.m_coordsX = coordsX + self.m_coordsY = coordsY + self.m_coordsZ = coordsZ + self.Modified() + + def setGenerateCellsGlobalIds( self: Self, generate: bool ) -> None: + """Set the generation of global cells ids to be True or False. + + Args: + self (Self) + generate (bool) + """ + self.m_generateCellsGlobalIds = generate + self.Modified() + + def setGeneratePointsGlobalIds( self: Self, generate: bool ) -> None: + """Set the generation of global points ids to be True or False. + + Args: + self (Self) + generate (bool) + """ + self.m_generatePointsGlobalIds = generate + self.Modified() + + def setFields( self: Self, fields: Iterable[ FieldInfo ] ) -> None: + """Specify the cells or points array to be added to the grid. + + Args: + self (Self) + fields (Iterable[ FieldInfo ]) + """ + self.m_fields = fields + self.Modified() + + def setNumberElements( self: Self, numberElementsX: Sequence[ int ], numberElementsY: Sequence[ int ], + numberElementsZ: Sequence[ int ] ) -> None: + """For each block that was defined in setCoordinates, specify the number of cells that they should contain. + + Args: + self (Self) + numberElementsX (Sequence[ int ]) + numberElementsY (Sequence[ int ]) + numberElementsZ (Sequence[ int ]) + """ + self.m_numberElementsX = numberElementsX + self.m_numberElementsY = numberElementsY + self.m_numberElementsZ = numberElementsZ + self.Modified() diff --git a/geos-mesh/src/geos/mesh/doctor/filters/MeshDoctorBase.py b/geos-mesh/src/geos/mesh/doctor/filters/MeshDoctorBase.py new file mode 100644 index 00000000..d85213ab --- /dev/null +++ b/geos-mesh/src/geos/mesh/doctor/filters/MeshDoctorBase.py @@ -0,0 +1,190 @@ +from typing_extensions import Self +from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase +from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from geos.mesh.doctor.parsing.cli_parsing import setup_logger +from geos.mesh.io.vtkIO import VtkOutput, write_mesh + +__doc__ = """ +MeshDoctorBase module provides base classes for all mesh doctor VTK filters. + +MeshDoctorBase serves as the foundation class for filters that process existing meshes, +while MeshDoctorGenerator is for filters that generate new meshes from scratch. + +These base classes provide common functionality including: +- Logger management and setup +- Grid access and manipulation methods +- File I/O operations for writing VTK unstructured grids +- Standard VTK filter interface implementation + +All mesh doctor filters should inherit from one of these base classes to ensure +consistent behavior and interface across the mesh doctor toolkit. + +Example usage patterns: + +.. code-block:: python + + # For filters that process existing meshes + from filters.MeshDoctorBase import MeshDoctorBase + + class MyProcessingFilter(MeshDoctorBase): + def __init__(self): + super().__init__(nInputPorts=1, nOutputPorts=1) + + def RequestData(self, request, inInfoVec, outInfo): + # Process input mesh and create output + pass + + # For filters that generate meshes from scratch + from filters.MeshDoctorBase import MeshDoctorGenerator + + class MyGeneratorFilter(MeshDoctorGenerator): + def __init__(self): + super().__init__(nOutputPorts=1) + + def RequestData(self, request, inInfo, outInfo): + # Generate new mesh + pass +""" + + +class MeshDoctorBase( VTKPythonAlgorithmBase ): + """Base class for all mesh doctor VTK filters. + + This class provides common functionality shared across all mesh doctor filters, + including logger management, grid access, and file writing capabilities. + """ + + def __init__( self: Self, + nInputPorts: int = 1, + nOutputPorts: int = 1, + inputType: str = 'vtkUnstructuredGrid', + outputType: str = 'vtkUnstructuredGrid' ) -> None: + """Initialize the base mesh doctor filter. + + Args: + nInputPorts (int): Number of input ports. Defaults to 1. + nOutputPorts (int): Number of output ports. Defaults to 1. + inputType (str): Input data type. Defaults to 'vtkUnstructuredGrid'. + outputType (str): Output data type. Defaults to 'vtkUnstructuredGrid'. + """ + super().__init__( nInputPorts=nInputPorts, + nOutputPorts=nOutputPorts, + inputType=inputType if nInputPorts > 0 else None, + outputType=outputType ) + self.m_logger = setup_logger + + def FillInputPortInformation( self: Self, port: int, info: vtkInformation ) -> int: + """Inherited from VTKPythonAlgorithmBase::FillInputPortInformation. + + Args: + port (int): input port + info (vtkInformationVector): info + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + if port == 0: + info.Set( self.INPUT_REQUIRED_DATA_TYPE(), "vtkUnstructuredGrid" ) + return 1 + + def RequestInformation( + self: Self, + request: vtkInformation, # noqa: F841 + inInfoVec: list[ vtkInformationVector ], # noqa: F841 + outInfoVec: vtkInformationVector, + ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestInformation. + + Args: + request (vtkInformation): request + inInfoVec (list[vtkInformationVector]): input objects + outInfoVec (vtkInformationVector): output objects + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + executive = self.GetExecutive() # noqa: F841 + outInfo = outInfoVec.GetInformationObject( 0 ) # noqa: F841 + return 1 + + def SetLogger( self: Self, logger ) -> None: + """Set the logger. + + Args: + logger: Logger instance to use + """ + self.m_logger = logger + self.Modified() + + def getGrid( self: Self ) -> vtkUnstructuredGrid: + """Returns the vtkUnstructuredGrid output. + + Args: + self (Self) + + Returns: + vtkUnstructuredGrid: The output grid + """ + self.Update() # triggers RequestData + return self.GetOutputDataObject( 0 ) + + def writeGrid( self: Self, filepath: str, is_data_mode_binary: bool = True, canOverwrite: bool = False ) -> None: + """Writes a .vtu file of the vtkUnstructuredGrid at the specified filepath. + + Args: + filepath (str): /path/to/your/file.vtu + is_data_mode_binary (bool, optional): Writes the file in binary format or ascii. Defaults to True. + canOverwrite (bool, optional): Allows or not to overwrite if the filepath already leads to an existing file. + Defaults to False. + """ + mesh: vtkUnstructuredGrid = self.getGrid() + if mesh: + vtk_output = VtkOutput( filepath, is_data_mode_binary ) + write_mesh( mesh, vtk_output, canOverwrite ) + else: + self.m_logger.error( f"No output grid was built. Cannot output vtkUnstructuredGrid at {filepath}." ) + + def copyInputToOutput( self: Self, input_mesh: vtkUnstructuredGrid ) -> vtkUnstructuredGrid: + """Helper method to copy input mesh structure and attributes to a new output mesh. + + Args: + input_mesh (vtkUnstructuredGrid): Input mesh to copy from + + Returns: + vtkUnstructuredGrid: New mesh with copied structure and attributes + """ + output_mesh: vtkUnstructuredGrid = input_mesh.NewInstance() + output_mesh.CopyStructure( input_mesh ) + output_mesh.CopyAttributes( input_mesh ) + return output_mesh + + +class MeshDoctorGenerator( MeshDoctorBase ): + """Base class for mesh doctor generator filters (no input required). + + This class extends MeshDoctorBase for filters that generate meshes + from scratch without requiring input meshes. + """ + + def __init__( self: Self, nOutputPorts: int = 1, outputType: str = 'vtkUnstructuredGrid' ) -> None: + """Initialize the base mesh doctor generator filter. + + Args: + nOutputPorts (int): Number of output ports. Defaults to 1. + outputType (str): Output data type. Defaults to 'vtkUnstructuredGrid'. + """ + super().__init__( nInputPorts=0, nOutputPorts=nOutputPorts, inputType=None, outputType=outputType ) + + def FillInputPortInformation( self: Self, port: int, info: vtkInformation ) -> int: + """Generator filters don't have input ports. + + Args: + port (int): input port (not used) + info (vtkInformationVector): info (not used) + + Returns: + int: Always returns 1 + """ + # Generator filters don't have input ports, so this method is not used + return 1 diff --git a/geos-mesh/src/geos/mesh/doctor/filters/NonConformal.py b/geos-mesh/src/geos/mesh/doctor/filters/NonConformal.py new file mode 100644 index 00000000..22d927ca --- /dev/null +++ b/geos-mesh/src/geos/mesh/doctor/filters/NonConformal.py @@ -0,0 +1,179 @@ +import numpy as np +import numpy.typing as npt +from typing_extensions import Self +from vtkmodules.util.numpy_support import numpy_to_vtk +from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector, vtkDataArray +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from geos.mesh.doctor.actions.non_conformal import Options, find_non_conformal_cells +from geos.mesh.doctor.filters.MeshDoctorBase import MeshDoctorBase + +__doc__ = """ +NonConformal module is a vtk filter that detects non-conformal mesh interfaces in a vtkUnstructuredGrid. +Non-conformal interfaces occur when adjacent cells do not share nodes or faces properly, which can indicate +mesh quality issues or intentional non-matching grid interfaces that require special handling. + +One filter input is vtkUnstructuredGrid, one filter output which is vtkUnstructuredGrid. + +To use the filter: + +.. code-block:: python + + from filters.NonConformal import NonConformal + + # instantiate the filter + nonConformalFilter: NonConformal = NonConformal() + + # set tolerance parameters + nonConformalFilter.setPointTolerance(1e-6) # tolerance for point matching + nonConformalFilter.setFaceTolerance(1e-6) # tolerance for face matching + nonConformalFilter.setAngleTolerance(10.0) # angle tolerance in degrees + + # optionally enable painting of non-conformal cells + nonConformalFilter.setPaintNonConformalCells(1) # 1 to enable, 0 to disable + + # set input mesh + nonConformalFilter.SetInputData(mesh) + + # execute the filter + output_mesh: vtkUnstructuredGrid = nonConformalFilter.getGrid() + + # get non-conformal cell pairs + non_conformal_cells = nonConformalFilter.getNonConformalCells() + # returns list of tuples with (cell1_id, cell2_id) for non-conformal interfaces + + # write the output mesh + nonConformalFilter.writeGrid("output/mesh_with_nonconformal_info.vtu") +""" + + +class NonConformal( MeshDoctorBase ): + + def __init__( self: Self ) -> None: + """Vtk filter to ... of a vtkUnstructuredGrid. + + Output mesh is vtkUnstructuredGrid. + """ + super().__init__( nInputPorts=1, + nOutputPorts=1, + inputType='vtkUnstructuredGrid', + outputType='vtkUnstructuredGrid' ) + self.m_angle_tolerance: float = 10.0 + self.m_face_tolerance: float = 0.0 + self.m_point_tolerance: float = 0.0 + self.m_non_conformal_cells: list[ tuple[ int, int ] ] = list() + self.m_paintNonConformalCells: int = 0 + + def RequestData( self: Self, request: vtkInformation, inInfoVec: list[ vtkInformationVector ], + outInfo: vtkInformationVector ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestData. + + Args: + request (vtkInformation): request + inInfoVec (list[vtkInformationVector]): input objects + outInfoVec (vtkInformationVector): output objects + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + input_mesh: vtkUnstructuredGrid = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] ) + output = vtkUnstructuredGrid.GetData( outInfo ) + + options = Options( self.m_angle_tolerance, self.m_point_tolerance, self.m_face_tolerance ) + non_conformal_cells = find_non_conformal_cells( input_mesh, options ) + self.m_non_conformal_cells = non_conformal_cells + + non_conformal_cells_extended = [ cell_id for pair in non_conformal_cells for cell_id in pair ] + unique_non_conformal_cells = frozenset( non_conformal_cells_extended ) + self.m_logger.info( f"You have {len( unique_non_conformal_cells )} non conformal cells.\n" + f"{', '.join( map( str, sorted( non_conformal_cells_extended ) ) )}" ) + + output_mesh: vtkUnstructuredGrid = input_mesh.NewInstance() + output_mesh.CopyStructure( input_mesh ) + output_mesh.CopyAttributes( input_mesh ) + + if self.m_paintNonConformalCells: + arrayNC: npt.NDArray = np.zeros( ( output_mesh.GetNumberOfCells(), 1 ), dtype=int ) + arrayNC[ unique_non_conformal_cells ] = 1 + vtkArrayNC: vtkDataArray = numpy_to_vtk( arrayNC ) + vtkArrayNC.SetName( "IsNonConformal" ) + output_mesh.GetCellData().AddArray( vtkArrayNC ) + + output.ShallowCopy( output_mesh ) + + return 1 + + def getAngleTolerance( self: Self ) -> float: + """Returns the angle tolerance. + + Args: + self (Self) + + Returns: + float + """ + return self.m_angle_tolerance + + def getfaceTolerance( self: Self ) -> float: + """Returns the face tolerance. + + Args: + self (Self) + + Returns: + float + """ + return self.m_face_tolerance + + def getPointTolerance( self: Self ) -> float: + """Returns the point tolerance. + + Args: + self (Self) + + Returns: + float + """ + return self.m_point_tolerance + + def setPaintNonConformalCells( self: Self, choice: int ) -> None: + """Set 0 or 1 to choose if you want to create a new "IsNonConformal" array in your output data. + + Args: + self (Self) + choice (int): 0 or 1 + """ + if choice not in [ 0, 1 ]: + self.m_logger.error( f"setPaintNonConformalCells: Please choose either 0 or 1 not '{choice}'." ) + else: + self.m_paintNonConformalCells = choice + self.Modified() + + def setAngleTolerance( self: Self, tolerance: float ) -> None: + """Set the angle tolerance parameter in degree. + + Args: + self (Self) + tolerance (float) + """ + self.m_angle_tolerance = tolerance + self.Modified() + + def setFaceTolerance( self: Self, tolerance: float ) -> None: + """Set the face tolerance parameter. + + Args: + self (Self) + tolerance (float) + """ + self.m_face_tolerance = tolerance + self.Modified() + + def setPointTolerance( self: Self, tolerance: float ) -> None: + """Set the point tolerance parameter. + + Args: + self (Self) + tolerance (float) + """ + self.m_point_tolerance = tolerance + self.Modified() diff --git a/geos-mesh/src/geos/mesh/doctor/filters/SelfIntersectingElements.py b/geos-mesh/src/geos/mesh/doctor/filters/SelfIntersectingElements.py new file mode 100644 index 00000000..5db9dc2b --- /dev/null +++ b/geos-mesh/src/geos/mesh/doctor/filters/SelfIntersectingElements.py @@ -0,0 +1,218 @@ +import numpy as np +import numpy.typing as npt +from typing_extensions import Self +from vtkmodules.util.numpy_support import numpy_to_vtk +from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector, vtkDataArray +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from geos.mesh.doctor.actions.self_intersecting_elements import get_invalid_cell_ids +from geos.mesh.doctor.filters.MeshDoctorBase import MeshDoctorBase + +__doc__ = """ +SelfIntersectingElements module is a vtk filter that identifies various types of invalid or problematic elements +in a vtkUnstructuredGrid. It detects elements with intersecting edges, intersecting faces, non-contiguous edges, +non-convex shapes, incorrectly oriented faces, and wrong number of points. + +One filter input is vtkUnstructuredGrid, one filter output which is vtkUnstructuredGrid. + +To use the filter: + +.. code-block:: python + + from filters.SelfIntersectingElements import SelfIntersectingElements + + # instantiate the filter + selfIntersectingElementsFilter: SelfIntersectingElements = SelfIntersectingElements() + + # set minimum distance parameter for intersection detection + selfIntersectingElementsFilter.setMinDistance(1e-6) + + # optionally enable painting of invalid elements + selfIntersectingElementsFilter.setPaintInvalidElements(1) # 1 to enable, 0 to disable + + # set input mesh + selfIntersectingElementsFilter.SetInputData(mesh) + + # execute the filter + output_mesh: vtkUnstructuredGrid = selfIntersectingElementsFilter.getGrid() + + # get different types of problematic elements + wrong_points_elements = selfIntersectingElementsFilter.getWrongNumberOfPointsElements() + intersecting_edges_elements = selfIntersectingElementsFilter.getIntersectingEdgesElements() + intersecting_faces_elements = selfIntersectingElementsFilter.getIntersectingFacesElements() + non_contiguous_edges_elements = selfIntersectingElementsFilter.getNonContiguousEdgesElements() + non_convex_elements = selfIntersectingElementsFilter.getNonConvexElements() + wrong_oriented_faces_elements = selfIntersectingElementsFilter.getFacesOrientedIncorrectlyElements() + + # write the output mesh + selfIntersectingElementsFilter.writeGrid("output/mesh_with_invalid_elements.vtu") +""" + + +class SelfIntersectingElements( MeshDoctorBase ): + + def __init__( self: Self ) -> None: + """Vtk filter to find invalid elements of a vtkUnstructuredGrid. + + Output mesh is vtkUnstructuredGrid. + """ + super().__init__( nInputPorts=1, + nOutputPorts=1, + inputType='vtkUnstructuredGrid', + outputType='vtkUnstructuredGrid' ) + self.m_min_distance: float = 0.0 + self.m_wrong_number_of_points_elements: list[ int ] = list() + self.m_intersecting_edges_elements: list[ int ] = list() + self.m_intersecting_faces_elements: list[ int ] = list() + self.m_non_contiguous_edges_elements: list[ int ] = list() + self.m_non_convex_elements: list[ int ] = list() + self.m_faces_oriented_incorrectly_elements: list[ int ] = list() + self.m_paintInvalidElements: int = 0 + + def RequestData( self: Self, request: vtkInformation, inInfoVec: list[ vtkInformationVector ], + outInfo: vtkInformationVector ) -> int: + """Inherited from VTKPythonAlgorithmBase::RequestData. + + Args: + request (vtkInformation): request + inInfoVec (list[vtkInformationVector]): input objects + outInfoVec (vtkInformationVector): output objects + + Returns: + int: 1 if calculation successfully ended, 0 otherwise. + """ + input_mesh: vtkUnstructuredGrid = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] ) + output = vtkUnstructuredGrid.GetData( outInfo ) + + invalid_cells = get_invalid_cell_ids( input_mesh, self.m_min_distance ) + + self.m_wrong_number_of_points_elements = invalid_cells.get( "wrong_number_of_points_elements", [] ) + self.m_intersecting_edges_elements = invalid_cells.get( "intersecting_edges_elements", [] ) + self.m_intersecting_faces_elements = invalid_cells.get( "intersecting_faces_elements", [] ) + self.m_non_contiguous_edges_elements = invalid_cells.get( "non_contiguous_edges_elements", [] ) + self.m_non_convex_elements = invalid_cells.get( "non_convex_elements", [] ) + self.m_faces_oriented_incorrectly_elements = invalid_cells.get( "faces_oriented_incorrectly_elements", [] ) + + # Log the results + total_invalid = sum( len( invalid_list ) for invalid_list in invalid_cells.values() ) + self.m_logger.info( f"Found {total_invalid} invalid elements:" ) + for criterion, cell_list in invalid_cells.items(): + if cell_list: + self.m_logger.info( f" {criterion}: {len( cell_list )} elements - {cell_list}" ) + + output_mesh: vtkUnstructuredGrid = input_mesh.NewInstance() + output_mesh.CopyStructure( input_mesh ) + output_mesh.CopyAttributes( input_mesh ) + + if self.m_paintInvalidElements: + # Create arrays to mark invalid elements + for criterion, cell_list in invalid_cells.items(): + if cell_list: + array: npt.NDArray = np.zeros( ( output_mesh.GetNumberOfCells(), 1 ), dtype=int ) + array[ cell_list ] = 1 + vtkArray: vtkDataArray = numpy_to_vtk( array ) + vtkArray.SetName( f"Is{criterion.replace('_', '').title()}" ) + output_mesh.GetCellData().AddArray( vtkArray ) + + output.ShallowCopy( output_mesh ) + + return 1 + + def getMinDistance( self: Self ) -> float: + """Returns the minimum distance. + + Args: + self (Self) + + Returns: + float + """ + return self.m_min_distance + + def getWrongNumberOfPointsElements( self: Self ) -> list[ int ]: + """Returns elements with wrong number of points. + + Args: + self (Self) + + Returns: + list[int] + """ + return self.m_wrong_number_of_points_elements + + def getIntersectingEdgesElements( self: Self ) -> list[ int ]: + """Returns elements with intersecting edges. + + Args: + self (Self) + + Returns: + list[int] + """ + return self.m_intersecting_edges_elements + + def getIntersectingFacesElements( self: Self ) -> list[ int ]: + """Returns elements with intersecting faces. + + Args: + self (Self) + + Returns: + list[int] + """ + return self.m_intersecting_faces_elements + + def getNonContiguousEdgesElements( self: Self ) -> list[ int ]: + """Returns elements with non-contiguous edges. + + Args: + self (Self) + + Returns: + list[int] + """ + return self.m_non_contiguous_edges_elements + + def getNonConvexElements( self: Self ) -> list[ int ]: + """Returns non-convex elements. + + Args: + self (Self) + + Returns: + list[int] + """ + return self.m_non_convex_elements + + def getFacesOrientedIncorrectlyElements( self: Self ) -> list[ int ]: + """Returns elements with incorrectly oriented faces. + + Args: + self (Self) + + Returns: + list[int] + """ + return self.m_faces_oriented_incorrectly_elements + + def setPaintInvalidElements( self: Self, choice: int ) -> None: + """Set 0 or 1 to choose if you want to create arrays marking invalid elements in your output data. + + Args: + self (Self) + choice (int): 0 or 1 + """ + if choice not in [ 0, 1 ]: + self.m_logger.error( f"setPaintInvalidElements: Please choose either 0 or 1 not '{choice}'." ) + else: + self.m_paintInvalidElements = choice + self.Modified() + + def setMinDistance( self: Self, distance: float ) -> None: + """Set the minimum distance parameter. + + Args: + self (Self) + distance (float) + """ + self.m_min_distance = distance + self.Modified() diff --git a/geos-mesh/src/geos/mesh/doctor/filters/SupportedElements.py b/geos-mesh/src/geos/mesh/doctor/filters/SupportedElements.py new file mode 100644 index 00000000..4a3e26ea --- /dev/null +++ b/geos-mesh/src/geos/mesh/doctor/filters/SupportedElements.py @@ -0,0 +1,229 @@ +# TODO Find an implementation to keep multiprocessing while using vtkFilter + +# import numpy as np +# import numpy.typing as npt +# from typing_extensions import Self +# from vtkmodules.util.numpy_support import numpy_to_vtk +# from vtkmodules.util.vtkAlgorithm import VTKPythonAlgorithmBase +# from vtkmodules.vtkCommonCore import vtkInformation, vtkInformationVector, vtkDataArray, VTK_INT +# from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +# from geos.mesh.doctor.actions.supported_elements import ( Options, find_unsupported_std_elements_types, +# find_unsupported_polyhedron_elements ) +# from geos.mesh.io.vtkIO import VtkOutput, write_mesh +# from geos.utils.Logger import Logger, getLogger + +# __doc__ = """ +# SupportedElements module is a vtk filter that identifies unsupported element types and problematic polyhedron +# elements in a vtkUnstructuredGrid. It checks for element types that are not supported by GEOS and validates +# polyhedron elements for geometric correctness. +# +# One filter input is vtkUnstructuredGrid, one filter output which is vtkUnstructuredGrid. +# +# To use the filter: +# +# .. code-block:: python +# +# from filters.SupportedElements import SupportedElements +# +# # instantiate the filter +# supportedElementsFilter: SupportedElements = SupportedElements() +# +# # optionally enable painting of unsupported element types +# supportedElementsFilter.setPaintUnsupportedElementTypes(1) # 1 to enable, 0 to disable +# +# # set input mesh +# supportedElementsFilter.SetInputData(mesh) +# +# # execute the filter +# output_mesh: vtkUnstructuredGrid = supportedElementsFilter.getGrid() +# +# # get unsupported elements +# unsupported_elements = supportedElementsFilter.getUnsupportedElements() +# +# # write the output mesh +# supportedElementsFilter.writeGrid("output/mesh_with_support_info.vtu") +# +# Note: This filter is currently disabled due to multiprocessing requirements. +# """ + +# class SupportedElements( VTKPythonAlgorithmBase ): + +# def __init__( self: Self ) -> None: +# """Vtk filter to ... a vtkUnstructuredGrid. + +# Output mesh is vtkUnstructuredGrid. +# """ +# super().__init__( nInputPorts=1, +# nOutputPorts=1, +# inputType='vtkUnstructuredGrid', +# outputType='vtkUnstructuredGrid' ) +# self.m_paintUnsupportedElementTypes: int = 0 +# # TODO Needs parallelism to work +# # self.m_paintUnsupportedPolyhedrons: int = 0 +# # self.m_chunk_size: int = 1 +# # self.m_num_proc: int = 1 +# self.m_logger: Logger = getLogger( "Element Volumes Filter" ) + +# def FillInputPortInformation( self: Self, port: int, info: vtkInformation ) -> int: +# """Inherited from VTKPythonAlgorithmBase::RequestInformation. + +# Args: +# port (int): input port +# info (vtkInformationVector): info + +# Returns: +# int: 1 if calculation successfully ended, 0 otherwise. +# """ +# if port == 0: +# info.Set( self.INPUT_REQUIRED_DATA_TYPE(), "vtkUnstructuredGrid" ) +# return 1 + +# def RequestInformation( +# self: Self, +# request: vtkInformation, # noqa: F841 +# inInfoVec: list[ vtkInformationVector ], # noqa: F841 +# outInfoVec: vtkInformationVector, +# ) -> int: +# """Inherited from VTKPythonAlgorithmBase::RequestInformation. + +# Args: +# request (vtkInformation): request +# inInfoVec (list[vtkInformationVector]): input objects +# outInfoVec (vtkInformationVector): output objects + +# Returns: +# int: 1 if calculation successfully ended, 0 otherwise. +# """ +# executive = self.GetExecutive() # noqa: F841 +# outInfo = outInfoVec.GetInformationObject( 0 ) # noqa: F841 +# return 1 + +# def RequestData( +# self: Self, +# request: vtkInformation, +# inInfoVec: list[ vtkInformationVector ], +# outInfo: vtkInformationVector +# ) -> int: +# """Inherited from VTKPythonAlgorithmBase::RequestData. + +# Args: +# request (vtkInformation): request +# inInfoVec (list[vtkInformationVector]): input objects +# outInfoVec (vtkInformationVector): output objects + +# Returns: +# int: 1 if calculation successfully ended, 0 otherwise. +# """ +# input_mesh: vtkUnstructuredGrid = vtkUnstructuredGrid.GetData( inInfoVec[ 0 ] ) +# output = vtkUnstructuredGrid.GetData( outInfo ) + +# output_mesh: vtkUnstructuredGrid = input_mesh.NewInstance() +# output_mesh.CopyStructure( input_mesh ) +# output_mesh.CopyAttributes( input_mesh ) + +# unsupported_std_elt_types: set[ int ] = find_unsupported_std_elements_types( input_mesh ) +# if len( unsupported_std_elt_types ) > 0: +# self.m_logger.info( "The following vtk element types in your mesh are not supported by GEOS:" ) +# self.m_logger.info( unsupported_std_elt_types ) + +# if self.m_paintUnsupportedElementTypes: +# nbr_cells: int = output_mesh.GetNumberOfCells() +# arrayCellTypes: npt.NDArray = np.zeros( nbr_cells, dtype=int ) +# for i in range( nbr_cells ): +# arrayCellTypes[ i ] = output_mesh.GetCellType(i) + +# arrayUET: npt.NDArray = np.zeros( nbr_cells, dtype=int ) +# arrayUET[ np.isin( arrayCellTypes, list( unsupported_std_elt_types ) ) ] = 1 +# vtkArrayWSP: vtkDataArray = numpy_to_vtk( arrayUET ) +# vtkArrayWSP.SetName( "HasUnsupportedType" ) +# output_mesh.GetCellData().AddArray( vtkArrayWSP ) + +# # TODO Needs parallelism to work +# # options = Options( self.m_num_proc, self.m_chunk_size ) +# # unsupported_polyhedron_elts: list[ int ] = find_unsupported_polyhedron_elements( input_mesh, options ) +# # if len( unsupported_polyhedron_elts ) > 0: +# # self.m_logger.info( "These vtk polyhedron cell indexes in your mesh are not supported by GEOS:" ) +# # self.m_logger.info( unsupported_polyhedron_elts ) + +# # if self.m_paintUnsupportedPolyhedrons: +# # arrayUP: npt.NDArray = np.zeros( output_mesh.GetNumberOfCells(), dtype=int ) +# # arrayUP[ unsupported_polyhedron_elts ] = 1 +# # self.m_logger.info( f"arrayUP: {arrayUP}" ) +# # vtkArrayWSP: vtkDataArray = numpy_to_vtk( arrayUP ) +# # vtkArrayWSP.SetName( "IsUnsupportedPolyhedron" ) +# # output_mesh.GetCellData().AddArray( vtkArrayWSP ) + +# output.ShallowCopy( output_mesh ) + +# return 1 + +# def SetLogger( self: Self, logger: Logger ) -> None: +# """Set the logger. + +# Args: +# logger (Logger): logger +# """ +# self.m_logger = logger +# self.Modified() + +# def getGrid( self: Self ) -> vtkUnstructuredGrid: +# """Returns the vtkUnstructuredGrid with volumes. + +# Args: +# self (Self) + +# Returns: +# vtkUnstructuredGrid +# """ +# self.Update() # triggers RequestData +# return self.GetOutputDataObject( 0 ) + +# def setPaintUnsupportedElementTypes( self: Self, choice: int ) -> None: +# """Set 0 or 1 to choose if you want to create a new "HasUnsupportedType" array in your output data. + +# Args: +# self (Self) +# choice (int): 0 or 1 +# """ +# if choice not in [ 0, 1 ]: +# self.m_logger.error( f"setPaintUnsupportedElementTypes: Please choose either 0 or 1 not '{choice}'." ) +# else: +# self.m_paintUnsupportedElementTypes = choice +# self.Modified() + +# # TODO Needs parallelism to work +# # def setPaintUnsupportedPolyhedrons( self: Self, choice: int ) -> None: +# # """Set 0 or 1 to choose if you want to create a new "IsUnsupportedPolyhedron" array in your output data. + +# # Args: +# # self (Self) +# # choice (int): 0 or 1 +# # """ +# # if choice not in [ 0, 1 ]: +# # self.m_logger.error( f"setPaintUnsupportedPolyhedrons: Please choose either 0 or 1 not '{choice}'." ) +# # else: +# # self.m_paintUnsupportedPolyhedrons = choice +# # self.Modified() + +# # def setChunkSize( self: Self, new_chunk_size: int ) -> None: +# # self.m_chunk_size = new_chunk_size +# # self.Modified() + +# # def setNumProc( self: Self, new_num_proc: int ) -> None: +# # self.m_num_proc = new_num_proc +# # self.Modified() + +# def writeGrid( self: Self, filepath: str, is_data_mode_binary: bool = True, canOverwrite: bool = False ) -> None: +# """Writes a .vtu file of the vtkUnstructuredGrid at the specified filepath with volumes. + +# Args: +# filepath (str): /path/to/your/file.vtu +# is_data_mode_binary (bool, optional): Writes the file in binary format or ascii. Defaults to True. +# canOverwrite (bool, optional): Allows or not to overwrite if the filepath already leads to an existing +# file. Defaults to False. +# """ +# mesh: vtkUnstructuredGrid = self.getGrid() +# if mesh: +# write_mesh( filepath, VtkOutput( filepath, is_data_mode_binary ), canOverwrite ) +# else: +# self.m_logger.error( f"No output grid was built. Cannot output vtkUnstructuredGrid at {filepath}." ) diff --git a/geos-mesh/src/geos/mesh/doctor/filters/__init__.py b/geos-mesh/src/geos/mesh/doctor/filters/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/geos-mesh/src/geos/mesh/io/vtkIO.py b/geos-mesh/src/geos/mesh/io/vtkIO.py index 1b93648a..dfc8951b 100644 --- a/geos-mesh/src/geos/mesh/io/vtkIO.py +++ b/geos-mesh/src/geos/mesh/io/vtkIO.py @@ -1,16 +1,17 @@ # SPDX-License-Identifier: Apache-2.0 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. # SPDX-FileContributor: Alexandre Benedicto - -import os.path -import logging from dataclasses import dataclass +from enum import Enum +import os.path from typing import Optional -from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkStructuredGrid, vtkPointSet -from vtkmodules.vtkIOLegacy import vtkUnstructuredGridWriter, vtkUnstructuredGridReader -from vtkmodules.vtkIOXML import ( vtkXMLUnstructuredGridReader, vtkXMLUnstructuredGridWriter, - vtkXMLStructuredGridReader, vtkXMLPUnstructuredGridReader, +from vtkmodules.vtkCommonDataModel import vtkPointSet, vtkUnstructuredGrid +from vtkmodules.vtkIOCore import vtkWriter +from vtkmodules.vtkIOLegacy import vtkDataReader, vtkUnstructuredGridWriter, vtkUnstructuredGridReader +from vtkmodules.vtkIOXML import ( vtkXMLDataReader, vtkXMLUnstructuredGridReader, vtkXMLUnstructuredGridWriter, + vtkXMLWriter, vtkXMLStructuredGridReader, vtkXMLPUnstructuredGridReader, vtkXMLPStructuredGridReader, vtkXMLStructuredGridWriter ) +from geos.utils.Logger import getLogger __doc__ = """ Input and Ouput methods for VTK meshes: @@ -18,175 +19,202 @@ - VTK, VTS, VTU writers """ +io_logger = getLogger( "IO for geos-mesh" ) +io_logger.propagate = False + + +class VtkFormat( Enum ): + """Enumeration for supported VTK file formats and their extensions.""" + VTK = ".vtk" + VTS = ".vts" + VTU = ".vtu" + PVTU = ".pvtu" + PVTS = ".pvts" + + +# Centralized mapping of formats to their corresponding reader classes +READER_MAP: dict[ VtkFormat, vtkDataReader | vtkXMLDataReader ] = { + VtkFormat.VTK: vtkUnstructuredGridReader, + VtkFormat.VTS: vtkXMLStructuredGridReader, + VtkFormat.VTU: vtkXMLUnstructuredGridReader, + VtkFormat.PVTU: vtkXMLPUnstructuredGridReader, + VtkFormat.PVTS: vtkXMLPStructuredGridReader +} + +# Centralized mapping of formats to their corresponding writer classes +WRITER_MAP: dict[ VtkFormat, vtkWriter | vtkXMLWriter ] = { + VtkFormat.VTK: vtkUnstructuredGridWriter, + VtkFormat.VTS: vtkXMLStructuredGridWriter, + VtkFormat.VTU: vtkXMLUnstructuredGridWriter, +} + @dataclass( frozen=True ) class VtkOutput: + """Configuration for writing a VTK file.""" output: str - is_data_mode_binary: bool + is_data_mode_binary: bool = 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..." ) - reader.SetFileName( vtk_input_file ) - if reader.IsFileUnstructuredGrid(): - logging.info( f"Reader matches. Reading file \"{vtk_input_file}\" using legacy format reader." ) - reader.Update() - return reader.GetOutput() - else: - logging.info( "Reader did not match the input file format." ) - return None +def _read_data( filepath: str, reader_class: vtkDataReader | vtkXMLDataReader ) -> Optional[ vtkPointSet ]: + """Generic helper to read a VTK file using a specific reader class.""" + reader: vtkDataReader | vtkXMLDataReader = reader_class() + io_logger.info( f"Attempting to read '{filepath}' with {reader_class.__name__}..." ) + # VTK readers have different methods to check file compatibility + can_read: bool = False + if hasattr( reader, 'CanReadFile' ): + can_read = reader.CanReadFile( filepath ) + elif hasattr( reader, 'IsFileUnstructuredGrid' ): # Legacy reader + can_read = reader.IsFileUnstructuredGrid() -def __read_vts( vtk_input_file: str ) -> Optional[ vtkStructuredGrid ]: - reader = vtkXMLStructuredGridReader() - logging.info( f"Testing file format \"{vtk_input_file}\" using XML format reader..." ) - if reader.CanReadFile( vtk_input_file ): - reader.SetFileName( vtk_input_file ) - logging.info( f"Reader matches. Reading file \"{vtk_input_file}\" using XML format reader." ) + if can_read: + reader.SetFileName( filepath ) reader.Update() + io_logger.info( "Read successful." ) return reader.GetOutput() - else: - logging.info( "Reader did not match the input file format." ) - return None + io_logger.info( "Reader did not match the file format." ) + return None -def __read_vtu( vtk_input_file: str ) -> Optional[ vtkUnstructuredGrid ]: - reader = vtkXMLUnstructuredGridReader() - logging.info( f"Testing file format \"{vtk_input_file}\" using XML format reader..." ) - if reader.CanReadFile( vtk_input_file ): - reader.SetFileName( vtk_input_file ) - logging.info( f"Reader matches. Reading file \"{vtk_input_file}\" using XML format reader." ) - reader.Update() - return reader.GetOutput() - else: - logging.info( "Reader did not match the input file format." ) - return None +def _write_data( mesh: vtkPointSet, writer_class: vtkWriter | vtkXMLWriter, output_path: str, is_binary: bool ) -> int: + """Generic helper to write a VTK file using a specific writer class.""" + io_logger.info( f"Writing mesh to '{output_path}' using {writer_class.__name__}..." ) + writer: vtkWriter | vtkXMLWriter = writer_class() + writer.SetFileName( output_path ) + writer.SetInputData( mesh ) -def __read_pvts( vtk_input_file: str ) -> Optional[ vtkStructuredGrid ]: - reader = vtkXMLPStructuredGridReader() - logging.info( f"Testing file format \"{vtk_input_file}\" using XML format reader..." ) - if reader.CanReadFile( vtk_input_file ): - reader.SetFileName( vtk_input_file ) - logging.info( f"Reader matches. Reading file \"{vtk_input_file}\" using XML format reader." ) - reader.Update() - return reader.GetOutput() - else: - logging.info( "Reader did not match the input file format." ) - return None + # Set data mode only for XML writers that support it + if hasattr( writer, 'SetDataModeToBinary' ): + if is_binary: + writer.SetDataModeToBinary() + io_logger.info( "Data mode set to Binary." ) + else: + writer.SetDataModeToAscii() + io_logger.info( "Data mode set to ASCII." ) + return writer.Write() -def __read_pvtu( vtk_input_file: str ) -> Optional[ vtkUnstructuredGrid ]: - reader = vtkXMLPUnstructuredGridReader() - logging.info( f"Testing file format \"{vtk_input_file}\" using XML format reader..." ) - if reader.CanReadFile( vtk_input_file ): - reader.SetFileName( vtk_input_file ) - logging.info( f"Reader matches. Reading file \"{vtk_input_file}\" using XML format reader." ) - reader.Update() - return reader.GetOutput() - else: - logging.info( "Reader did not match the input file format." ) - return None +def read_mesh( filepath: str ) -> vtkPointSet: + """ + Reads a VTK file, automatically detecting the format. -def read_mesh( vtk_input_file: str ) -> vtkPointSet: - """Read vtk file and build either an unstructured grid or a structured grid from it. + It first tries the reader associated with the file extension, then falls + back to trying all available readers if the first attempt fails. Args: - vtk_input_file (str): The file name. Extension will be used to guess file format\ - If first guess fails, other available readers will be tried. + filepath (str): The path to the VTK file. Raises: - ValueError: Invalid file path error - ValueError: No appropriate reader available for the file format + FileNotFoundError: If the input file does not exist. + ValueError: If no suitable reader can be found for the file. Returns: - vtkPointSet: Mesh read + vtkPointSet: The resulting mesh data. """ - if not os.path.exists( 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 ] - extension_to_reader = { - ".vtk": __read_vtk, - ".vts": __read_vts, - ".vtu": __read_vtu, - ".pvtu": __read_pvtu, - ".pvts": __read_pvts - } - # Testing first the reader that should match - if file_extension in extension_to_reader: - output_mesh = extension_to_reader.pop( file_extension )( vtk_input_file ) - if output_mesh: - return output_mesh - # If it does not match, then test all the others. - for reader in extension_to_reader.values(): - output_mesh = reader( vtk_input_file ) - if output_mesh: - return output_mesh - # No reader did work. - err_msg = f"Could not find the appropriate VTK reader for file \"{vtk_input_file}\"." - logging.error( err_msg ) - raise ValueError( err_msg ) - - -def __write_vtk( mesh: vtkUnstructuredGrid, output: str ) -> int: - logging.info( f"Writing mesh into file \"{output}\" using legacy format." ) - writer = vtkUnstructuredGridWriter() - writer.SetFileName( output ) - writer.SetInputData( mesh ) - return writer.Write() + if not os.path.exists( filepath ): + raise FileNotFoundError( f"Invalid file path: '{filepath}' does not exist." ) + _, extension = os.path.splitext( filepath ) + output_mesh: Optional[ vtkPointSet ] = None -def __write_vts( mesh: vtkStructuredGrid, output: str, toBinary: bool = False ) -> int: - logging.info( f"Writing mesh into file \"{output}\" using XML format." ) - writer = vtkXMLStructuredGridWriter() - writer.SetFileName( output ) - writer.SetInputData( mesh ) - writer.SetDataModeToBinary() if toBinary else writer.SetDataModeToAscii() - return writer.Write() + # 1. Try the reader associated with the file extension first + try: + file_format = VtkFormat( extension ) + if file_format in READER_MAP: + reader_class = READER_MAP[ file_format ] + output_mesh = _read_data( filepath, reader_class ) + except ValueError: + io_logger.warning( f"Unknown file extension '{extension}'. Trying all readers." ) + # 2. If the first attempt failed or extension was unknown, try all readers + if not output_mesh: + for reader_class in set( READER_MAP.values() ): # Use set to avoid duplicates + output_mesh = _read_data( filepath, reader_class ) + if output_mesh: + break -def __write_vtu( mesh: vtkUnstructuredGrid, output: str, toBinary: bool = False ) -> int: - logging.info( f"Writing mesh into file \"{output}\" using XML format." ) - writer = vtkXMLUnstructuredGridWriter() - writer.SetFileName( output ) - writer.SetInputData( mesh ) - writer.SetDataModeToBinary() if toBinary else writer.SetDataModeToAscii() - return writer.Write() + if not output_mesh: + raise ValueError( f"Could not find a suitable reader for '{filepath}'." ) + + return output_mesh + + +def read_unstructured_grid( filepath: str ) -> vtkUnstructuredGrid: + """ + Reads a VTK file and ensures it is a vtkUnstructuredGrid. + + This function uses the general `read_mesh` to load the data and then + validates its type. + + Args: + filepath (str): The path to the VTK file. + + Raises: + FileNotFoundError: If the input file does not exist. + ValueError: If no suitable reader can be found for the file. + TypeError: If the file is read successfully but is not a vtkUnstructuredGrid. + + Returns: + vtkUnstructuredGrid: The resulting unstructured grid data. + """ + io_logger.info( f"Reading file '{filepath}' and expecting vtkUnstructuredGrid." ) + + # Reuse the generic mesh reader + mesh = read_mesh( filepath ) + # Check the type of the resulting mesh + if not isinstance( mesh, vtkUnstructuredGrid ): + error_msg = ( f"File '{filepath}' was read successfully, but it is of type " + f"'{type(mesh).__name__}', not the expected vtkUnstructuredGrid." ) + io_logger.error( error_msg ) + raise TypeError( error_msg ) -def write_mesh( mesh: vtkPointSet, vtk_output: VtkOutput, canOverwrite: bool = False ) -> int: - """Write mesh to disk. + io_logger.info( "Validation successful. Mesh is a vtkUnstructuredGrid." ) + return mesh + + +def write_mesh( mesh: vtkPointSet, vtk_output: VtkOutput, can_overwrite: bool = False ) -> int: + """ + Writes a vtkPointSet to a file. - Nothing is done if file already exists. + The format is determined by the file extension in `VtkOutput.output`. Args: - mesh (vtkPointSet): Grid to write - vtk_output (VtkOutput): File path. File extension will be used to select VTK file format - canOverwrite (bool, optional): Authorize overwriting the file. Defaults to False. + mesh (vtkPointSet): The grid data to write. + vtk_output (VtkOutput): Configuration for the output file. + can_overwrite (bool, optional): If False, raises an error if the file + already exists. Defaults to False. Raises: - ValueError: Invalid VTK format. + FileExistsError: If the output file exists and `can_overwrite` is False. + ValueError: If the file extension is not a supported write format. + RuntimeError: If the VTK writer fails to write the file. Returns: - int: 0 if success + int: Returns 1 on success, consistent with the VTK writer's return code. """ - if os.path.exists( vtk_output.output ) and canOverwrite: - logging.error( f"File \"{vtk_output.output}\" already exists, nothing done." ) - return 1 - file_extension = os.path.splitext( vtk_output.output )[ -1 ] - if file_extension == ".vtk": - success_code = __write_vtk( mesh, vtk_output.output ) - elif file_extension == ".vts": - success_code = __write_vts( mesh, vtk_output.output, vtk_output.is_data_mode_binary ) - 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. - 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. + if os.path.exists( vtk_output.output ) and not can_overwrite: + raise FileExistsError( f"File '{vtk_output.output}' already exists. Set can_overwrite=True to replace it." ) + + _, extension = os.path.splitext( vtk_output.output ) + + try: + file_format = VtkFormat( extension ) + if file_format not in WRITER_MAP: + raise ValueError( f"Writing to extension '{extension}' is not supported." ) + + writer_class = WRITER_MAP[ file_format ] + success_code = _write_data( mesh, writer_class, vtk_output.output, vtk_output.is_data_mode_binary ) + + if not success_code: + raise RuntimeError( f"VTK writer failed to write file '{vtk_output.output}'." ) + + io_logger.info( f"Successfully wrote mesh to '{vtk_output.output}'." ) + return success_code # VTK writers return 1 for success + + except ValueError as e: + io_logger.error( e ) + raise diff --git a/geos-mesh/tests/test_Checks.py b/geos-mesh/tests/test_Checks.py new file mode 100644 index 00000000..5661c068 --- /dev/null +++ b/geos-mesh/tests/test_Checks.py @@ -0,0 +1,361 @@ +import pytest +from vtkmodules.vtkCommonCore import vtkPoints +from vtkmodules.vtkCommonDataModel import ( vtkCellArray, vtkTetra, vtkUnstructuredGrid, VTK_TETRA ) +from geos.mesh.doctor.filters.Checks import AllChecks, MainChecks + + +@pytest.fixture +def simple_mesh_with_issues() -> vtkUnstructuredGrid: + """Create a simple test mesh with known issues for testing checks. + + This mesh includes: + - Collocated nodes (points 0 and 3 are at the same location) + - A very small volume element + - Wrong support elements (duplicate node indices in cells) + + Returns: + vtkUnstructuredGrid: Test mesh with various issues + """ + mesh = vtkUnstructuredGrid() + + # Create points with some collocated nodes + points = vtkPoints() + points.InsertNextPoint( 0.0, 0.0, 0.0 ) # Point 0 + points.InsertNextPoint( 1.0, 0.0, 0.0 ) # Point 1 + points.InsertNextPoint( 0.0, 1.0, 0.0 ) # Point 2 + points.InsertNextPoint( 0.0, 0.0, 0.0 ) # Point 3 - duplicate of Point 0 + points.InsertNextPoint( 0.0, 0.0, 1.0 ) # Point 4 + points.InsertNextPoint( 2.0, 0.0, 0.0 ) # Point 5 + points.InsertNextPoint( 2.01, 0.0, 0.0 ) # Point 6 - very close to Point 5 (small volume) + points.InsertNextPoint( 2.0, 0.01, 0.0 ) # Point 7 - creates tiny element + points.InsertNextPoint( 2.0, 0.0, 0.01 ) # Point 8 - creates tiny element + mesh.SetPoints( points ) + + # Create cells + cells = vtkCellArray() + cell_types = [] + + # Normal tetrahedron + tet1 = vtkTetra() + tet1.GetPointIds().SetId( 0, 0 ) + tet1.GetPointIds().SetId( 1, 1 ) + tet1.GetPointIds().SetId( 2, 2 ) + tet1.GetPointIds().SetId( 3, 4 ) + cells.InsertNextCell( tet1 ) + cell_types.append( VTK_TETRA ) + + # Tetrahedron with duplicate node indices (wrong support) + tet2 = vtkTetra() + tet2.GetPointIds().SetId( 0, 3 ) # This is collocated with point 0 + tet2.GetPointIds().SetId( 1, 1 ) + tet2.GetPointIds().SetId( 2, 2 ) + tet2.GetPointIds().SetId( 3, 0 ) # Duplicate reference to same location + cells.InsertNextCell( tet2 ) + cell_types.append( VTK_TETRA ) + + # Very small volume tetrahedron + tet3 = vtkTetra() + tet3.GetPointIds().SetId( 0, 5 ) + tet3.GetPointIds().SetId( 1, 6 ) + tet3.GetPointIds().SetId( 2, 7 ) + tet3.GetPointIds().SetId( 3, 8 ) + cells.InsertNextCell( tet3 ) + cell_types.append( VTK_TETRA ) + + mesh.SetCells( cell_types, cells ) + return mesh + + +@pytest.fixture +def clean_mesh() -> vtkUnstructuredGrid: + """Create a clean test mesh without issues. + + Returns: + vtkUnstructuredGrid: Clean test mesh + """ + mesh = vtkUnstructuredGrid() + + # Create well-separated points + points = vtkPoints() + points.InsertNextPoint( 0.0, 0.0, 0.0 ) # Point 0 + points.InsertNextPoint( 1.0, 0.0, 0.0 ) # Point 1 + points.InsertNextPoint( 0.0, 1.0, 0.0 ) # Point 2 + points.InsertNextPoint( 0.0, 0.0, 1.0 ) # Point 3 + mesh.SetPoints( points ) + + # Create a single clean tetrahedron + cells = vtkCellArray() + cell_types = [] + + tet = vtkTetra() + tet.GetPointIds().SetId( 0, 0 ) + tet.GetPointIds().SetId( 1, 1 ) + tet.GetPointIds().SetId( 2, 2 ) + tet.GetPointIds().SetId( 3, 3 ) + cells.InsertNextCell( tet ) + cell_types.append( VTK_TETRA ) + + mesh.SetCells( cell_types, cells ) + return mesh + + +@pytest.fixture +def all_checks_filter() -> AllChecks: + """Create a fresh AllChecks filter for each test.""" + return AllChecks() + + +@pytest.fixture +def main_checks_filter() -> MainChecks: + """Create a fresh MainChecks filter for each test.""" + return MainChecks() + + +class TestAllChecks: + """Test class for AllChecks filter functionality.""" + + def test_filter_creation( self, all_checks_filter: AllChecks ): + """Test that AllChecks filter can be created successfully.""" + assert all_checks_filter is not None + assert hasattr( all_checks_filter, 'getAvailableChecks' ) + assert hasattr( all_checks_filter, 'setChecksToPerform' ) + assert hasattr( all_checks_filter, 'setCheckParameter' ) + + def test_available_checks( self, all_checks_filter: AllChecks ): + """Test that all expected checks are available.""" + available_checks = all_checks_filter.getAvailableChecks() + + # Check that we have the expected checks for AllChecks + expected_checks = [ + 'collocated_nodes', 'element_volumes', 'non_conformal', 'self_intersecting_elements', 'supported_elements' + ] + + for check in expected_checks: + assert check in available_checks, f"Check '{check}' should be available" + + def test_default_parameters( self, all_checks_filter: AllChecks ): + """Test that default parameters are correctly retrieved.""" + available_checks = all_checks_filter.getAvailableChecks() + + for check_name in available_checks: + defaults = all_checks_filter.getDefaultParameters( check_name ) + assert isinstance( defaults, dict ), f"Default parameters for '{check_name}' should be a dict" + + # Test specific known defaults + collocated_defaults = all_checks_filter.getDefaultParameters( 'collocated_nodes' ) + assert 'tolerance' in collocated_defaults + + volume_defaults = all_checks_filter.getDefaultParameters( 'element_volumes' ) + assert 'min_volume' in volume_defaults + + def test_set_checks_to_perform( self, all_checks_filter: AllChecks ): + """Test setting specific checks to perform.""" + # Set specific checks + checks_to_perform = [ 'collocated_nodes', 'element_volumes' ] + all_checks_filter.setChecksToPerform( checks_to_perform ) + + # Verify by checking if the filter state changed + assert hasattr( all_checks_filter, 'm_checks_to_perform' ) + assert all_checks_filter.m_checks_to_perform == checks_to_perform + + def test_set_check_parameter( self, all_checks_filter: AllChecks ): + """Test setting parameters for specific checks.""" + # Set a tolerance parameter for collocated nodes + all_checks_filter.setCheckParameter( 'collocated_nodes', 'tolerance', 1e-6 ) + + # Set minimum volume for element volumes + all_checks_filter.setCheckParameter( 'element_volumes', 'min_volume', 0.001 ) + + # Verify parameters are stored + assert 'collocated_nodes' in all_checks_filter.m_check_parameters + assert all_checks_filter.m_check_parameters[ 'collocated_nodes' ][ 'tolerance' ] == 1e-6 + assert all_checks_filter.m_check_parameters[ 'element_volumes' ][ 'min_volume' ] == 0.001 + + def test_set_all_checks_parameter( self, all_checks_filter: AllChecks ): + """Test setting a parameter that applies to all compatible checks.""" + # Set tolerance for all checks that support it + all_checks_filter.setAllChecksParameter( 'tolerance', 1e-8 ) + + # Check that tolerance was set for checks that support it + if 'collocated_nodes' in all_checks_filter.m_check_parameters: + assert all_checks_filter.m_check_parameters[ 'collocated_nodes' ][ 'tolerance' ] == 1e-8 + + def test_process_mesh_with_issues( self, all_checks_filter: AllChecks, + simple_mesh_with_issues: vtkUnstructuredGrid ): + """Test processing a mesh with known issues.""" + # Configure for specific checks + all_checks_filter.setChecksToPerform( [ 'collocated_nodes', 'element_volumes' ] ) + all_checks_filter.setCheckParameter( 'collocated_nodes', 'tolerance', 1e-12 ) + all_checks_filter.setCheckParameter( 'element_volumes', 'min_volume', 1e-3 ) + + # Process the mesh + all_checks_filter.SetInputDataObject( 0, simple_mesh_with_issues ) + all_checks_filter.Update() + + # Check results + results = all_checks_filter.getCheckResults() + + assert 'collocated_nodes' in results + assert 'element_volumes' in results + + # Check that collocated nodes were found + collocated_result = results[ 'collocated_nodes' ] + assert hasattr( collocated_result, 'nodes_buckets' ) + # We expect to find collocated nodes (points 0 and 3) + assert len( collocated_result.nodes_buckets ) > 0 + + # Check that volume issues were detected + volume_result = results[ 'element_volumes' ] + assert hasattr( volume_result, 'element_volumes' ) + + def test_process_clean_mesh( self, all_checks_filter: AllChecks, clean_mesh: vtkUnstructuredGrid ): + """Test processing a clean mesh without issues.""" + # Configure checks + all_checks_filter.setChecksToPerform( [ 'collocated_nodes', 'element_volumes' ] ) + all_checks_filter.setCheckParameter( 'collocated_nodes', 'tolerance', 1e-12 ) + all_checks_filter.setCheckParameter( 'element_volumes', 'min_volume', 1e-6 ) + + # Process the mesh + all_checks_filter.SetInputDataObject( 0, clean_mesh ) + all_checks_filter.Update() + + # Check results + results = all_checks_filter.getCheckResults() + + assert 'collocated_nodes' in results + assert 'element_volumes' in results + + # Check that no issues were found + collocated_result = results[ 'collocated_nodes' ] + assert len( collocated_result.nodes_buckets ) == 0 + + volume_result = results[ 'element_volumes' ] + assert len( volume_result.element_volumes ) == 0 + + def test_output_mesh_unchanged( self, all_checks_filter: AllChecks, clean_mesh: vtkUnstructuredGrid ): + """Test that the output mesh is unchanged from the input (checks don't modify geometry).""" + original_num_points = clean_mesh.GetNumberOfPoints() + original_num_cells = clean_mesh.GetNumberOfCells() + + # Process the mesh + all_checks_filter.SetInputDataObject( 0, clean_mesh ) + all_checks_filter.Update() + + # Get output mesh + output_mesh = all_checks_filter.getGrid() + + # Verify structure is unchanged + assert output_mesh.GetNumberOfPoints() == original_num_points + assert output_mesh.GetNumberOfCells() == original_num_cells + + # Verify points are the same + for i in range( original_num_points ): + original_point = clean_mesh.GetPoint( i ) + output_point = output_mesh.GetPoint( i ) + assert original_point == output_point + + +class TestMainChecks: + """Test class for MainChecks filter functionality.""" + + def test_filter_creation( self, main_checks_filter: MainChecks ): + """Test that MainChecks filter can be created successfully.""" + assert main_checks_filter is not None + assert hasattr( main_checks_filter, 'getAvailableChecks' ) + assert hasattr( main_checks_filter, 'setChecksToPerform' ) + assert hasattr( main_checks_filter, 'setCheckParameter' ) + + def test_available_checks( self, main_checks_filter: MainChecks ): + """Test that main checks are available (subset of all checks).""" + available_checks = main_checks_filter.getAvailableChecks() + + # MainChecks should have a subset of checks + expected_main_checks = [ 'collocated_nodes', 'element_volumes', 'self_intersecting_elements' ] + + for check in expected_main_checks: + assert check in available_checks, f"Main check '{check}' should be available" + + def test_process_mesh( self, main_checks_filter: MainChecks, simple_mesh_with_issues: vtkUnstructuredGrid ): + """Test processing a mesh with MainChecks.""" + # Process the mesh with default configuration + main_checks_filter.SetInputDataObject( 0, simple_mesh_with_issues ) + main_checks_filter.Update() + + # Check that results are obtained + results = main_checks_filter.getCheckResults() + assert isinstance( results, dict ) + assert len( results ) > 0 + + # Check that main checks were performed + available_checks = main_checks_filter.getAvailableChecks() + for check_name in available_checks: + if check_name in results: + result = results[ check_name ] + assert result is not None + + +class TestFilterComparison: + """Test class for comparing AllChecks and MainChecks filters.""" + + def test_all_checks_vs_main_checks_availability( self, all_checks_filter: AllChecks, + main_checks_filter: MainChecks ): + """Test that MainChecks is a subset of AllChecks.""" + all_checks = set( all_checks_filter.getAvailableChecks() ) + main_checks = set( main_checks_filter.getAvailableChecks() ) + + # MainChecks should be a subset of AllChecks + assert main_checks.issubset( all_checks ), "MainChecks should be a subset of AllChecks" + + # AllChecks should have more checks than MainChecks + assert len( all_checks ) >= len( main_checks ), "AllChecks should have at least as many checks as MainChecks" + + def test_parameter_consistency( self, all_checks_filter: AllChecks, main_checks_filter: MainChecks ): + """Test that parameter handling is consistent between filters.""" + # Get common checks + all_checks = set( all_checks_filter.getAvailableChecks() ) + main_checks = set( main_checks_filter.getAvailableChecks() ) + common_checks = all_checks.intersection( main_checks ) + + # Test that default parameters are the same for common checks + for check_name in common_checks: + all_defaults = all_checks_filter.getDefaultParameters( check_name ) + main_defaults = main_checks_filter.getDefaultParameters( check_name ) + assert all_defaults == main_defaults, f"Default parameters should be the same for '{check_name}'" + + +class TestErrorHandling: + """Test class for error handling and edge cases.""" + + def test_invalid_check_name( self, all_checks_filter: AllChecks ): + """Test handling of invalid check names.""" + # Try to set an invalid check + invalid_checks = [ 'nonexistent_check' ] + all_checks_filter.setChecksToPerform( invalid_checks ) + + # The filter should handle this gracefully + # (The actual behavior depends on implementation - it might warn or ignore) + assert all_checks_filter.m_checks_to_perform == invalid_checks + + def test_invalid_parameter_name( self, all_checks_filter: AllChecks ): + """Test handling of invalid parameter names.""" + # Try to set an invalid parameter + all_checks_filter.setCheckParameter( 'collocated_nodes', 'invalid_param', 123 ) + + # This should not crash the filter + assert 'collocated_nodes' in all_checks_filter.m_check_parameters + assert 'invalid_param' in all_checks_filter.m_check_parameters[ 'collocated_nodes' ] + + def test_empty_mesh( self, all_checks_filter: AllChecks ): + """Test handling of empty mesh.""" + # Create an empty mesh + empty_mesh = vtkUnstructuredGrid() + empty_mesh.SetPoints( vtkPoints() ) + + # Process the empty mesh + all_checks_filter.setChecksToPerform( [ 'collocated_nodes' ] ) + all_checks_filter.SetInputDataObject( 0, empty_mesh ) + all_checks_filter.Update() + + # Should complete without error + results = all_checks_filter.getCheckResults() + assert isinstance( results, dict ) diff --git a/geos-mesh/tests/test_collocated_nodes.py b/geos-mesh/tests/test_collocated_nodes.py index 86f798f7..c4de479b 100644 --- a/geos-mesh/tests/test_collocated_nodes.py +++ b/geos-mesh/tests/test_collocated_nodes.py @@ -2,7 +2,8 @@ from typing import Iterator, Tuple from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkTetra, vtkUnstructuredGrid, VTK_TETRA -from geos.mesh.doctor.actions.collocated_nodes import Options, __action +from geos.mesh.doctor.actions.collocated_nodes import Options, mesh_action +from geos.mesh.doctor.filters.CollocatedNodes import CollocatedNodes def get_points() -> Iterator[ Tuple[ vtkPoints, int ] ]: @@ -27,7 +28,7 @@ def test_simple_collocated_points( data: Tuple[ vtkPoints, int ] ): mesh = vtkUnstructuredGrid() mesh.SetPoints( points ) - result = __action( mesh, Options( tolerance=1.e-12 ) ) + result = mesh_action( mesh, Options( tolerance=1.e-12 ) ) assert len( result.wrong_support_elements ) == 0 assert len( result.nodes_buckets ) == num_nodes_bucket @@ -58,8 +59,152 @@ def test_wrong_support_elements(): mesh.SetPoints( points ) mesh.SetCells( cell_types, cells ) - result = __action( mesh, Options( tolerance=1.e-12 ) ) + result = mesh_action( mesh, Options( tolerance=1.e-12 ) ) assert len( result.nodes_buckets ) == 0 assert len( result.wrong_support_elements ) == 1 assert result.wrong_support_elements[ 0 ] == 0 + + +# Create a test mesh with collocated nodes +@pytest.fixture +def sample_mesh() -> vtkUnstructuredGrid: + # Create a simple mesh with duplicate points + mesh = vtkUnstructuredGrid() + + # Create points + points = vtkPoints() + points.InsertNextPoint( 0.0, 0.0, 0.0 ) # Point 0 + points.InsertNextPoint( 1.0, 0.0, 0.0 ) # Point 1 + points.InsertNextPoint( 0.0, 1.0, 0.0 ) # Point 2 + points.InsertNextPoint( 0.0, 0.0, 0.0 ) # Point 3 - duplicate of Point 0 + points.InsertNextPoint( 2.0, 0.0, 0.0 ) # Point 4 + mesh.SetPoints( points ) + + # Create cells + cells = vtkCellArray() + # Create a triangular cell with normal connectivity + cells.InsertNextCell( 3 ) + cells.InsertCellPoint( 0 ) + cells.InsertCellPoint( 1 ) + cells.InsertCellPoint( 2 ) + + # Create a cell with duplicate point indices + cells.InsertNextCell( 3 ) + cells.InsertCellPoint( 3 ) # This is a duplicate of point 0 + cells.InsertCellPoint( 1 ) + cells.InsertCellPoint( 4 ) + mesh.SetCells( 5, cells ) # 5 is VTK_TRIANGLE + return mesh + + +@pytest.fixture +def collocated_nodes_filter() -> CollocatedNodes: + return CollocatedNodes() + + +def test_init( collocated_nodes_filter: CollocatedNodes ): + """Test initialization of the CollocatedNodes filter.""" + assert collocated_nodes_filter.m_collocatedNodesBuckets == list() + assert collocated_nodes_filter.m_paintWrongSupportElements == 0 + assert collocated_nodes_filter.m_tolerance == 0.0 + assert collocated_nodes_filter.m_wrongSupportElements == list() + + +def test_collocated_nodes_detection( sample_mesh: vtkUnstructuredGrid, collocated_nodes_filter: CollocatedNodes ): + """Test the filter's ability to detect collocated nodes.""" + # Set input mesh + collocated_nodes_filter.SetInputDataObject( sample_mesh ) + + # Set tolerance + collocated_nodes_filter.setTolerance( 1e-6 ) + + # Run filter + collocated_nodes_filter.Update() + + # Check results + buckets = collocated_nodes_filter.getCollocatedNodeBuckets() + assert len( buckets ) > 0 + + # We expect points 0 and 3 to be in the same bucket + bucket_with_duplicates = None + for bucket in buckets: + if 0 in bucket and 3 in bucket: + bucket_with_duplicates = bucket + break + + assert bucket_with_duplicates is not None, "Failed to detect collocated nodes 0 and 3" + + +def test_wrong_support_elements2( sample_mesh: vtkUnstructuredGrid, collocated_nodes_filter: CollocatedNodes ): + """Test the filter's ability to detect elements with wrong support.""" + # Set input mesh + collocated_nodes_filter.SetInputDataObject( sample_mesh ) + + # Run filter + collocated_nodes_filter.Update() + + # Check results + wrong_elements = collocated_nodes_filter.getWrongSupportElements() + + # In our test mesh, we don't have cells with duplicate point indices within the same cell + # So this should be empty + assert isinstance( wrong_elements, list ) + + +def test_paint_wrong_support_elements( sample_mesh: vtkUnstructuredGrid, collocated_nodes_filter: CollocatedNodes ): + """Test the painting of wrong support elements.""" + # Set input mesh + collocated_nodes_filter.SetInputDataObject( sample_mesh ) + + # Enable painting + collocated_nodes_filter.setPaintWrongSupportElements( 1 ) + + # Run filter + collocated_nodes_filter.Update() + + # Get output mesh + output_mesh = collocated_nodes_filter.getGrid() + + # Check if the array was added + cell_data = output_mesh.GetCellData() + has_array = cell_data.HasArray( "HasDuplicatedNodes" ) + assert has_array, "The HasDuplicatedNodes array wasn't added to cell data" + + +def test_set_paint_wrong_support_elements( collocated_nodes_filter: CollocatedNodes ): + """Test setPaintWrongSupportElements method.""" + # Valid input + collocated_nodes_filter.setPaintWrongSupportElements( 1 ) + assert collocated_nodes_filter.m_paintWrongSupportElements == 1 + + # Valid input + collocated_nodes_filter.setPaintWrongSupportElements( 0 ) + assert collocated_nodes_filter.m_paintWrongSupportElements == 0 + + # Invalid input + collocated_nodes_filter.setPaintWrongSupportElements( 2 ) + # Should remain unchanged + assert collocated_nodes_filter.m_paintWrongSupportElements == 0 + + +def test_set_tolerance( collocated_nodes_filter: CollocatedNodes ): + """Test setTolerance method.""" + collocated_nodes_filter.setTolerance( 0.001 ) + assert collocated_nodes_filter.m_tolerance == 0.001 + + +def test_get_collocated_node_buckets( collocated_nodes_filter: CollocatedNodes ): + """Test getCollocatedNodeBuckets method.""" + # Set a value manually + collocated_nodes_filter.m_collocatedNodesBuckets = [ ( 0, 1 ), ( 2, 3 ) ] + result = collocated_nodes_filter.getCollocatedNodeBuckets() + assert result == [ ( 0, 1 ), ( 2, 3 ) ] + + +def test_get_wrong_support_elements( collocated_nodes_filter: CollocatedNodes ): + """Test getWrongSupportElements method.""" + # Set a value manually + collocated_nodes_filter.m_wrongSupportElements = [ 0, 3, 5 ] + result = collocated_nodes_filter.getWrongSupportElements() + assert result == [ 0, 3, 5 ] diff --git a/geos-mesh/tests/test_element_volumes.py b/geos-mesh/tests/test_element_volumes.py index dccbda93..64b903bd 100644 --- a/geos-mesh/tests/test_element_volumes.py +++ b/geos-mesh/tests/test_element_volumes.py @@ -1,7 +1,228 @@ -import numpy -from vtkmodules.vtkCommonCore import vtkPoints -from vtkmodules.vtkCommonDataModel import VTK_TETRA, vtkCellArray, vtkTetra, vtkUnstructuredGrid -from geos.mesh.doctor.actions.element_volumes import Options, __action +import numpy as np +import numpy.typing as npt +import pytest +from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkHexahedron, vtkTetra, vtkUnstructuredGrid, VTK_TETRA +from vtkmodules.vtkCommonCore import vtkPoints, vtkIdList +from vtkmodules.util.numpy_support import vtk_to_numpy +from geos.mesh.doctor.actions.element_volumes import Options, mesh_action +from geos.mesh.doctor.filters.ElementVolumes import ElementVolumes + + +@pytest.fixture +def tetra_mesh() -> vtkUnstructuredGrid: + """Create a simple tetrahedron with known volume (1/6)""" + points = vtkPoints() + points.InsertNextPoint( 0, 0, 0 ) # Point 0 + points.InsertNextPoint( 1, 0, 0 ) # Point 1 + points.InsertNextPoint( 0, 1, 0 ) # Point 2 + points.InsertNextPoint( 0, 0, 1 ) # Point 3 + + tetra = vtkTetra() + tetra.GetPointIds().SetId( 0, 0 ) + tetra.GetPointIds().SetId( 1, 1 ) + tetra.GetPointIds().SetId( 2, 2 ) + tetra.GetPointIds().SetId( 3, 3 ) + + ug = vtkUnstructuredGrid() + ug.SetPoints( points ) + ug.InsertNextCell( tetra.GetCellType(), tetra.GetPointIds() ) + return ug + + +@pytest.fixture +def hexa_mesh() -> vtkUnstructuredGrid: + """Create a simple hexahedron with known volume (1.0)""" + points = vtkPoints() + points.InsertNextPoint( 0, 0, 0 ) # Point 0 + points.InsertNextPoint( 1, 0, 0 ) # Point 1 + points.InsertNextPoint( 1, 1, 0 ) # Point 2 + points.InsertNextPoint( 0, 1, 0 ) # Point 3 + points.InsertNextPoint( 0, 0, 1 ) # Point 4 + points.InsertNextPoint( 1, 0, 1 ) # Point 5 + points.InsertNextPoint( 1, 1, 1 ) # Point 6 + points.InsertNextPoint( 0, 1, 1 ) # Point 7 + + hexa = vtkHexahedron() + for i in range( 8 ): + hexa.GetPointIds().SetId( i, i ) + + ug = vtkUnstructuredGrid() + ug.SetPoints( points ) + ug.InsertNextCell( hexa.GetCellType(), hexa.GetPointIds() ) + return ug + + +@pytest.fixture +def negative_vol_mesh() -> vtkUnstructuredGrid: + """Create a tetrahedron with negative volume (wrong winding)""" + points = vtkPoints() + points.InsertNextPoint( 0, 0, 0 ) # Point 0 + points.InsertNextPoint( 1, 0, 0 ) # Point 1 + points.InsertNextPoint( 0, 1, 0 ) # Point 2 + points.InsertNextPoint( 0, 0, 1 ) # Point 3 + + tetra = vtkTetra() + # Switch two points to create negative volume + tetra.GetPointIds().SetId( 0, 0 ) + tetra.GetPointIds().SetId( 1, 2 ) # Swapped from normal order + tetra.GetPointIds().SetId( 2, 1 ) # Swapped from normal order + tetra.GetPointIds().SetId( 3, 3 ) + + ug = vtkUnstructuredGrid() + ug.SetPoints( points ) + ug.InsertNextCell( tetra.GetCellType(), tetra.GetPointIds() ) + return ug + + +@pytest.fixture +def zero_vol_mesh() -> vtkUnstructuredGrid: + """Create a tetrahedron with zero volume (coplanar points)""" + points = vtkPoints() + points.InsertNextPoint( 0, 0, 0 ) # Point 0 + points.InsertNextPoint( 1, 0, 0 ) # Point 1 + points.InsertNextPoint( 0, 1, 0 ) # Point 2 + points.InsertNextPoint( 1, 1, 0 ) # Point 3 (coplanar with others) + + tetra = vtkTetra() + tetra.GetPointIds().SetId( 0, 0 ) + tetra.GetPointIds().SetId( 1, 1 ) + tetra.GetPointIds().SetId( 2, 2 ) + tetra.GetPointIds().SetId( 3, 3 ) + + ug = vtkUnstructuredGrid() + ug.SetPoints( points ) + ug.InsertNextCell( tetra.GetCellType(), tetra.GetPointIds() ) + return ug + + +@pytest.fixture +def volume_filter() -> ElementVolumes: + """Create a fresh ElementVolumes filter for each test""" + return ElementVolumes() + + +def test_tetrahedron_volume( tetra_mesh: vtkUnstructuredGrid, volume_filter: ElementVolumes ) -> None: + """Test volume calculation for a regular tetrahedron""" + volume_filter.SetInputDataObject( 0, tetra_mesh ) + volume_filter.Update() + output: vtkUnstructuredGrid = volume_filter.getGrid() + + volumes: npt.NDArray = vtk_to_numpy( output.GetCellData().GetArray( "MESH_DOCTOR_VOLUME" ) ) + expected_volume: float = 1 / 6 # Tetrahedron volume + + assert len( volumes ) == 1 + assert volumes[ 0 ] == pytest.approx( expected_volume, abs=1e-6 ) + + +def test_hexahedron_volume( hexa_mesh: vtkUnstructuredGrid, volume_filter: ElementVolumes ) -> None: + """Test volume calculation for a regular hexahedron""" + volume_filter.SetInputDataObject( 0, hexa_mesh ) + volume_filter.Update() + output: vtkUnstructuredGrid = volume_filter.getGrid() + + volumes: npt.NDArray = vtk_to_numpy( output.GetCellData().GetArray( "MESH_DOCTOR_VOLUME" ) ) + expected_volume: float = 1.0 # Unit cube volume + + assert len( volumes ) == 1 + assert volumes[ 0 ] == pytest.approx( expected_volume, abs=1e-6 ) + + +def test_negative_volume_detection( negative_vol_mesh: vtkUnstructuredGrid, volume_filter: ElementVolumes ) -> None: + """Test detection of negative volumes""" + volume_filter.SetInputDataObject( 0, negative_vol_mesh ) + volume_filter.setReturnNegativeZeroVolumes( True ) + volume_filter.Update() + + output: vtkUnstructuredGrid = volume_filter.getGrid() + volumes: npt.NDArray = vtk_to_numpy( output.GetCellData().GetArray( "MESH_DOCTOR_VOLUME" ) ) + + assert len( volumes ) == 1 + assert volumes[ 0 ] < 0 + + # Test getNegativeZeroVolumes method + negative_zero_volumes: npt.NDArray = volume_filter.getNegativeZeroVolumes() + assert len( negative_zero_volumes ) == 1 + assert negative_zero_volumes[ 0, 0 ] == 0 # First cell index + assert negative_zero_volumes[ 0, 1 ] == volumes[ 0 ] # Volume value + + +def test_zero_volume_detection( zero_vol_mesh: vtkUnstructuredGrid, volume_filter: ElementVolumes ) -> None: + """Test detection of zero volumes""" + volume_filter.SetInputDataObject( 0, zero_vol_mesh ) + volume_filter.setReturnNegativeZeroVolumes( True ) + volume_filter.Update() + + output: vtkUnstructuredGrid = volume_filter.getGrid() + volumes: npt.NDArray = vtk_to_numpy( output.GetCellData().GetArray( "MESH_DOCTOR_VOLUME" ) ) + + assert len( volumes ) == 1 + assert volumes[ 0 ] == pytest.approx( 0.0, abs=1e-6 ) + + # Test getNegativeZeroVolumes method + negative_zero_volumes: npt.NDArray = volume_filter.getNegativeZeroVolumes() + assert len( negative_zero_volumes ) == 1 + assert negative_zero_volumes[ 0, 0 ] == 0 # First cell index + assert negative_zero_volumes[ 0, 1 ] == pytest.approx( 0.0, abs=1e-6 ) # Volume value + + +def test_return_negative_zero_volumes_flag( volume_filter: ElementVolumes ) -> None: + """Test setting and getting the returnNegativeZeroVolumes flag""" + # Default should be False + assert not volume_filter.m_returnNegativeZeroVolumes + + # Set to True and verify + volume_filter.setReturnNegativeZeroVolumes( True ) + assert volume_filter.m_returnNegativeZeroVolumes + + # Set to False and verify + volume_filter.setReturnNegativeZeroVolumes( False ) + assert not volume_filter.m_returnNegativeZeroVolumes + + +def test_mixed_mesh( tetra_mesh: vtkUnstructuredGrid, hexa_mesh: vtkUnstructuredGrid, + volume_filter: ElementVolumes ) -> None: + """Test with a combined mesh containing multiple element types""" + # Create a mixed mesh with both tet and hex + mixed_mesh = vtkUnstructuredGrid() + + # Copy points from tetra_mesh + tetra_points: vtkPoints = tetra_mesh.GetPoints() + points = vtkPoints() + for i in range( tetra_points.GetNumberOfPoints() ): + points.InsertNextPoint( tetra_points.GetPoint( i ) ) + + # Add points from hexa_mesh with offset + hexa_points: vtkPoints = hexa_mesh.GetPoints() + offset: int = points.GetNumberOfPoints() + for i in range( hexa_points.GetNumberOfPoints() ): + x, y, z = hexa_points.GetPoint( i ) + points.InsertNextPoint( x + 2, y, z ) # Offset in x-direction + + mixed_mesh.SetPoints( points ) + + # Add tetra cell + tetra_cell: vtkTetra = tetra_mesh.GetCell( 0 ) + ids: vtkIdList = tetra_cell.GetPointIds() + mixed_mesh.InsertNextCell( tetra_cell.GetCellType(), ids.GetNumberOfIds(), + [ ids.GetId( i ) for i in range( ids.GetNumberOfIds() ) ] ) + + # Add hexa cell with offset ids + hexa_cell: vtkHexahedron = hexa_mesh.GetCell( 0 ) + ids: vtkIdList = hexa_cell.GetPointIds() + hexa_ids: list[ int ] = [ ids.GetId( i ) + offset for i in range( ids.GetNumberOfIds() ) ] + mixed_mesh.InsertNextCell( hexa_cell.GetCellType(), len( hexa_ids ), hexa_ids ) + + # Apply filter + volume_filter.SetInputDataObject( 0, mixed_mesh ) + volume_filter.Update() + output: vtkUnstructuredGrid = volume_filter.getGrid() + + # Check volumes + volumes: npt.NDArray = vtk_to_numpy( output.GetCellData().GetArray( "MESH_DOCTOR_VOLUME" ) ) + + assert len( volumes ) == 2 + assert volumes[ 0 ] == pytest.approx( 1 / 6, abs=1e-6 ) # Tetrahedron volume + assert volumes[ 1 ] == pytest.approx( 1.0, abs=1e-6 ) # Hexahedron volume def test_simple_tet(): @@ -28,12 +249,12 @@ def test_simple_tet(): mesh.SetPoints( points ) mesh.SetCells( cell_types, cells ) - result = __action( mesh, Options( min_volume=1. ) ) + result = mesh_action( mesh, Options( min_volume=1. ) ) assert len( result.element_volumes ) == 1 assert result.element_volumes[ 0 ][ 0 ] == 0 - assert abs( result.element_volumes[ 0 ][ 1 ] - 1. / 6. ) < 10 * numpy.finfo( float ).eps + assert abs( result.element_volumes[ 0 ][ 1 ] - 1. / 6. ) < 10 * np.finfo( float ).eps - result = __action( mesh, Options( min_volume=0. ) ) + result = mesh_action( mesh, Options( min_volume=0. ) ) assert len( result.element_volumes ) == 0 diff --git a/geos-mesh/tests/test_generate_cube.py b/geos-mesh/tests/test_generate_cube.py index d02ef68b..22d0ef06 100644 --- a/geos-mesh/tests/test_generate_cube.py +++ b/geos-mesh/tests/test_generate_cube.py @@ -1,4 +1,47 @@ -from geos.mesh.doctor.actions.generate_cube import __build, Options, FieldInfo +import pytest +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, vtkPointData, vtkCellData +from vtkmodules.vtkCommonCore import vtkDataArray +from geos.mesh.doctor.actions.generate_cube import FieldInfo, Options, __build +from geos.mesh.doctor.filters.GenerateRectilinearGrid import GenerateRectilinearGrid + + +@pytest.fixture +def generate_rectilinear_grid_filter() -> GenerateRectilinearGrid: + filter = GenerateRectilinearGrid() + filter.setCoordinates( [ 0.0, 5.0, 10.0 ], [ 0.0, 10.0, 20.0 ], [ 0.0, 50.0 ] ) + filter.setNumberElements( [ 5, 5 ], [ 5, 5 ], [ 10 ] ) # 10 cells along X, Y, Z axis + filter.setGenerateCellsGlobalIds( True ) + filter.setGeneratePointsGlobalIds( True ) + + cells_dim1 = FieldInfo( "cell1", 1, "CELLS" ) + cells_dim3 = FieldInfo( "cell3", 3, "CELLS" ) + points_dim1 = FieldInfo( "point1", 1, "POINTS" ) + points_dim3 = FieldInfo( "point3", 3, "POINTS" ) + filter.setFields( [ cells_dim1, cells_dim3, points_dim1, points_dim3 ] ) + + return filter + + +def test_generate_rectilinear_grid( generate_rectilinear_grid_filter: GenerateRectilinearGrid ) -> None: + generate_rectilinear_grid_filter.Update() + mesh = generate_rectilinear_grid_filter.GetOutputDataObject( 0 ) + + assert isinstance( mesh, vtkUnstructuredGrid ) + assert mesh.GetNumberOfCells() == 1000 + assert mesh.GetNumberOfPoints() == 1331 + assert mesh.GetBounds() == ( 0.0, 10.0, 0.0, 20.0, 0.0, 50.0 ) + + pointData: vtkPointData = mesh.GetPointData() + ptArray1: vtkDataArray = pointData.GetArray( "point1" ) + ptArray3: vtkDataArray = pointData.GetArray( "point3" ) + assert ptArray1.GetNumberOfComponents() == 1 + assert ptArray3.GetNumberOfComponents() == 3 + + cellData: vtkCellData = mesh.GetCellData() + cellArray1: vtkDataArray = cellData.GetArray( "cell1" ) + cellArray3: vtkDataArray = cellData.GetArray( "cell3" ) + assert cellArray1.GetNumberOfComponents() == 1 + assert cellArray3.GetNumberOfComponents() == 3 def test_generate_cube(): diff --git a/geos-mesh/tests/test_generate_fractures.py b/geos-mesh/tests/test_generate_fractures.py index 66c9496f..fb6fd978 100644 --- a/geos-mesh/tests/test_generate_fractures.py +++ b/geos-mesh/tests/test_generate_fractures.py @@ -6,8 +6,9 @@ from vtkmodules.util.numpy_support import numpy_to_vtk, vtk_to_numpy from geos.mesh.doctor.actions.check_fractures import format_collocated_nodes from geos.mesh.doctor.actions.generate_cube import build_rectilinear_blocks_mesh, XYZ -from geos.mesh.doctor.actions.generate_fractures import ( __split_mesh_on_fractures, Options, FracturePolicy, +from geos.mesh.doctor.actions.generate_fractures import ( split_mesh_on_fractures, Options, FracturePolicy, Coordinates3D, IDMapping ) +from geos.mesh.doctor.filters.GenerateFractures import GenerateFractures from geos.mesh.utils.genericHelpers import to_vtk_id_list FaceNodesCoords = tuple[ tuple[ float ] ] @@ -202,7 +203,7 @@ def __generate_test_data() -> Iterator[ TestCase ]: @pytest.mark.parametrize( "test_case", __generate_test_data() ) def test_generate_fracture( test_case: TestCase ): - main_mesh, fracture_meshes = __split_mesh_on_fractures( test_case.input_mesh, test_case.options ) + main_mesh, fracture_meshes = split_mesh_on_fractures( test_case.input_mesh, test_case.options ) fracture_mesh: vtkUnstructuredGrid = fracture_meshes[ 0 ] assert main_mesh.GetNumberOfPoints() == test_case.result.main_mesh_num_points assert main_mesh.GetNumberOfCells() == test_case.result.main_mesh_num_cells @@ -214,6 +215,32 @@ def test_generate_fracture( test_case: TestCase ): assert len( res ) == test_case.result.fracture_mesh_num_points +@pytest.mark.parametrize( "test_case_filter", __generate_test_data() ) +def test_GenerateFracture( test_case_filter: TestCase ): + genFracFilter = GenerateFractures() + genFracFilter.SetInputDataObject( 0, test_case_filter.input_mesh ) + genFracFilter.setFieldName( test_case_filter.options.field ) + field_values: str = ','.join( map( str, test_case_filter.options.field_values_combined ) ) + genFracFilter.setFieldValues( field_values ) + genFracFilter.setFracturesOutputDirectory( "." ) + if test_case_filter.options.policy == FracturePolicy.FIELD: + genFracFilter.setPolicy( 0 ) + else: + genFracFilter.setPolicy( 1 ) + genFracFilter.Update() + + main_mesh, fracture_meshes = genFracFilter.getAllGrids() + fracture_mesh: vtkUnstructuredGrid = fracture_meshes[ 0 ] + assert main_mesh.GetNumberOfPoints() == test_case_filter.result.main_mesh_num_points + assert main_mesh.GetNumberOfCells() == test_case_filter.result.main_mesh_num_cells + assert fracture_mesh.GetNumberOfPoints() == test_case_filter.result.fracture_mesh_num_points + assert fracture_mesh.GetNumberOfCells() == test_case_filter.result.fracture_mesh_num_cells + + res = format_collocated_nodes( fracture_mesh ) + assert res == test_case_filter.collocated_nodes + assert len( res ) == test_case_filter.result.fracture_mesh_num_points + + def add_simplified_field_for_cells( mesh: vtkUnstructuredGrid, field_name: str, field_dimension: int ): """Reduce functionality obtained from src.geos.mesh.doctor.actions.generate_fracture.__add_fields where the goal is to add a cell data array with incrementing values. @@ -299,7 +326,7 @@ def add_quad( mesh: vtkUnstructuredGrid, face: FaceNodesCoords ): @pytest.mark.skip( "Test to be fixed" ) def test_copy_fields_when_splitting_mesh(): """This test is designed to check the __copy_fields method from generate_fractures, - that will be called when using __split_mesh_on_fractures method from generate_fractures. + that will be called when using split_mesh_on_fractures method from generate_fractures. """ # Generating the rectilinear grid and its quads on all borders x: numpy.array = numpy.array( [ 0, 1, 2 ] ) @@ -330,7 +357,7 @@ def test_copy_fields_when_splitting_mesh(): field_values_per_fracture=[ frozenset( map( int, [ "9" ] ) ) ], mesh_VtkOutput=None, all_fractures_VtkOutput=None ) - main_mesh, fracture_meshes = __split_mesh_on_fractures( mesh, options ) + main_mesh, fracture_meshes = split_mesh_on_fractures( mesh, options ) fracture_mesh: vtkUnstructuredGrid = fracture_meshes[ 0 ] assert main_mesh.GetCellData().GetNumberOfArrays() == 1 assert fracture_mesh.GetCellData().GetNumberOfArrays() == 1 @@ -344,7 +371,7 @@ def test_copy_fields_when_splitting_mesh(): # Test for invalid point field name add_simplified_field_for_cells( mesh, "GLOBAL_IDS_POINTS", 1 ) with pytest.raises( ValueError ) as pytest_wrapped_e: - main_mesh, fracture_meshes = __split_mesh_on_fractures( mesh, options ) + main_mesh, fracture_meshes = split_mesh_on_fractures( mesh, options ) assert pytest_wrapped_e.type == ValueError # Test for invalid cell field name mesh: vtkUnstructuredGrid = build_rectilinear_blocks_mesh( [ xyzs ] ) @@ -356,5 +383,5 @@ def test_copy_fields_when_splitting_mesh(): add_simplified_field_for_cells( mesh, "GLOBAL_IDS_CELLS", 1 ) assert mesh.GetCellData().GetNumberOfArrays() == 2 with pytest.raises( ValueError ) as pytest_wrapped_e: - main_mesh, fracture_meshes = __split_mesh_on_fractures( mesh, options ) + main_mesh, fracture_meshes = split_mesh_on_fractures( mesh, options ) assert pytest_wrapped_e.type == ValueError diff --git a/geos-mesh/tests/test_generate_global_ids.py b/geos-mesh/tests/test_generate_global_ids.py index 614f771c..8f1ad25a 100644 --- a/geos-mesh/tests/test_generate_global_ids.py +++ b/geos-mesh/tests/test_generate_global_ids.py @@ -1,6 +1,6 @@ from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkUnstructuredGrid, vtkVertex, VTK_VERTEX -from geos.mesh.doctor.actions.generate_global_ids import __build_global_ids +from geos.mesh.doctor.actions.generate_global_ids import build_global_ids def test_generate_global_ids(): @@ -17,7 +17,7 @@ def test_generate_global_ids(): mesh.SetPoints( points ) mesh.SetCells( [ VTK_VERTEX ], vertices ) - __build_global_ids( mesh, True, True ) + build_global_ids( mesh, True, True ) global_cell_ids = mesh.GetCellData().GetGlobalIds() global_point_ids = mesh.GetPointData().GetGlobalIds() diff --git a/geos-mesh/tests/test_non_conformal.py b/geos-mesh/tests/test_non_conformal.py index 9f6da41a..e95c4697 100644 --- a/geos-mesh/tests/test_non_conformal.py +++ b/geos-mesh/tests/test_non_conformal.py @@ -1,6 +1,6 @@ import numpy from geos.mesh.doctor.actions.generate_cube import build_rectilinear_blocks_mesh, XYZ -from geos.mesh.doctor.actions.non_conformal import Options, __action +from geos.mesh.doctor.actions.non_conformal import Options, mesh_action def test_two_close_hexs(): @@ -12,13 +12,13 @@ def test_two_close_hexs(): # Close enough, but points tolerance is too strict to consider the faces matching. options = Options( angle_tolerance=1., point_tolerance=delta / 2, face_tolerance=delta * 2 ) - results = __action( mesh, options ) + results = mesh_action( mesh, options ) assert len( results.non_conformal_cells ) == 1 assert set( results.non_conformal_cells[ 0 ] ) == { 0, 1 } # Close enough, and points tolerance is loose enough to consider the faces matching. options = Options( angle_tolerance=1., point_tolerance=delta * 2, face_tolerance=delta * 2 ) - results = __action( mesh, options ) + results = mesh_action( mesh, options ) assert len( results.non_conformal_cells ) == 0 @@ -31,7 +31,7 @@ def test_two_distant_hexs(): options = Options( angle_tolerance=1., point_tolerance=delta / 2., face_tolerance=delta / 2. ) - results = __action( mesh, options ) + results = mesh_action( mesh, options ) assert len( results.non_conformal_cells ) == 0 @@ -44,7 +44,7 @@ def test_two_close_shifted_hexs(): options = Options( angle_tolerance=1., point_tolerance=delta_x * 2, face_tolerance=delta_x * 2 ) - results = __action( mesh, options ) + results = mesh_action( mesh, options ) assert len( results.non_conformal_cells ) == 1 assert set( results.non_conformal_cells[ 0 ] ) == { 0, 1 } @@ -58,6 +58,6 @@ def test_big_elem_next_to_small_elem(): options = Options( angle_tolerance=1., point_tolerance=delta * 2, face_tolerance=delta * 2 ) - results = __action( mesh, options ) + results = mesh_action( mesh, options ) assert len( results.non_conformal_cells ) == 1 assert set( results.non_conformal_cells[ 0 ] ) == { 0, 1 } diff --git a/geos-mesh/tests/test_self_intersecting_elements.py b/geos-mesh/tests/test_self_intersecting_elements.py index 45216f01..99299222 100644 --- a/geos-mesh/tests/test_self_intersecting_elements.py +++ b/geos-mesh/tests/test_self_intersecting_elements.py @@ -1,6 +1,8 @@ from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkHexahedron, vtkUnstructuredGrid, VTK_HEXAHEDRON -from geos.mesh.doctor.actions.self_intersecting_elements import Options, __action +from geos.mesh.doctor.actions.self_intersecting_elements import Options, mesh_action +from geos.mesh.doctor.filters.SelfIntersectingElements import SelfIntersectingElements +import pytest def test_jumbled_hex(): @@ -35,7 +37,144 @@ def test_jumbled_hex(): mesh.SetPoints( points ) mesh.SetCells( cell_types, cells ) - result = __action( mesh, Options( min_distance=0. ) ) + result = mesh_action( mesh, Options( min_distance=0. ) ) assert len( result.intersecting_faces_elements ) == 1 assert result.intersecting_faces_elements[ 0 ] == 0 + + +@pytest.fixture +def jumbled_hex_mesh(): + """Create a hexahedron with intentionally swapped nodes to create self-intersecting faces.""" + points = vtkPoints() + points.SetNumberOfPoints( 8 ) + points.SetPoint( 0, ( 0, 0, 0 ) ) + points.SetPoint( 1, ( 1, 0, 0 ) ) + points.SetPoint( 2, ( 1, 1, 0 ) ) + points.SetPoint( 3, ( 0, 1, 0 ) ) + points.SetPoint( 4, ( 0, 0, 1 ) ) + points.SetPoint( 5, ( 1, 0, 1 ) ) + points.SetPoint( 6, ( 1, 1, 1 ) ) + points.SetPoint( 7, ( 0, 1, 1 ) ) + + cell_types = [ VTK_HEXAHEDRON ] + cells = vtkCellArray() + cells.AllocateExact( 1, 8 ) + + hex = vtkHexahedron() + hex.GetPointIds().SetId( 0, 0 ) + hex.GetPointIds().SetId( 1, 1 ) + hex.GetPointIds().SetId( 2, 3 ) # Intentionally wrong + hex.GetPointIds().SetId( 3, 2 ) # Intentionally wrong + hex.GetPointIds().SetId( 4, 4 ) + hex.GetPointIds().SetId( 5, 5 ) + hex.GetPointIds().SetId( 6, 6 ) + hex.GetPointIds().SetId( 7, 7 ) + cells.InsertNextCell( hex ) + + mesh = vtkUnstructuredGrid() + mesh.SetPoints( points ) + mesh.SetCells( cell_types, cells ) + return mesh + + +@pytest.fixture +def valid_hex_mesh(): + """Create a properly ordered hexahedron with no self-intersecting faces.""" + points = vtkPoints() + points.SetNumberOfPoints( 8 ) + points.SetPoint( 0, ( 0, 0, 0 ) ) + points.SetPoint( 1, ( 1, 0, 0 ) ) + points.SetPoint( 2, ( 1, 1, 0 ) ) + points.SetPoint( 3, ( 0, 1, 0 ) ) + points.SetPoint( 4, ( 0, 0, 1 ) ) + points.SetPoint( 5, ( 1, 0, 1 ) ) + points.SetPoint( 6, ( 1, 1, 1 ) ) + points.SetPoint( 7, ( 0, 1, 1 ) ) + + cell_types = [ VTK_HEXAHEDRON ] + cells = vtkCellArray() + cells.AllocateExact( 1, 8 ) + + hex = vtkHexahedron() + for i in range( 8 ): + hex.GetPointIds().SetId( i, i ) + cells.InsertNextCell( hex ) + + mesh = vtkUnstructuredGrid() + mesh.SetPoints( points ) + mesh.SetCells( cell_types, cells ) + return mesh + + +def test_self_intersecting_elements_filter_detects_invalid_elements( jumbled_hex_mesh ): + """Test that the SelfIntersectingElements filter correctly detects invalid elements.""" + filter = SelfIntersectingElements() + filter.setMinDistance( 0.0 ) + filter.SetInputDataObject( 0, jumbled_hex_mesh ) + filter.Update() + + output = filter.getGrid() + # Check that the filter detected the invalid element + intersecting_faces = filter.getIntersectingFacesElements() + assert len( intersecting_faces ) == 1 + assert intersecting_faces[ 0 ] == 0 + + # Check that output mesh has same structure + assert output.GetNumberOfCells() == 1 + assert output.GetNumberOfPoints() == 8 + + +def test_self_intersecting_elements_filter_valid_mesh( valid_hex_mesh ): + """Test that the SelfIntersectingElements filter finds no issues in a valid mesh.""" + filter = SelfIntersectingElements() + filter.setMinDistance( 1e-12 ) # Use small tolerance instead of 0.0 + filter.SetInputDataObject( 0, valid_hex_mesh ) + filter.Update() + + output = filter.getGrid() + # Check that no invalid elements were found + assert len( filter.getIntersectingFacesElements() ) == 0 + assert len( filter.getWrongNumberOfPointsElements() ) == 0 + assert len( filter.getIntersectingEdgesElements() ) == 0 + assert len( filter.getNonContiguousEdgesElements() ) == 0 + assert len( filter.getNonConvexElements() ) == 0 + assert len( filter.getFacesOrientedIncorrectlyElements() ) == 0 + + # Check that output mesh has same structure + assert output.GetNumberOfCells() == 1 + assert output.GetNumberOfPoints() == 8 + + +def test_self_intersecting_elements_filter_paint_invalid_elements( jumbled_hex_mesh ): + """Test that the SelfIntersectingElements filter can paint invalid elements.""" + filter = SelfIntersectingElements() + filter.setMinDistance( 0.0 ) + filter.setPaintInvalidElements( 1 ) # Enable painting + filter.SetInputDataObject( 0, jumbled_hex_mesh ) + filter.Update() + + output = filter.getGrid() + # Check that painting arrays were added to the output + cell_data = output.GetCellData() + + # Should have arrays marking the invalid elements + # The exact array names depend on the implementation + assert cell_data.GetNumberOfArrays() > 0 + + # Check that invalid elements were detected + intersecting_faces = filter.getIntersectingFacesElements() + assert len( intersecting_faces ) == 1 + + +def test_self_intersecting_elements_filter_getters_setters(): + """Test getter and setter methods of the SelfIntersectingElements filter.""" + filter = SelfIntersectingElements() + + # Test min distance getter/setter + filter.setMinDistance( 0.5 ) + assert filter.getMinDistance() == 0.5 + + # Test paint invalid elements setter (no getter available) + filter.setPaintInvalidElements( 1 ) + # No exception should be raised diff --git a/geos-mesh/tests/test_supported_elements.py b/geos-mesh/tests/test_supported_elements.py index 07321abc..297c3899 100644 --- a/geos-mesh/tests/test_supported_elements.py +++ b/geos-mesh/tests/test_supported_elements.py @@ -1,10 +1,13 @@ # import os import pytest +import numpy as np from typing import Tuple from vtkmodules.vtkCommonCore import vtkIdList, vtkPoints -from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid, VTK_POLYHEDRON +from vtkmodules.vtkCommonDataModel import ( vtkUnstructuredGrid, vtkQuad, vtkTetra, vtkHexahedron, vtkPolyhedron, + vtkCellArray, VTK_POLYHEDRON, VTK_QUAD, VTK_TETRA, VTK_HEXAHEDRON ) # from geos.mesh.doctor.actions.supported_elements import Options, action, __action from geos.mesh.doctor.actions.vtk_polyhedron import parse_face_stream, FaceStream +# from geos.mesh.doctor.filters.SupportedElements import SupportedElements from geos.mesh.utils.genericHelpers import to_vtk_id_list @@ -116,3 +119,277 @@ def test_parse_face_stream() -> None: face_stream = FaceStream.build_from_vtk_id_list( faces ) assert face_stream.num_faces == 12 assert face_stream.num_support_points == 20 + + +def create_simple_tetra_grid(): + """Create a simple tetrahedral grid for testing""" + # Create an unstructured grid + points_tetras: vtkPoints = vtkPoints() + points_tetras_coords: list[ tuple[ float ] ] = [ + ( 1.0, 0.5, 0.0 ), # point0 + ( 1.0, 0.0, 1.0 ), + ( 1.0, 1.0, 1.0 ), + ( 0.0, 0.5, 0.5 ), + ( 2.0, 0.5, 0.5 ), + ( 1.0, 0.5, 2.0 ), # point5 + ( 0.0, 0.5, 1.5 ), + ( 2.0, 0.5, 1.5 ) + ] + for point_tetra in points_tetras_coords: + points_tetras.InsertNextPoint( point_tetra ) + + tetra1: vtkTetra = vtkTetra() + tetra1.GetPointIds().SetId( 0, 0 ) + tetra1.GetPointIds().SetId( 1, 1 ) + tetra1.GetPointIds().SetId( 2, 2 ) + tetra1.GetPointIds().SetId( 3, 3 ) + + tetra2: vtkTetra = vtkTetra() + tetra2.GetPointIds().SetId( 0, 0 ) + tetra2.GetPointIds().SetId( 1, 2 ) + tetra2.GetPointIds().SetId( 2, 1 ) + tetra2.GetPointIds().SetId( 3, 4 ) + + tetra3: vtkTetra = vtkTetra() + tetra3.GetPointIds().SetId( 0, 1 ) + tetra3.GetPointIds().SetId( 1, 5 ) + tetra3.GetPointIds().SetId( 2, 2 ) + tetra3.GetPointIds().SetId( 3, 6 ) + + tetra4: vtkTetra = vtkTetra() + tetra4.GetPointIds().SetId( 0, 1 ) + tetra4.GetPointIds().SetId( 1, 2 ) + tetra4.GetPointIds().SetId( 2, 5 ) + tetra4.GetPointIds().SetId( 3, 7 ) + + 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 ) + return tetras_grid + + +def create_mixed_grid(): + """Create a grid with supported and unsupported cell types, 4 Hexahedrons with 2 quad fracs vertical""" + # Create an unstructured grid + four_hexs_points: vtkPoints = vtkPoints() + four_hexs_points_coords: list[ tuple[ float ] ] = [ + ( 0.0, 0.0, 0.0 ), # point0 + ( 1.0, 0.0, 0.0 ), # point1 + ( 2.0, 0.0, 0.0 ), # point2 + ( 0.0, 1.0, 0.0 ), # point3 + ( 1.0, 1.0, 0.0 ), # point4 + ( 2.0, 1.0, 0.0 ), # point5 + ( 0.0, 0.0, 1.0 ), # point6 + ( 1.0, 0.0, 1.0 ), # point7 + ( 2.0, 0.0, 1.0 ), # point8 + ( 0.0, 1.0, 1.0 ), # point9 + ( 1.0, 1.0, 1.0 ), # point10 + ( 2.0, 1.0, 1.0 ), # point11 + ( 0.0, 0.0, 2.0 ), # point12 + ( 1.0, 0.0, 2.0 ), # point13 + ( 2.0, 0.0, 2.0 ), # point14 + ( 0.0, 1.0, 2.0 ), # point15 + ( 1.0, 1.0, 2.0 ), # point16 + ( 2.0, 1.0, 2.0 ) + ] + for four_hexs_point in four_hexs_points_coords: + four_hexs_points.InsertNextPoint( four_hexs_point ) + + # hex1 + four_hex1: vtkHexahedron = vtkHexahedron() + four_hex1.GetPointIds().SetId( 0, 0 ) + four_hex1.GetPointIds().SetId( 1, 1 ) + four_hex1.GetPointIds().SetId( 2, 4 ) + four_hex1.GetPointIds().SetId( 3, 3 ) + four_hex1.GetPointIds().SetId( 4, 6 ) + four_hex1.GetPointIds().SetId( 5, 7 ) + four_hex1.GetPointIds().SetId( 6, 10 ) + four_hex1.GetPointIds().SetId( 7, 9 ) + + # hex2 + four_hex2: vtkHexahedron = vtkHexahedron() + four_hex2.GetPointIds().SetId( 0, 0 + 1 ) + four_hex2.GetPointIds().SetId( 1, 1 + 1 ) + four_hex2.GetPointIds().SetId( 2, 4 + 1 ) + four_hex2.GetPointIds().SetId( 3, 3 + 1 ) + four_hex2.GetPointIds().SetId( 4, 6 + 1 ) + four_hex2.GetPointIds().SetId( 5, 7 + 1 ) + four_hex2.GetPointIds().SetId( 6, 10 + 1 ) + four_hex2.GetPointIds().SetId( 7, 9 + 1 ) + + # hex3 + four_hex3: vtkHexahedron = vtkHexahedron() + four_hex3.GetPointIds().SetId( 0, 0 + 6 ) + four_hex3.GetPointIds().SetId( 1, 1 + 6 ) + four_hex3.GetPointIds().SetId( 2, 4 + 6 ) + four_hex3.GetPointIds().SetId( 3, 3 + 6 ) + four_hex3.GetPointIds().SetId( 4, 6 + 6 ) + four_hex3.GetPointIds().SetId( 5, 7 + 6 ) + four_hex3.GetPointIds().SetId( 6, 10 + 6 ) + four_hex3.GetPointIds().SetId( 7, 9 + 6 ) + + # hex4 + four_hex4: vtkHexahedron = vtkHexahedron() + four_hex4.GetPointIds().SetId( 0, 0 + 7 ) + four_hex4.GetPointIds().SetId( 1, 1 + 7 ) + four_hex4.GetPointIds().SetId( 2, 4 + 7 ) + four_hex4.GetPointIds().SetId( 3, 3 + 7 ) + four_hex4.GetPointIds().SetId( 4, 6 + 7 ) + four_hex4.GetPointIds().SetId( 5, 7 + 7 ) + four_hex4.GetPointIds().SetId( 6, 10 + 7 ) + four_hex4.GetPointIds().SetId( 7, 9 + 7 ) + + # quad1 + four_hex_quad1: vtkQuad = vtkQuad() + four_hex_quad1.GetPointIds().SetId( 0, 1 ) + four_hex_quad1.GetPointIds().SetId( 1, 4 ) + four_hex_quad1.GetPointIds().SetId( 2, 10 ) + four_hex_quad1.GetPointIds().SetId( 3, 7 ) + + # quad2 + four_hex_quad2: vtkQuad = vtkQuad() + four_hex_quad2.GetPointIds().SetId( 0, 1 + 6 ) + four_hex_quad2.GetPointIds().SetId( 1, 4 + 6 ) + four_hex_quad2.GetPointIds().SetId( 2, 10 + 6 ) + four_hex_quad2.GetPointIds().SetId( 3, 7 + 6 ) + + four_hex_grid_2_quads = vtkUnstructuredGrid() + four_hex_grid_2_quads.SetPoints( four_hexs_points ) + all_cell_types_four_hex_grid_2_quads = [ VTK_HEXAHEDRON ] * 4 + [ VTK_QUAD ] * 2 + all_cells_four_hex_grid_2_quads = [ four_hex1, four_hex2, four_hex3, four_hex4, four_hex_quad1, four_hex_quad2 ] + for cell_type, cell in zip( all_cell_types_four_hex_grid_2_quads, all_cells_four_hex_grid_2_quads ): + four_hex_grid_2_quads.InsertNextCell( cell_type, cell.GetPointIds() ) + return four_hex_grid_2_quads + + +def create_unsupported_polyhedron_grid(): + """Create a grid with an unsupported polyhedron (non-convex)""" + grid = vtkUnstructuredGrid() + # Create points for the grid + points = vtkPoints() # Need to import vtkPoints + # Create points for a non-convex polyhedron + point_coords = np.array( [ + [ 0.0, 0.0, 0.0 ], # 0 + [ 1.0, 0.0, 0.0 ], # 1 + [ 1.0, 1.0, 0.0 ], # 2 + [ 0.0, 1.0, 0.0 ], # 3 + [ 0.0, 0.0, 1.0 ], # 4 + [ 1.0, 0.0, 1.0 ], # 5 + [ 1.0, 1.0, 1.0 ], # 6 + [ 0.0, 1.0, 1.0 ], # 7 + [ 0.5, 0.5, -0.5 ] # 8 (point makes it non-convex) + ] ) + # Add points to the points array + for point in point_coords: + points.InsertNextPoint( point ) + # Set the points in the grid + grid.SetPoints( points ) + # Create a polyhedron + polyhedron = vtkPolyhedron() + # For simplicity, we'll create a polyhedron that would be recognized as unsupported + # This is a simplified example - you may need to adjust based on your actual implementation + polyhedron.GetPointIds().SetNumberOfIds( 9 ) + for i in range( 9 ): + polyhedron.GetPointIds().SetId( i, i ) + # Add the polyhedron to the grid + grid.InsertNextCell( polyhedron.GetCellType(), polyhedron.GetPointIds() ) + return grid + + +# TODO reimplement once SupportedElements can handle multiprocessing +# class TestSupportedElements: + +# def test_only_supported_elements( self ): +# """Test a grid with only supported element types""" +# # Create grid with only supported elements (tetra) +# grid = create_simple_tetra_grid() +# # Apply the filter +# filter = SupportedElements() +# filter.SetInputDataObject( grid ) +# filter.Update() +# result = filter.getGrid() +# assert result is not None +# # Verify no arrays were added (since all elements are supported) +# assert result.GetCellData().GetArray( "HasUnsupportedType" ) is None +# assert result.GetCellData().GetArray( "IsUnsupportedPolyhedron" ) is None + +# def test_unsupported_element_types( self ): +# """Test a grid with unsupported element types""" +# # Create grid with unsupported elements +# grid = create_mixed_grid() +# # Apply the filter with painting enabled +# filter = SupportedElements() +# filter.m_logger.critical( "test_unsupported_element_types" ) +# filter.SetInputDataObject( grid ) +# filter.setPaintUnsupportedElementTypes( 1 ) +# filter.Update() +# result = filter.getGrid() +# assert result is not None +# # Verify the array was added +# unsupported_array = result.GetCellData().GetArray( "HasUnsupportedType" ) +# assert unsupported_array is not None +# for i in range( 0, 4 ): +# assert unsupported_array.GetValue( i ) == 0 # Hexahedron should be supported +# for j in range( 4, 6 ): +# assert unsupported_array.GetValue( j ) == 1 # Quad should not be supported + +# TODO Needs parallelism to work +# def test_unsupported_polyhedron( self ): +# """Test a grid with unsupported polyhedron""" +# # Create grid with unsupported polyhedron +# grid = create_unsupported_polyhedron_grid() +# # Apply the filter with painting enabled +# filter = SupportedElements() +# filter.m_logger.critical( "test_unsupported_polyhedron" ) +# filter.SetInputDataObject( grid ) +# filter.setPaintUnsupportedPolyhedrons( 1 ) +# filter.Update() +# result = filter.getGrid() +# assert result is not None +# # Verify the array was added +# polyhedron_array = result.GetCellData().GetArray( "IsUnsupportedPolyhedron" ) +# assert polyhedron_array is None +# # Since we created an unsupported polyhedron, it should be marked +# assert polyhedron_array.GetValue( 0 ) == 1 + +# def test_paint_flags( self ): +# """Test setting invalid paint flags""" +# filter = SupportedElements() +# # Should log an error but not raise an exception +# filter.setPaintUnsupportedElementTypes( 2 ) # Invalid value +# filter.setPaintUnsupportedPolyhedrons( 2 ) # Invalid value +# # Values should remain unchanged +# assert filter.m_paintUnsupportedElementTypes == 0 +# assert filter.m_paintUnsupportedPolyhedrons == 0 + +# def test_set_chunk_size( self ): +# """Test that setChunkSize properly updates the chunk size""" +# # Create filter instance +# filter = SupportedElements() +# # Note the initial value +# initial_chunk_size = filter.m_chunk_size +# # Set a new chunk size +# new_chunk_size = 100 +# filter.setChunkSize( new_chunk_size ) +# # Verify the chunk size was updated +# assert filter.m_chunk_size == new_chunk_size +# assert filter.m_chunk_size != initial_chunk_size + +# def test_set_num_proc( self ): +# """Test that setNumProc properly updates the number of processors""" +# # Create filter instance +# filter = SupportedElements() +# # Note the initial value +# initial_num_proc = filter.m_num_proc +# # Set a new number of processors +# new_num_proc = 4 +# filter.setNumProc( new_num_proc ) +# # Verify the number of processors was updated +# assert filter.m_num_proc == new_num_proc +# assert filter.m_num_proc != initial_num_proc