Skip to content

Checks: Add failed mesh elements/sides to debug mesh #160

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Open
kopperp wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Checks: Add failed mesh elements/sides to debug mesh #160

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
1192845
Checks: Align CheckOrient logic better
kopperp Aug 16, 2026
247c413
Checks: Add failed sides to debug mesh
kopperp Aug 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 15 additions & 9 deletions pyhope/basis/basis_connect.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,9 +27,9 @@
# ----------------------------------------------------------------------------------------------------------------------------------
from __future__ import annotations
import re
import sys
from typing import Final, Optional, cast
from collections.abc import Iterable
from operator import itemgetter
# ----------------------------------------------------------------------------------------------------------------------------------
# Third-party libraries
# ----------------------------------------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -142,10 +142,12 @@ def CheckConnect() -> None:
""" Check if the mesh is correctly connected
"""
# Local imports ----------------------------------------
import pyhope.io.io_vars as io_vars
import pyhope.output.output as hopout
import pyhope.mesh.mesh_vars as mesh_vars
from pyhope.common.common_parallel import run_in_parallel
from pyhope.common.common_vars import np_mtp
from pyhope.io.io_debug import DebugIO
from pyhope.readintools.readintools import GetLogical
# ------------------------------------------------------

Expand Down Expand Up @@ -178,13 +180,13 @@ def CheckConnect() -> None:
ordering = False, # noqa: E251
)
else:
res = [elem for elem in elems if check_sides(elem, failed_only=True)]
res = [check_sides(elem, failed_only=True) for elem in elems]

if len(res) > 0: # pragma: no cover
# Flatten per-element results (skip None placeholders)
results = tuple(result for elem_results in res if isinstance(elem_results, Iterable) and elem_results is not None
for result in elem_results) # noqa: E272
# Flatten per-element results (skip None placeholders)
results = tuple(result for elem_results in res if isinstance(elem_results, Iterable) and elem_results is not None
for result in elem_results) # noqa: E272

if len(results): # pragma: no cover
nGeo: Final[int] = mesh_vars.nGeo
sides: Final[list] = mesh_vars.sides
points: Final[np.ndarray] = mesh_vars.mesh.points
Expand Down Expand Up @@ -216,7 +218,7 @@ def CheckConnect() -> None:
nodes = elem.nodes[face_to_nodes( side.face, elem.type, nGeo)]
nbnodes = nbelem.nodes[face_to_nodes(nbside.face, nbelem.type, nGeo)]

print()
hopout.info('')
# Check if side is oriented inwards
errStr = 'Side connectivity does not match the calculated neighbour side'
print(hopout.warn(errStr, length=len(errStr)+16))
Expand All @@ -235,5 +237,9 @@ def CheckConnect() -> None:
print(hopout.warn('- Coordinates : [' + ' '.join(f'{s:12.3f}' for s in points[nbnodes[-1, 0]]) + ']')) # noqa: E271
print(hopout.warn('- Coordinates : [' + ' '.join(f'{s:12.3f}' for s in points[nbnodes[-1, -1]]) + ']')) # noqa: E271

hopout.warning(f'Connectivity check failed for {len(results)} / {nconn} connections!')
sys.exit(1)
# Write a low-order debug mesh if requested
if io_vars.debugmesh:
hopout.info('')
DebugIO(errSides=list(map(itemgetter(1), results)))

hopout.error(f'Connectivity check failed for {len(results)} / {nconn} connections!')
183 changes: 91 additions & 92 deletions pyhope/basis/basis_orient.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,9 @@
from __future__ import annotations
import gc
import re
import sys
from typing import Final, Optional
from typing import Final, Optional, cast
from collections.abc import Iterable
from operator import itemgetter
# ----------------------------------------------------------------------------------------------------------------------------------
# Third-party libraries
# ----------------------------------------------------------------------------------------------------------------------------------
Expand All @@ -38,84 +39,100 @@
# Typing libraries
# ----------------------------------------------------------------------------------------------------------------------------------
import typing
if typing.TYPE_CHECKING:
from pyhope.common.common_numba import NUMBA_AVAILABLE
if typing.TYPE_CHECKING or NUMBA_AVAILABLE:
import numpy.typing as npt
from pyhope.mesh.mesh_vars import ELEM
# ----------------------------------------------------------------------------------------------------------------------------------
# Local imports
# ----------------------------------------------------------------------------------------------------------------------------------
import pyhope.mesh.mesh_vars as mesh_vars
from pyhope.basis.basis_watertight import eval_nsurf
from pyhope.mesh.mesh_common import LINMAP
from pyhope.mesh.mesh_common import dir_to_nodes, faces
# ----------------------------------------------------------------------------------------------------------------------------------
# Local definitions
# ----------------------------------------------------------------------------------------------------------------------------------
from pyhope.mesh.mesh_common import face_to_nodes, faces
# ==================================================================================================================================


