pymrm.particles.GridParticle

class GridParticle(Particle):
    """Particle from level-function samples on its own body-frame grid.

    The samples are interpolated with a cubic B-spline
    (:class:`scipy.interpolate.RegularGridInterpolator`), so the particle can
    be translated and rotated for free.  The local grid **must extend beyond
    the particle surface** (positive samples all around); queries outside the
    local grid return the clamped boundary value plus the clamping distance,
    keeping the sign correct far away.

    Parameters
    ----------
    values : ndarray
        Level-function samples, negative inside.
    x_local : list of 1-D arrays
        Body-frame cell coordinates of the sample grid, one per axis.
    position, orientation : see :class:`Particle`.
    method : str, optional
        Interpolation method (default ``"cubic"``).
    """

    def __init__(self, values, x_local, position, orientation=None,
                 method="cubic"):
        position = np.asarray(position, dtype=float)
        super().__init__(position, orientation)
        from scipy.interpolate import RegularGridInterpolator
        values = np.asarray(values, dtype=float)
        x_local = [np.asarray(x, dtype=float) for x in x_local]
        if values.ndim != self.ndim or len(x_local) != self.ndim:
            raise ValueError("values/x_local dimensionality mismatch")
        self._lo = np.array([x[0] for x in x_local])
        self._hi = np.array([x[-1] for x in x_local])
        self._interp = RegularGridInterpolator(
            x_local, values, method=method, bounds_error=False, fill_value=None)
        self._eps = 1e-5 * float(np.max(self._hi - self._lo))

    def level_body(self, coords):
        coords = np.asarray(coords, dtype=float)
        clipped = np.clip(coords, self._lo, self._hi)
        excess = np.linalg.norm(coords - clipped, axis=-1)
        return self._interp(clipped) + excess

    def normal_body(self, coords):
        coords = np.asarray(coords, dtype=float)
        eps = self._eps
        g = np.empty(coords.shape)
        for a in range(self.ndim):
            dp = coords.copy(); dp[..., a] += eps
            dm = coords.copy(); dm[..., a] -= eps
            g[..., a] = (self.level_body(dp) - self.level_body(dm)) / (2 * eps)
        return g

    def bounding_box_body(self):
        return tuple((lo, hi) for lo, hi in zip(self._lo, self._hi))