From 186f0d03cff97bd7eb0176f8a491d8b1ed3c4094 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Tue, 3 Feb 2026 16:07:49 -0600 Subject: [PATCH 1/3] Draft --- .../src/geos/mesh_doctor/actions/orphan2d.py | 252 ++++++++++++++++++ .../src/geos/mesh_doctor/parsing/__init__.py | 1 + .../mesh_doctor/parsing/orphan2dParsing.py | 99 +++++++ mesh-doctor/src/geos/mesh_doctor/register.py | 2 +- 4 files changed, 353 insertions(+), 1 deletion(-) create mode 100644 mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py create mode 100644 mesh-doctor/src/geos/mesh_doctor/parsing/orphan2dParsing.py diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py b/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py new file mode 100644 index 000000000..0a21a2bf7 --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py @@ -0,0 +1,252 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from dataclasses import dataclass +from typing import Optional +import vtk +from tqdm import tqdm +from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid +from geos.mesh_doctor.parsing.cliParsing import setupLogger +from geos.mesh.io.vtkIO import readUnstructuredGrid, writeMesh, VtkOutput + + +@dataclass( frozen=True ) +class Options: + orphanVtkOutput: Optional[ VtkOutput ] + cleanVtkOutput: Optional[ VtkOutput ] + + +@dataclass( frozen=True ) +class Result: + total2dCells: int + total3dCells: int + matched2dCells: int + orphaned2dCells: int + orphaned2dIndices: list[ int ] + orphanMesh: Optional[ vtkUnstructuredGrid ] + cleanMesh: Optional[ vtkUnstructuredGrid ] + + +def getCellFacePoints( cell: vtk.vtkCell ) -> list[ tuple[ int, ...] ]: + """Extract all face point sets from a 3D cell. + + Args: + cell: A 3D VTK cell. + + Returns: + List of tuples containing sorted global point IDs for each face. + """ + faces = [] + cellType = cell.GetCellType() + + if cellType == vtk.VTK_TETRA: + # 4 triangular faces + faces = [ [ 0, 1, 2 ], [ 0, 1, 3 ], [ 0, 2, 3 ], [ 1, 2, 3 ] ] + elif cellType == vtk.VTK_HEXAHEDRON: + # 6 quadrilateral faces + faces = [ [ 0, 1, 2, 3 ], [ 4, 5, 6, 7 ], [ 0, 1, 5, 4 ], [ 1, 2, 6, 5 ], [ 2, 3, 7, 6 ], [ 3, 0, 4, 7 ] ] + elif cellType == vtk.VTK_WEDGE: + # 2 triangular + 3 quadrilateral faces + faces = [ [ 0, 1, 2 ], [ 3, 4, 5 ], [ 0, 1, 4, 3 ], [ 1, 2, 5, 4 ], [ 2, 0, 3, 5 ] ] + elif cellType == vtk.VTK_PYRAMID: + # 1 quadrilateral + 4 triangular faces + faces = [ [ 0, 1, 2, 3 ], [ 0, 1, 4 ], [ 1, 2, 4 ], [ 2, 3, 4 ], [ 3, 0, 4 ] ] + + # Convert local indices to global point IDs + pointIds = cell.GetPointIds() + globalFaces = [] + for face in faces: + globalFace = tuple( sorted( [ pointIds.GetId( i ) for i in face ] ) ) + globalFaces.append( globalFace ) + + return globalFaces + + +def buildMeshSubset( mesh: vtkUnstructuredGrid, cellIndices: list[ int ], description: str ) -> vtkUnstructuredGrid: + """Build a new mesh containing only the specified cells. + + Args: + mesh: Input unstructured grid. + cellIndices: List of cell indices to include. + description: Description for progress messages. + + Returns: + A new vtkUnstructuredGrid containing only the specified cells. + """ + if not cellIndices: + setupLogger.info( f"No {description} to create." ) + return None + + setupLogger.info( f"Creating mesh with {len(cellIndices)} {description}..." ) + + # Create new unstructured grid + newMesh = vtkUnstructuredGrid() + newPoints = vtk.vtkPoints() + + # Preserve point data type (precision) from input mesh + oldPoints = mesh.GetPoints() + newPoints.SetDataType( oldPoints.GetDataType() ) + + # Map old point IDs to new point IDs + pointMap = {} + newPointId = 0 + + # Add cells + for oldCellId in tqdm( cellIndices, desc=f"Building {description} mesh" ): + cell = mesh.GetCell( oldCellId ) + pointIds = cell.GetPointIds() + + # Create new point IDs if needed + newCellPointIds = vtk.vtkIdList() + for i in range( pointIds.GetNumberOfIds() ): + oldId = pointIds.GetId( i ) + if oldId not in pointMap: + pointMap[ oldId ] = newPointId + newPoints.InsertNextPoint( mesh.GetPoint( oldId ) ) + newPointId += 1 + newCellPointIds.InsertNextId( pointMap[ oldId ] ) + + newMesh.InsertNextCell( cell.GetCellType(), newCellPointIds ) + + newMesh.SetPoints( newPoints ) + + # Copy point data + setupLogger.debug( "Copying point data..." ) + for i in range( mesh.GetPointData().GetNumberOfArrays() ): + oldArray = mesh.GetPointData().GetArray( i ) + + # Create new array with same type and number of components + newArray = oldArray.NewInstance() + newArray.SetName( oldArray.GetName() ) + newArray.SetNumberOfComponents( oldArray.GetNumberOfComponents() ) + + for oldId in sorted( pointMap.keys() ): + newArray.InsertNextTuple( oldArray.GetTuple( oldId ) ) + + newMesh.GetPointData().AddArray( newArray ) + + # Copy cell data + setupLogger.debug( "Copying cell data..." ) + for i in range( mesh.GetCellData().GetNumberOfArrays() ): + oldArray = mesh.GetCellData().GetArray( i ) + + # Create new array with same type and number of components + newArray = oldArray.NewInstance() + newArray.SetName( oldArray.GetName() ) + newArray.SetNumberOfComponents( oldArray.GetNumberOfComponents() ) + + for oldCellId in cellIndices: + newArray.InsertNextTuple( oldArray.GetTuple( oldCellId ) ) + + newMesh.GetCellData().AddArray( newArray ) + + return newMesh + + +def meshAction( mesh: vtkUnstructuredGrid, options: Options ) -> Result: + """Check if 2D cells are faces of 3D cells. + + Args: + mesh: The input mesh to analyze. + options: The options for processing. + + Returns: + Result: The result of the orphan check. + """ + setupLogger.info( "Starting orphan 2D cell check..." ) + + # Separate 2D and 3D cells + cell2dIndices = [] + cell3dIndices = [] + + setupLogger.info( "Classifying cells by dimension..." ) + for i in tqdm( range( mesh.GetNumberOfCells() ), desc="Scanning cells" ): + cell = mesh.GetCell( i ) + dim = cell.GetCellDimension() + if dim == 2: + cell2dIndices.append( i ) + elif dim == 3: + cell3dIndices.append( i ) + + setupLogger.info( f"Found {len(cell2dIndices)} 2D cells" ) + setupLogger.info( f"Found {len(cell3dIndices)} 3D cells" ) + + # Build set of all faces from 3D cells + setupLogger.info( "Extracting faces from 3D cells..." ) + all3dFaces = set() + for idx in tqdm( cell3dIndices, desc="Processing 3D cells" ): + cell = mesh.GetCell( idx ) + faces = getCellFacePoints( cell ) + all3dFaces.update( faces ) + + setupLogger.info( f"Total unique faces from 3D cells: {len(all3dFaces)}" ) + + # Check each 2D cell + setupLogger.info( "Checking 2D cells against 3D faces..." ) + orphaned2dIndices = [] + matched2dCount = 0 + + for idx in tqdm( cell2dIndices, desc="Checking 2D cells" ): + cell = mesh.GetCell( idx ) + pointIds = cell.GetPointIds() + + # Get sorted point IDs for the 2D cell + cellPoints = tuple( sorted( [ pointIds.GetId( i ) for i in range( pointIds.GetNumberOfIds() ) ] ) ) + + if cellPoints in all3dFaces: + matched2dCount += 1 + else: + orphaned2dIndices.append( idx ) + + setupLogger.info( f"2D cells that ARE faces of 3D cells: {matched2dCount}" ) + setupLogger.info( f"2D cells that are NOT faces of 3D cells: {len(orphaned2dIndices)}" ) + + # Build orphan mesh if requested + orphanMesh = None + if options.orphanVtkOutput and orphaned2dIndices: + orphanMesh = buildMeshSubset( mesh, orphaned2dIndices, "orphaned 2D cells" ) + + # Build clean mesh if requested + cleanMesh = None + if options.cleanVtkOutput: + if orphaned2dIndices: + allCellIndices = list( range( mesh.GetNumberOfCells() ) ) + orphanedSet = set( orphaned2dIndices ) + cleanCellIndices = [ i for i in allCellIndices if i not in orphanedSet ] + cleanMesh = buildMeshSubset( mesh, cleanCellIndices, "clean cells" ) + else: + # No orphans, return the original mesh + cleanMesh = mesh + + return Result( total2dCells=len( cell2dIndices ), + total3dCells=len( cell3dIndices ), + matched2dCells=matched2dCount, + orphaned2dCells=len( orphaned2dIndices ), + orphaned2dIndices=orphaned2dIndices, + orphanMesh=orphanMesh, + cleanMesh=cleanMesh ) + + +def action( vtuInputFile: str, options: Options ) -> Result: + """Read a VTU file and check for orphaned 2D cells. + + Args: + vtuInputFile: The path to the input VTU file. + options: The options for processing. + + Returns: + Result: The result of the orphan check. + """ + mesh = readUnstructuredGrid( vtuInputFile ) + result = meshAction( mesh, options ) + + # Write output files if requested + if options.orphanVtkOutput and result.orphanMesh: + setupLogger.info( f"Writing orphaned cells to {options.orphanVtkOutput.output}" ) + writeMesh( result.orphanMesh, options.orphanVtkOutput ) + + if options.cleanVtkOutput and result.cleanMesh: + setupLogger.info( f"Writing clean mesh to {options.cleanVtkOutput.output}" ) + writeMesh( result.cleanMesh, options.cleanVtkOutput ) + + return result diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py b/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py index 7169fad39..2aca9fcdb 100644 --- a/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py @@ -16,6 +16,7 @@ NON_CONFORMAL = "nonConformal" SELF_INTERSECTING_ELEMENTS = "selfIntersectingElements" SUPPORTED_ELEMENTS = "supportedElements" +ORPHAN_2D = "orphan2d" @dataclass( frozen=True ) diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/orphan2dParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/orphan2dParsing.py new file mode 100644 index 000000000..f8777aecd --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/orphan2dParsing.py @@ -0,0 +1,99 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +from __future__ import annotations +from argparse import _SubParsersAction +from typing import Any, Optional +from geos.mesh_doctor.actions.orphan2d import Options, Result +from geos.mesh_doctor.parsing import ORPHAN_2D +from geos.mesh_doctor.parsing.cliParsing import setupLogger, addVtuInputFileArgument +from geos.mesh.io.vtkIO import VtkOutput + +__ORPHAN_OUTPUT = "orphanOutput" +__CLEAN_OUTPUT = "cleanOutput" +__DATA_MODE = "dataMode" +__DATA_MODE_VALUES = "binary", "ascii" +__DATA_MODE_DEFAULT = __DATA_MODE_VALUES[ 0 ] + + +def convert( parsedOptions: dict[ str, Any ] ) -> Options: + """Convert parsed command-line options to Options object. + + Args: + parsedOptions: Dictionary of parsed command-line options. + + Returns: + Options: Configuration options for orphan check. + """ + orphanOutput: Optional[ str ] = parsedOptions.get( __ORPHAN_OUTPUT ) + cleanOutput: Optional[ str ] = parsedOptions.get( __CLEAN_OUTPUT ) + dataMode: str = parsedOptions.get( __DATA_MODE, __DATA_MODE_DEFAULT ) + isDataModeBinary: bool = dataMode == __DATA_MODE_DEFAULT + + # Create VtkOutput for orphan file if specified + orphanVtkOutput = None + if orphanOutput: + orphanVtkOutput = VtkOutput( output=orphanOutput, isDataModeBinary=isDataModeBinary ) + + # Create VtkOutput for clean file if specified + cleanVtkOutput = None + if cleanOutput: + cleanVtkOutput = VtkOutput( output=cleanOutput, isDataModeBinary=isDataModeBinary ) + + return Options( orphanVtkOutput=orphanVtkOutput, cleanVtkOutput=cleanVtkOutput ) + + +def fillSubparser( subparsers: _SubParsersAction[ Any ] ) -> None: + """Add orphan 2D cell check subparser with its arguments. + + Args: + subparsers: The subparsers action to add the parser to. + """ + p = subparsers.add_parser( ORPHAN_2D, + help="Check if 2D cells are faces of 3D cells and identify orphaned 2D elements." ) + addVtuInputFileArgument( p ) + p.add_argument( '--' + __ORPHAN_OUTPUT, + type=str, + metavar='FILE', + help="[string]: Output VTU file for orphaned 2D cells (cells not matching any 3D face)." ) + p.add_argument( '--' + __CLEAN_OUTPUT, + type=str, + metavar='FILE', + help="[string]: Output VTU file with orphaned 2D cells removed." ) + p.add_argument( + '--' + __DATA_MODE, + type=str, + metavar=", ".join( __DATA_MODE_VALUES ), + default=__DATA_MODE_DEFAULT, + help='[string]: For ".vtu" output format, the data mode can be binary or ascii. Defaults to binary.' ) + + +def displayResults( options: Options, result: Result ) -> None: + """Display the results of the orphan 2D cell check. + + Args: + options: The options used for the check. + result: The result of the orphan check. + """ + setupLogger.results( "=" * 80 ) + setupLogger.results( "ORPHAN 2D CELL CHECK RESULTS" ) + setupLogger.results( "=" * 80 ) + setupLogger.results( f"Total 2D cells: {result.total2dCells}" ) + setupLogger.results( f"Total 3D cells: {result.total3dCells}" ) + setupLogger.results( f"2D cells matching 3D faces: {result.matched2dCells}" ) + setupLogger.results( f"Orphaned 2D cells: {result.orphaned2dCells}" ) + + if result.orphaned2dCells == 0: + setupLogger.results( "Validation PASSED: All 2D cells are faces of 3D cells" ) + else: + setupLogger.results( f"Validation FAILED: {result.orphaned2dCells} 2D cells are not faces of 3D cells" ) + + # Show first few orphaned cell indices for debugging + if result.orphaned2dIndices: + numToShow = min( 10, len( result.orphaned2dIndices ) ) + indices = ", ".join( map( str, result.orphaned2dIndices[ :numToShow ] ) ) + if len( result.orphaned2dIndices ) > numToShow: + indices += ", ..." + setupLogger.results( f" First orphaned cell indices: {indices}" ) + + setupLogger.results( "=" * 80 ) diff --git a/mesh-doctor/src/geos/mesh_doctor/register.py b/mesh-doctor/src/geos/mesh_doctor/register.py index bb677f4ee..3a39b705a 100644 --- a/mesh-doctor/src/geos/mesh_doctor/register.py +++ b/mesh-doctor/src/geos/mesh_doctor/register.py @@ -58,7 +58,7 @@ def registerParsingActions( for actionName in ( parsing.ALL_CHECKS, parsing.COLLOCATES_NODES, parsing.ELEMENT_VOLUMES, parsing.FIX_ELEMENTS_ORDERINGS, parsing.GENERATE_CUBE, parsing.GENERATE_FRACTURES, parsing.GENERATE_GLOBAL_IDS, parsing.MAIN_CHECKS, parsing.NON_CONFORMAL, - parsing.SELF_INTERSECTING_ELEMENTS, parsing.SUPPORTED_ELEMENTS ): + parsing.SELF_INTERSECTING_ELEMENTS, parsing.SUPPORTED_ELEMENTS, parsing.ORPHAN_2D ): __HELPERS[ actionName ] = actionName __ACTIONS[ actionName ] = actionName From 6a1c607f685df9ae8a977d3b22d4b4ef618f5a52 Mon Sep 17 00:00:00 2001 From: Bertrand Denel <120652669+bd713@users.noreply.github.com> Date: Wed, 4 Feb 2026 10:12:45 -0800 Subject: [PATCH 2/3] Apply suggestions from code review Co-authored-by: Jacques Franc <49998870+jafranc@users.noreply.github.com> --- mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py b/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py index 0a21a2bf7..32502ebfb 100644 --- a/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py @@ -1,6 +1,6 @@ # SPDX-License-Identifier: Apache-2.0 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. -# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +# SPDX-FileContributor: Bertrand Denel from dataclasses import dataclass from typing import Optional import vtk @@ -51,11 +51,14 @@ def getCellFacePoints( cell: vtk.vtkCell ) -> list[ tuple[ int, ...] ]: elif cellType == vtk.VTK_PYRAMID: # 1 quadrilateral + 4 triangular faces faces = [ [ 0, 1, 2, 3 ], [ 0, 1, 4 ], [ 1, 2, 4 ], [ 2, 3, 4 ], [ 3, 0, 4 ] ] + else: + raise NotImplementedError(f"Orphaned2d is not implemented for cell type {cellType}. It is supported for TETRAHEDRA ({vtk.VTK_TETRA}), HEXA {(vtk.VTK_HEXAHEDRON}), WEDGE ({vtk.VTK_WEDGE}) and PYRAMID ({vtk.VTK_PYRAMID})") # Convert local indices to global point IDs pointIds = cell.GetPointIds() globalFaces = [] for face in faces: + #sort to make comparison permutation invariant globalFace = tuple( sorted( [ pointIds.GetId( i ) for i in face ] ) ) globalFaces.append( globalFace ) From d2c2e18c6c664c465c665339fe0729598eff404e Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Thu, 5 Feb 2026 00:20:54 -0600 Subject: [PATCH 3/3] fix suggestion --- mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py b/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py index 32502ebfb..dfaef1182 100644 --- a/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py @@ -51,8 +51,10 @@ def getCellFacePoints( cell: vtk.vtkCell ) -> list[ tuple[ int, ...] ]: elif cellType == vtk.VTK_PYRAMID: # 1 quadrilateral + 4 triangular faces faces = [ [ 0, 1, 2, 3 ], [ 0, 1, 4 ], [ 1, 2, 4 ], [ 2, 3, 4 ], [ 3, 0, 4 ] ] - else: - raise NotImplementedError(f"Orphaned2d is not implemented for cell type {cellType}. It is supported for TETRAHEDRA ({vtk.VTK_TETRA}), HEXA {(vtk.VTK_HEXAHEDRON}), WEDGE ({vtk.VTK_WEDGE}) and PYRAMID ({vtk.VTK_PYRAMID})") + else: + raise NotImplementedError( + f"Orphan2d is not implemented for cell type {cellType}. It is supported for TETRAHEDRA ({vtk.VTK_TETRA}), HEXA ({vtk.VTK_HEXAHEDRON}), WEDGE ({vtk.VTK_WEDGE}) and PYRAMID ({vtk.VTK_PYRAMID})" + ) # Convert local indices to global point IDs pointIds = cell.GetPointIds()