def check_orientation(ionodes : npt.NDArray,
elemType : int,
VdmEqToGP: npt.NDArray[np.float64],
DGP : npt.NDArray[np.float64],
weights : npt.NDArray[np.float64],
) -> tuple[bool, Optional[str]]:
def check_orientation(elem : ELEM,
VdmEqToGP : npt.NDArray[np.float64],
DGP : npt.NDArray[np.float64],
weights : npt.NDArray[np.float64],
failed_only: bool = False,
) -> Optional[list[tuple]]:
""" Check the orientation of the surface normals
"""
mapLin = LINMAP(elemType, order=mesh_vars.nGeo)
iopoints = mesh_vars.mesh.points
mapnodes = ionodes[mapLin]
points = iopoints[mapnodes]
points = mesh_vars.mesh.points
nGeo = mesh_vars.nGeo
elemType = elem.type

if elemType % 10 != 8:
return None

# Center of element
cElem = points.reshape(-1, 3).mean(axis=0)
cElem = points[elem.nodes].reshape(-1, 3).mean(axis=0)
results = None

success = True
sface = None
# Fall back to geometry topology faces if sides array is not initialized
sides = getattr(mesh_vars, 'sides', None)
sideIter = ((sides[s].face, s) for s in elem.sides) if (elem.sides is not None and sides is not None) else \
((f, None) for f in faces(elemType))

for face in faces(elemType):
indices, doTransp = dir_to_nodes(face, elemType, mesh_vars.nGeo)
fnodes = mapnodes[indices] if not doTransp else mapnodes[indices].transpose()
fpoints = iopoints[fnodes]
for face, SideID in sideIter:
idx = cast(np.ndarray, elem.nodes)[face_to_nodes(face, elemType, nGeo)]
fpoints = points[idx]

# Calculate high-order surface normal vector via Gauss integration
nSurf = eval_nsurf(fpoints.transpose(2, 0, 1), VdmEqToGP, DGP, weights)
nSurf = eval_nsurf(fpoints.transpose(2, 0, 1), VdmEqToGP, DGP, weights)

# Vector pointing from face center to element center
fCenter = fpoints.reshape(-1, 3).mean(axis=0)
vCenter = cElem - fCenter
fCenter = fpoints.reshape(-1, 3).mean(axis=0)
vCenter = cElem - fCenter

success = np.dot(nSurf, vCenter) <= 0.

# If requested, only return errors
if failed_only and success:
continue

# Lazily initialize results on first failure
if results is None:
results = []
results.append((success, elem.elemID, face, SideID))

# Avoid creating empty lists on elem_results
if results is None:
return None if failed_only else []

if np.dot(nSurf, vCenter) > 0.:
success = False
sface = face
break
return success, sface
return results


def process_chunk(chunk: tuple) -> list:
"""Process a chunk of elements by checking surface normal orientation
"""
# Only keep failures to reduce memory and avoid building large arrays of successes
chunk_results = []
for elemChunk in chunk:
iElem, ionodes, elemType = elemChunk
elem_result = check_orientation(ionodes, elemType,
for elem in chunk:
elem_result = check_orientation(elem,
process_chunk.VdmEqToGP, # ty: ignore[unresolved-attribute]
process_chunk.DGP, # ty: ignore[unresolved-attribute]
process_chunk.weights) # ty: ignore[unresolved-attribute]
# Append a lightweight sentinel (None) for successes, actual failure list otherwise
chunk_results.append((elem_result, iElem) if not elem_result[0] else None)
process_chunk.weights, # ty: ignore[unresolved-attribute]
failed_only=True)
chunk_results.append(elem_result)
return chunk_results


