diff --git a/pyhope/basis/basis_connect.py b/pyhope/basis/basis_connect.py index 9598bfe..f85bfea 100644 --- a/pyhope/basis/basis_connect.py +++ b/pyhope/basis/basis_connect.py @@ -26,6 +26,7 @@ # Standard libraries # ---------------------------------------------------------------------------------------------------------------------------------- from __future__ import annotations +import gc import re from typing import Final, Optional, cast from collections.abc import Iterable @@ -159,6 +160,9 @@ def CheckConnect() -> None: if not checkConnectivity: return None + # Temporally disable the garbage collector + gc.disable() + # Check all sides elems: Final[list] = mesh_vars.elems @@ -243,3 +247,7 @@ def CheckConnect() -> None: DebugIO(errSides=list(map(itemgetter(1), results))) hopout.error(f'Connectivity check failed for {len(results)} / {nconn} connections!') + + # Run garbage collector to release memory + gc.collect() + gc.enable() diff --git a/pyhope/basis/basis_orient.py b/pyhope/basis/basis_orient.py index 9850f16..9e9053e 100644 --- a/pyhope/basis/basis_orient.py +++ b/pyhope/basis/basis_orient.py @@ -144,6 +144,9 @@ def CheckOrient() -> None: if not checkSurfaceNormals: return None + # Temporally disable the garbage collector + gc.disable() + nGeo: Final[int] = mesh_vars.nGeo elems: Final[list] = mesh_vars.elems @@ -193,3 +196,4 @@ def CheckOrient() -> None: # Run garbage collector to release memory gc.collect() + gc.enable() diff --git a/pyhope/basis/basis_watertight.py b/pyhope/basis/basis_watertight.py index c49089f..b755853 100644 --- a/pyhope/basis/basis_watertight.py +++ b/pyhope/basis/basis_watertight.py @@ -276,6 +276,9 @@ def CheckWatertight() -> None: if not checkWatertightness: return None + # Temporally disable the garbage collector + gc.disable() + nGeo: Final[int] = mesh_vars.nGeo # Compute the equidistant point set used by meshIO @@ -384,3 +387,4 @@ def CheckWatertight() -> None: # Run garbage collector to release memory gc.collect() + gc.enable() diff --git a/pyhope/common/common_tools.py b/pyhope/common/common_tools.py index cbf9f94..755df02 100644 --- a/pyhope/common/common_tools.py +++ b/pyhope/common/common_tools.py @@ -26,6 +26,7 @@ # Standard libraries # ---------------------------------------------------------------------------------------------------------------------------------- from __future__ import annotations +# import gc import time # from sortedcontainers import SortedDict from collections import defaultdict @@ -118,6 +119,21 @@ def temporary_assign(obj: Any, attr: Any, value: Any) -> GeneratorType: # noqa: setattr(obj, attr, orig) +# @contextmanager +# def pause_gc(collect: bool = False) -> GeneratorType: +# """ Temporally pause the cyclic garbace collector +# """ +# wasEnabled = gc.isenabled() +# gc.disable() +# try: +# yield +# finally: +# if wasEnabled: +# if collect: +# gc.collect() +# gc.enable() + + class IndexedLists: def __init__(self) -> None: # PERF: The unsorted dict is faster for adding and removing elements diff --git a/pyhope/mesh/connect/connect.py b/pyhope/mesh/connect/connect.py index 88bed33..b07172d 100644 --- a/pyhope/mesh/connect/connect.py +++ b/pyhope/mesh/connect/connect.py @@ -30,6 +30,7 @@ import gc import sys from collections import defaultdict +from itertools import chain from typing import Final, Optional, cast # from multiprocessing import Pool # ---------------------------------------------------------------------------------------------------------------------------------- @@ -68,6 +69,34 @@ def flip_analytic(side: int, nbside: npt.NDArray) -> int: return find_index(nbside, side) +def stack_corners(sides: list) -> npt.NDArray: + """ Corner nodes of all sides (nSides, 4), triangles are padded with -1 + """ + nSides = len(sides) + try: + corners = np.array([s.corners for s in sides], dtype=np.int64) + except ValueError: # Mixed triangles and quadrilaterals + corners = np.full((nSides, 4), -1, dtype=np.int64) + lengths = np.fromiter((len(s.corners) for s in sides), dtype=np.int64, count=nSides) + for nCorners in (3, 4): + ids = np.flatnonzero(lengths == nCorners) + if ids.size > 0: + corners[ids, :nCorners] = np.array([sides[i].corners for i in ids], dtype=np.int64) + return corners + + if corners.shape[1] == 3: # Only triangles, pad at the end so that the first column is always a corner + corners = np.pad(corners, ((0, 0), (0, 1)), constant_values=-1) + + return corners + + +def side_keys(corners: npt.NDArray) -> list[int]: + """ Build side keys from sorted corner nodes + """ + # Two sides share a key if they share all corners, triangles are padded with -1 + return [hash(tuple(c) if c[0] >= 0 else tuple(c[1:])) for c in np.sort(corners, axis=1).tolist()] + + def connect_sides(sideIDs: list[int], sides: list, flipID: int) -> None: """ Connect the master and slave sides """ @@ -99,6 +128,30 @@ def connect_sides(sideIDs: list[int], sides: list, flipID: int) -> None: # slaveSide.nbLocSide = sides[sideIDs[0]].locSide # noqa: E251 +def connect_plain_sides(sides: list, corners: npt.NDArray, pairs: npt.NDArray) -> None: + """ Connect all pairs of sides + > Only valid for sides without boundary condition + """ + first, second = pairs[:, 0], pairs[:, 1] + hit = corners[second] == corners[first, :1] + flip = np.where(hit.any(axis=1), hit.argmax(axis=1), -1) + 1 + + # Periodic sides must be sorted along the SFC, the master side is the one with the lower ID + master = np.minimum(first, second).tolist() + slave = np.maximum(first, second).tolist() + + for iMaster, iSlave, flipID in zip(master, slave, flip.tolist(), strict=True): + masterSide = sides[iMaster] + slaveSide = sides[iSlave] + + masterSide.MS = 1 # noqa: E251 + masterSide.connection = iSlave # noqa: E251 + masterSide.flip = flipID # noqa: E251 + slaveSide.MS = 0 # noqa: E251 + slaveSide.connection = iMaster # noqa: E251 + slaveSide.flip = flipID # noqa: E251 + + def find_bc_index(bcs: list | tuple, key: str) -> Optional[int]: """ Find the index of a BC from its name in the list of BCs """ @@ -266,6 +319,9 @@ def ConnectMesh() -> None: doPeriodicCorrect = mesh_vars.doPeriodicCorrect doMortars = mesh_vars.doMortars + # Temporally disable the garbage collector + gc.disable() + mesh = mesh_vars.mesh elems = mesh_vars.elems sides = mesh_vars.sides @@ -308,11 +364,13 @@ def ConnectMesh() -> None: # side_corners = dict(results) # else: # side_corners = {side: hash(np.sort(sides[side].corners).tobytes()) for elem in elems for side in elem.sides} + # Build the reverse dictionary + sideCorners = stack_corners(sides) + sideKeys = side_keys(sideCorners) corner_side = defaultdict(list) - for elem in elems: - for side in cast(list | np.ndarray, elem.sides): - corner_side[hash(tuple(sorted(sides[side].corners)))].append(side) + for side in chain.from_iterable(elem.sides for elem in elems): + corner_side[sideKeys[side]].append(side) # > Create a dict containing only the periodic corners peri_corners = {} @@ -402,87 +460,77 @@ def ConnectMesh() -> None: # Try to connect the inner / periodic sides passedTypes = {} + internal = [] for val in corner_side.values(): match len(val): case 1: # BC side continue case 2: # Internal side - sideIDs = val - # Flip pyramids - # if elems[sides[sideIDs[0]].elemID].type % 10 != 5 and \ - # elems[sides[sideIDs[1]].elemID].type % 10 == 5: - # sideIDs = sideIDs[::-1] - - side0 = sides[sideIDs[0]] - side1 = sides[sideIDs[1]] - - # Ignore true BCs on inner faces for connect - if side0.bcid is not None: - if side1.bcid is None: - hopout.error('Found internal face with one BC, exiting...') - if bcs[side0.bcid].type[0] not in (0, 1, 100): - continue - - corners = side0.corners - nbcorners = side1.corners - - # Translate to periodic nodes if required - if side0.bcid is not None and side1.bcid is not None and bcs[side1.bcid].type[0] == 1: - nbcorners = np.fromiter((mesh_vars.periNodes[(s, bcs[side1.bcid].name)] for s in side1.corners), dtype=int) - - flipID = flip_analytic(corners[0], nbcorners) + 1 - - # Sanity check the flip with the other nodes - # > INFO: MOVED TO OWN CHECKCONNECT ROUTINE - # elem = (elems[side0.elemID], elems[side1.elemID]) - # if elem[0].type % 10 == 8 and elem[1].type % 10 == 8: - # # Map the meshio nodes to the tensor-product nodes - # elemType = elem[0].type - # nodes = elem[0].nodes[sidetovol2(nGeo, 0 , side0.face, elemType)] - # nbNodes = elem[1].nodes[sidetovol2(nGeo, flipID, side1.face, elemType)] - # - # # Translate to periodic nodes if required - # if side0.bcid is not None and side1.bcid is not None and bcs[side1.bcid].type[0] == 1: - # nbNodes = np.vectorize(lambda s: mesh_vars.periNodes[(s, bcs[side1.bcid].name)], otypes=[int])(nbNodes) - # - # # Check if the node IDs match - # if not np.array_equal(nodes, nbNodes): - # # Print the node IDs - # print(hopout.warn(f'NodeIDs side[-]: {nodes }')) - # print(hopout.warn(f'NodeIDs side[+]: {nbNodes}')) - # hopout.error('Error in connectivity check, node IDs do not match!') - - # Connect the sides - connect_sides(sideIDs, sides, flipID) - # Update the progress bar - bar.step(2) - - # Use guard clauses to avoid unnecessary checks - if side0.bcid is None: - continue # No boundary condition on first side - if bcs[side0.bcid].type[0] != 1: - continue # Not a periodic BC on first side - if bcs[side1.bcid].type[0] != 1: - hopout.error('Found internal side with inconsistent BC types, exiting...') - if not doPeriodicCorrect: - continue # Periodic correction not enabled - - # At this point, we know both sides have periodic BCs - iVV = bcs[side0.bcid].type[3] - VV = vvs[np.abs(iVV) - 1]['Dir'] * np.sign(iVV) - locSides = tuple(sides[s] for s in sideIDs) # noqa: E272 - locElems = tuple(elems[s.elemID] for s in locSides) - - # Only update hexahedral elements - if any(e.type % 10 != 8 for e in locElems): - for e in locElems: - passedTypes[e.type] = passedTypes.get(e.type, 0) + 1 - else: - periodic_update(locSides, locElems, VV) - + internal.append(val) case _: # Zero or more than 2 sides hopout.error('Found internal side with more than two adjacent elements, exiting...', traceback=True) + # Connect sides without any boundary condition + internal = np.array(internal, dtype=np.int64).reshape(-1, 2) + hasBC = np.fromiter((s.bcid is not None for s in sides), dtype=bool, count=len(sides)) + plain = ~hasBC[internal].any(axis=1) + connect_plain_sides(sides, sideCorners, internal[plain]) + + bar.step(2*int(np.count_nonzero(plain))) + + # The remaining sides carry an inner / periodic BC + for val in internal[~plain].tolist(): + sideIDs = val + # Flip pyramids + # if elems[sides[sideIDs[0]].elemID].type % 10 != 5 and \ + # elems[sides[sideIDs[1]].elemID].type % 10 == 5: + # sideIDs = sideIDs[::-1] + + side0 = sides[sideIDs[0]] + side1 = sides[sideIDs[1]] + + # Ignore true BCs on inner faces for connect + if side0.bcid is not None: + if side1.bcid is None: + hopout.error('Found internal face with one BC, exiting...') + if bcs[side0.bcid].type[0] not in (0, 1, 100): + continue + + corners = side0.corners + nbcorners = side1.corners + + # Translate to periodic nodes if required + if side0.bcid is not None and side1.bcid is not None and bcs[side1.bcid].type[0] == 1: + nbcorners = np.fromiter((mesh_vars.periNodes[(s, bcs[side1.bcid].name)] for s in side1.corners), dtype=int) + + # Connect the sides + connect_sides(sideIDs, sides, flip_analytic(corners[0], nbcorners) + 1) + # Update the progress bar + bar.step(2) + + # Use guard clauses to avoid unnecessary checks + if side0.bcid is None: + continue # No boundary condition on first side + if bcs[side0.bcid].type[0] != 1: + continue # Not a periodic BC on first side + if bcs[side1.bcid].type[0] != 1: + hopout.error('Found internal side with inconsistent BC types, exiting...') + if not doPeriodicCorrect: + continue # Periodic correction not enabled + + # At this point, we know both sides have periodic BCs + iVV = bcs[side0.bcid].type[3] + VV = vvs[np.abs(iVV) - 1]['Dir'] * np.sign(iVV) + locSides = tuple(sides[s] for s in sideIDs) # noqa: E272 + locElems = tuple(elems[s.elemID] for s in locSides) + + # Only update hexahedral elements + if any(e.type % 10 != 8 for e in locElems): + for e in locElems: + passedTypes[e.type] = passedTypes.get(e.type, 0) + 1 + else: + periodic_update(locSides, locElems, VV) + if passedTypes: print(hopout.warn(hopout.Colors.WARN + '─'*(46-16) + hopout.Colors.END)) print(hopout.warn('Periodic correction skipped for elements:')) @@ -572,6 +620,7 @@ def ConnectMesh() -> None: # Run garbage collector to release memory gc.collect() + gc.enable() # hopout.sep() # hopout.info('CONNECT MESH DONE!') diff --git a/pyhope/mesh/extrude/mesh_extrude.py b/pyhope/mesh/extrude/mesh_extrude.py index c4486d0..6679a61 100644 --- a/pyhope/mesh/extrude/mesh_extrude.py +++ b/pyhope/mesh/extrude/mesh_extrude.py @@ -91,6 +91,9 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: hopout.sep() hopout.routine(f' Template: {extrTemplate}') + # Temporally disable the garbage collector + gc.disable() + # Setup the extrusion extrShifts = LoadTemplate(extrTemplate.strip().lower(), __file__, 'Extrusion').ExtrudeTemplate() @@ -384,6 +387,7 @@ def MeshExtrude(mesh: meshio.Mesh) -> meshio.Mesh: # Run garbage collector to release memory gc.collect() + gc.enable() hopout.sep() return mesh diff --git a/pyhope/mesh/mesh_builtin.py b/pyhope/mesh/mesh_builtin.py index 3308eed..f84a3b3 100644 --- a/pyhope/mesh/mesh_builtin.py +++ b/pyhope/mesh/mesh_builtin.py @@ -113,6 +113,9 @@ def MeshCartesian() -> meshio.Mesh: offsetp = 0 offsets = 0 + # Temporally disable the garbage collector + gc.disable() + # GMSH only supports mesh elements within a single model # > https://gitlab.onelab.info/gmsh/gmsh/-/issues/2836 gmsh.model.add('Domain') @@ -410,6 +413,7 @@ def MeshCartesian() -> meshio.Mesh: # Run garbage collector to release memory gc.collect() + gc.enable() # Store the element types mesh_vars.nZones = nZones diff --git a/pyhope/mesh/mesh_duplicates.py b/pyhope/mesh/mesh_duplicates.py index a7b97ad..2294490 100644 --- a/pyhope/mesh/mesh_duplicates.py +++ b/pyhope/mesh/mesh_duplicates.py @@ -163,8 +163,12 @@ def EliminateDuplicates() -> None: from pyhope.common.common_unique import unique from pyhope.mesh.connect.connect import find_bc_index # ------------------------------------------------------ + hopout.routine('Removing duplicate points') + # Temporally disable the garbage collector + gc.disable() + bcs: Final[list] = mesh_vars.bcs vvs: Final[list] = mesh_vars.vvs @@ -316,3 +320,4 @@ def EliminateDuplicates() -> None: # Run garbage collector to release memory gc.collect() + gc.enable() diff --git a/pyhope/mesh/mesh_external.py b/pyhope/mesh/mesh_external.py index 0be53d6..dae75fd 100644 --- a/pyhope/mesh/mesh_external.py +++ b/pyhope/mesh/mesh_external.py @@ -62,6 +62,9 @@ def MeshExternal() -> meshio.Mesh: hopout.routine('Setting boundary conditions') hopout.sep() + # Temporally disable the garbage collector + gc.disable() + # Load the boundary conditions nBCs = CountOption('BoundaryName') mesh_vars.bcs = [BC() for _ in range(nBCs)] @@ -162,6 +165,7 @@ def MeshExternal() -> meshio.Mesh: # Run garbage collector to release memory gc.collect() + gc.enable() return mesh diff --git a/pyhope/mesh/mesh_sides.py b/pyhope/mesh/mesh_sides.py index a2921d3..a9f8115 100644 --- a/pyhope/mesh/mesh_sides.py +++ b/pyhope/mesh/mesh_sides.py @@ -25,7 +25,7 @@ # ---------------------------------------------------------------------------------------------------------------------------------- # Standard libraries # ---------------------------------------------------------------------------------------------------------------------------------- -from typing import cast +import gc # ---------------------------------------------------------------------------------------------------------------------------------- # Third-party libraries # ---------------------------------------------------------------------------------------------------------------------------------- @@ -54,7 +54,6 @@ def GenerateSides() -> None: mesh = mesh_vars.mesh nElems = 0 nSides = 0 - sCount = 0 # INFO: Tuples should be faster than lists but cumbersome to add elements # mesh_vars.elems = () # mesh_vars.sides = () @@ -72,9 +71,12 @@ def GenerateSides() -> None: totalElems += mesh.get_cells_type(elemType).shape[0] # Use a moderate chunk size to bound intermediate progress updates - chunk = max(1, min(1000, max(10, int(len(elems)/(400))))) + chunk = max(1, min(1000, max(10, int(totalElems/(400))))) bar = ProgressBar(value=totalElems, title='│ Processing Elements', length=33, chunk=chunk) + # Temporally disable the garbage collector + gc.disable() + # Loop over all element types for elemType in mesh.cells_dict: # Only consider three-dimensional types @@ -82,13 +84,14 @@ def GenerateSides() -> None: continue # Get the elements - ioelems = mesh.get_cells_type(elemType) - elemMap = mesh_vars.ELEMMAP(elemType) - nIOElems = ioelems.shape[0] - nIOSides = len(faces(elemType)) + ioelems = mesh.get_cells_type(elemType) + elemMap = mesh_vars.ELEMMAP(elemType) + nIOElems = ioelems.shape[0] + elemFaces = faces(elemType) + nIOSides = len(elemFaces) # Map volume cell sets to elements - iocsets = mesh.cell_sets_dict + iocsets = mesh.cell_sets_dict elemSet: list[int | None] = [None for _ in range(nIOElems)] for key, val in iocsets.items(): @@ -116,70 +119,55 @@ def GenerateSides() -> None: elems.extend([ELEM() for _ in range(nIOElems )]) # pyright: ignore[reportArgumentType] sides.extend([SIDE() for _ in range(nIOElems*nIOSides)]) # pyright: ignore[reportArgumentType] + sideIDs = nSides + np.arange(nIOElems*nIOSides).reshape(nIOElems, nIOSides) + # Create the corner faces - corner_faces = tuple(face_to_cgns(s, elemType) for s in faces(elemType)) # noqa: E272 + corner_faces = tuple(face_to_cgns(s, elemType) for s in elemFaces) # noqa: E272 corner_length = tuple(len(c) for c in corner_faces) # noqa: E272 corner_index = tuple(np.array(c, dtype=int) for c in corner_faces) # noqa: E272 + corner_coll = tuple(ioelems[:, c] for c in corner_index) # noqa: E272 # Create dictionaries - for iElem in range(nElems, nElems+nIOElems): - nodes = ioelems[iElem-nElems] - # elems[iElem].update(type = elemMap, # noqa: E251 - # elemID = iElem, # noqa: E251 - # sides = [], # noqa: E251 - # nodes = ioelems[iElem]) # noqa: E251 - elems[iElem].type = elemMap # noqa: E251 - elems[iElem].elemID = iElem # noqa: E251 - elems[iElem].sides = [] # noqa: E251 - elems[iElem].nodes = nodes # noqa: E251 + for e in range(nIOElems): + iElem = nElems + e + elem = elems[iElem] + + # elem.update(type = elemMap, # noqa: E251 + # elemID = iElem, # noqa: E251 + # sides = [], # noqa: E251 + # nodes = ioelems[iElem]) # noqa: E251 + elem.type = elemMap # noqa: E251 + elem.elemID = iElem # noqa: E251 + elem.sides = sideIDs[e] # noqa: E251 + elem.nodes = ioelems[e] # noqa: E251 # Create the zone # > Account for different elemTypes in the index - if elemSet[iElem-nElems] is not None: - elems[iElem].zone = elemSet[iElem-nElems] - - # Create the sides - for key, val in enumerate(corner_length): - # sides[iSide].update(sideType=4) - sides[nSides + key].sideType = val - - # Assign corners to sides, CGNS format - for index, face in enumerate(faces(elemType)): - # PERF: It is faster to use a list comprehension than numpy.fromiter - # corners = np.fromiter((nodes[s] for s in corner_faces[index]), dtype=int) - # corners = nodes[corner_index[index]] - # sides[sCount].update(face = face, # noqa: E251 - # elemID = iElem, # noqa: E251 - # sideID = sCount, # noqa: E251 - # locSide = index+1, # noqa: E251 - # corners = np.array(corners)) # noqa: E251 - sides[sCount].face = face # noqa: E251 - sides[sCount].elemID = iElem # noqa: E251 - sides[sCount].sideID = sCount # noqa: E251 - sides[sCount].locSide = index+1 # noqa: E251 - # PERF: It is faster to use a list comprehension than numpy.fromiter - # sides[sCount].corners = corners # noqa: E251 - sides[sCount].corners = nodes[corner_index[index]] # noqa: E251 - sCount += 1 - - # Add to nSides - nSides += nIOSides + zone = elemSet[e] + if zone is not None: + elem.zone = zone + + # Create the sides and assign the corners to the sides, CGNS format + sCount = nSides + e*nIOSides + for index in range(nIOSides): + side = sides[sCount + index] + side.sideType = corner_length[index] # noqa: E251 + side.face = elemFaces[index] # noqa: E251 + side.elemID = iElem # noqa: E251 + side.sideID = sCount + index # noqa: E251 + side.locSide = index + 1 # noqa: E251 + side.corners = corner_coll[index][e] # noqa: E251 # Update progress bar bar.step() - # Add to nElems + # Add to nSides and nElems + nSides += nIOElems*nIOSides nElems += nIOElems # Close the progress bar bar.close() - # Append sides to elem - for side in sides: - elemID = cast(int, side.elemID) - sideID = cast(int, side.sideID) - cast(list, elems[elemID].sides).append(sideID) - - # Convert lists to numpy arrays - for elem in elems: - elem.sides = np.array(elem.sides) + # Re-enable the garbace collector + gc.collect() + gc.enable() diff --git a/pyhope/mesh/mesh_sort.py b/pyhope/mesh/mesh_sort.py index 1e76665..89bcd2f 100644 --- a/pyhope/mesh/mesh_sort.py +++ b/pyhope/mesh/mesh_sort.py @@ -76,47 +76,39 @@ def UpdateElemID(elems : list, import pyhope.io.io_vars as io_vars # ------------------------------------------------------ - totalElems = len(elems) - totalSides = len(sides) - - # Initialize sorted cells - sorted_elems = [None] * totalElems - sorted_sides = [None] * totalSides + sorted_elems = [elems[i] for i in sorted_indices.tolist()] # Python ints, not np.int64 + nElems = len(elems) + nSides = [len(e.sides) for e in sorted_elems] + + # Read ALL old side IDs before anything is overwritten + sorted_sides = [sides[s] for e in sorted_elems for s in e.sides] + offsetSides = np.concatenate(([0], np.cumsum(nSides))).tolist() + offsetESide = np.repeat(np.arange(nElems), nSides).tolist() + + for sideID, (s, elemSide) in enumerate(zip(sorted_sides, offsetESide, strict=True)): + s.sideID = sideID + s.elemID = elemSide + + # Calculate IJK position for this element + if nElemsIJK is not None: + nI, nJ = int(nElemsIJK[0]), int(nElemsIJK[1]) + ids = np.arange(nElems, dtype=np.int64) + k = ids // (nI*nJ) + rem = ids - k*nI*nJ + j = rem // nI + i = rem - j*nI + ijk = (np.stack((i, j, k), axis=1) + 1).astype(io_vars.outputbytes) bar.title('│ Processing Elements') - # Initialize the sideID and offset - offsetSide = 0 - sideID = 0 - # Overwrite the elem/side IDs - for newElemID, oldElemID in enumerate(sorted_indices): - elem = elems[oldElemID] - elem.elemID = newElemID - - # Calculate IJK position for this element - if nElemsIJK is not None: - k = newElemID // (nElemsIJK[0] * nElemsIJK[1]) - j = (newElemID - k * nElemsIJK[0] * nElemsIJK[1]) // nElemsIJK[0] - i = newElemID - k * nElemsIJK[0] * nElemsIJK[1] - j * nElemsIJK[0] - elem.elemIJK = np.array([i+1, j+1, k+1], dtype=io_vars.outputbytes) - - sorted_elems[newElemID] = elem - - # Correct the sideID - nSides = 0 - for key, val in enumerate(elem.sides): - side = sides[val] - side.sideID = offsetSide + key - side.elemID = newElemID - sorted_sides[sideID] = side - sideID += 1 - nSides += 1 + for newElemID, oldElemID in enumerate(sorted_elems): + oldElemID.elemID = newElemID # Correct the sideID - # nSides = len(elem.sides) - elem.sides = list(range(offsetSide, offsetSide + nSides)) - offsetSide += nSides + oldElemID.sides = list(range(offsetSides[newElemID], offsetSides[newElemID+1])) + if nElemsIJK is not None: + oldElemID.elemIJK = ijk[newElemID] bar.step() @@ -464,6 +456,9 @@ def SortMesh() -> None: meshsort = GetIntFromStr('MeshSorting', default=default) + # Temporally disable the garbage collector + gc.disable() + # Sort the mesh match meshsort: case MeshSort.NONE.value: @@ -480,3 +475,4 @@ def SortMesh() -> None: # Run garbage collector to release memory gc.collect() + gc.enable() diff --git a/pyhope/mesh/reader/reader_gambit.py b/pyhope/mesh/reader/reader_gambit.py index 7da94bc..9375490 100644 --- a/pyhope/mesh/reader/reader_gambit.py +++ b/pyhope/mesh/reader/reader_gambit.py @@ -80,6 +80,9 @@ def ReadGambit(fnames: list, mesh: meshio.Mesh) -> meshio.Mesh: hopout.sep() + # Temporally disable the garbage collector + gc.disable() + # Create an empty meshio object points = mesh.points if len(mesh.points.shape)>1 else np.zeros((0, 3), dtype=np.float64) pointl = cast(list, points.tolist()) @@ -200,7 +203,7 @@ def ReadGambit(fnames: list, mesh: meshio.Mesh) -> meshio.Mesh: try: _, zoneNum, _, zoneNElems, _, zoneMat, _, zoneFlags = tokens - zoneNum, zoneNElems, zoneMat, zoneFlags = int(zoneNum), int(zoneNElems), int(zoneMat), int(zoneFlags) + zoneNum, zoneNElems, _zoneMat, _zoneFlags = int(zoneNum), int(zoneNElems), int(zoneMat), int(zoneFlags) except ValueError: continue @@ -355,5 +358,6 @@ def ReadGambit(fnames: list, mesh: meshio.Mesh) -> meshio.Mesh: # Run garbage collector to release memory gc.collect() + gc.enable() return mesh diff --git a/pyhope/mesh/reader/reader_gmsh.py b/pyhope/mesh/reader/reader_gmsh.py index 3cf3737..30f723b 100644 --- a/pyhope/mesh/reader/reader_gmsh.py +++ b/pyhope/mesh/reader/reader_gmsh.py @@ -111,6 +111,9 @@ def ReadGMSH(fnames: list) -> meshio.Mesh: # Setup stacksize resource.setrlimit(resource.RLIMIT_STACK, (resource.RLIM_INFINITY, resource.RLIM_INFINITY)) + # Temporally disable the garbage collector + gc.disable() + hopout.sep() gmsh.initialize() @@ -269,6 +272,7 @@ def ReadGMSH(fnames: list) -> meshio.Mesh: # Run garbage collector to release memory gc.collect() + gc.enable() return mesh @@ -288,6 +292,9 @@ def BCCGNS(mesh: meshio.Mesh, fnames: list) -> meshio.Mesh: hopout.routine('Applying boundary conditions') # hopout.sep() + # Temporally disable the garbage collector + gc.disable() + points = mesh.points cells = mesh.cells # elems = mesh_vars.elems @@ -401,6 +408,7 @@ def BCCGNS(mesh: meshio.Mesh, fnames: list) -> meshio.Mesh: # Run garbage collector to release memory gc.collect() + gc.enable() hopout.sep() return mesh diff --git a/pyhope/mesh/reader/reader_hopr.py b/pyhope/mesh/reader/reader_hopr.py index 019a31b..44e56c3 100644 --- a/pyhope/mesh/reader/reader_hopr.py +++ b/pyhope/mesh/reader/reader_hopr.py @@ -64,6 +64,9 @@ def ReadHOPR(fnames: list, mesh: meshio.Mesh) -> meshio.Mesh: hopout.sep() + # Temporally disable the garbage collector + gc.disable() + # Create or re-use the meshio object points = mesh.points if len(mesh.points.shape)>1 else np.zeros((0, 3), dtype=np.float64) pointl = cast(list, points.tolist()) @@ -331,5 +334,6 @@ def ReadHOPR(fnames: list, mesh: meshio.Mesh) -> meshio.Mesh: # Run garbage collector to release memory gc.collect() + gc.enable() return mesh diff --git a/pyhope/mesh/topology/mesh_serendipity.py b/pyhope/mesh/topology/mesh_serendipity.py index 08e82ce..1364cfc 100644 --- a/pyhope/mesh/topology/mesh_serendipity.py +++ b/pyhope/mesh/topology/mesh_serendipity.py @@ -58,6 +58,9 @@ def convertSerendipityToFullLagrange(mesh: meshio.Mesh) -> meshio.Mesh: hopout.routine('Converting serendipity to full langrange mesh') + # Temporally disable the garbage collector + gc.disable() + # Copy original points points = mesh.points.copy() @@ -250,5 +253,6 @@ def convertSerendipityToFullLagrange(mesh: meshio.Mesh) -> meshio.Mesh: # Run garbage collector to release memory gc.collect() + gc.enable() return mesh diff --git a/pyhope/mesh/topology/mesh_topology.py b/pyhope/mesh/topology/mesh_topology.py index c4398ec..1f0d82e 100644 --- a/pyhope/mesh/topology/mesh_topology.py +++ b/pyhope/mesh/topology/mesh_topology.py @@ -85,6 +85,9 @@ def MeshChangeElemType(mesh: meshio.Mesh) -> meshio.Mesh: hopout.info('Converting hexahedral elements to simplex elements') + # Temporally disable the garbage collector + gc.disable() + # Instantiate ELEMTYPE elemTypeInam = ELEMTYPE().inam @@ -300,6 +303,7 @@ def MeshChangeElemType(mesh: meshio.Mesh) -> meshio.Mesh: # Run garbage collector to release memory gc.collect() + gc.enable() hopout.sep() return mesh