diff --git a/pyhope/basis/basis_connect.py b/pyhope/basis/basis_connect.py index 078264e..45c7354 100644 --- a/pyhope/basis/basis_connect.py +++ b/pyhope/basis/basis_connect.py @@ -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 # ---------------------------------------------------------------------------------------------------------------------------------- @@ -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 # ------------------------------------------------------ @@ -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 @@ -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)) @@ -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!') diff --git a/pyhope/basis/basis_orient.py b/pyhope/basis/basis_orient.py index 860ae7e..c1cbb40 100644 --- a/pyhope/basis/basis_orient.py +++ b/pyhope/basis/basis_orient.py @@ -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 # ---------------------------------------------------------------------------------------------------------------------------------- @@ -38,57 +39,70 @@ # 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: @@ -96,26 +110,29 @@ def process_chunk(chunk: tuple) -> list: """ # 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 # ------------------------------------------------------ @@ -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) @@ -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() diff --git a/pyhope/basis/basis_watertight.py b/pyhope/basis/basis_watertight.py index 3952256..72fa551 100644 --- a/pyhope/basis/basis_watertight.py +++ b/pyhope/basis/basis_watertight.py @@ -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 # ---------------------------------------------------------------------------------------------------------------------------------- @@ -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 # ------------------------------------------------------ @@ -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): @@ -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}' @@ -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 diff --git a/pyhope/io/io_debug.py b/pyhope/io/io_debug.py index 12f78fe..4268802 100644 --- a/pyhope/io/io_debug.py +++ b/pyhope/io/io_debug.py @@ -129,13 +129,14 @@ def FillSideData(msides: list, sInv : dict[str, int], pMap : npt.NDArray, bcs : list, + errOut: bool = False, ) -> tuple[dict, dict[str, list]]: # Create a defaultdict since we have optional members sidedata = defaultdict(lambda: [[] for _ in range(len(sypes))]) # Populate connectivity and data - for side in tuple(s for s in msides if s.bcid is not None): + for side in tuple(s for s in msides if s.bcid is not None or errOut): sideType = 'triangle' if side.sideType == 3 else 'quad' sidx = sInv[sideType] @@ -144,13 +145,16 @@ def FillSideData(msides: list, sides[sideType].append(sideNodes) # Add the sideData - bcID = side.bcid - bc = bcs[bcID] sidedata['ElemID' ][sidx].append(side.elemID + 1) - sidedata['BCID' ][sidx].append(bcID + 1) - sidedata['BCType' ][sidx].append(bc.type[0] ) - sidedata['BCState' ][sidx].append(bc.type[2] ) - sidedata['BCAlpha' ][sidx].append(bc.type[3] ) + + # Add the boundary dataa + if side.bcid is not None: + bcID = side.bcid + bc = bcs[bcID] + sidedata['BCID' ][sidx].append(bcID + 1) + sidedata['BCType' ][sidx].append(bc.type[0] ) + sidedata['BCState' ][sidx].append(bc.type[2] ) + sidedata['BCAlpha' ][sidx].append(bc.type[3] ) return sides, dict(sidedata) @@ -224,7 +228,8 @@ def FillNodeData(melems: list, return nodes, dict(nodedata) -def DebugIO() -> None: +def DebugIO(errElems: Optional[list] = None, + errSides: Optional[list] = None) -> None: """ Routine to output the debug mesh. Downcast the existing PyHOPE format to first order, enrich with debug information and output in XDMF format @@ -236,6 +241,10 @@ def DebugIO() -> None: from pyhope.mesh.mesh_vars import ELEMTYPE # ------------------------------------------------------ + # Write a low-order debug mesh if requested + if not io_vars.debugmesh: + return None + mesh : Final = mesh_vars.mesh mpoints: Final[np.ndarray] = mesh.points melems : Final[list] = mesh_vars.elems @@ -243,6 +252,10 @@ def DebugIO() -> None: bcs : Final[list] = mesh_vars.bcs pname : Final[str] = io_vars.projectname + # Convert error lists to sets + errElems = set(errElems) if errElems is not None else () + errSides = set(errSides) if errSides is not None else () + # Create empty meshio objects elems = {} elemtypes = set() @@ -269,13 +282,12 @@ def DebugIO() -> None: elemtypes.add(elemType) # Add the first-order sides to the sides set - for sideID in melem.sides: # ty: ignore [not-iterable] - # Only consider boundary sides - if msides[sideID].bcid is not None: - sideType = 'triangle' if msides[sideID].sideType == 3 else 'quad' - if sideType not in sInv: - sInv[sideType] = len(sidetypes) - sidetypes.add(sideType) + for sideID in (s for s in melem.sides if msides[s].bcid is not None or s in (errSides if errSides is not None else ())): # ty: ignore [not-iterable] + # Only consider boundary/error sides + sideType = 'triangle' if msides[sideID].sideType == 3 else 'quad' + if sideType not in sInv: + sInv[sideType] = len(sidetypes) + sidetypes.add(sideType) # Create ordered mapping from first-order points to high-order points points = np.concatenate([np.asarray(melem.nodes)[:melem.type % 10] for melem in melems]) @@ -300,12 +312,15 @@ def DebugIO() -> None: # Fill elem data elems, elemdata = FillElemData(melems, elems, hasIJK, types, tInv, pMap) # Fill the side data (optional) - sides, sidedata = FillSideData(msides, sides, sypes, sInv, pMap, bcs) + sides, sidedata = FillSideData(msides, sides, sypes, sInv, pMap, bcs, bool(errSides) and not any(s.bcid is not None for s in msides)) # ruff: ignore[line-too-long] # Fill the edge data (optional) edges, edgedata = FillEdgeData(melems, edges, hasFEM, pMap) # Fill the node data (optional) nodes, nodedata = FillNodeData(melems, nodes, hasFEM, pMap) + # Create the output list + debugOut = [] + # Update points to unique first-order coords coords = mpoints[pMap] # Find the mapping from the cell keys to the elemtypes @@ -314,11 +329,6 @@ def DebugIO() -> None: # Ensure cell_data lists are aligned to the actual cell block order used by meshio.Mesh elemdata = {k: [np.asarray(v[idx]) for idx in elemOrder] for k, v in elemdata.items()} eleminfo = {'name': 'Volume'} - # Clean-up for memory safety - del elemOrder - - # Create the output list - debugOut = [] # Create the final debugElem with first-order elements debugElem = meshio.Mesh(points = coords, # noqa: E251 @@ -328,14 +338,38 @@ def DebugIO() -> None: ) debugOut.append(debugElem) + if len(errElems): # pragma: no cover + # Create empty meshio objects + errorElems = {} + + # Prepare element and side containers + for t in elemtypes: + errorElems.setdefault(t, []) + + # Fill error elem data + errorElems, errorElemdata = FillElemData((e for e in melems if e.elemID in errElems), errorElems, hasIJK, types, tInv, pMap) # ruff: ignore[line-too-long] + + # Ensure cell_data lists are aligned to the actual cell block order used by meshio.Mesh + errorElemdata = {k: [np.asarray(v[idx]) for idx in elemOrder] for k, v in errorElemdata.items()} + errorEleminfo = {'name': 'Volume [Error]'} + + # Create the final debugElem with first-order elements + debugElem = meshio.Mesh(points = coords, # noqa: E251 + cells = errorElems, # noqa: E251 + cell_data = errorElemdata, # noqa: E251 + info = errorEleminfo, # noqa: E251 + ) + debugOut.append(debugElem) + + # Clean-up for memory safety + del elemOrder + # Find the mapping from the side keys to the elemtypes sideOrder = [sInv[cb] for cb in sides] # Ensure cell_data lists are aligned to the actual cell block order used by meshio.Mesh sidedata = {k: [np.asarray(v[idx]) for idx in sideOrder] for k, v in sidedata.items()} sideinfo = {'name': 'Surface'} - # Clean-up for memory safety - del sideOrder # Create the final debugSide with first-order elements debugSide = meshio.Mesh(points = coords, # noqa: E251 @@ -345,6 +379,32 @@ def DebugIO() -> None: ) debugOut.append(debugSide) + if len(errSides): # pragma: no cover + # Create empty meshio objects + errorSides = {} + + # Prepare element and side containers + for st in sidetypes: + errorSides.setdefault(st, []) + + # Fill the error side data + errorSides, errorSidedata = FillSideData((s for s in msides if s.sideID in errSides), errorSides, sypes, sInv, pMap, bcs, True) # ruff: ignore[line-too-long] + + # Ensure cell_data lists are aligned to the actual cell block order used by meshio.Mesh + errorSidedata = {k: [np.asarray(v[idx]) for idx in sideOrder] for k, v in errorSidedata.items()} + errorSideinfo = {'name': 'Surface [Error]'} + + # Create the final debugSide with first-order elements + debugSide = meshio.Mesh(points = coords, # noqa: E251 + cells = errorSides, # noqa: E251 + cell_data = errorSidedata, # noqa: E251 + info = errorSideinfo, # noqa: E251 + ) + debugOut.append(debugSide) + + # Clean-up for memory safety + del sideOrder + if edgedata is not None: # Ensure cell_data lists are aligned to the actual cell block order used by meshio.Mesh edgedata = {k: [np.asarray(v)] for k, v in edgedata.items()} @@ -372,3 +432,10 @@ def DebugIO() -> None: fname = f'{pname}_DebugMesh.xdmf' hopout.routine(f'Writing XDMF mesh to "{fname}"') meshio.xdmf.main.XdmfWriter(fname, debugOut) + + if len(errElems): # pragma: no cover + print('│' + hopout.Symbols.INFO[:3] + + f'Reason: Detected {len(errElems)} / {len(melems)} erroneous elements, written to "Volume [Error]"') + if len(errSides): # pragma: no cover + print('│' + hopout.Symbols.INFO[:3] + + f'Reason: Detected {len(errSides)} / {len(msides)} erroneous sides, written to "Surface [Error]"') diff --git a/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index 626eded..8682685 100644 --- a/pyhope/mesh/extrude/mesh_extrude.py +++ b/pyhope/mesh/extrude/mesh_extrude.py @@ -61,7 +61,7 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: 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_vars import nGeo + from pyhope.mesh.mesh_vars import ELEM, nGeo from pyhope.mesh.topology.mesh_topology import appendBCSet from pyhope.readintools.readintools import GetInt, GetStr # ------------------------------------------------------ @@ -377,8 +377,9 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: if isinstance(elemType, str): elemType = elemNames[elemType] - result = check_orientation(ionodes, elemType, VdmEqToGP, DGP, weights) - if result is not None and not result[0]: + elem = ELEM(type=elemType, nodes=ionodes, elemID=0, sides=None) + results = check_orientation(elem, VdmEqToGP, DGP, weights) + if results and any(not res[0] for res in results): hopout.error('Extruded element has inward pointing normal vector. Wrong extrusion direction?') # Run garbage collector to release memory diff --git a/pyproject.toml b/pyproject.toml index 82a0a6f..94fe5b6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -243,7 +243,9 @@ ignore_names = [ "__init__", "data_counter", "h5_file", "_meshio_to_gmsh_order", # Exclude library functions "Basis", "Mapping", "Mesh", - "mesh_format_to_tensor_product" + "mesh_format_to_tensor_product", + # Exclude lookup functions + "dir_to_nodes", ] # Exclude unit tests exclude = ["check_unittest.py"]