def CheckOrient() -> None:
""" Check if element surface normals point outwards
"""
# Local imports ----------------------------------------
import pyhope.mesh.mesh_vars as mesh_vars
import pyhope.io.io_vars as io_vars
import pyhope.output.output as hopout
import pyhope.mesh.mesh_vars as mesh_vars
from pyhope.basis.basis_basis import barycentric_weights, legendre_gauss_nodes
from pyhope.basis.basis_basis import calc_vandermonde, polynomial_derivative_matrix
from pyhope.basis.basis_watertight import init_worker
from pyhope.common.common_parallel import run_in_parallel
from pyhope.common.common_vars import np_mtp
from pyhope.io.io_debug import DebugIO
from pyhope.readintools.readintools import GetLogical
# ------------------------------------------------------

Expand All @@ -127,8 +144,8 @@ def CheckOrient() -> None:
if not checkSurfaceNormals:
return None

mesh = mesh_vars.mesh
nGeo: Final[int] = mesh_vars.nGeo
nGeo: Final[int] = mesh_vars.nGeo
elems: Final[list] = mesh_vars.elems

# Setup mathematical structures identical to CheckWatertight
xEq: Final[npt.NDArray[np.float64]] = np.linspace(-1., 1., nGeo+1)
Expand All @@ -139,58 +156,40 @@ def CheckOrient() -> None:
VdmEqToGP: Final[npt.NDArray[np.float64]] = calc_vandermonde(nGeo+1, nGeo+1, wBaryEq, xEq, xGP)
weights: Final[npt.NDArray[np.float64]] = np.outer(wGP, wGP)

elemNames: Final[dict] = mesh_vars.ELEMTYPE.name
elemKeys : Final = mesh_vars.ELEMTYPE.type.keys()
nElems = 0
passedTypes = []

for elemType in mesh.cells_dict:
# Only consider three-dimensional types
if not any(s in elemType for s in elemKeys):
continue

# Only consider hexahedrons
if 'hexahedron' not in elemType:
passedTypes.append(elemType)
continue

# Get the elements
ioelems = mesh.get_cells_type(elemType)
nIOElems = ioelems.shape[0]

if isinstance(elemType, str):
elemType = elemNames[elemType]

# Prepare elements for parallel processing
if np_mtp > 0:
tasks = tuple((iElem, ioelems[iElem - nElems], elemType)
for iElem in range(nElems, nElems + nIOElems))
# Run in parallel with a chunk size
# > Dispatch the tasks to the workers, minimum 10 tasks per worker, maximum 1000 tasks per worker
res = run_in_parallel(process_chunk, # noqa: E251
tasks, # noqa: E251
chunk_size = max(1, min(1000, max(10, int(len(tasks)/(40.*np_mtp))))), # noqa: E251
initializer = init_worker, # noqa: E251
init_args = (process_chunk, VdmEqToGP, DGP, weights), # noqa: E251
ordering = False, # noqa: E251
)
else:
res = np.fromiter(((check_orientation(ioelems[iElem - nElems], elemType, VdmEqToGP, DGP, weights), iElem)
for iElem in range(nElems, nElems + nIOElems)), dtype=object)

if len(res) > 0 and not np.all([success for (success, _), _ in res]):
failed_elems = [(iElem + 1, face) for (success, face), iElem in res if not success]
for iElem, face in failed_elems:
print(hopout.warn(f'> Element {iElem}, Side {face}'))
sys.exit(1)

# Add to nElems
nElems += nIOElems
# Prepare elements for parallel processing
if np_mtp > 0:
# Run in parallel with a chunk size
# > Dispatch the tasks to the workers, minimum 10 tasks per worker, maximum 1000 tasks per worker
res = run_in_parallel(process_chunk,
elems,
chunk_size = max(1, min(1000, max(10, int(len(elems)/(40.*np_mtp))))), # noqa: E251
initializer = init_worker, # noqa: E251
init_args = (process_chunk, VdmEqToGP, DGP, weights), # noqa: E251
ordering = False, # noqa: E251
)
else:
res = [check_orientation(elem, VdmEqToGP, DGP, weights, failed_only=True) for elem in elems]

# Flatten results while filtering empty items
results = tuple(result for elem_results in res if isinstance(elem_results, Iterable) and elem_results is not None
for result in elem_results) # noqa: E272

