diff --git a/docs/user-guide/parameter-file.md b/docs/user-guide/parameter-file.md index f6867ea0..b5cd5a0d 100644 --- a/docs/user-guide/parameter-file.md +++ b/docs/user-guide/parameter-file.md @@ -114,10 +114,11 @@ Mesh parameters controlling the respective mesh generation mode. |
Parameter
|
Type
|
Default
|
Allowed
| Explanation | | ---------------------------------------- | ---------------------------------- | ------------------------------------- | -------------------------------------- | ---------------------------------------- | -| `Mode` | int | string | — | `1` (internal) | `3` (external) | Mesh generation mode. `1` builds meshes using the internal mesh generator. `3` reads and converts meshes from external mesh generators. | +| `Mode` | int | string | — | `1` (internal) | `3` (external) | Mesh generation mode. `1` builds meshes using the internal mesh generator. `3` reads and converts meshes from external mesh generators. | | `NGeo` | int | `1` | ≥ 1 | Polynomial order used for representation of curved elements. Higher orders improve geometric fidelity. | | `BoundaryOrder` | int | `2` | ≥ 2 | Legacy parameter for polynomial order used for representation of curved elements. Prefer `NGeo` where applicable; kept for backward compatibility. | -| `MeshSorting` | sting | `SFC` | `SFC`, `IJK`, `LEX`, `Snake`, `None` | Mesh sorting mode to reorder the elements. | +| `MeshSorting` | string | `SFC` | `SFC`, `IJK`, `LEX`, `Snake`, `None` | Mesh sorting mode to reorder the elements. | +| `MeshSortingSFC` | string | `default` | `default` (PYPI hilbertcurve), `hilbert` (HOPR), `hilbertZ` (HOPR), `morton` (HOPR), `mortonZ` (HOPR) | The "default" space-filling curve is from the Python PYPI package and the others are C-ported from HOPR. | | `doSortIJK` | bool | `F` | `T` | `F` | Reorder elements primarily along I, then J, then K instead of the default space-filling Hilbert curve (legacy). | ### Internal Mesh Generator (Mode = `1` | `internal`) diff --git a/pyhope/mesh/mesh_sort.py b/pyhope/mesh/mesh_sort.py index 52904105..8c5faa82 100644 --- a/pyhope/mesh/mesh_sort.py +++ b/pyhope/mesh/mesh_sort.py @@ -140,57 +140,97 @@ def _set_bounding_box(self, mini, maxi): def SortMeshBySFC() -> None: # Local imports ---------------------------------------- from hilbertcurve.hilbertcurve import HilbertCurve + from pyhope.common.common_progress import ProgressBar from pyhope.common.common_vars import np_mtp from pyhope.mesh.mesh_common import calc_elem_bary - from pyhope.common.common_progress import ProgressBar + from pyhope.mesh.mesh_vars import MeshSortSFC + from pyhope.readintools.readintools import CreateIntFromString, CreateIntOption, GetIntFromStr import pyhope.mesh.mesh_vars as mesh_vars import pyhope.output.output as hopout + from pyhope.mesh.sort.sort_hilbert import hilbert, morton # Monkey-patching HilbertCurve from pyhope.mesh.sort.sort_hilbert import HilbertCurveNumpy # INFO: Alternative Hilbert curve sorting (not on PyPI) # from hilsort import hilbert_sort # ------------------------------------------------------ - # Monkey-patching HilbertCurve - HilbertCurveNumpy() + CreateIntFromString('MeshSortingSFC', default=MeshSortSFC.default.name, + help=f'Mesh sorting mode for SFC [{", ".join([s.name for s in MeshSortSFC])}]') # noqa: E501 + CreateIntOption( 'MeshSortingSFC', number=MeshSortSFC.default .value, name=MeshSortSFC.default .name) + CreateIntOption( 'MeshSortingSFC', number=MeshSortSFC.hilbert .value, name=MeshSortSFC.hilbert .name) + CreateIntOption( 'MeshSortingSFC', number=MeshSortSFC.hilbertZ.value, name=MeshSortSFC.hilbertZ.name) + CreateIntOption( 'MeshSortingSFC', number=MeshSortSFC.morton .value, name=MeshSortSFC.morton .name) + CreateIntOption( 'MeshSortingSFC', number=MeshSortSFC.mortonZ .value, name=MeshSortSFC.mortonZ .name) + + sfc_type: Final[int] = GetIntFromStr('MeshSortingSFC') + + hopout.sep() hopout.routine('Sorting elements along space-filling curve') - mesh = mesh_vars.mesh - elems = mesh_vars.elems - sides = mesh_vars.sides + mesh = mesh_vars.mesh + elems = mesh_vars.elems + sides = mesh_vars.sides + points = mesh.points # Use a moderate chunk size to bound intermediate progress updates chunk = max(1, min(1000, max(10, int(len(elems)/(400))))) bar = ProgressBar(value=len(elems), title='│ Preparing Elements', length=33, chunk=chunk) - # Global bounding box - points = mesh.points - xmin = points.min(axis=0) - xmax = points.max(axis=0) + match sfc_type: + case MeshSortSFC.default.value: + # Monkey-patching HilbertCurve + HilbertCurveNumpy() - # Calculate the element barycenters and associated element offsets - elem_bary = calc_elem_bary(elems) + # Global bounding box + xmin = points.min(axis=0) + xmax = points.max(axis=0) - # Calculate the space-filling curve resolution for the given KIND - kind: Final[int] = 4 - nbits, spacing = SFCResolution(kind, xmin, xmax) + # Calculate the element barycenters and associated element offsets + elem_bary = calc_elem_bary(elems) - # Discretize the element positions according to the chosen resolution - elem_disc = Coords2Int(elem_bary, spacing, xmin) + kind: Final[int] = 4 + nbits, spacing = SFCResolution(kind, xmin, xmax) + elem_disc = Coords2Int(elem_bary, spacing, xmin) - # Generate the space-filling curve and order elements along it - hc = HilbertCurve(p=nbits, n=3, n_procs=np_mtp) + hc = HilbertCurve(p=nbits, n=3, n_procs=np_mtp) + distances = cast(npt.ArrayLike, hc.distances_from_points(elem_disc)) - distances = cast(npt.ArrayLike, hc.distances_from_points(elem_disc)) # bottleneck - sorted_indices = np.argsort(distances) + case MeshSortSFC.hilbert .value | \ + MeshSortSFC.hilbertZ.value | \ + MeshSortSFC.morton .value | \ + MeshSortSFC.mortonZ .value: - # INFO: Alternative Hilbert curve sorting (not on PyPI) - # distances = np.array(hilbert_sort(8, elem_bary)) - # Find the new sorting with the old elem_bary - # value_to_index = {tuple(value.tolist()): idx for idx, value in enumerate(distances)} - # Now, create an array that maps each element to the new sorting - # sorted_indices = np.array([value_to_index[tuple(val.tolist())] for val in elem_bary]) + # Calculate the element barycenters and associated element offsets + elem_bary = calc_elem_bary(elems) + + gmin = elem_bary.min() # scalar global min + gmax = elem_bary.max() # scalar global max + + # nbits = (64-1)//3 = 21, matching HOPR setBoundingBox with INTEGER(KIND=8) + nbits = (np.iinfo(np.int64).bits - 1) // 3 + intfact = (1 << nbits) - 1 + spacing = np.float64(intfact) / (gmax - gmin) + elem_disc = Coords2Int(elem_bary, spacing, gmin) + + match sfc_type: + case MeshSortSFC.hilbert.value: + distances = hilbert(3, nbits, elem_disc) + + case MeshSortSFC.hilbertZ.value: + # Hilbert on (x,y), linear stride along z + distances = hilbert(2, nbits, elem_disc[:, :2]) + elem_disc[:, 2] * intfact * intfact + + case MeshSortSFC.morton.value: + distances = morton( 3, nbits, elem_disc) + + case MeshSortSFC.mortonZ.value: + # Morton on (x,y), linear stride along z + distances = morton( 2, nbits, elem_disc[:, :2]) + elem_disc[:, 2] * intfact * intfact + case _: + raise ValueError(f'Unknown sfc_type "{sfc_type}" (valid: [{", ".join([s.name for s in MeshSortSFC])}])') + + sorted_indices = np.argsort(distances) sorted_elems, sorted_sides = UpdateElemID(elems, sides, sorted_indices, bar) mesh_vars.elems = sorted_elems @@ -208,6 +248,7 @@ def SortMeshByIJK() -> None: from pyhope.common.common_progress import ProgressBar # ------------------------------------------------------ + hopout.sep() hopout.routine('Sorting elements along I,J,K direction') mesh = mesh_vars.mesh @@ -420,8 +461,6 @@ def SortMesh() -> None: meshsort = GetIntFromStr('MeshSorting', default=default) - hopout.sep() - # Sort the mesh match meshsort: case MeshSort.NONE.value: diff --git a/pyhope/mesh/mesh_vars.py b/pyhope/mesh/mesh_vars.py index 934e4204..f6e9496b 100644 --- a/pyhope/mesh/mesh_vars.py +++ b/pyhope/mesh/mesh_vars.py @@ -103,6 +103,15 @@ class MeshSort(Enum): Snake = 4 +@unique +class MeshSortSFC(Enum): + default = 0 + hilbert = 1 + hilbertZ = 2 + morton = 3 + mortonZ = 4 + + @dataclass(init=False, repr=False, eq=False, slots=False) class CGNS: regenerate_BCs: bool = False # Flag if CGNS needs BC regeneration diff --git a/pyhope/mesh/sort/sort_hilbert.py b/pyhope/mesh/sort/sort_hilbert.py index 9a6367b2..9b7b65e6 100644 --- a/pyhope/mesh/sort/sort_hilbert.py +++ b/pyhope/mesh/sort/sort_hilbert.py @@ -8,6 +8,23 @@ # Copyright (c) 2024 Numerics Research Group, University of Stuttgart, Prof. Andrea Beck # Copyright (c) 2022 Gabriel Altay (Original Version) # +# Hilbert(Z) curve adapted from +# hilbert.c - Computes Hilbert space-filling curve coordinates, without recursion, from integer index, and vice versa, and other +# Hilbert-related calculations. Also known as Pi-order or Peano scan. +# +# Author: Doug Moore +# Dept. of Computational and Applied Math +# Rice University +# http://www.caam.rice.edu/~dougm +# Date: Sun Feb 20 2000 +# Copyright (c) 1998-2000, Rice University +# +# Acknowledgement: +# This implementation is based on the work of A. R. Butz ("Alternative Algorithm for Hilbert's Space-Filling Curve", IEEE Trans. +# Comp., April, 1971, pp 424-426) and its interpretation by Spencer W. Thomas, University of Michigan +# (http://www-personal.umich.edu/~spencer/Home.html) in his widely available C software. While the implementation here differs +# considerably from his, the first two interfaces and the style of some comments are very much derived from his work. +# # PyHOPE is free software: you can redistribute it and/or modify it under the # terms of the GNU General Public License as published by the Free Software # Foundation, either version 3 of the License, or (at your option) any later @@ -96,10 +113,13 @@ def _distances_from_points_numpy(self, q = np.uint64(1) << np.uint64(pbits - 1) # m while q > 1: pmask = q - 1 + for i in range(n): mask = (upts[:, i] & q) != 0 + if mask.any(): upts[mask, 0] ^= pmask + if (~mask).any(): t = (upts[~mask, 0] ^ upts[~mask, i]) & pmask upts[~mask, 0] ^= t @@ -111,12 +131,13 @@ def _distances_from_points_numpy(self, upts[:, i] ^= upts[:, i - 1] tmask = np.zeros(M, dtype=np.uint64) - q = np.uint64(1) << np.uint64(pbits - 1) + q = np.uint64(1) << np.uint64(pbits - 1) while q > 1: mask = (upts[:, n - 1] & q) != 0 if mask.any(): tmask[mask] ^= (q - 1) q >>= 1 + upts ^= tmask[:, None] # Interleave bit-planes into a big (potentially >64-bit) integer per row @@ -168,3 +189,144 @@ def _dfp_patched(self, points: Iterable[Iterable[int]], match_type: bool = False HilbertCurve.distances_from_points = _dfp_patched HilbertCurve._numpy_patch_applied = True + + +def _bit_transpose(n_dims: int, n_bits: int, coords: npt.NDArray) -> npt.NDArray: + """ Transpose the bit-planes of a packed coordinate word + """ + in_b = n_bits + in_field_ends = np.uint64(1) + in_mask = np.uint64((1 << in_b) - 1) + result = np.zeros_like(coords) + + while (ut_b := in_b >> 1): + shift_amt = np.uint64((n_dims - 1) * ut_b) + ut_field_ends = in_field_ends | (in_field_ends << (shift_amt + ut_b)) + ut_mask = (ut_field_ends << np.uint64(ut_b)) - ut_field_ends + ut_coords = np.zeros_like(coords) + + if in_b & 1: + # Odd number of bits: peel off the top bit of each field separately + in_field_starts = in_field_ends << np.uint64(in_b - 1) + odd_shift = np.uint64(2 * shift_amt) + for d in range(n_dims): + chunk = coords & in_mask + coords >>= np.uint64(in_b) + result |= (chunk & in_field_starts) << odd_shift + odd_shift += np.uint64(1) + chunk &= ~in_field_starts + chunk = (chunk | (chunk << shift_amt)) & ut_mask + ut_coords |= chunk << np.uint64(d * ut_b) + else: + for d in range(n_dims): + chunk = coords & in_mask + coords >>= np.uint64(in_b) + chunk = (chunk | (chunk << shift_amt)) & ut_mask + ut_coords |= chunk << np.uint64(d * ut_b) + + coords = ut_coords + in_b = ut_b + in_field_ends = ut_field_ends + in_mask = ut_mask + + return result | coords + + +def _adjust_rotation(n_dims : int , bits : npt.NDArray, + rotation: npt.NDArray, nd1_ones: int) -> npt.NDArray: + """ rotation = (rotation + 1 + ffs(bits)) % nDims + """ + nd1 = np.uint64(nd1_ones) + b = bits & (-bits.astype(np.int64)).astype(np.uint64) & nd1 + # Count trailing zeros of b per element: each '>>1' step increments rotation + while np.any(b): + active = b != 0 + rotation = np.where(active, rotation + np.uint64(1), rotation) + b = np.where(active, b >> np.uint64(1), b) + rotation += np.uint64(1) + rotation = np.where(rotation >= np.uint64(n_dims), rotation - np.uint64(n_dims), rotation) + + return rotation + + +# ---------------------------------------------------------------------------------------------------------------------------------- +# Hilbert (Z-order) curve +# ---------------------------------------------------------------------------------------------------------------------------------- +def hilbert(n_dims: int, n_bits: int, + coords: npt.NDArray[np.int64]) -> npt.NDArray[np.int64]: + """ Convert integer coordinates to Hilbert curve indices + """ + M = coords.shape[0] + u = coords.astype(np.uint64) + + if n_dims == 1: # pragma: no cover + return u[:, 0].astype(np.int64) + + n_dims_bits = n_dims * n_bits + nd_ones = np.uint64(( 1 << n_dims ) - 1) + nth_bits = np.uint64(((1 << n_dims_bits) - 1) // int(nd_ones)) + nd1_ones = int(nd_ones >> np.uint64(1)) + + # Pack coordinates into a single word: packed[m] = u[m,0] | u[m,1]< 1: + # Bit-transpose then Gray-decode + packed = _bit_transpose(n_dims, n_bits, packed.copy()) + packed ^= packed >> np.uint64(n_dims) + + rotation = np.zeros(M, dtype=np.uint64) + flip_bit = np.zeros(M, dtype=np.uint64) + result = np.zeros(M, dtype=np.uint64) + + b = n_dims_bits + while b > 0: + b -= n_dims + bits = (packed >> np.uint64(b)) & nd_ones + bits = (flip_bit ^ bits) + + # rotateRight(bits, rotation, n_dims) + rot_r = rotation % np.uint64(n_dims) + bits = ((bits >> rot_r) | (bits << (np.uint64(n_dims) - rot_r))) & nd_ones + + result <<= np.uint64(n_dims) + result |= bits + flip_bit = np.uint64(1) << rotation + + rotation = _adjust_rotation( n_dims, bits.copy(), rotation, nd1_ones) + + result ^= nth_bits >> np.uint64(1) + + else: + # n_bits == 1: trivial Gray decode + result = packed + + # Final Gray decode of the index itself + d = 1 + while d < n_dims_bits: + result ^= result >> np.uint64(d) + d *= 2 + + return result.astype(np.int64) + + +# ---------------------------------------------------------------------------------------------------------------------------------- +# Morton (Z-order) curve +# ---------------------------------------------------------------------------------------------------------------------------------- +def morton(n_dims: int, n_bits: int, + coords: npt.NDArray[np.int64]) -> npt.NDArray[np.int64]: + """ Convert integer coordinates to Morton (Z-order) indices + """ + M = coords.shape[0] + result = np.zeros(M, dtype=np.int64) + + for d in range(n_dims): + col = coords[:, d].astype(np.int64) + for i in range(n_bits): + # Bit i of dimension d maps to position n_dims*i + (n_dims-1-d) + bit_pos = n_dims * i + (n_dims - 1 - d) + result |= ((col >> np.int64(i)) & np.int64(1)) << np.int64(bit_pos) + + return result diff --git a/tutorials/2-06-external_mesh_Gmsh_v4_3D/analyze.toml b/tutorials/2-06-external_mesh_Gmsh_v4_3D/analyze.toml index 5b9160a5..bf142f49 100644 --- a/tutorials/2-06-external_mesh_Gmsh_v4_3D/analyze.toml +++ b/tutorials/2-06-external_mesh_Gmsh_v4_3D/analyze.toml @@ -17,7 +17,7 @@ mean = 0.00659394280873397 stddev = 0.015885236562431546 [SideInfo] -min = -62317.0 +min = -62318.0 max = 62319.0 -mean = 2174.244025661588 -stddev = 16744.4540292264 +mean = 2192.275235565357 +stddev = 16691.18003373148 diff --git a/tutorials/2-06-external_mesh_Gmsh_v4_3D/parameter.ini b/tutorials/2-06-external_mesh_Gmsh_v4_3D/parameter.ini index 33ec01f1..b75d726a 100644 --- a/tutorials/2-06-external_mesh_Gmsh_v4_3D/parameter.ini +++ b/tutorials/2-06-external_mesh_Gmsh_v4_3D/parameter.ini @@ -9,9 +9,11 @@ OutputFormat = HDF5 ! Mesh output format (HDF5 V !=============================================================================== ! ! MESH !=============================================================================== ! -FileName = 2-06-external_mesh_Gmsh_v4_3D.msh ! Name of external mesh file -Mode = 3 ! 1 Cartesian 3 External -MeshScale = 0.001 +FileName = 2-06-external_mesh_Gmsh_v4_3D.msh ! Name of external mesh file +Mode = 3 ! 1 Cartesian 3 External +MeshScale = 0.001 +MeshSorting = SFC +MeshSortingSFC = hilbert !=============================================================================== ! ! BOUNDARY CONDITIONS diff --git a/tutorials/2-14-external_mesh_Gmsh_curved_mixed/analyze.toml b/tutorials/2-14-external_mesh_Gmsh_curved_mixed/analyze.toml index 658d19f9..619c964a 100644 --- a/tutorials/2-14-external_mesh_Gmsh_curved_mixed/analyze.toml +++ b/tutorials/2-14-external_mesh_Gmsh_curved_mixed/analyze.toml @@ -1,23 +1,23 @@ [ElemInfo] min = 0.0 max = 129114.0 -mean = 26323.34916683952 -stddev = 35350.49825262569 +mean = 26354.91592339072 +stddev = 35389.12703130935 [GlobalNodeIDs] min = 1.0 max = 57639.0 mean = 28734.790309339034 -stddev = 16597.10203326578 +stddev = 16597.102033265783 [NodeCoords] min = -10.0 max = 10.0 -mean = 0.8225186730641619 -stddev = 1.7610539123328763 +mean = 0.8225186730641623 +stddev = 1.7610539123328766 [SideInfo] -min = -14484.0 +min = -14485.0 max = 14486.0 -mean = 496.05026207088065 -stddev = 3910.450564423621 +mean = 495.9934326078656 +stddev = 3909.2393955719863 diff --git a/tutorials/2-14-external_mesh_Gmsh_curved_mixed/parameter.ini b/tutorials/2-14-external_mesh_Gmsh_curved_mixed/parameter.ini index 27627f29..a3dac594 100644 --- a/tutorials/2-14-external_mesh_Gmsh_curved_mixed/parameter.ini +++ b/tutorials/2-14-external_mesh_Gmsh_curved_mixed/parameter.ini @@ -14,7 +14,8 @@ Mode = 3 ! 1 Cartesian 3 External NGeo = 2 ! Polynomial order of the mapping MeshIsAlreadyCurved = T ! Enable curving FileName = ./2-14-external_mesh_Gmsh_curved_mixed.msh - +MeshSorting = SFC +MeshSortingSFC = hilbertZ CheckWatertightness = T ! Check watertighness of the mesh !=============================================================================== ! diff --git a/tutorials/2-15-external_mesh_Gmsh_extrude/analyze.toml b/tutorials/2-15-external_mesh_Gmsh_extrude/analyze.toml index 8b6669e2..7281296c 100644 --- a/tutorials/2-15-external_mesh_Gmsh_extrude/analyze.toml +++ b/tutorials/2-15-external_mesh_Gmsh_extrude/analyze.toml @@ -1,8 +1,8 @@ [ElemInfo] min = 0.0 max = 208.0 -mean = 65.16666666666667 -stddev = 73.51587130227956 +mean = 63.5 +stddev = 73.43572397379599 [GlobalNodeIDs] min = 1.0 @@ -13,11 +13,11 @@ stddev = 22.240503591420765 [NodeCoords] min = 0.0 max = 2.0 -mean = 0.5888888888888288 +mean = 0.588888888888829 stddev = 0.47218954135280516 [SideInfo] -min = -22.0 +min = -20.0 max = 61.0 -mean = 6.848484848484849 -stddev = 13.52107493512574 +mean = 6.872727272727273 +stddev = 13.465180072326001 diff --git a/tutorials/2-15-external_mesh_Gmsh_extrude/parameter.ini b/tutorials/2-15-external_mesh_Gmsh_extrude/parameter.ini index 79974c0f..d3def5e0 100644 --- a/tutorials/2-15-external_mesh_Gmsh_extrude/parameter.ini +++ b/tutorials/2-15-external_mesh_Gmsh_extrude/parameter.ini @@ -12,6 +12,8 @@ OutputFormat = HDF5 ! Mesh output format (HDF5 Mode = 3 ! 1 Cartesian 3 External NGeo = 2 ! Polynomial order of the mapping FileName = ./2-15-external_mesh_Gmsh_extrude.geo +MeshSorting = SFC +MeshSortingSFC = mortonZ !=============================================================================== ! ! EXTRUSION diff --git a/tutorials/3-01-agglomeration_CGNS_NACA/analyze.toml b/tutorials/3-01-agglomeration_CGNS_NACA/analyze.toml index c7768380..34875897 100644 --- a/tutorials/3-01-agglomeration_CGNS_NACA/analyze.toml +++ b/tutorials/3-01-agglomeration_CGNS_NACA/analyze.toml @@ -17,7 +17,7 @@ mean = 748.6569646374462 stddev = 2397.4014609017318 [SideInfo] -min = -7209.0 +min = -7208.0 max = 7210.0 -mean = 253.97320011337868 -stddev = 1940.0479806901267 +mean = 254.41757369614513 +stddev = 1936.9666264837401 diff --git a/tutorials/3-01-agglomeration_CGNS_NACA/parameter.ini b/tutorials/3-01-agglomeration_CGNS_NACA/parameter.ini index ce2f420b..4bc88e02 100644 --- a/tutorials/3-01-agglomeration_CGNS_NACA/parameter.ini +++ b/tutorials/3-01-agglomeration_CGNS_NACA/parameter.ini @@ -9,9 +9,10 @@ OutputFormat = HDF5 ! Mesh output format (HDF5 V !=============================================================================== ! ! MESH !=============================================================================== ! -FileName = 3-01-agglomeration_CGNS_NACA.cgns ! Name of mesh file -Mode = 3 ! 1 Cartesian 3 External -! MeshScale = 0.001 ! Scale all input meshes by factor +FileName = 3-01-agglomeration_CGNS_NACA.cgns ! Name of mesh file +Mode = 3 ! 1 Cartesian 3 External +MeshSorting = SFC +MeshSortingSFC = morton !=============================================================================== ! ! CURVING BY AGGLOMERATION