diff --git a/.mypy.ini b/.mypy.ini index 50a8effb..66c805f0 100644 --- a/.mypy.ini +++ b/.mypy.ini @@ -1,5 +1,6 @@ [mypy] python_version = 3.10 +exclude = /(build|dist|\.eggs)/ # strict = true warn_unreachable = True diff --git a/mesh-doctor/tests/helpers.py b/mesh-doctor/tests/helpers.py new file mode 100644 index 00000000..85ef588c --- /dev/null +++ b/mesh-doctor/tests/helpers.py @@ -0,0 +1,644 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2026 TotalEnergies. +# SPDX-FileContributor: jafranc +"""Mesh-building helpers for mesh-doctor unit tests. + +Generators extracted from misc/mesh_gen.py (--gen branch) plus +synthetic defect configurations for testing each mesh-doctor check. + +Public API +---------- +Mesh generators (from gen branch): + build_hex_mesh clean z-major hex grid + build_tetra_mesh Delaunay3D tetra on z-major grid + build_tetra_mesh_patterned Delaunay3D on checkerboard-filtered points ("shifted") + build_hex_mesh_with_fractures hex + quad fractures, 'attribute' cell array + build_tetra_mesh_with_fractures tetra + tri fractures, 'attribute' cell array + +Non-manifold configurations: + build_surface_with_non_manifold_edge 3 quads sharing one edge (flap) + build_volume_with_non_manifold_edge 2 hexes sharing only one edge (bowtie) + +Mesh-doctor defect configurations: + build_mesh_with_collocated_nodes → collocatedNodes check + build_mesh_with_orphan_2d_cells → orphan2d check + build_mesh_with_negative_volume → elementVolumes check + build_mesh_with_self_intersecting_elements → selfIntersectingElements check +""" + +import numpy as np +import numpy.typing as npt +import vtk +from vtkmodules.vtkCommonCore import vtkIdList, vtkIntArray, vtkPoints +from vtkmodules.vtkCommonDataModel import ( vtkPolyData, vtkUnstructuredGrid, VTK_HEXAHEDRON, VTK_QUAD ) +from vtkmodules.vtkFiltersCore import vtkAppendFilter +from vtkmodules.util.numpy_support import numpy_to_vtk + +# ────────────────────────────────────────────────────────────── +# Private geometry primitives +# ────────────────────────────────────────────────────────────── + + +def _gen_box_z_major( Lx: float, Ly: float, Lz: float, nx: int, ny: int, nz: int ) -> npt.NDArray[ np.float64 ]: + """Grid points with x slowest, y medium, z fastest (z-major ordering). + + Point at (i, j, k) has linear index k + j*nz + i*nz*ny. + """ + xs = np.linspace( 0.0, Lx, nx ) + ys = np.linspace( 0.0, Ly, ny ) + zs = np.linspace( 0.0, Lz, nz ) + return np.array( [ [ xs[ i ], ys[ j ], zs[ k ] ] for i in range( nx ) for j in range( ny ) for k in range( nz ) ], + dtype=float ) + + +def _gen_rect_yz_z_major( Ly: float, Lz: float, ny: int, nz: int, fixx: float, offy: int, + offz: int ) -> npt.NDArray[ np.float64 ]: + """Rectangular grid in the Y–Z plane at x = fixx; z fastest, then y. + + Trims offy rows on each y-side and offz rows on each z-side. + Row stride in the returned array is nz - 2*offz. + """ + ys = np.linspace( 0.0, Ly, ny ) + zs = np.linspace( 0.0, Lz, nz ) + return np.array( [ [ fixx, ys[ j ], zs[ k ] ] for j in range( offy, ny - offy ) for k in range( offz, nz - offz ) ], + dtype=float ) + + +def _patterned_indices_z_major( nx: int, ny: int, nz: int ) -> npt.NDArray[ np.intp ]: + """Checkerboard-patterned index selection for the z-major point grid. + + Returns every other point per row with alternating start (parity), + producing a non-uniform "shifted" distribution suitable for + an irregular Delaunay tetra mesh. + """ + blocks = [] + for row in range( nx * ny ): + x = row // ny + y = row % ny + parity = ( x % 2 ) ^ ( y % 2 ) + start = row * nz + parity + blocks.append( np.arange( start, ( row + 1 ) * nz, 2 ) ) + return np.concatenate( blocks ) + + +def _paint_and_append( main: vtkUnstructuredGrid, fracs: list[ vtkUnstructuredGrid ] ) -> vtkUnstructuredGrid: + """Combine volume mesh and fracture meshes, attaching an 'attribute' array. + + The 'attribute' value is 1 for matrix cells and 2, 3, … for each fracture + in the order they appear in *fracs*. This is the layout expected by + mesh-doctor's generateFractures action. + """ + append = vtkAppendFilter() + append.MergePointsOff() + append.AddInputData( main ) + for frac in fracs: + append.AddInputData( frac ) + append.Update() + output: vtkUnstructuredGrid = append.GetOutput() + + n_main = main.GetNumberOfCells() + n_total = output.GetNumberOfCells() + + attr = vtkIntArray() + attr.SetName( "attribute" ) + attr.SetNumberOfComponents( 1 ) + attr.SetNumberOfTuples( n_total ) + + for c in range( n_main ): + attr.SetValue( c, 1 ) + offset = n_main + for fi, frac in enumerate( fracs ): + for c in range( frac.GetNumberOfCells() ): + attr.SetValue( offset + c, fi + 2 ) + offset += frac.GetNumberOfCells() + + output.GetCellData().AddArray( attr ) + return output + + +# ────────────────────────────────────────────────────────────── +# Public mesh generators (from mesh_gen.py --gen branch) +# ────────────────────────────────────────────────────────────── + + +def build_hex_mesh( nx: int = 5, + ny: int = 5, + nz: int = 5, + Lx: float = 1.0, + Ly: float = 1.0, + Lz: float = 1.0 ) -> vtkUnstructuredGrid: + """Regular hexahedral mesh with z-major point ordering (the "shifted" variant). + + Points are stored with x slowest, z fastest. VTK hex nodes follow the + standard ordering: bottom face (k, k+1) then top face (i+1 slab), matching + the z-major CELL_ORDERING from mesh_gen.py. + + Args: + nx: Number of nodes in x (cells = nx-1). + ny: Number of nodes in y (cells = ny-1). + nz: Number of nodes in z (cells = nz-1). + Lx: Domain length in x. + Ly: Domain length in y. + Lz: Domain length in z. + + Returns: + vtkUnstructuredGrid with (nx-1)*(ny-1)*(nz-1) VTK_HEXAHEDRON cells. + """ + pts = _gen_box_z_major( Lx, Ly, Lz, nx, ny, nz ) + + mesh = vtkUnstructuredGrid() + mesh.Allocate( ( nx - 1 ) * ( ny - 1 ) * ( nz - 1 ) ) + vpts = vtkPoints() + vpts.SetData( numpy_to_vtk( pts ) ) + mesh.SetPoints( vpts ) + + for i in range( nx - 1 ): + for j in range( ny - 1 ): + for k in range( nz - 1 ): + # z-major base index: k + j*nz + i*nz*ny + base = k + j * nz + i * nz * ny + ids = vtkIdList() + ids.InsertNextId( base ) # VTK node 0: (i , j , k ) + ids.InsertNextId( base + nz * ny ) # VTK node 1: (i+1, j , k ) + ids.InsertNextId( base + nz * ny + nz ) # VTK node 2: (i+1, j+1, k ) + ids.InsertNextId( base + nz ) # VTK node 3: (i , j+1, k ) + ids.InsertNextId( base + 1 ) # VTK node 4: (i , j , k+1) + ids.InsertNextId( base + nz * ny + 1 ) # VTK node 5: (i+1, j , k+1) + ids.InsertNextId( base + nz * ny + nz + 1 ) # VTK node 6: (i+1, j+1, k+1) + ids.InsertNextId( base + nz + 1 ) # VTK node 7: (i , j+1, k+1) + mesh.InsertNextCell( VTK_HEXAHEDRON, ids ) + + return mesh + + +def build_tetra_mesh( nx: int = 5, + ny: int = 5, + nz: int = 5, + Lx: float = 1.0, + Ly: float = 1.0, + Lz: float = 1.0 ) -> vtkUnstructuredGrid: + """Tetrahedral mesh built via Delaunay3D on a full z-major grid. + + Uses all (nx*ny*nz) grid points as Delaunay input, yielding a structured + but tetrahedral discretisation. Equivalent to mesh_gen.py --tri --structured. + + Args: + nx: Node count along x. + ny: Node count along y. + nz: Node count along z. + Lx: Domain extent along x. + Ly: Domain extent along y. + Lz: Domain extent along z. + + Returns: + vtkUnstructuredGrid of VTK_TETRA cells. + """ + pts = _gen_box_z_major( Lx, Ly, Lz, nx, ny, nz ) + cloud = vtkPolyData() + vpts = vtkPoints() + vpts.SetData( numpy_to_vtk( pts ) ) + cloud.SetPoints( vpts ) + + delaunay = vtk.vtkDelaunay3D() + delaunay.SetInputData( cloud ) + delaunay.Update() + return delaunay.GetOutput() + + +def build_tetra_mesh_patterned( nx: int = 5, + ny: int = 5, + nz: int = 5, + Lx: float = 1.0, + Ly: float = 1.0, + Lz: float = 1.0 ) -> vtkUnstructuredGrid: + """Irregular ("shifted") tetrahedral mesh via Delaunay3D on checkerboard-filtered points. + + Every other grid point is dropped in a checkerboard pattern before + Delaunay triangulation, producing a more irregular connectivity that + exercises more code paths in mesh-doctor checks. Equivalent to + mesh_gen.py --tri (unstructured variant, odd nx forced if even). + + Args: + nx: Node count in x; incremented to odd if even (required by patterning). + ny: Node count along y. + nz: Node count along z. + Lx: Domain extent along x. + Ly: Domain extent along y. + Lz: Domain extent along z. + + Returns: + vtkUnstructuredGrid of VTK_TETRA cells. + """ + if nx % 2 == 0: + nx += 1 + + pts = _gen_box_z_major( Lx, Ly, Lz, nx, ny, nz ) + pts = pts[ _patterned_indices_z_major( nx, ny, nz ), : ] + + cloud = vtkPolyData() + vpts = vtkPoints() + vpts.SetData( numpy_to_vtk( pts ) ) + cloud.SetPoints( vpts ) + + delaunay = vtk.vtkDelaunay3D() + delaunay.SetInputData( cloud ) + delaunay.Update() + return delaunay.GetOutput() + + +def build_hex_mesh_with_fractures( nx: int = 5, + ny: int = 5, + nz: int = 5, + Lx: float = 1.0, + Ly: float = 1.0, + Lz: float = 1.0, + offx: int = 1, + offy: int = 1, + nfrac: int = 1 ) -> vtkUnstructuredGrid: + """Hex mesh with embedded quad fracture surfaces. + + Fractures are Y–Z planes located at x = Lx / (2^nfrac) * (i+1) for + i = 0 … nfrac-1. Each fracture surface is trimmed by *offx* nodes on + each z-side and *offy* nodes on each y-side so it doesn't touch the domain + boundary. + + The returned mesh has an 'attribute' cell data array: + 1 → matrix cells + 2, 3, … → fracture cells (in insertion order). + + This layout is consumed by mesh-doctor's generateFractures action. + + Args: + nx: Node count along x. + ny: Node count along y. + nz: Node count along z. + Lx: Domain extent along x. + Ly: Domain extent along y. + Lz: Domain extent along z. + offx: Trim margin in z (nodes to skip on each z-side of fracture). + offy: Trim margin in y (nodes to skip on each y-side of fracture). + nfrac: Number of fracture planes to embed. + + Returns: + vtkUnstructuredGrid with VTK_HEXAHEDRON + VTK_QUAD cells and 'attribute'. + """ + main = build_hex_mesh( nx, ny, nz, Lx, Ly, Lz ) + fracs: list[ vtkUnstructuredGrid ] = [] + stride = nz - 2 * offx # points per row on the fracture surface (z-direction) + + for fi in range( nfrac ): + x_pos = Lx / ( 2**nfrac ) * ( fi + 1 ) + pts = _gen_rect_yz_z_major( Ly, Lz, ny, nz, x_pos, offy, offx ) + + frac = vtkUnstructuredGrid() + vpts = vtkPoints() + vpts.SetData( numpy_to_vtk( pts ) ) + frac.SetPoints( vpts ) + frac.Allocate() + + for j in range( ny - 1 - 2 * offy ): + for k in range( nz - 1 - 2 * offx ): + base = k + j * stride + ids = vtkIdList() + ids.InsertNextId( base ) + ids.InsertNextId( base + 1 ) + ids.InsertNextId( base + stride + 1 ) + ids.InsertNextId( base + stride ) + frac.InsertNextCell( VTK_QUAD, ids ) + fracs.append( frac ) + + return _paint_and_append( main, fracs ) + + +def build_tetra_mesh_with_fractures( nx: int = 5, + ny: int = 5, + nz: int = 5, + Lx: float = 1.0, + Ly: float = 1.0, + Lz: float = 1.0, + offx: int = 1, + offy: int = 1, + nfrac: int = 1 ) -> vtkUnstructuredGrid: + """Tetra mesh with embedded triangulated fracture surfaces. + + Fractures are Y–Z planes triangulated via vtkDelaunay2D on a trimmed + rectangular point cloud. Equivalent to mesh_gen.py --tri --quad=False. + The returned mesh has an 'attribute' cell data array (same convention as + build_hex_mesh_with_fractures). + + Args: + nx: Node count along x. + ny: Node count along y. + nz: Node count along z. + Lx: Domain extent along x. + Ly: Domain extent along y. + Lz: Domain extent along z. + offx: Trim margin in z for fracture surfaces. + offy: Trim margin in y for fracture surfaces. + nfrac: Number of fracture planes. + + Returns: + vtkUnstructuredGrid with VTK_TETRA + VTK_TRIANGLE cells and 'attribute'. + """ + main = build_tetra_mesh( nx, ny, nz, Lx, Ly, Lz ) + fracs: list[ vtkUnstructuredGrid ] = [] + + for fi in range( nfrac ): + x_pos = Lx / ( 2**nfrac ) * ( fi + 1 ) + pts = _gen_rect_yz_z_major( Ly, Lz, ny, nz, x_pos, offy, offx ) + + cloud = vtkPolyData() + vpts = vtkPoints() + vpts.SetData( numpy_to_vtk( pts ) ) + cloud.SetPoints( vpts ) + + delaunay = vtk.vtkDelaunay2D() + delaunay.SetInputData( cloud ) + delaunay.SetProjectionPlaneMode( vtk.VTK_BEST_FITTING_PLANE ) + delaunay.Update() + tri_poly = delaunay.GetOutput() + + frac = vtkUnstructuredGrid() + frac.SetPoints( tri_poly.GetPoints() ) + frac.Allocate() + for c in range( tri_poly.GetNumberOfCells() ): + cell = tri_poly.GetCell( c ) + frac.InsertNextCell( cell.GetCellType(), cell.GetPointIds() ) + fracs.append( frac ) + + return _paint_and_append( main, fracs ) + + +# ────────────────────────────────────────────────────────────── +# Non-manifold edge configurations +# ────────────────────────────────────────────────────────────── + + +def build_surface_with_non_manifold_edge() -> vtkUnstructuredGrid: + """Surface mesh with one edge shared by three quad faces (non-manifold). + + Layout (edge 1–2 is the non-manifold edge):: + + Quad A z=0 plane x ∈ [0,1] + Quad B z=0 plane x ∈ [1,2] + Quad C x=1 plane perpendicular "flap" + + Points: + 0:(0,0,0) 1:(1,0,0) 2:(1,1,0) 3:(0,1,0) + 4:(2,0,0) 5:(2,1,0) + 6:(1,0,1) 7:(1,1,1) ← flap top edge + + Edge 1–2 (x=1, y∈[0,1], z=0) appears in A, B, and C → non-manifold. + + Detected by ``euler.meshAction`` as ``numNonManifoldEdges > 0`` + (applied to the surface mesh directly, not via volume extraction). + + Returns: + vtkUnstructuredGrid of 3 VTK_QUAD cells. + """ + coords = [ + ( 0.0, 0.0, 0.0 ), # 0 + ( 1.0, 0.0, 0.0 ), # 1 ← shared-edge endpoint + ( 1.0, 1.0, 0.0 ), # 2 ← shared-edge endpoint + ( 0.0, 1.0, 0.0 ), # 3 + ( 2.0, 0.0, 0.0 ), # 4 + ( 2.0, 1.0, 0.0 ), # 5 + ( 1.0, 0.0, 1.0 ), # 6 + ( 1.0, 1.0, 1.0 ), # 7 + ] + mesh = _make_ugrid( coords ) + for quad_pts in [ ( 0, 1, 2, 3 ), ( 1, 4, 5, 2 ), ( 1, 6, 7, 2 ) ]: + _add_cell( mesh, VTK_QUAD, quad_pts ) + return mesh + + +def build_volume_with_non_manifold_edge() -> vtkUnstructuredGrid: + """Two hexahedra sharing only one edge (bowtie / pinched configuration). + + Hex A occupies [0,1]³. Hex B occupies [1,2]×[1,2]×[0,1]. + They share exactly the edge from (1,1,0) to (1,1,1). + + On the extracted boundary surface that edge is surrounded by four faces + (two from each hex), making it a non-manifold boundary edge:: + + Boundary faces containing edge (1,1,0)→(1,1,1): + Hex A x-max face: (1,0,0)-(1,1,0)-(1,1,1)-(1,0,1) + Hex A y-max face: (0,1,0)-(1,1,0)-(1,1,1)-(0,1,1) + Hex B x-min face: (1,1,0)-(1,2,0)-(1,2,1)-(1,1,1) + Hex B y-min face: (1,1,0)-(2,1,0)-(2,1,1)-(1,1,1) + + Detected by ``euler.meshAction`` as ``numNonManifoldEdges > 0``. + + Returns: + vtkUnstructuredGrid of 2 VTK_HEXAHEDRON cells (14 points). + """ + coords = [ + # Hex A: [0,1]³ + ( 0.0, 0.0, 0.0 ), + ( 1.0, 0.0, 0.0 ), + ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), # 0-3 + ( 0.0, 0.0, 1.0 ), + ( 1.0, 0.0, 1.0 ), + ( 1.0, 1.0, 1.0 ), + ( 0.0, 1.0, 1.0 ), # 4-7 + # Hex B: [1,2]×[1,2]×[0,1] — points 2 and 6 from Hex A are reused + ( 2.0, 1.0, 0.0 ), + ( 2.0, 2.0, 0.0 ), + ( 1.0, 2.0, 0.0 ), # 8-10 + ( 2.0, 1.0, 1.0 ), + ( 2.0, 2.0, 1.0 ), + ( 1.0, 2.0, 1.0 ), # 11-13 + ] + mesh = _make_ugrid( coords ) + _add_cell( mesh, VTK_HEXAHEDRON, ( 0, 1, 2, 3, 4, 5, 6, 7 ) ) # Hex A + _add_cell( mesh, VTK_HEXAHEDRON, ( 2, 8, 9, 10, 6, 11, 12, 13 ) ) # Hex B (reuses 2,6) + return mesh + + +# ────────────────────────────────────────────────────────────── +# Mesh-doctor targeted defect configurations +# ────────────────────────────────────────────────────────────── + + +def build_mesh_with_collocated_nodes() -> vtkUnstructuredGrid: + """Two adjacent hexahedra whose shared-face nodes are duplicated in place. + + The mesh has 16 points instead of 12: the four nodes on the shared face + (x=1 plane) are represented twice at identical coordinates. + + Expected result from ``collocatedNodes.meshAction``:: + + len(result.nodesBuckets) == 4 + each bucket contains 2 node IDs + + Returns: + vtkUnstructuredGrid of 2 VTK_HEXAHEDRON cells with 16 points. + """ + coords = [ + # Hex A [0,1]³ + ( 0.0, 0.0, 0.0 ), + ( 1.0, 0.0, 0.0 ), + ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), # 0-3 + ( 0.0, 0.0, 1.0 ), + ( 1.0, 0.0, 1.0 ), + ( 1.0, 1.0, 1.0 ), + ( 0.0, 1.0, 1.0 ), # 4-7 + # Hex B [1,2]³ — points 8,11,12,15 are collocated duplicates of 1,2,5,6 + ( 1.0, 0.0, 0.0 ), + ( 2.0, 0.0, 0.0 ), + ( 2.0, 1.0, 0.0 ), + ( 1.0, 1.0, 0.0 ), # 8-11 + ( 1.0, 0.0, 1.0 ), + ( 2.0, 0.0, 1.0 ), + ( 2.0, 1.0, 1.0 ), + ( 1.0, 1.0, 1.0 ), # 12-15 + ] + mesh = _make_ugrid( coords ) + _add_cell( mesh, VTK_HEXAHEDRON, ( 0, 1, 2, 3, 4, 5, 6, 7 ) ) + _add_cell( mesh, VTK_HEXAHEDRON, ( 8, 9, 10, 11, 12, 13, 14, 15 ) ) + return mesh + + +def build_mesh_with_orphan_2d_cells() -> vtkUnstructuredGrid: + """One hex with one matching quad face and one orphan quad not matching any 3D face. + + Layout:: + + Hex: unit cube [0,1]³ (nodes 0-7) + Quad M: z=0 face of the hex (nodes 0,1,2,3) → matched + Quad O: z=2 plane (nodes 8,9,10,11) → orphan + + Expected result from ``orphan2d.meshAction``:: + + result.matched2dCells == 1 + result.orphaned2dCells == 1 + result.orphaned2dIndices == [2] (third cell, 0-indexed) + + Returns: + vtkUnstructuredGrid of 1 VTK_HEXAHEDRON + 2 VTK_QUAD cells. + """ + coords = [ + ( 0.0, 0.0, 0.0 ), + ( 1.0, 0.0, 0.0 ), + ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), # 0-3 hex base + ( 0.0, 0.0, 1.0 ), + ( 1.0, 0.0, 1.0 ), + ( 1.0, 1.0, 1.0 ), + ( 0.0, 1.0, 1.0 ), # 4-7 hex top + ( 0.0, 0.0, 2.0 ), + ( 1.0, 0.0, 2.0 ), + ( 1.0, 1.0, 2.0 ), + ( 0.0, 1.0, 2.0 ), # 8-11 orphan + ] + mesh = _make_ugrid( coords ) + _add_cell( mesh, VTK_HEXAHEDRON, ( 0, 1, 2, 3, 4, 5, 6, 7 ) ) # hex + _add_cell( mesh, VTK_QUAD, ( 0, 1, 2, 3 ) ) # matched (face [0,1,2,3]) + _add_cell( mesh, VTK_QUAD, ( 8, 9, 10, 11 ) ) # orphan + _add_cell( mesh, VTK_QUAD, ( 0, 3, 6, 5 ) ) # orphan + return mesh + + +def build_mesh_with_negative_volume() -> vtkUnstructuredGrid: + """Two hexahedra: one correctly oriented, one with top/bottom faces swapped (inverted). + + Swapping the bottom and top node groups of a VTK hex reverses the + Jacobian sign, yielding a negative volume. + + Expected result from ``elementVolumes.meshAction(mesh, Options(minVolume=0))``:: + + len(result.elementVolumes) == 1 + result.elementVolumes[0][0] == 1 (second cell) + + Returns: + vtkUnstructuredGrid of 2 VTK_HEXAHEDRON cells. + """ + coords = [ + # Hex A [0,1]³ correct ordering + ( 0.0, 0.0, 0.0 ), + ( 1.0, 0.0, 0.0 ), + ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), # 0-3 + ( 0.0, 0.0, 1.0 ), + ( 1.0, 0.0, 1.0 ), + ( 1.0, 1.0, 1.0 ), + ( 0.0, 1.0, 1.0 ), # 4-7 + # Hex B [2,3]³ inverted: bottom/top groups swapped + ( 2.0, 0.0, 0.0 ), + ( 3.0, 0.0, 0.0 ), + ( 3.0, 1.0, 0.0 ), + ( 2.0, 1.0, 0.0 ), # 8-11 ← physical bottom + ( 2.0, 0.0, 1.0 ), + ( 3.0, 0.0, 1.0 ), + ( 3.0, 1.0, 1.0 ), + ( 2.0, 1.0, 1.0 ), # 12-15 ← physical top + ] + mesh = _make_ugrid( coords ) + _add_cell( mesh, VTK_HEXAHEDRON, ( 0, 1, 2, 3, 4, 5, 6, 7 ) ) # Hex A correct + _add_cell( mesh, VTK_HEXAHEDRON, ( 12, 13, 14, 15, 8, 9, 10, 11 ) ) # Hex B top→bottom first + return mesh + + +def build_mesh_with_self_intersecting_elements() -> vtkUnstructuredGrid: + """Two hexahedra that geometrically overlap in the x ∈ [0.5, 1.0] band. + + Hex A spans [0,1]³; Hex B spans [0.5,1.5]×[0,1]×[0,1]. They share no + mesh nodes but their geometric extents intersect. + + Expected result from ``selfIntersectingElements.meshAction``:: + + The pair (0, 1) is reported as intersecting. + + Returns: + vtkUnstructuredGrid of 2 VTK_HEXAHEDRON cells (16 points, no shared nodes). + """ + coords = [ + # Hex A [0,1]³ + ( 0.0, 0.0, 0.0 ), + ( 1.0, 0.0, 0.0 ), + ( 1.0, 1.0, 0.0 ), + ( 0.0, 1.0, 0.0 ), # 0-3 + ( 0.0, 0.0, 1.0 ), + ( 1.0, 0.0, 1.0 ), + ( 1.0, 1.0, 1.0 ), + ( 0.0, 1.0, 1.0 ), # 4-7 + # Hex B [0.5,1.5]×[0,1]×[0,1] + ( 0.5, 0.0, 0.0 ), + ( 1.5, 0.0, 0.0 ), + ( 1.5, 1.0, 0.0 ), + ( 0.5, 1.0, 0.0 ), # 8-11 + ( 0.5, 0.0, 1.0 ), + ( 1.5, 0.0, 1.0 ), + ( 1.5, 1.0, 1.0 ), + ( 0.5, 1.0, 1.0 ), # 12-15 + ] + mesh = _make_ugrid( coords ) + _add_cell( mesh, VTK_HEXAHEDRON, ( 0, 1, 2, 3, 4, 5, 6, 7 ) ) + _add_cell( mesh, VTK_HEXAHEDRON, ( 8, 9, 10, 11, 12, 13, 14, 15 ) ) + return mesh + + +# ────────────────────────────────────────────────────────────── +# Internal construction helpers +# ────────────────────────────────────────────────────────────── + + +def _make_ugrid( coords: list[ tuple[ float, float, float ] ] ) -> vtkUnstructuredGrid: + """Create an empty vtkUnstructuredGrid with the given point coordinates.""" + points = vtkPoints() + for c in coords: + points.InsertNextPoint( c ) + mesh = vtkUnstructuredGrid() + mesh.SetPoints( points ) + mesh.Allocate() + return mesh + + +def _add_cell( mesh: vtkUnstructuredGrid, cell_type: int, point_ids: tuple[ int, ...] ) -> None: + """Append one cell to *mesh* in-place.""" + ids = vtkIdList() + for pid in point_ids: + ids.InsertNextId( pid ) + mesh.InsertNextCell( cell_type, ids ) diff --git a/mesh-doctor/tests/test_collocatedNodes.py b/mesh-doctor/tests/test_collocatedNodes.py index 281f7df2..05dabc8a 100644 --- a/mesh-doctor/tests/test_collocatedNodes.py +++ b/mesh-doctor/tests/test_collocatedNodes.py @@ -1,11 +1,12 @@ # SPDX-License-Identifier: Apache-2.0 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. -# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto, Jacques Franc import pytest from typing import Iterator, Tuple from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkTetra, vtkUnstructuredGrid, VTK_TETRA from geos.mesh_doctor.actions.collocatedNodes import Options, meshAction +from .helpers import build_mesh_with_collocated_nodes def getPoints() -> Iterator[ Tuple[ vtkPoints, int ] ]: @@ -70,3 +71,28 @@ def test_wrongSupportElements() -> None: assert len( result.nodesBuckets ) == 0 assert len( result.wrongSupportElements ) == 1 assert result.wrongSupportElements[ 0 ] == 0 + + +def test_collocated_shared_face_nodes() -> None: + """Two hexes with duplicated shared-face nodes: 4 buckets with exact pair indices. + + ``build_mesh_with_collocated_nodes()`` has 16 points instead of 12. + The four shared-face nodes on the x=1 plane are duplicated in place: + + original │ duplicate + ──────────┼─────────── + 1 │ 8 + 2 │ 11 + 5 │ 12 + 6 │ 15 + + Each bucket is a tuple ``(first_inserted_id, rejected_duplicate_id)``. + """ + mesh = build_mesh_with_collocated_nodes() + result = meshAction( mesh, Options( tolerance=1.e-6 ) ) + + assert len( result.wrongSupportElements ) == 0 + assert len( result.nodesBuckets ) == 4 + + buckets = sorted( result.nodesBuckets ) + assert buckets == [ ( 1, 8 ), ( 2, 11 ), ( 5, 12 ), ( 6, 15 ) ] diff --git a/mesh-doctor/tests/test_elementVolumes.py b/mesh-doctor/tests/test_elementVolumes.py index 669a5e92..ef1e623c 100644 --- a/mesh-doctor/tests/test_elementVolumes.py +++ b/mesh-doctor/tests/test_elementVolumes.py @@ -1,10 +1,11 @@ # SPDX-License-Identifier: Apache-2.0 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. -# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto, Jacques Franc import numpy from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkCommonDataModel import VTK_TETRA, vtkCellArray, vtkTetra, vtkUnstructuredGrid from geos.mesh_doctor.actions.elementVolumes import Options, meshAction +from .helpers import build_mesh_with_negative_volume def test_simpleTet() -> None: @@ -41,3 +42,32 @@ def test_simpleTet() -> None: result = meshAction( mesh, Options( minVolume=0. ) ) assert len( result.elementVolumes ) == 0 + + +def test_negative_volume_hex_detected() -> None: + """Hexahedron with top/bottom node groups swapped has negative signed volume. + + ``build_mesh_with_negative_volume()`` contains two hexes: + - cell 0: correctly oriented unit cube → positive volume (+1.0) + - cell 1: top and bottom face groups swapped → negative volume (−1.0) + + With ``minVolume=0`` only cells whose volume ≤ minVolume are reported. + Cell 1 (volume ≈ −1.0) must appear; cell 0 must not. + """ + mesh = build_mesh_with_negative_volume() + result = meshAction( mesh, Options( minVolume=0.0 ) ) + + assert len( result.elementVolumes ) == 1 + assert result.elementVolumes[ 0 ][ 0 ] == 1 # second cell + assert result.elementVolumes[ 0 ][ 1 ] < 0.0 # negative volume + + +def test_positive_volume_hex_not_reported() -> None: + """The correctly oriented hex in the same mesh is not flagged.""" + mesh = build_mesh_with_negative_volume() + + # With threshold just below the correct cell volume (+1.0) nothing is reported + result = meshAction( mesh, Options( minVolume=-0.5 ) ) + + flagged_ids = { entry[ 0 ] for entry in result.elementVolumes } + assert 0 not in flagged_ids # cell 0 is valid, volume ≈ +1.0 diff --git a/mesh-doctor/tests/test_euler.py b/mesh-doctor/tests/test_euler.py new file mode 100644 index 00000000..a82ba732 --- /dev/null +++ b/mesh-doctor/tests/test_euler.py @@ -0,0 +1,74 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2026 TotalEnergies. +# SPDX-FileContributor: Jacques Franc +import pytest +from vtkmodules.vtkFiltersCore import vtkFeatureEdges +from vtkmodules.vtkFiltersGeometry import vtkGeometryFilter +from geos.mesh_doctor.actions.euler import Options, meshAction +from .helpers import build_volume_with_non_manifold_edge, build_surface_with_non_manifold_edge + +_OPTS = Options() + + +def test_volume_non_manifold_edge_detected() -> None: + """Bowtie of two hexes sharing only one edge has a non-manifold boundary. + + ``build_volume_with_non_manifold_edge()`` produces two hexahedra that share + only the edge (1,1,0)→(1,1,1). On the extracted boundary surface that edge + is touched by four faces (two from each hex), making it non-manifold. + + Expected topology (V=14, E=23, F=12, C=2): + solidEulerCharacteristic = 14 - 23 + 12 - 2 = 1 + numNonManifoldEdges = 1 + numConnectedComponents = 1 (hexes share 2 points → same component) + """ + result = meshAction( build_volume_with_non_manifold_edge(), _OPTS ) + + assert result.numNonManifoldEdges == 1 + assert result.solidEulerCharacteristic == 1 + assert result.numConnectedComponents == 1 + assert result.numBoundaryEdges == 0 + + +def test_volume_non_manifold_topology_counts() -> None: + """Exact V/E/F/C counts for the bowtie configuration.""" + result = meshAction( build_volume_with_non_manifold_edge(), _OPTS ) + + assert result.numVertices == 14 + assert result.numEdges == 23 + assert result.numFaces == 12 + assert result.numCells == 2 + + +def test_surface_non_manifold_edge_detected() -> None: + """Three quads sharing one edge form a non-manifold surface. + + ``build_surface_with_non_manifold_edge()`` returns a vtkUnstructuredGrid of + 3 VTK_QUAD cells. Edge (1,0,0)→(1,1,0) is shared by all three faces. + + ``euler.meshAction`` requires 3D cells and raises ``RuntimeError`` on a + pure surface mesh. The non-manifold edge is verified directly via + ``vtkFeatureEdges`` — the same filter ``euler.py`` uses internally. + """ + ugrid = build_surface_with_non_manifold_edge() + + # euler.meshAction cannot handle a mesh with no 3D cells + with pytest.raises( RuntimeError ): + meshAction( ugrid, _OPTS ) + + # Convert to polydata and run the same vtkFeatureEdges query euler.py uses + gf = vtkGeometryFilter() + gf.SetInputData( ugrid ) + gf.Update() + surface = gf.GetOutput() + + fe = vtkFeatureEdges() + fe.SetInputData( surface ) + fe.BoundaryEdgesOff() + fe.ManifoldEdgesOff() + fe.NonManifoldEdgesOn() + fe.FeatureEdgesOff() + fe.Update() + + # The edge (1,0,0)→(1,1,0) is shared by 3 quads → 1 non-manifold edge segment + assert fe.GetOutput().GetNumberOfCells() > 0 diff --git a/mesh-doctor/tests/test_orphan.py b/mesh-doctor/tests/test_orphan.py new file mode 100644 index 00000000..4f0b515c --- /dev/null +++ b/mesh-doctor/tests/test_orphan.py @@ -0,0 +1,69 @@ +# SPDX-License-Identifier: Apache-2.0 +# SPDX-FileCopyrightText: Copyright 2023-2026 TotalEnergies. +# SPDX-FileContributor: Jacques Franc +from geos.mesh_doctor.actions.orphan2d import Options, meshAction +from .helpers import ( build_hex_mesh, build_hex_mesh_with_fractures, build_tetra_mesh_with_fractures, + build_mesh_with_orphan_2d_cells ) + +_OPTS = Options( orphanVtkOutput=None, cleanVtkOutput=None ) + + +def test_dedicated_orphan_cell() -> None: + """Mesh with one matched and one orphaned 2D quad is correctly classified. + + ``build_mesh_with_orphan_2d_cells()`` contains: + - 1 hex (3D) + - 1 quad that IS the bottom face of the hex → matched + - 1 quad at z=2 that is NOT any face of the hex → orphan, cell index 2 + - 1 quad at z=[0,1] that is diagonal to hex and is NOT any face of the hex → orphan, cell index 3 + """ + result = meshAction( build_mesh_with_orphan_2d_cells(), _OPTS ) + + assert result.total2dCells == 3 + assert result.total3dCells == 1 + assert result.matched2dCells == 1 + assert result.orphaned2dCells == 2 + assert result.orphaned2dIndices == [ 2, 3 ] + + +def test_clean_hex_has_no_orphans() -> None: + """A pure volumetric hex mesh (no 2D cells) reports zero orphans.""" + result = meshAction( build_hex_mesh( 4, 4, 4 ), _OPTS ) + + assert result.total2dCells == 0 + assert result.orphaned2dCells == 0 + + +def test_hex_fracture_quads_are_orphans() -> None: + """Fracture quads appended with MergePointsOff have disjoint point indices. + + ``build_hex_mesh_with_fractures`` combines the volume and fracture sub-meshes + via ``vtkAppendFilter(MergePointsOff)``. Because the fracture quad node IDs + differ from those of the hex faces (even at the same coordinates), + ``orphan2d.meshAction`` cannot match them → all fracture quads are orphans. + + 5×5×5 grid with 1 fracture plane and default trim margins (offx=offy=1) + produces 4 fracture quads, all orphaned. + """ + mesh = build_hex_mesh_with_fractures( 5, 5, 5, nfrac=1 ) + result = meshAction( mesh, _OPTS ) + + assert result.total3dCells == 64 + assert result.total2dCells == 4 + assert result.matched2dCells == 0 + assert result.orphaned2dCells == 4 + + +def test_tetra_fracture_triangles_are_orphans() -> None: + """Same MergePointsOff behaviour applies to tetra meshes with tri fractures. + + 5×5×5 Delaunay tetra mesh with 1 fracture plane produces 8 triangles, + all orphaned. + """ + mesh = build_tetra_mesh_with_fractures( 5, 5, 5, nfrac=1 ) + result = meshAction( mesh, _OPTS ) + + assert result.total3dCells > 0 + assert result.total2dCells == 8 + assert result.matched2dCells == 0 + assert result.orphaned2dCells == 8 diff --git a/mesh-doctor/tests/test_selfIntersectingElements.py b/mesh-doctor/tests/test_selfIntersectingElements.py index a477e912..30f7c2ed 100644 --- a/mesh-doctor/tests/test_selfIntersectingElements.py +++ b/mesh-doctor/tests/test_selfIntersectingElements.py @@ -1,9 +1,10 @@ # SPDX-License-Identifier: Apache-2.0 # SPDX-FileCopyrightText: Copyright 2023-2024 TotalEnergies. -# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto +# SPDX-FileContributor: Thomas Gazolla, Alexandre Benedicto, Jacques Franc from vtkmodules.vtkCommonCore import vtkPoints from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkHexahedron, vtkUnstructuredGrid, VTK_HEXAHEDRON from geos.mesh_doctor.actions.selfIntersectingElements import Options, meshAction +from .helpers import build_mesh_with_negative_volume def test_jumbledHex() -> None: @@ -43,3 +44,23 @@ def test_jumbledHex() -> None: assert len( result.invalidCellIds[ "intersectingFacesElements" ] ) == 1 assert result.invalidCellIds[ "intersectingFacesElements" ][ 0 ] == 0 + + +def test_inverted_hex_orientation_from_gen_mesh() -> None: + """Inverted hex (top/bottom node groups swapped) is caught by orientation checks. + + ``build_mesh_with_negative_volume()`` contains two hexes: + - cell 0: correctly oriented unit cube + - cell 1: top and bottom face groups swapped → non-convex + wrong face orientation + + ``vtkCellValidator`` flags cell 1 via: + - ``nonConvexElements`` (concavity caused by the inversion) + - ``facesOrientedIncorrectlyElements`` (inward-pointing normals) + + Cell 0 must NOT appear in either of those categories. + """ + mesh = build_mesh_with_negative_volume() + result = meshAction( mesh, Options( minDistance=0. ) ) + + assert result.invalidCellIds[ "nonConvexElements" ] == [ 1 ] + assert result.invalidCellIds[ "facesOrientedIncorrectlyElements" ] == [ 1 ]