Mesh: Port the original Hilbert(Z)/Morton(Z) from HOPR

docs/user-guide/parameter-file.md

@@ -114,10 +114,11 @@ Mesh parameters controlling the respective mesh generation mode.
 
 | <div style="width:200px">Parameter</div> | <div style="width:90px">Type</div> | <div style="width:50px">Default</div> | <div style="width:200px">Allowed</div> | Explanation                            |
 | ---------------------------------------- | ---------------------------------- | ------------------------------------- | -------------------------------------- | ---------------------------------------- |
-| `Mode`                                   | int &#124; string                  | —                                    | `1` (internal) &#124; `3` (external)    | Mesh generation mode. `1` builds meshes using the internal mesh generator. `3` reads and converts meshes from external mesh generators. |
+| `Mode`                                   | int &#124; string                  | —                                     | `1` (internal) &#124; `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` &#124; `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`)

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:

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

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]<<n_bits | ...
+    packed = np.zeros(M, dtype=np.uint64)
+    for d in range(n_dims):
+        packed |= u[:, d] << np.uint64(d * n_bits)
+
+    if n_bits > 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

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

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