if len(results):
for result in cast(tuple[tuple], results):
_, elemID, face, sideID = result
hopout.info('')
print(hopout.warn(f'Side is oriented inwards! Element {elemID + 1}, Face {face}, Side {sideID + 1}'))

if io_vars.debugmesh:
hopout.info('')
DebugIO(errSides=list(map(itemgetter(3), results)))

hopout.error(f'Surface normals check failed for {len(results)} / {len(elems)} elements!')

# Warn if we passed any element types
if len(passedTypes) > 0:
print(hopout.warn('Ignored element type{}: {}'.format('s' if len(passedTypes) > 1 else '',
[re.sub(r"\d+$", "", s) for s in passedTypes])))
if any(e.type % 10 != 8 for e in elems):
elemTypes = list({e.type for e in elems if e.type % 10 != 8})
print(hopout.warn('Ignored element type: {}'.format([re.sub(r"\d+$", "", mesh_vars.ELEMTYPE.inam[e][0]) for e in elemTypes]))) # ruff: ignore[line-too-long]

# Run garbage collector to release memory
gc.collect()
22 changes: 14 additions & 8 deletions pyhope/basis/basis_watertight.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@
import re
from typing import Final, Optional, cast
from collections.abc import Callable, Iterable
from operator import itemgetter
# ----------------------------------------------------------------------------------------------------------------------------------
# Third-party libraries
# ----------------------------------------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -256,12 +257,14 @@ def CheckWatertight() -> None:
""" Check if the mesh is watertight
"""
# Local imports ----------------------------------------
import pyhope.io.io_vars as io_vars
import pyhope.output.output as hopout
import pyhope.mesh.mesh_vars as mesh_vars
from pyhope.basis.basis_basis import barycentric_weights, legendre_gauss_nodes
from pyhope.basis.basis_basis import calc_vandermonde, polynomial_derivative_matrix
from pyhope.common.common_parallel import run_in_parallel
from pyhope.common.common_vars import np_mtp
from pyhope.io.io_debug import DebugIO
from pyhope.readintools.readintools import GetLogical
# ------------------------------------------------------

Expand Down Expand Up @@ -314,15 +317,13 @@ def CheckWatertight() -> None:
ordering = False, # noqa: E251
)
else:
res = [elem for elem in elems if check_sides(elem,
VdmEqToGP, DGP, weights,
failed_only=True)]
res = [check_sides(elem, VdmEqToGP, DGP, weights, failed_only=True) for elem in elems]

if len(res) > 0: # pragma: no cover
# Flatten per-element results (skip None placeholders)
results = tuple(result for elem_results in res if isinstance(elem_results, Iterable) and elem_results is not None
for result in elem_results) # noqa: E272
# Flatten per-element results (skip None placeholders)
results = tuple(result for elem_results in res if isinstance(elem_results, Iterable) and elem_results is not None
for result in elem_results) # noqa: E272

if len(results): # pragma: no cover
# Compute total number of checked connections without materializing all results
nconn = 0
for SideID, side in enumerate(sides):
Expand Down Expand Up @@ -352,7 +353,7 @@ def CheckWatertight() -> None:
nodes = elem.nodes[face_to_nodes( side.face, elem.type, nGeo)]
nbnodes = nbelem.nodes[face_to_nodes(nbside.face, nbelem.type, nGeo)]

print()
hopout.info('')
# Check if side is oriented inwards
errStr = 'Side is oriented inwards!' if nSurfErr < 0 \
else f'Surface normals are not within tolerance {nSurfErr:9.6e} > {tol:9.6e}'
Expand All @@ -374,6 +375,11 @@ def CheckWatertight() -> None:
print(hopout.warn('- Coordinates : [' + ' '.join(f'{s:12.3f}' for s in points[nbnodes[-1, 0]]) + ']')) # noqa: E271
print(hopout.warn('- Coordinates : [' + ' '.join(f'{s:12.3f}' for s in points[nbnodes[-1, -1]]) + ']')) # noqa: E271

# Write a low-order debug mesh if requested
if io_vars.debugmesh:
hopout.info('')
DebugIO(errSides=list(map(itemgetter(1), results)))

hopout.error(f'Watertightness check failed for {len(results)} / {nconn} connections!')

# Run garbage collector to release memory
Expand Down
Loading
Loading