Skip to content

Checks: Perform normal calculation using full integration #156

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

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

Checks: Perform normal calculation using full integration #156

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
2e49f91
Checks: Perform normal calculation using full integration
kopperp Jun 26, 2026
ff7c488
Checks: Re-use normal calculation, update README.md
kopperp Jun 26, 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
12 changes: 8 additions & 4 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,7 @@ $ pyhope tutorials/1-01-cartbox/parameter.ini
│ INIT OUTPUT...
│ ProjectName │ 1-01-cartbox │ *CUSTOM │
│ OutputFormat │ 0 [HDF5] │ *CUSTOM │
│ OutputBytes │ 0 [int32] │ DEFAULT │
│ DebugMesh │ T │ *CUSTOM │
│ DebugVisu │ F │ *CUSTOM │
├─────────────────────────────────────────────
Expand Down Expand Up @@ -122,23 +123,26 @@ $ pyhope tutorials/1-01-cartbox/parameter.ini
├── BUILD DATA STRUCTURE...
├────
├── Removing duplicate points
├── Ensuring normals point outward
├────
│ CheckSurfaceNormals │ True │ DEFAULT │
│ Processing Elements |█████████████████████████████████| 512/512 [100%] in 0.0s (24000.00/s)
├────
├── Generating sides
├─────────────────────────────────────────────
│ SORT MESH...
├────
│ MeshSorting │ 1 [SFC] │ DEFAULT │
│ MeshSortingSFC │ 0 [default] │ DEFAULT │
├────
├── Sorting elements along space-filling curve
├─────────────────────────────────────────────
│ CHECK NORMALS POINTING OUTWARDS...
├────
│ CheckSurfaceNormals │ True │ DEFAULT │
│ Processing Elements |█████████████████████████████████| 512/512 [100%] in 0.1s (2559.71/s)
├─────────────────────────────────────────────
│ CONNECT MESH...
├────
│ doPeriodicCorrect │ False │ DEFAULT │
│ doMortars │ True │ DEFAULT │
│ doMortarRebuild │ 1 [auto] │ DEFAULT │
│ Processing Sides |█████████████████████████████████| 3072/3072 [100%] in 0.0s (24000.00/s)
├────
│ Number of sides : 3072
Expand Down
72 changes: 36 additions & 36 deletions pyhope/mesh/mesh_orient.py → pyhope/basis/basis_orient.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@
# Local imports
# ----------------------------------------------------------------------------------------------------------------------------------
import pyhope.mesh.mesh_vars as mesh_vars
from pyhope.common.common_numba import jit, types
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
# ----------------------------------------------------------------------------------------------------------------------------------
Expand All @@ -53,29 +53,11 @@
# ==================================================================================================================================


@jit((types.float64)(types.float64[:, :, ::1], types.float64[::1]), nopython=True, cache=True, nogil=True)
def eval_dotprod(fpoints: npt.NDArray[np.float64], nVecFace: npt.NDArray[np.float64]) -> np.float64:
# nVecFace = nVecFace / np.linalg.norm(nVecFace)
# > Manually compute norm
nVecFace = nVecFace / np.sqrt(nVecFace[0]*nVecFace[0] + nVecFace[1]*nVecFace[1] + nVecFace[2]*nVecFace[2])

vec1 = fpoints[-1, 0, :] - fpoints[0, 0, :]
vec2 = fpoints[0, -1, :] - fpoints[0, 0, :]

# normal = np.cross(vec1, vec2)
# > Manually compute cross product
normal = np.empty_like(vec1)
normal[0] = vec1[1] * vec2[2] - vec1[2] * vec2[1]
normal[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2]
normal[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0]

# Dot product and check if normal points outwards
# > Manually compute dot product
return nVecFace[0]*normal[0] + nVecFace[1]*normal[1] + nVecFace[2]*normal[2]


