From 400d48cb518afd05fb0a30b8740eb7a4415db476 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Mon, 2 Mar 2026 17:27:13 -0600 Subject: [PATCH 1/3] Draft --- .../src/geos/mesh_doctor/actions/euler.py | 291 ++++++++++++++++++ .../src/geos/mesh_doctor/parsing/__init__.py | 2 + .../geos/mesh_doctor/parsing/eulerParsing.py | 107 +++++++ mesh-doctor/src/geos/mesh_doctor/register.py | 3 +- 4 files changed, 402 insertions(+), 1 deletion(-) create mode 100644 mesh-doctor/src/geos/mesh_doctor/actions/euler.py create mode 100644 mesh-doctor/src/geos/mesh_doctor/parsing/eulerParsing.py diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/euler.py b/mesh-doctor/src/geos/mesh_doctor/actions/euler.py new file mode 100644 index 00000000..4f3c9f5d --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/actions/euler.py @@ -0,0 +1,291 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +"""Compute Euler Characteristic for mesh files (3D elements only).""" + +from dataclasses import dataclass +import vtk +from tqdm import tqdm + +from geos.mesh_doctor.parsing.cliParsing import setupLogger +from geos.mesh.io.vtkIO import readUnstructuredGrid + + +@dataclass( frozen=True ) +class Options: + """Options for Euler characteristic computation.""" + pass + + +@dataclass( frozen=True ) +class Result: + """Result of Euler characteristic computation. + + Attributes: + numVertices: Number of surface vertices (V) + numEdges: Number of surface edges (E) + numFaces: Number of surface faces (F) + eulerCharacteristic: Euler characteristic (χ = V - E + F) + num3dCells: Number of 3D volumetric cells in input + num2dCells: Number of 2D surface cells in input (ignored) + numOtherCells: Number of other cells in input + numBoundaryEdges: Number of boundary edges + numNonManifoldEdges: Number of non-manifold edges + """ + numVertices: int + numEdges: int + numFaces: int + eulerCharacteristic: int + num3dCells: int + num2dCells: int + numOtherCells: int + numBoundaryEdges: int + numNonManifoldEdges: int + + +def __filter3dElements( mesh: vtk.vtkUnstructuredGrid ) -> tuple[ vtk.vtkUnstructuredGrid, int, int, int, bool ]: + """Filter only 3D volumetric elements from unstructured grid. + + Removes 2D faces, 1D edges, and 0D vertices. + + Args: + mesh: Input unstructured grid + + Returns: + Tuple of (filtered_mesh, n_3d, n_2d, n_other, has_3d) + """ + # Cell types that are 3D volumes + volumeTypes = [ + vtk.VTK_TETRA, # 10 + vtk.VTK_HEXAHEDRON, # 12 + vtk.VTK_WEDGE, # 13 + vtk.VTK_PYRAMID, # 14 + vtk.VTK_VOXEL, # 11 + vtk.VTK_PENTAGONAL_PRISM, # 15 + vtk.VTK_HEXAGONAL_PRISM, # 16 + vtk.VTK_QUADRATIC_TETRA, # 24 + vtk.VTK_QUADRATIC_HEXAHEDRON, # 25 + vtk.VTK_QUADRATIC_WEDGE, # 26 + vtk.VTK_QUADRATIC_PYRAMID, # 27 + ] + + surfaceTypes = [ + vtk.VTK_TRIANGLE, vtk.VTK_QUAD, vtk.VTK_POLYGON, vtk.VTK_QUADRATIC_TRIANGLE, vtk.VTK_QUADRATIC_QUAD + ] + + # Count cell types + n3d = 0 + n2d = 0 + nOther = 0 + + setupLogger.info( "Classifying cell types..." ) + for i in tqdm( range( mesh.GetNumberOfCells() ), desc="Scanning cells" ): + cellType = mesh.GetCellType( i ) + + if cellType in volumeTypes: + n3d += 1 + elif cellType in surfaceTypes: + n2d += 1 + else: + nOther += 1 + + setupLogger.info( f"Cell type breakdown:" ) + setupLogger.info( f" 3D volumetric cells: {n3d}" ) + setupLogger.info( f" 2D surface cells: {n2d}" ) + setupLogger.info( f" Other cells: {nOther}" ) + + # Check if we have 3D cells + has3d = n3d > 0 + + if not has3d: + setupLogger.warning( "No 3D volumetric elements found!" ) + setupLogger.warning( "This appears to be a pure surface mesh." ) + setupLogger.warning( "Computing Euler characteristic from all cells..." ) + return mesh, n3d, n2d, nOther, False + + if n2d > 0: + setupLogger.info( f"Filtering out {n2d} 2D boundary cells..." ) + setupLogger.info( f"Using only {n3d} volumetric elements for surface extraction" ) + + # Add cell type array for filtering + cellTypes = vtk.vtkIntArray() + cellTypes.SetName( "CellType" ) + cellTypes.SetNumberOfTuples( mesh.GetNumberOfCells() ) + + for i in range( mesh.GetNumberOfCells() ): + cellTypes.SetValue( i, mesh.GetCellType( i ) ) + + mesh.GetCellData().AddArray( cellTypes ) + + # Threshold filter to extract only 3D cells + threshold = vtk.vtkThreshold() + threshold.SetInputData( mesh ) + threshold.SetInputArrayToProcess( 0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS, "CellType" ) + threshold.SetThresholdFunction( vtk.vtkThreshold.THRESHOLD_BETWEEN ) + threshold.SetLowerThreshold( 10 ) # Minimum 3D cell type + threshold.SetUpperThreshold( 100 ) # Maximum reasonable cell type + threshold.Update() + + filteredMesh = threshold.GetOutput() + + return filteredMesh, n3d, n2d, nOther, has3d + + +def __extractSurface( mesh: vtk.vtkUnstructuredGrid ) -> vtk.vtkPolyData: + """Extract surface from unstructured grid (3D elements only). + + Args: + mesh: Input unstructured grid (3D elements) + + Returns: + Surface as polydata + """ + setupLogger.info( "Extracting surface from 3D elements..." ) + surfaceFilter = vtk.vtkDataSetSurfaceFilter() + surfaceFilter.SetInputData( mesh ) + surfaceFilter.Update() + return surfaceFilter.GetOutput() + + +def __countUniqueEdges( surface: vtk.vtkPolyData ) -> int: + """Count unique edges in the surface mesh. + + Args: + surface: Surface mesh + + Returns: + Number of unique edges + """ + setupLogger.info( "Counting unique edges..." ) + edges = set() + nFaces = surface.GetNumberOfCells() + + for i in tqdm( range( nFaces ), desc="Processing faces" ): + cell = surface.GetCell( i ) + nEdges = cell.GetNumberOfEdges() + for j in range( nEdges ): + edge = cell.GetEdge( j ) + pt1 = edge.GetPointId( 0 ) + pt2 = edge.GetPointId( 1 ) + edges.add( tuple( sorted( [ pt1, pt2 ] ) ) ) + + return len( edges ) + + +def __checkMeshQuality( surface: vtk.vtkPolyData ) -> tuple[ int, int ]: + """Check for boundary edges and non-manifold features. + + Args: + surface: Surface mesh + + Returns: + Tuple of (boundary_edges, non_manifold_edges) + """ + setupLogger.info( "Checking mesh quality..." ) + + # Count boundary edges + featureEdgesBoundary = vtk.vtkFeatureEdges() + featureEdgesBoundary.SetInputData( surface ) + featureEdgesBoundary.BoundaryEdgesOn() + featureEdgesBoundary.ManifoldEdgesOff() + featureEdgesBoundary.NonManifoldEdgesOff() + featureEdgesBoundary.FeatureEdgesOff() + featureEdgesBoundary.Update() + boundaryEdges = featureEdgesBoundary.GetOutput().GetNumberOfCells() + + # Count non-manifold edges + featureEdgesNm = vtk.vtkFeatureEdges() + featureEdgesNm.SetInputData( surface ) + featureEdgesNm.BoundaryEdgesOff() + featureEdgesNm.ManifoldEdgesOff() + featureEdgesNm.NonManifoldEdgesOn() + featureEdgesNm.FeatureEdgesOff() + featureEdgesNm.Update() + nonManifoldEdges = featureEdgesNm.GetOutput().GetNumberOfCells() + + return boundaryEdges, nonManifoldEdges + + +def meshAction( mesh: vtk.vtkUnstructuredGrid, options: Options ) -> Result: + """Compute Euler characteristic for a mesh. + + Only considers 3D volumetric elements. Extracts surface and computes V - E + F. + + Args: + mesh: Input unstructured grid + options: Computation options + + Returns: + Result with Euler characteristic and topology information + """ + setupLogger.info( "Starting Euler characteristic computation..." ) + setupLogger.info( f"Input mesh: {mesh.GetNumberOfPoints()} points, {mesh.GetNumberOfCells()} cells" ) + + # Filter to 3D elements only + mesh3d, n3d, n2d, nOther, has3d = __filter3dElements( mesh ) + + # Extract surface from 3D elements + surface = __extractSurface( mesh3d ) + + # Get counts + V = surface.GetNumberOfPoints() + F = surface.GetNumberOfCells() + E = __countUniqueEdges( surface ) + + setupLogger.info( f"Surface topology:" ) + setupLogger.info( f" Vertices (V): {V:,}" ) + setupLogger.info( f" Edges (E): {E:,}" ) + setupLogger.info( f" Faces (F): {F:,}" ) + + # Calculate Euler characteristic + euler = V - E + F + + setupLogger.info( f"Euler characteristic (χ = V - E + F): {euler}" ) + + # Interpret result + if euler == 2: + setupLogger.info( "Topology: Closed surface (sphere-like)" ) + elif euler == 0: + setupLogger.info( "Topology: Torus-like or open cylinder" ) + elif euler == 1: + setupLogger.info( "Topology: Disk-like (open surface)" ) + else: + genus = ( 2 - euler ) / 2 + setupLogger.info( f"Topology: Complex (genus g ≈ {genus:.1f})" ) + + # Check mesh quality + boundaryEdges, nonManifoldEdges = __checkMeshQuality( surface ) + + setupLogger.info( f"Mesh quality:" ) + setupLogger.info( f" Boundary edges: {boundaryEdges:,}" ) + setupLogger.info( f" Non-manifold edges: {nonManifoldEdges:,}" ) + + if boundaryEdges == 0 and nonManifoldEdges == 0: + setupLogger.info( " Perfect closed manifold mesh!" ) + elif boundaryEdges > 0: + setupLogger.warning( " Open surface detected (has boundaries)" ) + if nonManifoldEdges > 0: + setupLogger.warning( " Non-manifold edges detected (mesh has issues!)" ) + + return Result( numVertices=V, + numEdges=E, + numFaces=F, + eulerCharacteristic=euler, + num3dCells=n3d, + num2dCells=n2d, + numOtherCells=nOther, + numBoundaryEdges=boundaryEdges, + numNonManifoldEdges=nonManifoldEdges ) + + +def action( vtuInputFile: str, options: Options ) -> Result: + """Compute Euler characteristic for a VTU file. + + Args: + vtuInputFile: Path to input VTU file + options: Computation options + + Returns: + Result with Euler characteristic and topology information + """ + mesh = readUnstructuredGrid( vtuInputFile ) + return meshAction( mesh, options ) diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py b/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py index 2aca9fcd..ecdb80d9 100644 --- a/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/__init__.py @@ -17,6 +17,8 @@ SELF_INTERSECTING_ELEMENTS = "selfIntersectingElements" SUPPORTED_ELEMENTS = "supportedElements" ORPHAN_2D = "orphan2d" +CHECK_INTERNAL_TAGS = "checkInternalTags" +EULER = "euler" @dataclass( frozen=True ) diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/eulerParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/eulerParsing.py new file mode 100644 index 00000000..dbf9fa6f --- /dev/null +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/eulerParsing.py @@ -0,0 +1,107 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. +"""Command line parsing for Euler characteristic computation.""" + +from __future__ import annotations +import argparse +from typing import Any, Optional + +from geos.mesh_doctor.actions.euler import Options, Result, action +from geos.mesh_doctor.parsing import EULER +from geos.mesh_doctor.parsing.cliParsing import setupLogger, addVtuInputFileArgument + + +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 Euler computation. + """ + return Options() + + +def fillSubparser( subparsers: argparse._SubParsersAction[ Any ] ) -> None: + """Fill the argument parser for the Euler characteristic action. + + Args: + subparsers: Subparsers from the main argument parser + """ + p = subparsers.add_parser( + EULER, + help="Compute Euler characteristic (χ = V - E + F) for mesh topology analysis.", + description="""\ +Computes the Euler characteristic for a mesh by: + 1. Filtering to 3D volumetric elements only (ignores 2D boundary cells) + 2. Extracting the surface from these 3D elements + 3. Counting vertices (V), edges (E), and faces (F) + 4. Computing χ = V - E + F + +Expected Euler values: + χ = 2 → Closed surface (sphere-like, e.g., cube, tetrahedron) + χ = 0 → Torus or cylinder topology + χ = 1 → Disk or open surface + Other → Complex topology with genus g = (2 - χ)/2 + +This tool helps verify that meshes form proper closed volumes for simulation. +""" ) + + addVtuInputFileArgument( p ) + + +def displayResults( options: Options, result: Result ) -> None: + """Display the results of the Euler characteristic computation. + + Args: + options: The options used for the computation. + result: The result of the computation. + """ + setupLogger.results( "=" * 80 ) + setupLogger.results( "EULER CHARACTERISTIC RESULTS" ) + setupLogger.results( "=" * 80 ) + setupLogger.results( f"Mesh Topology (from {result.num3dCells} 3D elements):" ) + setupLogger.results( f" Surface Vertices (V): {result.numVertices:,}" ) + setupLogger.results( f" Surface Edges (E): {result.numEdges:,}" ) + setupLogger.results( f" Surface Faces (F): {result.numFaces:,}" ) + setupLogger.results( "-" * 80 ) + setupLogger.results( f" Euler Characteristic (χ = V - E + F): {result.eulerCharacteristic}" ) + setupLogger.results( "=" * 80 ) + + # Interpret topology + if result.eulerCharacteristic == 2: + setupLogger.results( "Topology: Closed surface (sphere-like)" ) + elif result.eulerCharacteristic == 0: + setupLogger.results( "Topology: Torus-like or open cylinder" ) + elif result.eulerCharacteristic == 1: + setupLogger.results( "Topology: Disk-like (open surface)" ) + else: + genus = ( 2 - result.eulerCharacteristic ) / 2 + setupLogger.results( f"Topology: Complex (genus g ≈ {genus:.1f})" ) + + # Show input cell breakdown if there were 2D cells filtered out + if result.num2dCells > 0: + setupLogger.results( "" ) + setupLogger.results( f"Input cell breakdown:" ) + setupLogger.results( f" 3D volumetric cells: {result.num3dCells} (used for surface extraction)" ) + setupLogger.results( f" 2D surface cells: {result.num2dCells} (filtered out)" ) + if result.numOtherCells > 0: + setupLogger.results( f" Other cells: {result.numOtherCells} (filtered out)" ) + + # Quality check results (always performed) + setupLogger.results( "" ) + setupLogger.results( "Mesh Quality Check:" ) + setupLogger.results( "-" * 80 ) + setupLogger.results( f" Boundary edges: {result.numBoundaryEdges:,}" ) + setupLogger.results( f" Non-manifold edges: {result.numNonManifoldEdges:,}" ) + + if result.numBoundaryEdges == 0 and result.numNonManifoldEdges == 0: + setupLogger.results( " Status: Perfect closed manifold mesh!" ) + else: + if result.numBoundaryEdges > 0: + setupLogger.results( " Status: Open surface detected (has boundaries)" ) + if result.numNonManifoldEdges > 0: + setupLogger.results( " Status: Non-manifold edges detected (mesh has issues!)" ) + + setupLogger.results( "=" * 80 ) diff --git a/mesh-doctor/src/geos/mesh_doctor/register.py b/mesh-doctor/src/geos/mesh_doctor/register.py index 3a39b705..99675376 100644 --- a/mesh-doctor/src/geos/mesh_doctor/register.py +++ b/mesh-doctor/src/geos/mesh_doctor/register.py @@ -58,7 +58,8 @@ 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.ORPHAN_2D ): + parsing.SELF_INTERSECTING_ELEMENTS, parsing.SUPPORTED_ELEMENTS, parsing.ORPHAN_2D, + parsing.CHECK_INTERNAL_TAGS, parsing.EULER ): __HELPERS[ actionName ] = actionName __ACTIONS[ actionName ] = actionName From 7724055f11c8fcc3c866346dd8732227a6b87316 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Mon, 2 Mar 2026 19:44:44 -0600 Subject: [PATCH 2/3] clean --- .../src/geos/mesh_doctor/actions/euler.py | 202 ++++++++++-------- .../src/geos/mesh_doctor/actions/orphan2d.py | 5 +- .../geos/mesh_doctor/parsing/eulerParsing.py | 93 +++++--- 3 files changed, 184 insertions(+), 116 deletions(-) diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/euler.py b/mesh-doctor/src/geos/mesh_doctor/actions/euler.py index 4f3c9f5d..6c91ac23 100644 --- a/mesh-doctor/src/geos/mesh_doctor/actions/euler.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/euler.py @@ -19,17 +19,18 @@ class Options: @dataclass( frozen=True ) class Result: """Result of Euler characteristic computation. - + Attributes: numVertices: Number of surface vertices (V) numEdges: Number of surface edges (E) numFaces: Number of surface faces (F) - eulerCharacteristic: Euler characteristic (χ = V - E + F) + eulerCharacteristic: Euler characteristic (X = V - E + F) num3dCells: Number of 3D volumetric cells in input num2dCells: Number of 2D surface cells in input (ignored) numOtherCells: Number of other cells in input numBoundaryEdges: Number of boundary edges numNonManifoldEdges: Number of non-manifold edges + numConnectedComponents: Number of disconnected mesh regions """ numVertices: int numEdges: int @@ -40,55 +41,64 @@ class Result: numOtherCells: int numBoundaryEdges: int numNonManifoldEdges: int + numConnectedComponents: int + + +def __countConnectedComponents( mesh: vtk.vtkUnstructuredGrid ) -> int: + """Count number of disconnected mesh components. + + Args: + mesh: Input unstructured grid + + Returns: + Number of disconnected regions + """ + setupLogger.info( "Checking for disconnected components..." ) + + connectivity = vtk.vtkConnectivityFilter() + connectivity.SetInputData( mesh ) + connectivity.SetExtractionModeToAllRegions() + connectivity.ColorRegionsOn() + connectivity.Update() + + numRegions = connectivity.GetNumberOfExtractedRegions() + + setupLogger.info( f"Found {numRegions} disconnected component(s)" ) + + return numRegions def __filter3dElements( mesh: vtk.vtkUnstructuredGrid ) -> tuple[ vtk.vtkUnstructuredGrid, int, int, int, bool ]: """Filter only 3D volumetric elements from unstructured grid. - + Removes 2D faces, 1D edges, and 0D vertices. - + Args: mesh: Input unstructured grid - + Returns: Tuple of (filtered_mesh, n_3d, n_2d, n_other, has_3d) """ - # Cell types that are 3D volumes - volumeTypes = [ - vtk.VTK_TETRA, # 10 - vtk.VTK_HEXAHEDRON, # 12 - vtk.VTK_WEDGE, # 13 - vtk.VTK_PYRAMID, # 14 - vtk.VTK_VOXEL, # 11 - vtk.VTK_PENTAGONAL_PRISM, # 15 - vtk.VTK_HEXAGONAL_PRISM, # 16 - vtk.VTK_QUADRATIC_TETRA, # 24 - vtk.VTK_QUADRATIC_HEXAHEDRON, # 25 - vtk.VTK_QUADRATIC_WEDGE, # 26 - vtk.VTK_QUADRATIC_PYRAMID, # 27 - ] - - surfaceTypes = [ - vtk.VTK_TRIANGLE, vtk.VTK_QUAD, vtk.VTK_POLYGON, vtk.VTK_QUADRATIC_TRIANGLE, vtk.VTK_QUADRATIC_QUAD - ] - - # Count cell types + # Classify cells by dimension + cell3dIds = [] n3d = 0 n2d = 0 nOther = 0 setupLogger.info( "Classifying cell types..." ) for i in tqdm( range( mesh.GetNumberOfCells() ), desc="Scanning cells" ): - cellType = mesh.GetCellType( i ) + cell = mesh.GetCell( i ) + cellDim = cell.GetCellDimension() - if cellType in volumeTypes: + if cellDim == 3: + cell3dIds.append( i ) n3d += 1 - elif cellType in surfaceTypes: + elif cellDim == 2: n2d += 1 else: nOther += 1 - setupLogger.info( f"Cell type breakdown:" ) + setupLogger.info( "Cell type breakdown:" ) setupLogger.info( f" 3D volumetric cells: {n3d}" ) setupLogger.info( f" 2D surface cells: {n2d}" ) setupLogger.info( f" Other cells: {nOther}" ) @@ -106,36 +116,27 @@ def __filter3dElements( mesh: vtk.vtkUnstructuredGrid ) -> tuple[ vtk.vtkUnstruc setupLogger.info( f"Filtering out {n2d} 2D boundary cells..." ) setupLogger.info( f"Using only {n3d} volumetric elements for surface extraction" ) - # Add cell type array for filtering - cellTypes = vtk.vtkIntArray() - cellTypes.SetName( "CellType" ) - cellTypes.SetNumberOfTuples( mesh.GetNumberOfCells() ) - - for i in range( mesh.GetNumberOfCells() ): - cellTypes.SetValue( i, mesh.GetCellType( i ) ) - - mesh.GetCellData().AddArray( cellTypes ) + # Extract only 3D cells using vtkExtractCells + idList = vtk.vtkIdList() + for cellId in cell3dIds: + idList.InsertNextId( cellId ) - # Threshold filter to extract only 3D cells - threshold = vtk.vtkThreshold() - threshold.SetInputData( mesh ) - threshold.SetInputArrayToProcess( 0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS, "CellType" ) - threshold.SetThresholdFunction( vtk.vtkThreshold.THRESHOLD_BETWEEN ) - threshold.SetLowerThreshold( 10 ) # Minimum 3D cell type - threshold.SetUpperThreshold( 100 ) # Maximum reasonable cell type - threshold.Update() + extractor = vtk.vtkExtractCells() + extractor.SetInputData( mesh ) + extractor.SetCellList( idList ) + extractor.Update() - filteredMesh = threshold.GetOutput() + filteredMesh = extractor.GetOutput() return filteredMesh, n3d, n2d, nOther, has3d def __extractSurface( mesh: vtk.vtkUnstructuredGrid ) -> vtk.vtkPolyData: """Extract surface from unstructured grid (3D elements only). - + Args: mesh: Input unstructured grid (3D elements) - + Returns: Surface as polydata """ @@ -148,35 +149,29 @@ def __extractSurface( mesh: vtk.vtkUnstructuredGrid ) -> vtk.vtkPolyData: def __countUniqueEdges( surface: vtk.vtkPolyData ) -> int: """Count unique edges in the surface mesh. - + Args: surface: Surface mesh - + Returns: Number of unique edges """ setupLogger.info( "Counting unique edges..." ) - edges = set() - nFaces = surface.GetNumberOfCells() - for i in tqdm( range( nFaces ), desc="Processing faces" ): - cell = surface.GetCell( i ) - nEdges = cell.GetNumberOfEdges() - for j in range( nEdges ): - edge = cell.GetEdge( j ) - pt1 = edge.GetPointId( 0 ) - pt2 = edge.GetPointId( 1 ) - edges.add( tuple( sorted( [ pt1, pt2 ] ) ) ) + # Use VTK's edge extraction (faster than manual iteration) + edgeExtractor = vtk.vtkExtractEdges() + edgeExtractor.SetInputData( surface ) + edgeExtractor.Update() - return len( edges ) + return edgeExtractor.GetOutput().GetNumberOfCells() def __checkMeshQuality( surface: vtk.vtkPolyData ) -> tuple[ int, int ]: """Check for boundary edges and non-manifold features. - + Args: surface: Surface mesh - + Returns: Tuple of (boundary_edges, non_manifold_edges) """ @@ -207,13 +202,13 @@ def __checkMeshQuality( surface: vtk.vtkPolyData ) -> tuple[ int, int ]: def meshAction( mesh: vtk.vtkUnstructuredGrid, options: Options ) -> Result: """Compute Euler characteristic for a mesh. - + Only considers 3D volumetric elements. Extracts surface and computes V - E + F. - + Args: mesh: Input unstructured grid options: Computation options - + Returns: Result with Euler characteristic and topology information """ @@ -223,6 +218,9 @@ def meshAction( mesh: vtk.vtkUnstructuredGrid, options: Options ) -> Result: # Filter to 3D elements only mesh3d, n3d, n2d, nOther, has3d = __filter3dElements( mesh ) + # Count connected components BEFORE extracting surface + numComponents = __countConnectedComponents( mesh3d ) + # Extract surface from 3D elements surface = __extractSurface( mesh3d ) @@ -231,7 +229,7 @@ def meshAction( mesh: vtk.vtkUnstructuredGrid, options: Options ) -> Result: F = surface.GetNumberOfCells() E = __countUniqueEdges( surface ) - setupLogger.info( f"Surface topology:" ) + setupLogger.info( "Surface topology:" ) setupLogger.info( f" Vertices (V): {V:,}" ) setupLogger.info( f" Edges (E): {E:,}" ) setupLogger.info( f" Faces (F): {F:,}" ) @@ -239,32 +237,65 @@ def meshAction( mesh: vtk.vtkUnstructuredGrid, options: Options ) -> Result: # Calculate Euler characteristic euler = V - E + F - setupLogger.info( f"Euler characteristic (χ = V - E + F): {euler}" ) - - # Interpret result - if euler == 2: - setupLogger.info( "Topology: Closed surface (sphere-like)" ) - elif euler == 0: - setupLogger.info( "Topology: Torus-like or open cylinder" ) - elif euler == 1: - setupLogger.info( "Topology: Disk-like (open surface)" ) + setupLogger.info( f"Euler characteristic (X = V - E + F): {euler}" ) + setupLogger.info( f"Connected components: {numComponents}" ) + + # Interpret result with component information + if numComponents == 1: + # Single component - standard topology interpretation + if euler == 2: + setupLogger.info( "Topology: Single closed surface (sphere-like) - VALID for simulation!" ) + elif euler == 0: + setupLogger.warning( "Topology: Torus-like (genus 1)" ) + setupLogger.warning( " Expected X = 2 for standard simulation mesh" ) + elif euler == 1: + setupLogger.warning( "Topology: Disk-like (open surface)" ) + setupLogger.warning( " Mesh is not closed! Expected X = 2" ) + elif euler < 0: + genus = ( 2 - euler ) // 2 + setupLogger.warning( f"Topology: Complex surface with handles (genus g = {genus})" ) + setupLogger.warning( " Expected X = 2 for standard simulation mesh" ) + else: + setupLogger.warning( f"Topology: Unexpected (X = {euler})" ) else: - genus = ( 2 - euler ) / 2 - setupLogger.info( f"Topology: Complex (genus g ≈ {genus:.1f})" ) + # Multiple components detected + expectedEuler = 2 * numComponents # If all components are sphere-like + + setupLogger.error( f"Topology: Mesh has {numComponents} DISCONNECTED COMPONENTS!" ) + setupLogger.error( f" Euler characteristic: {euler}" ) + setupLogger.error( f" Expected for {numComponents} spheres: X = {expectedEuler}" ) + + if numComponents > 1: + # Internal cavities = additional components beyond the outer surface + numCavities = numComponents - 1 + setupLogger.error( f" Likely {numCavities} internal void(s)/cavity(ies)" ) + + if euler == expectedEuler: + setupLogger.error( " All components are sphere-like (X=2 each)" ) + elif euler < expectedEuler: + setupLogger.error( " Some components have genus > 0 (holes/handles)" ) + else: + setupLogger.error( " Unusual topology - verify mesh integrity" ) + + setupLogger.error( " This mesh is NOT suitable for simulation without fixing!" ) + setupLogger.error( " Expected: single closed volume (X = 2, components = 1)" ) # Check mesh quality boundaryEdges, nonManifoldEdges = __checkMeshQuality( surface ) - setupLogger.info( f"Mesh quality:" ) + setupLogger.info( "Mesh quality:" ) setupLogger.info( f" Boundary edges: {boundaryEdges:,}" ) setupLogger.info( f" Non-manifold edges: {nonManifoldEdges:,}" ) if boundaryEdges == 0 and nonManifoldEdges == 0: - setupLogger.info( " Perfect closed manifold mesh!" ) + if euler == 2 and numComponents == 1: + setupLogger.info( " Perfect closed manifold mesh - READY for simulation!" ) + else: + setupLogger.warning( " Perfect manifold but topology issues - verify before simulation" ) elif boundaryEdges > 0: setupLogger.warning( " Open surface detected (has boundaries)" ) if nonManifoldEdges > 0: - setupLogger.warning( " Non-manifold edges detected (mesh has issues!)" ) + setupLogger.error( " Non-manifold edges detected (mesh has issues!)" ) return Result( numVertices=V, numEdges=E, @@ -274,16 +305,17 @@ def meshAction( mesh: vtk.vtkUnstructuredGrid, options: Options ) -> Result: num2dCells=n2d, numOtherCells=nOther, numBoundaryEdges=boundaryEdges, - numNonManifoldEdges=nonManifoldEdges ) + numNonManifoldEdges=nonManifoldEdges, + numConnectedComponents=numComponents ) def action( vtuInputFile: str, options: Options ) -> Result: """Compute Euler characteristic for a VTU file. - + Args: vtuInputFile: Path to input VTU file options: Computation options - + Returns: Result with Euler characteristic and topology information """ diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py b/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py index dfaef118..1c3fbe4f 100644 --- a/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/orphan2d.py @@ -67,7 +67,8 @@ def getCellFacePoints( cell: vtk.vtkCell ) -> list[ tuple[ int, ...] ]: return globalFaces -def buildMeshSubset( mesh: vtkUnstructuredGrid, cellIndices: list[ int ], description: str ) -> vtkUnstructuredGrid: +def buildMeshSubset( mesh: vtkUnstructuredGrid, cellIndices: list[ int ], + description: str ) -> Optional[ vtkUnstructuredGrid ]: """Build a new mesh containing only the specified cells. Args: @@ -76,7 +77,7 @@ def buildMeshSubset( mesh: vtkUnstructuredGrid, cellIndices: list[ int ], descri description: Description for progress messages. Returns: - A new vtkUnstructuredGrid containing only the specified cells. + A new vtkUnstructuredGrid containing only the specified cells, or None if no cells. """ if not cellIndices: setupLogger.info( f"No {description} to create." ) diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/eulerParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/eulerParsing.py index dbf9fa6f..d061ba54 100644 --- a/mesh-doctor/src/geos/mesh_doctor/parsing/eulerParsing.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/eulerParsing.py @@ -4,9 +4,9 @@ from __future__ import annotations import argparse -from typing import Any, Optional +from typing import Any -from geos.mesh_doctor.actions.euler import Options, Result, action +from geos.mesh_doctor.actions.euler import Options, Result from geos.mesh_doctor.parsing import EULER from geos.mesh_doctor.parsing.cliParsing import setupLogger, addVtuInputFileArgument @@ -29,21 +29,20 @@ def fillSubparser( subparsers: argparse._SubParsersAction[ Any ] ) -> None: Args: subparsers: Subparsers from the main argument parser """ - p = subparsers.add_parser( - EULER, - help="Compute Euler characteristic (χ = V - E + F) for mesh topology analysis.", - description="""\ + p = subparsers.add_parser( EULER, + help="Compute Euler characteristic (X = V - E + F) for mesh topology analysis.", + description="""\ Computes the Euler characteristic for a mesh by: 1. Filtering to 3D volumetric elements only (ignores 2D boundary cells) 2. Extracting the surface from these 3D elements 3. Counting vertices (V), edges (E), and faces (F) - 4. Computing χ = V - E + F + 4. Computing X = V - E + F Expected Euler values: - χ = 2 → Closed surface (sphere-like, e.g., cube, tetrahedron) - χ = 0 → Torus or cylinder topology - χ = 1 → Disk or open surface - Other → Complex topology with genus g = (2 - χ)/2 + X = 2 Closed surface (sphere-like, e.g., cube, tetrahedron) + X = 0 Torus or cylinder topology + X = 1 Disk or open surface + Other Complex topology with genus g = (2 - X)/2 This tool helps verify that meshes form proper closed volumes for simulation. """ ) @@ -61,35 +60,68 @@ def displayResults( options: Options, result: Result ) -> None: setupLogger.results( "=" * 80 ) setupLogger.results( "EULER CHARACTERISTIC RESULTS" ) setupLogger.results( "=" * 80 ) - setupLogger.results( f"Mesh Topology (from {result.num3dCells} 3D elements):" ) + setupLogger.results( f"Mesh Topology (from {result.num3dCells:,} 3D elements):" ) setupLogger.results( f" Surface Vertices (V): {result.numVertices:,}" ) setupLogger.results( f" Surface Edges (E): {result.numEdges:,}" ) setupLogger.results( f" Surface Faces (F): {result.numFaces:,}" ) setupLogger.results( "-" * 80 ) - setupLogger.results( f" Euler Characteristic (χ = V - E + F): {result.eulerCharacteristic}" ) + setupLogger.results( f" Euler Characteristic (X = V - E + F): {result.eulerCharacteristic}" ) + setupLogger.results( f" Connected Components: {result.numConnectedComponents}" ) setupLogger.results( "=" * 80 ) - # Interpret topology - if result.eulerCharacteristic == 2: - setupLogger.results( "Topology: Closed surface (sphere-like)" ) - elif result.eulerCharacteristic == 0: - setupLogger.results( "Topology: Torus-like or open cylinder" ) - elif result.eulerCharacteristic == 1: - setupLogger.results( "Topology: Disk-like (open surface)" ) + # Interpret topology with validation + if result.numConnectedComponents == 1: + # Single component - standard interpretation + if result.eulerCharacteristic == 2: + setupLogger.results( "VALIDATION PASSED: Single closed surface (sphere-like)" ) + setupLogger.results( " Status: MESH IS VALID FOR SIMULATION" ) + elif result.eulerCharacteristic == 0: + setupLogger.results( "WARNING: Torus-like topology (genus 1)" ) + setupLogger.results( " Expected X = 2 for standard simulation mesh" ) + elif result.eulerCharacteristic == 1: + setupLogger.results( "WARNING: Open surface (not closed)" ) + setupLogger.results( " Mesh has boundary - expected X = 2" ) + elif result.eulerCharacteristic < 0: + genus = ( 2 - result.eulerCharacteristic ) // 2 + setupLogger.results( f"WARNING: Complex topology with handles (genus g = {genus})" ) + setupLogger.results( f" Surface has {genus} hole(s)/handle(s)" ) + setupLogger.results( " Expected X = 2 for standard simulation mesh" ) + else: + setupLogger.results( f"WARNING: Unexpected topology (X = {result.eulerCharacteristic})" ) else: - genus = ( 2 - result.eulerCharacteristic ) / 2 - setupLogger.results( f"Topology: Complex (genus g ≈ {genus:.1f})" ) + # Multiple components - indicates disconnected regions or internal cavities + expectedEuler = 2 * result.numConnectedComponents + numCavities = result.numConnectedComponents - 1 + + setupLogger.results( f"VALIDATION FAILED: {result.numConnectedComponents} DISCONNECTED COMPONENTS DETECTED!" ) + setupLogger.results( f" Euler characteristic: {result.eulerCharacteristic}" ) + setupLogger.results( f" Expected (if all spherical): X = {expectedEuler}" ) + setupLogger.results( f" Likely internal void(s)/cavity(ies): {numCavities}" ) + setupLogger.results( "" ) + + if result.eulerCharacteristic == expectedEuler: + setupLogger.results( " All components are sphere-like (X=2 each)" ) + setupLogger.results( " Issue: Mesh has internal hollow regions (outer + inner surfaces)" ) + elif result.eulerCharacteristic < expectedEuler: + avgGenus = ( expectedEuler - result.eulerCharacteristic ) / 2 + setupLogger.results( f" Some components have handles/holes (avg genus ~{avgGenus:.1f})" ) + else: + setupLogger.results( " Unusual topology - verify mesh integrity" ) + + setupLogger.results( "" ) + setupLogger.results( " NOT SUITABLE FOR SIMULATION without fixing!" ) + setupLogger.results( " Expected: Single closed volume (X = 2, components = 1)" ) # Show input cell breakdown if there were 2D cells filtered out if result.num2dCells > 0: setupLogger.results( "" ) - setupLogger.results( f"Input cell breakdown:" ) - setupLogger.results( f" 3D volumetric cells: {result.num3dCells} (used for surface extraction)" ) - setupLogger.results( f" 2D surface cells: {result.num2dCells} (filtered out)" ) + setupLogger.results( "Input cell breakdown:" ) + setupLogger.results( f" 3D volumetric cells: {result.num3dCells:,} (used for surface extraction)" ) + setupLogger.results( f" 2D surface cells: {result.num2dCells:,} (filtered out)" ) if result.numOtherCells > 0: - setupLogger.results( f" Other cells: {result.numOtherCells} (filtered out)" ) + setupLogger.results( f" Other cells: {result.numOtherCells:,} (filtered out)" ) - # Quality check results (always performed) + # Quality check results setupLogger.results( "" ) setupLogger.results( "Mesh Quality Check:" ) setupLogger.results( "-" * 80 ) @@ -97,11 +129,14 @@ def displayResults( options: Options, result: Result ) -> None: setupLogger.results( f" Non-manifold edges: {result.numNonManifoldEdges:,}" ) if result.numBoundaryEdges == 0 and result.numNonManifoldEdges == 0: - setupLogger.results( " Status: Perfect closed manifold mesh!" ) + if result.eulerCharacteristic == 2 and result.numConnectedComponents == 1: + setupLogger.results( " Status: PERFECT - Ready for simulation!" ) + else: + setupLogger.results( " Status: Manifold but topology issues (X != 2 or internal cavities)" ) else: if result.numBoundaryEdges > 0: setupLogger.results( " Status: Open surface detected (has boundaries)" ) if result.numNonManifoldEdges > 0: - setupLogger.results( " Status: Non-manifold edges detected (mesh has issues!)" ) + setupLogger.results( " Status: Non-manifold edges detected!" ) setupLogger.results( "=" * 80 ) From dc4e51b8bbdb6ed5d1a7b943539f7cd32fe429da Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Mon, 2 Mar 2026 23:41:58 -0600 Subject: [PATCH 3/3] changed from surface to solid euler --- .../src/geos/mesh_doctor/actions/euler.py | 251 +++++++++++------- .../geos/mesh_doctor/parsing/eulerParsing.py | 225 ++++++++++------ 2 files changed, 298 insertions(+), 178 deletions(-) diff --git a/mesh-doctor/src/geos/mesh_doctor/actions/euler.py b/mesh-doctor/src/geos/mesh_doctor/actions/euler.py index 6c91ac23..18394f70 100644 --- a/mesh-doctor/src/geos/mesh_doctor/actions/euler.py +++ b/mesh-doctor/src/geos/mesh_doctor/actions/euler.py @@ -1,6 +1,6 @@ # SPDX-License-Identifier: Apache-2.0 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. -"""Compute Euler Characteristic for mesh files (3D elements only).""" +"""Compute Solid Euler Characteristic for mesh files (3D elements only).""" from dataclasses import dataclass import vtk @@ -18,24 +18,26 @@ class Options: @dataclass( frozen=True ) class Result: - """Result of Euler characteristic computation. + """Result of solid Euler characteristic computation. Attributes: - numVertices: Number of surface vertices (V) - numEdges: Number of surface edges (E) - numFaces: Number of surface faces (F) - eulerCharacteristic: Euler characteristic (X = V - E + F) + numVertices: Number of vertices (V) in 3D mesh + numEdges: Number of unique edges (E) in 3D mesh + numFaces: Number of unique faces (F) in 3D mesh + numCells: Number of 3D cells (C) + solidEulerCharacteristic: Solid Euler characteristic (chi = V - E + F - C) num3dCells: Number of 3D volumetric cells in input num2dCells: Number of 2D surface cells in input (ignored) numOtherCells: Number of other cells in input - numBoundaryEdges: Number of boundary edges - numNonManifoldEdges: Number of non-manifold edges - numConnectedComponents: Number of disconnected mesh regions + numBoundaryEdges: Number of boundary edges on surface + numNonManifoldEdges: Number of non-manifold edges on surface + numConnectedComponents: Number of disconnected 3D mesh regions """ numVertices: int numEdges: int numFaces: int - eulerCharacteristic: int + numCells: int + solidEulerCharacteristic: int num3dCells: int num2dCells: int numOtherCells: int @@ -53,7 +55,7 @@ def __countConnectedComponents( mesh: vtk.vtkUnstructuredGrid ) -> int: Returns: Number of disconnected regions """ - setupLogger.info( "Checking for disconnected components..." ) + setupLogger.info( "Checking for disconnected 3D components..." ) connectivity = vtk.vtkConnectivityFilter() connectivity.SetInputData( mesh ) @@ -63,7 +65,7 @@ def __countConnectedComponents( mesh: vtk.vtkUnstructuredGrid ) -> int: numRegions = connectivity.GetNumberOfExtractedRegions() - setupLogger.info( f"Found {numRegions} disconnected component(s)" ) + setupLogger.info( f"Found {numRegions} disconnected 3D component(s)" ) return numRegions @@ -109,12 +111,12 @@ def __filter3dElements( mesh: vtk.vtkUnstructuredGrid ) -> tuple[ vtk.vtkUnstruc if not has3d: setupLogger.warning( "No 3D volumetric elements found!" ) setupLogger.warning( "This appears to be a pure surface mesh." ) - setupLogger.warning( "Computing Euler characteristic from all cells..." ) + setupLogger.warning( "Cannot compute solid Euler characteristic." ) return mesh, n3d, n2d, nOther, False if n2d > 0: setupLogger.info( f"Filtering out {n2d} 2D boundary cells..." ) - setupLogger.info( f"Using only {n3d} volumetric elements for surface extraction" ) + setupLogger.info( f"Using only {n3d} volumetric elements" ) # Extract only 3D cells using vtkExtractCells idList = vtk.vtkIdList() @@ -131,6 +133,91 @@ def __filter3dElements( mesh: vtk.vtkUnstructuredGrid ) -> tuple[ vtk.vtkUnstruc return filteredMesh, n3d, n2d, nOther, has3d +def __countUniqueEdgesAndFaces( mesh: vtk.vtkUnstructuredGrid ) -> tuple[ int, int ]: + """Count unique edges and faces in 3D mesh (NumPy optimized). + + Args: + mesh: 3D unstructured grid + + Returns: + Tuple of (num_edges, num_faces) + """ + setupLogger.info( "Counting unique edges and faces in 3D mesh..." ) + + try: + import numpy as np + use_numpy = True + except ImportError: + use_numpy = False + + numCells = mesh.GetNumberOfCells() + + if use_numpy: + # Use numpy for faster operations + edge_list = [] + face_list = [] + + for i in tqdm( range( numCells ), desc="Processing cells", mininterval=1.0 ): + + cell = mesh.GetCell( i ) + + # Edges + numEdges = cell.GetNumberOfEdges() + for edge_idx in range( numEdges ): + edge = cell.GetEdge( edge_idx ) + p0 = edge.GetPointId( 0 ) + p1 = edge.GetPointId( 1 ) + edge_list.append( ( min( p0, p1 ), max( p0, p1 ) ) ) + + # Faces + numFaces = cell.GetNumberOfFaces() + for face_idx in range( numFaces ): + face = cell.GetFace( face_idx ) + num_pts = face.GetNumberOfPoints() + point_ids = tuple( sorted( [ face.GetPointId( j ) for j in range( num_pts ) ] ) ) + face_list.append( point_ids ) + + # Use numpy unique for deduplication (faster than set) + setupLogger.info( " Deduplicating edges and faces..." ) + + # For edges: convert to array and use unique + edge_array = np.array( edge_list, dtype=np.int64 ) + unique_edges = np.unique( edge_array, axis=0 ) + num_edges = len( unique_edges ) + + # For faces: use set (numpy doesn't handle variable-length well) + num_faces = len( set( face_list ) ) + + else: + # Fallback to optimized set-based approach + edge_set = set() + face_set = set() + + for i in tqdm( range( numCells ), desc="Processing cells", mininterval=0.5 ): + cell = mesh.GetCell( i ) + + numEdges = cell.GetNumberOfEdges() + for edge_idx in range( numEdges ): + edge = cell.GetEdge( edge_idx ) + p0 = edge.GetPointId( 0 ) + p1 = edge.GetPointId( 1 ) + edge_set.add( ( min( p0, p1 ), max( p0, p1 ) ) ) + + numFaces = cell.GetNumberOfFaces() + for face_idx in range( numFaces ): + face = cell.GetFace( face_idx ) + num_pts = face.GetNumberOfPoints() + face_set.add( tuple( sorted( [ face.GetPointId( j ) for j in range( num_pts ) ] ) ) ) + + num_edges = len( edge_set ) + num_faces = len( face_set ) + + setupLogger.info( f" Unique edges: {num_edges:,}" ) + setupLogger.info( f" Unique faces: {num_faces:,}" ) + + return num_edges, num_faces + + def __extractSurface( mesh: vtk.vtkUnstructuredGrid ) -> vtk.vtkPolyData: """Extract surface from unstructured grid (3D elements only). @@ -140,32 +227,13 @@ def __extractSurface( mesh: vtk.vtkUnstructuredGrid ) -> vtk.vtkPolyData: Returns: Surface as polydata """ - setupLogger.info( "Extracting surface from 3D elements..." ) + setupLogger.info( "Extracting boundary surface from 3D elements..." ) surfaceFilter = vtk.vtkDataSetSurfaceFilter() surfaceFilter.SetInputData( mesh ) surfaceFilter.Update() return surfaceFilter.GetOutput() -def __countUniqueEdges( surface: vtk.vtkPolyData ) -> int: - """Count unique edges in the surface mesh. - - Args: - surface: Surface mesh - - Returns: - Number of unique edges - """ - setupLogger.info( "Counting unique edges..." ) - - # Use VTK's edge extraction (faster than manual iteration) - edgeExtractor = vtk.vtkExtractEdges() - edgeExtractor.SetInputData( surface ) - edgeExtractor.Update() - - return edgeExtractor.GetOutput().GetNumberOfCells() - - def __checkMeshQuality( surface: vtk.vtkPolyData ) -> tuple[ int, int ]: """Check for boundary edges and non-manifold features. @@ -201,106 +269,97 @@ def __checkMeshQuality( surface: vtk.vtkPolyData ) -> tuple[ int, int ]: def meshAction( mesh: vtk.vtkUnstructuredGrid, options: Options ) -> Result: - """Compute Euler characteristic for a mesh. + """Compute solid Euler characteristic for a mesh. - Only considers 3D volumetric elements. Extracts surface and computes V - E + F. + Only considers 3D volumetric elements. Computes chi_solid = V - E + F - C. Args: mesh: Input unstructured grid options: Computation options Returns: - Result with Euler characteristic and topology information + Result with solid Euler characteristic and topology information """ - setupLogger.info( "Starting Euler characteristic computation..." ) + setupLogger.info( "Starting solid Euler characteristic computation..." ) setupLogger.info( f"Input mesh: {mesh.GetNumberOfPoints()} points, {mesh.GetNumberOfCells()} cells" ) # Filter to 3D elements only mesh3d, n3d, n2d, nOther, has3d = __filter3dElements( mesh ) - # Count connected components BEFORE extracting surface + if not has3d: + raise RuntimeError( "Cannot compute solid Euler - no 3D cells found" ) + + # Count connected components numComponents = __countConnectedComponents( mesh3d ) - # Extract surface from 3D elements - surface = __extractSurface( mesh3d ) + # Get basic counts + V = mesh3d.GetNumberOfPoints() + C = mesh3d.GetNumberOfCells() - # Get counts - V = surface.GetNumberOfPoints() - F = surface.GetNumberOfCells() - E = __countUniqueEdges( surface ) + # Count unique edges and faces in 3D mesh + E, F = __countUniqueEdgesAndFaces( mesh3d ) - setupLogger.info( "Surface topology:" ) + setupLogger.info( "Solid mesh topology:" ) setupLogger.info( f" Vertices (V): {V:,}" ) setupLogger.info( f" Edges (E): {E:,}" ) setupLogger.info( f" Faces (F): {F:,}" ) + setupLogger.info( f" Cells (C): {C:,}" ) + + # Calculate solid Euler characteristic + solidEuler = V - E + F - C - # Calculate Euler characteristic - euler = V - E + F + setupLogger.info( f"Solid Euler characteristic (chi = V - E + F - C): {solidEuler}" ) - setupLogger.info( f"Euler characteristic (X = V - E + F): {euler}" ) - setupLogger.info( f"Connected components: {numComponents}" ) + # Interpret result + setupLogger.info( "Topology interpretation:" ) + setupLogger.info( f" 3D connected components: {numComponents}" ) - # Interpret result with component information if numComponents == 1: - # Single component - standard topology interpretation - if euler == 2: - setupLogger.info( "Topology: Single closed surface (sphere-like) - VALID for simulation!" ) - elif euler == 0: - setupLogger.warning( "Topology: Torus-like (genus 1)" ) - setupLogger.warning( " Expected X = 2 for standard simulation mesh" ) - elif euler == 1: - setupLogger.warning( "Topology: Disk-like (open surface)" ) - setupLogger.warning( " Mesh is not closed! Expected X = 2" ) - elif euler < 0: - genus = ( 2 - euler ) // 2 - setupLogger.warning( f"Topology: Complex surface with handles (genus g = {genus})" ) - setupLogger.warning( " Expected X = 2 for standard simulation mesh" ) + if solidEuler == 1: + setupLogger.info( " chi = 1: Contractible (solid ball topology)" ) + setupLogger.info( " VALID simple 3D region for simulation" ) + elif solidEuler == 0: + setupLogger.warning( " chi = 0: Solid torus (has through-hole)" ) + setupLogger.warning( " Verify this matches your domain geometry" ) + elif solidEuler == 2: + setupLogger.warning( " chi = 2: Hollow shell or internal cavity topology" ) + setupLogger.warning( " Expected chi = 1 for simple solid ball" ) + setupLogger.warning( " This suggests internal void or nested structure" ) + setupLogger.warning( " 3D cells ARE connected (verified above) - verify intended" ) else: - setupLogger.warning( f"Topology: Unexpected (X = {euler})" ) + setupLogger.warning( f" chi = {solidEuler}: Unusual topology" ) + setupLogger.warning( " Verify mesh integrity" ) else: - # Multiple components detected - expectedEuler = 2 * numComponents # If all components are sphere-like - - setupLogger.error( f"Topology: Mesh has {numComponents} DISCONNECTED COMPONENTS!" ) - setupLogger.error( f" Euler characteristic: {euler}" ) - setupLogger.error( f" Expected for {numComponents} spheres: X = {expectedEuler}" ) - - if numComponents > 1: - # Internal cavities = additional components beyond the outer surface - numCavities = numComponents - 1 - setupLogger.error( f" Likely {numCavities} internal void(s)/cavity(ies)" ) - - if euler == expectedEuler: - setupLogger.error( " All components are sphere-like (X=2 each)" ) - elif euler < expectedEuler: - setupLogger.error( " Some components have genus > 0 (holes/handles)" ) - else: - setupLogger.error( " Unusual topology - verify mesh integrity" ) - - setupLogger.error( " This mesh is NOT suitable for simulation without fixing!" ) - setupLogger.error( " Expected: single closed volume (X = 2, components = 1)" ) + setupLogger.error( f" Mesh has {numComponents} disconnected 3D components!" ) + setupLogger.error( " This is NOT suitable for simulation without fixing" ) # Check mesh quality + surface = __extractSurface( mesh3d ) boundaryEdges, nonManifoldEdges = __checkMeshQuality( surface ) setupLogger.info( "Mesh quality:" ) setupLogger.info( f" Boundary edges: {boundaryEdges:,}" ) setupLogger.info( f" Non-manifold edges: {nonManifoldEdges:,}" ) - if boundaryEdges == 0 and nonManifoldEdges == 0: - if euler == 2 and numComponents == 1: - setupLogger.info( " Perfect closed manifold mesh - READY for simulation!" ) + # Final validation + if numComponents == 1 and boundaryEdges == 0 and nonManifoldEdges == 0: + if solidEuler == 1: + setupLogger.info( " Perfect: single connected volume, simple topology - READY!" ) else: - setupLogger.warning( " Perfect manifold but topology issues - verify before simulation" ) + setupLogger.warning( f" Connected volume but chi = {solidEuler}" ) + setupLogger.warning( " Verify internal features are intended" ) + elif numComponents > 1: + setupLogger.error( " Multiple disconnected regions - INVALID!" ) elif boundaryEdges > 0: - setupLogger.warning( " Open surface detected (has boundaries)" ) - if nonManifoldEdges > 0: - setupLogger.error( " Non-manifold edges detected (mesh has issues!)" ) + setupLogger.error( " Open surface detected - INVALID!" ) + elif nonManifoldEdges > 0: + setupLogger.error( " Non-manifold geometry detected - INVALID!" ) return Result( numVertices=V, numEdges=E, numFaces=F, - eulerCharacteristic=euler, + numCells=C, + solidEulerCharacteristic=solidEuler, num3dCells=n3d, num2dCells=n2d, numOtherCells=nOther, @@ -310,14 +369,14 @@ def meshAction( mesh: vtk.vtkUnstructuredGrid, options: Options ) -> Result: def action( vtuInputFile: str, options: Options ) -> Result: - """Compute Euler characteristic for a VTU file. + """Compute solid Euler characteristic for a VTU file. Args: vtuInputFile: Path to input VTU file options: Computation options Returns: - Result with Euler characteristic and topology information + Result with solid Euler characteristic and topology information """ mesh = readUnstructuredGrid( vtuInputFile ) return meshAction( mesh, options ) diff --git a/mesh-doctor/src/geos/mesh_doctor/parsing/eulerParsing.py b/mesh-doctor/src/geos/mesh_doctor/parsing/eulerParsing.py index d061ba54..a55a298c 100644 --- a/mesh-doctor/src/geos/mesh_doctor/parsing/eulerParsing.py +++ b/mesh-doctor/src/geos/mesh_doctor/parsing/eulerParsing.py @@ -1,6 +1,6 @@ # SPDX-License-Identifier: Apache-2.0 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. -"""Command line parsing for Euler characteristic computation.""" +"""Command line parsing for solid Euler characteristic computation.""" from __future__ import annotations import argparse @@ -24,119 +24,180 @@ def convert( parsedOptions: dict[ str, Any ] ) -> Options: def fillSubparser( subparsers: argparse._SubParsersAction[ Any ] ) -> None: - """Fill the argument parser for the Euler characteristic action. + """Fill the argument parser for the solid Euler characteristic action. Args: subparsers: Subparsers from the main argument parser """ - p = subparsers.add_parser( EULER, - help="Compute Euler characteristic (X = V - E + F) for mesh topology analysis.", - description="""\ -Computes the Euler characteristic for a mesh by: + p = subparsers.add_parser( + EULER, + help="Compute solid Euler characteristic (chi = V - E + F - C) for 3D mesh topology validation.", + description="""\ +Computes the solid Euler characteristic for a 3D mesh by: 1. Filtering to 3D volumetric elements only (ignores 2D boundary cells) - 2. Extracting the surface from these 3D elements - 3. Counting vertices (V), edges (E), and faces (F) - 4. Computing X = V - E + F - -Expected Euler values: - X = 2 Closed surface (sphere-like, e.g., cube, tetrahedron) - X = 0 Torus or cylinder topology - X = 1 Disk or open surface - Other Complex topology with genus g = (2 - X)/2 - -This tool helps verify that meshes form proper closed volumes for simulation. + 2. Counting vertices (V), edges (E), faces (F), and cells (C) in the 3D mesh + 3. Computing chi = V - E + F - C + +Solid Euler characteristic values: + chi = 1 Simple solid region (contractible, ball-like) - standard + chi = 0 Solid torus (has through-hole) + chi = 2 Complex internal topology (may indicate internal features) + chi > 2 Multiple components or complex structure + +Primary validation criteria (what matters for simulation): + - 3D connectivity = 1 region (cells properly connected via shared faces) + - Boundary edges = 0 (closed manifold surface) + - Non-manifold edges = 0 (no overlapping cells) + +The Euler characteristic provides topological information but is secondary. +3D connectivity and manifold properties are the definitive validation checks. """ ) addVtuInputFileArgument( p ) def displayResults( options: Options, result: Result ) -> None: - """Display the results of the Euler characteristic computation. + """Display the results of the solid Euler characteristic computation. Args: options: The options used for the computation. result: The result of the computation. """ setupLogger.results( "=" * 80 ) - setupLogger.results( "EULER CHARACTERISTIC RESULTS" ) + setupLogger.results( "SOLID EULER CHARACTERISTIC RESULTS" ) setupLogger.results( "=" * 80 ) - setupLogger.results( f"Mesh Topology (from {result.num3dCells:,} 3D elements):" ) - setupLogger.results( f" Surface Vertices (V): {result.numVertices:,}" ) - setupLogger.results( f" Surface Edges (E): {result.numEdges:,}" ) - setupLogger.results( f" Surface Faces (F): {result.numFaces:,}" ) + setupLogger.results( f"3D Mesh Topology (from {result.num3dCells:,} 3D elements):" ) + setupLogger.results( f" Vertices (V): {result.numVertices:,}" ) + setupLogger.results( f" Edges (E): {result.numEdges:,}" ) + setupLogger.results( f" Faces (F): {result.numFaces:,}" ) + setupLogger.results( f" Cells (C): {result.numCells:,}" ) setupLogger.results( "-" * 80 ) - setupLogger.results( f" Euler Characteristic (X = V - E + F): {result.eulerCharacteristic}" ) - setupLogger.results( f" Connected Components: {result.numConnectedComponents}" ) + setupLogger.results( f" Euler Characteristic (chi = V - E + F - C): {result.solidEulerCharacteristic}" ) + setupLogger.results( f" 3D Connected Components: {result.numConnectedComponents}" ) setupLogger.results( "=" * 80 ) - # Interpret topology with validation + # PRIMARY VALIDATION: 3D connectivity + setupLogger.results( "" ) + setupLogger.results( "3D CONNECTIVITY CHECK:" ) + setupLogger.results( "-" * 80 ) + if result.numConnectedComponents == 1: - # Single component - standard interpretation - if result.eulerCharacteristic == 2: - setupLogger.results( "VALIDATION PASSED: Single closed surface (sphere-like)" ) - setupLogger.results( " Status: MESH IS VALID FOR SIMULATION" ) - elif result.eulerCharacteristic == 0: - setupLogger.results( "WARNING: Torus-like topology (genus 1)" ) - setupLogger.results( " Expected X = 2 for standard simulation mesh" ) - elif result.eulerCharacteristic == 1: - setupLogger.results( "WARNING: Open surface (not closed)" ) - setupLogger.results( " Mesh has boundary - expected X = 2" ) - elif result.eulerCharacteristic < 0: - genus = ( 2 - result.eulerCharacteristic ) // 2 - setupLogger.results( f"WARNING: Complex topology with handles (genus g = {genus})" ) - setupLogger.results( f" Surface has {genus} hole(s)/handle(s)" ) - setupLogger.results( " Expected X = 2 for standard simulation mesh" ) - else: - setupLogger.results( f"WARNING: Unexpected topology (X = {result.eulerCharacteristic})" ) + setupLogger.results( "PASS: Single connected 3D region" ) + setupLogger.results( " All 3D cells form one continuous volume via shared faces" ) else: - # Multiple components - indicates disconnected regions or internal cavities - expectedEuler = 2 * result.numConnectedComponents - numCavities = result.numConnectedComponents - 1 - - setupLogger.results( f"VALIDATION FAILED: {result.numConnectedComponents} DISCONNECTED COMPONENTS DETECTED!" ) - setupLogger.results( f" Euler characteristic: {result.eulerCharacteristic}" ) - setupLogger.results( f" Expected (if all spherical): X = {expectedEuler}" ) - setupLogger.results( f" Likely internal void(s)/cavity(ies): {numCavities}" ) - setupLogger.results( "" ) + setupLogger.results( f"FAIL: {result.numConnectedComponents} disconnected 3D regions!" ) + setupLogger.results( " Mesh has separate volumes not connected by shared faces" ) + setupLogger.results( " NOT suitable for simulation - requires fixing" ) - if result.eulerCharacteristic == expectedEuler: - setupLogger.results( " All components are sphere-like (X=2 each)" ) - setupLogger.results( " Issue: Mesh has internal hollow regions (outer + inner surfaces)" ) - elif result.eulerCharacteristic < expectedEuler: - avgGenus = ( expectedEuler - result.eulerCharacteristic ) / 2 - setupLogger.results( f" Some components have handles/holes (avg genus ~{avgGenus:.1f})" ) - else: - setupLogger.results( " Unusual topology - verify mesh integrity" ) - - setupLogger.results( "" ) - setupLogger.results( " NOT SUITABLE FOR SIMULATION without fixing!" ) - setupLogger.results( " Expected: Single closed volume (X = 2, components = 1)" ) + # TOPOLOGY INTERPRETATION + setupLogger.results( "" ) + setupLogger.results( "TOPOLOGY INTERPRETATION:" ) + setupLogger.results( "-" * 80 ) - # Show input cell breakdown if there were 2D cells filtered out - if result.num2dCells > 0: - setupLogger.results( "" ) - setupLogger.results( "Input cell breakdown:" ) - setupLogger.results( f" 3D volumetric cells: {result.num3dCells:,} (used for surface extraction)" ) - setupLogger.results( f" 2D surface cells: {result.num2dCells:,} (filtered out)" ) - if result.numOtherCells > 0: - setupLogger.results( f" Other cells: {result.numOtherCells:,} (filtered out)" ) + if result.numConnectedComponents == 1: + # Single connected region - interpret Euler + if result.solidEulerCharacteristic == 1: + setupLogger.results( "chi = 1: Simple solid ball (contractible)" ) + setupLogger.results( " Standard topology for simulation meshes" ) + + elif result.solidEulerCharacteristic == 0: + setupLogger.results( "chi = 0: Solid torus topology" ) + setupLogger.results( " Domain has through-hole or tunnel" ) + setupLogger.results( " Verify this matches your expected geometry" ) + + elif result.solidEulerCharacteristic == 2: + setupLogger.results( "chi = 2: Hollow shell or internal cavity" ) + setupLogger.results( " Expected chi = 1 for simple solid ball" ) + setupLogger.results( " Suggests internal void or nested structure" ) + setupLogger.results( " 3D cells ARE connected (verified above)" ) + setupLogger.results( " ACCEPTABLE if internal structure is intentional" ) - # Quality check results + else: + setupLogger.results( f"chi = {result.solidEulerCharacteristic}: Unusual value" ) + if result.solidEulerCharacteristic > 2: + setupLogger.results( " May indicate complex internal structure" ) + setupLogger.results( " Verify mesh structure and topology are correct" ) + else: + # Multiple disconnected regions + expectedEuler = result.numConnectedComponents # Each component contributes + setupLogger.results( "Multiple regions: topology analysis not meaningful" ) + setupLogger.results( + f"Expected chi approx {expectedEuler} for {result.numConnectedComponents} disconnected balls" ) + setupLogger.results( f"Actual chi = {result.solidEulerCharacteristic}" ) + + # MESH QUALITY setupLogger.results( "" ) - setupLogger.results( "Mesh Quality Check:" ) + setupLogger.results( "MESH QUALITY:" ) setupLogger.results( "-" * 80 ) setupLogger.results( f" Boundary edges: {result.numBoundaryEdges:,}" ) setupLogger.results( f" Non-manifold edges: {result.numNonManifoldEdges:,}" ) - if result.numBoundaryEdges == 0 and result.numNonManifoldEdges == 0: - if result.eulerCharacteristic == 2 and result.numConnectedComponents == 1: - setupLogger.results( " Status: PERFECT - Ready for simulation!" ) + if result.numBoundaryEdges == 0: + setupLogger.results( " Closed manifold (no open boundaries)" ) + else: + setupLogger.results( " Open boundaries detected" ) + + if result.numNonManifoldEdges == 0: + setupLogger.results( " Manifold geometry (no overlapping cells)" ) + else: + setupLogger.results( " Non-manifold edges (overlapping/duplicate cells)" ) + + # FINAL VALIDATION + setupLogger.results( "" ) + setupLogger.results( "FINAL VALIDATION:" ) + setupLogger.results( "=" * 80 ) + + # Check all three critical criteria + has_single_component = result.numConnectedComponents == 1 + is_closed = result.numBoundaryEdges == 0 + is_manifold = result.numNonManifoldEdges == 0 + + if has_single_component and is_closed and is_manifold: + # All critical checks pass + if result.solidEulerCharacteristic == 1: + setupLogger.results( "STATUS: PERFECT" ) + setupLogger.results( " Single connected region: YES" ) + setupLogger.results( " Closed manifold: YES" ) + setupLogger.results( " Simple topology (chi = 1): YES" ) + setupLogger.results( " READY for simulation" ) else: - setupLogger.results( " Status: Manifold but topology issues (X != 2 or internal cavities)" ) + setupLogger.results( "STATUS: VALID (with complex topology)" ) + setupLogger.results( " Single connected region: YES" ) + setupLogger.results( " Closed manifold: YES" ) + setupLogger.results( f" Complex topology (chi = {result.solidEulerCharacteristic}): WARNING" ) + setupLogger.results( " ACCEPTABLE for simulation" ) + setupLogger.results( " Verify internal structure is intentional" ) else: - if result.numBoundaryEdges > 0: - setupLogger.results( " Status: Open surface detected (has boundaries)" ) - if result.numNonManifoldEdges > 0: - setupLogger.results( " Status: Non-manifold edges detected!" ) + # Failed critical checks + setupLogger.results( "STATUS: INVALID" ) + setupLogger.results( " Mesh has critical issues:" ) + + if not has_single_component: + setupLogger.results( f" Multiple disconnected regions ({result.numConnectedComponents}): FAIL" ) + else: + setupLogger.results( " Single connected region: PASS" ) + + if not is_closed: + setupLogger.results( f" Open boundaries ({result.numBoundaryEdges:,} boundary edges): FAIL" ) + else: + setupLogger.results( " Closed manifold: PASS" ) + + if not is_manifold: + setupLogger.results( f" Non-manifold geometry ({result.numNonManifoldEdges:,} edges): FAIL" ) + else: + setupLogger.results( " Manifold geometry: PASS" ) + + setupLogger.results( "" ) + setupLogger.results( " NOT SUITABLE for simulation without fixing" ) + + # Show input cell breakdown + if result.num2dCells > 0 or result.numOtherCells > 0: + setupLogger.results( "" ) + setupLogger.results( "Input cell breakdown:" ) + setupLogger.results( f" 3D volumetric cells: {result.num3dCells:,} (used for analysis)" ) + if result.num2dCells > 0: + setupLogger.results( f" 2D surface cells: {result.num2dCells:,} (filtered out)" ) + if result.numOtherCells > 0: + setupLogger.results( f" Other cells: {result.numOtherCells:,} (filtered out)" ) setupLogger.results( "=" * 80 )