diff --git a/README.md b/README.md index d6f46fdf..0ff2a65d 100644 --- a/README.md +++ b/README.md @@ -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 │ ├───────────────────────────────────────────── @@ -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 diff --git a/pyhope/mesh/mesh_orient.py b/pyhope/basis/basis_orient.py similarity index 76% rename from pyhope/mesh/mesh_orient.py rename to pyhope/basis/basis_orient.py index d9906ed6..860ae7e0 100644 --- a/pyhope/mesh/mesh_orient.py +++ b/pyhope/basis/basis_orient.py @@ -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 # ---------------------------------------------------------------------------------------------------------------------------------- @@ -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 """ @@ -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 @@ -117,22 +98,29 @@ 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') @@ -140,6 +128,16 @@ def OrientMesh() -> None: 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() @@ -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]): diff --git a/pyhope/basis/basis_watertight.py b/pyhope/basis/basis_watertight.py index 76992d72..3952256e 100644 --- a/pyhope/basis/basis_watertight.py +++ b/pyhope/basis/basis_watertight.py @@ -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 @@ -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) @@ -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) diff --git a/pyhope/io/io_debug.py b/pyhope/io/io_debug.py index 55c69307..30193e87 100644 --- a/pyhope/io/io_debug.py +++ b/pyhope/io/io_debug.py @@ -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: diff --git a/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index 5f818175..14efc1a1 100644 --- a/pyhope/mesh/extrude/mesh_extrude.py +++ b/pyhope/mesh/extrude/mesh_extrude.py @@ -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 @@ -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 @@ -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 diff --git a/pyhope/script/pyhope_cli.py b/pyhope/script/pyhope_cli.py index c1a37b1f..4691bd9f 100755 --- a/pyhope/script/pyhope_cli.py +++ b/pyhope/script/pyhope_cli.py @@ -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 @@ -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 @@ -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()