def check_orientation(ionodes : npt.NDArray,
elemType: int,
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]]:
""" Check the orientation of the surface normals
"""
Expand All @@ -91,19 +73,18 @@ def check_orientation(ionodes : npt.NDArray,
sface = None

for face in faces(elemType):
# Center of face
indices, doTransp = dir_to_nodes(face, elemType, mesh_vars.nGeo)
fnodes = mapnodes[indices] if not doTransp else mapnodes[indices].transpose()
fpoints = iopoints[fnodes]

# Tangent and normal vectors
nVecFace = cElem - fpoints.reshape(-1, 3).mean(axis=0)
# nVecFace = nVecFace / np.linalg.norm(nVecFace)
# nVecFace = nVecFace / np.sqrt(np.dot(nVecFace, nVecFace))
# > Manually compute dot product
dotprod = eval_dotprod(fpoints, nVecFace)
# Calculate high-order surface normal vector via Gauss integration
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

if dotprod < 0.:
if np.dot(nSurf, vCenter) > 0.:
success = False
sface = face
break
Expand All @@ -117,29 +98,46 @@ def process_chunk(chunk: tuple) -> list:
chunk_results = []
for elemChunk in chunk:
iElem, ionodes, elemType = elemChunk
elem_result = check_orientation(ionodes, elemType)
elem_result = check_orientation(ionodes, elemType,
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)
return chunk_results


def OrientMesh() -> None:
def CheckOrient() -> None:
# Local imports ----------------------------------------
import pyhope.mesh.mesh_vars as mesh_vars
import pyhope.output.output as hopout
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.readintools.readintools import GetLogical
# ------------------------------------------------------

hopout.routine('Ensuring normals point outward')
hopout.separator()
hopout.info('CHECK NORMALS POINTING OUTWARDS...')
hopout.sep()

checkSurfaceNormals = GetLogical('CheckSurfaceNormals')
if not checkSurfaceNormals:
return None

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

# Setup mathematical structures identical to CheckWatertight
xEq: Final[npt.NDArray[np.float64]] = np.linspace(-1., 1., nGeo+1)
wBaryEq: Final[npt.NDArray[np.float64]] = barycentric_weights(nGeo+1, xEq)

xGP, wGP = legendre_gauss_nodes(nGeo+1)
DGP: Final[npt.NDArray[np.float64]] = polynomial_derivative_matrix(nGeo+1, xGP)
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()
Expand Down Expand Up @@ -172,10 +170,12 @@ def OrientMesh() -> None:
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), iElem)
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]):
Expand Down
10 changes: 5 additions & 5 deletions pyhope/basis/basis_watertight.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,7 @@ def eval_nsurf(XGeo: npt.NDArray[np.float64],
# # dXdetaGP = np.moveaxis(dXdetaGP, 0, 1).reshape(3, -1) # Flatten for cross computation (slower)
# dXdetaGP = dXdetaGP.reshape(3, -1) # Flatten for cross computation

# Compute derivatives at all Gauss points using matrix multiplications to avoid extra copies
# Compute derivatives at all Gauss points
dXdxiGP = np.empty_like(xGP)
dXdetaGP = np.empty_like(xGP)
DT = DGP.T
Expand All @@ -120,7 +120,7 @@ def eval_nsurf(XGeo: npt.NDArray[np.float64],
# Integrate over the Gauss points
# return -np.sum(nVecW, axis=(1, 2)) # Sum over the last two axes

# Compute the weighted cross product integral directly
# Compute the weighted cross product integral
NSurf0, NSurf1, NSurf2 = eval_dotprod(dXdetaGP, dXdxiGP, weights)

return np.array((NSurf0, NSurf1, NSurf2), dtype=xGP.dtype)
Expand Down Expand Up @@ -243,9 +243,9 @@ def process_chunk(chunk: tuple) -> list:
chunk_results = []
for elem in chunk:
elem_result = check_sides(elem,
process_chunk.VdmEqToGP, # pyright: ignore[reportFunctionMemberAccess] # ty: ignore[unresolved-attribute]
process_chunk.DGP, # pyright: ignore[reportFunctionMemberAccess] # ty: ignore[unresolved-attribute]
process_chunk.weights, # pyright: ignore[reportFunctionMemberAccess] # ty: ignore[unresolved-attribute]
process_chunk.VdmEqToGP, # ty: ignore[unresolved-attribute]
process_chunk.DGP, # ty: ignore[unresolved-attribute]
process_chunk.weights, # ty: ignore[unresolved-attribute]
failed_only=True)
# Append a lightweight sentinel (None) for successes, actual failure list otherwise
chunk_results.append(elem_result)
Expand Down
2 changes: 2 additions & 0 deletions pyhope/io/io_debug.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,8 @@
#
# tree = ET.ElementTree(vtkfile)
# tree.write(filename, encoding='utf-8', xml_declaration=True)


@cache
def isValidInt(s: Any) -> bool: # noqa: ANN401
try:
Expand Down
17 changes: 15 additions & 2 deletions pyhope/mesh/extrude/mesh_extrude.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,12 +53,14 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh:
# Local imports ----------------------------------------
import pyhope.mesh.mesh_vars as mesh_vars
import pyhope.output.output as hopout
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_orient import check_orientation
from pyhope.common.common_progress import ProgressBar
from pyhope.common.common_template import LoadTemplate
from pyhope.common.common_tools import temporary_assign
from pyhope.io.io_gmsh import GMSHCELLTYPES
from pyhope.mesh.mesh_common import NDOFperElemType
from pyhope.mesh.mesh_orient import check_orientation
from pyhope.mesh.mesh_vars import nGeo
from pyhope.mesh.topology.mesh_topology import appendBCSet
from pyhope.readintools.readintools import GetInt, GetStr
Expand Down Expand Up @@ -354,6 +356,17 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh:

# Temporarily assign mesh_vars.mesh
with temporary_assign(mesh_vars, 'mesh', mesh):
# Compute the equidistant point set used by meshIO
xEq: Final[npt.NDArray[np.float64]] = np.linspace(-1., 1., nGeo+1)
wBaryEq: Final[npt.NDArray[np.float64]] = barycentric_weights(nGeo+1, xEq)

xGP, wGP = legendre_gauss_nodes(nGeo+1)
DGP: Final[npt.NDArray[np.float64]] = polynomial_derivative_matrix(nGeo+1, xGP)
VdmEqToGP: Final[npt.NDArray[np.float64]] = calc_vandermonde(nGeo+1, nGeo+1, wBaryEq, xEq, xGP)

# Compute the weights
weights: Final[npt.NDArray[np.float64]] = np.outer(wGP, wGP) # Shape: (N_GP+1, N_GP+1)

# Check if the surface normal of the first cell points outwards
for elemType in tuple(s for s in mesh.cells_dict if 'hexahedron' in s):
# Only check the first element, other elements get covered in OrientMesh
Expand All @@ -364,7 +377,7 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh:
if isinstance(elemType, str):
elemType = elemNames[elemType]

if not check_orientation(ionodes, elemType)[0]:
if not check_orientation(ionodes, elemType, VdmEqToGP, DGP, weights)[0]:
hopout.error('Extruded element has inward pointing normal vector. Wrong extrusion direction?')

# Run garbage collector to release memory
Expand Down
9 changes: 6 additions & 3 deletions pyhope/script/pyhope_cli.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ def main() -> None:
import pyhope.output.output as hopout
from pyhope.basis.basis_connect import CheckConnect
from pyhope.basis.basis_jacobian import CheckJacobians
from pyhope.basis.basis_orient import CheckOrient
from pyhope.basis.basis_watertight import CheckWatertight
from pyhope.check.check import Check
from pyhope.common.common import DefineCommon, InitCommon
Expand All @@ -54,7 +55,6 @@ def main() -> None:
from pyhope.mesh.fem.fem import FEMConnect
from pyhope.mesh.mesh import DefineMesh, InitMesh, GenerateMesh
from pyhope.mesh.mesh_duplicates import EliminateDuplicates
from pyhope.mesh.mesh_orient import OrientMesh
from pyhope.mesh.mesh_sides import GenerateSides
from pyhope.mesh.mesh_sort import SortMesh
from pyhope.mesh.mesh_mortar import RebuildMortarGeometry
Expand Down Expand Up @@ -123,12 +123,15 @@ def main() -> None:
hopout.sep()

EliminateDuplicates()
if not args.skip_checks:
OrientMesh()

# Build our data structures
GenerateSides()
SortMesh()

# Perform the mesh checks not depending on connectivity
if not args.skip_checks:
CheckOrient()

ConnectMesh()
TransformMesh()

Expand Down
Loading