[ENH] Introduce FiniteFault model with validation, serialization, and spline taper support

docs/implementation/finite_faults.md

@@ -0,0 +1,135 @@
+# Finite Faults Implementation
+
+## Goal
+
+Define finite faults as validated, serializable input data and use that
+definition to taper fault displacement on the interpolated fault surface.
+
+The implementation is based on the projection and UV taper prototype in
+`gempy_engine/modules/faults/finite_faults.py`. The prototype is useful, but it
+is not yet connected to the engine's older callable-based finite-fault hook.
+
+## Status
+
+### Phase 1: Serializable input definition
+
+- [x] Introduce a Pydantic dataclass for the finite-fault definition.
+- [x] Keep persisted fields JSON-native rather than storing NumPy arrays.
+- [x] Serialize taper types as stable string values.
+- [x] Validate centers, radii, rotation, and spline control points.
+- [x] Add JSON and Python round-trip tests using `pydantic.TypeAdapter`.
+- [x] Preserve the prototype import path and numerical API.
+
+### Phase 2: Projection correctness
+
+- [ ] Remove the incorrect fixed half-step from plane projection.
+- [ ] Define explicit behavior for near-zero gradients.
+- [ ] Test exact projection onto a plane and projection idempotency.
+- [ ] Decide whether nonlinear fields need iterative re-evaluation.
+- [ ] Fix the dense-grid gradient accessor before using it for projection.
+
+### Phase 3: Stack wiring
+
+- [ ] Replace the callable in `FiniteFaultData` with the declarative definition.
+- [ ] Associate at most one finite-fault definition with each fault stack.
+- [ ] Validate that definitions are attached only to `StackRelationType.FAULT` stacks.
+- [ ] Make scalar gradients available when a finite-fault stack is evaluated.
+- [ ] Pass the fault scalar field, gradients, and surface isovalue to the taper operation.
+- [ ] Apply the taper to the fault drift before dependent stacks are interpolated.
+- [ ] Cover both sequential and flat-stack interpolation paths.
+
+### Phase 4: Integration and backend support
+
+- [ ] Add a numerical integration test for a dependent stratigraphic stack.
+- [ ] Verify that displacement reaches zero at the finite-fault tips.
+- [ ] Document the approximately planar local-frame limitation.
+- [ ] Define and test NumPy and PyTorch backend behavior.
+- [ ] Add the finite-fault definition to the server payload when stack data is exposed there.
+
+## Input Contract
+
+`FiniteFault` is a frozen Pydantic dataclass. Its persisted representation
+contains only JSON-compatible values:
+
+| Field | Type | Meaning |
+| --- | --- | --- |
+| `center` | 3-tuple of floats | Point at the center of the finite-fault footprint |
+| `strike_radius` | positive float or 2-tuple | Positive and negative strike radii |
+| `dip_radius` | positive float or 2-tuple | Positive and negative dip radii |
+| `taper` | `cubic`, `quadratic`, or `spline` | Slip taper profile |
+| `rotation_deg` | float | In-plane rotation in degrees |
+| `spline_control_points` | optional sequence of 2-tuples | Distance-to-slip profile for a spline taper |
+
+For an asymmetric radius, tuple order is `(positive_direction,
+negative_direction)`. All radii must be finite and greater than zero.
+
+Spline control points use `(normalized_distance, slip_multiplier)`. Distances
+must be strictly increasing from `0` to `1`, and multipliers must remain in the
+range `[0, 1]`. Omitting the points for a spline taper selects the engine's
+default profile. Supplying spline points for another taper is invalid.
+
+`normal_radius` is intentionally not part of this contract. The UV workflow
+projects points onto the fault surface, so normal distance is not part of its
+two-dimensional footprint. A volumetric ellipsoid would be a separate model.
+
+## Serialization
+
+Pydantic dataclasses use `TypeAdapter` for serialization and deserialization:
+
+```python
+from pydantic import TypeAdapter
+
+from gempy_engine.core.data.finite_fault import FiniteFault, TaperType
+
+adapter = TypeAdapter(FiniteFault)
+
+finite_fault = FiniteFault(
+    center=(0.0, 0.0, 0.0),
+    strike_radius=(2.0, 1.0),
+    dip_radius=0.75,
+    taper=TaperType.SPLINE,
+    rotation_deg=15.0,
+    spline_control_points=(
+        (0.0, 1.0),
+        (0.5, 0.8),
+        (1.0, 0.0),
+    ),
+)
+
+payload: bytes = adapter.dump_json(finite_fault)
+restored: FiniteFault = adapter.validate_json(payload)
+```
+
+Equivalent JSON:
+
+```json
+{
+  "center": [0.0, 0.0, 0.0],
+  "strike_radius": [2.0, 1.0],
+  "dip_radius": 0.75,
+  "taper": "spline",
+  "rotation_deg": 15.0,
+  "spline_control_points": [
+    [0.0, 1.0],
+    [0.5, 0.8],
+    [1.0, 0.0]
+  ]
+}
+```
+
+## Known Prototype Issues
+
+- `project_points_onto_surface` currently moves points only halfway to a linear plane.
+- Near-zero gradients are silently replaced with a denominator of one.
+- `FiniteFault.calculate_slip` expects callers to project points separately.
+- The local strike/dip frame is constant and therefore approximates curved faults.
+- The engine-integrated `FiniteFaultData` loses its callable when serialized.
+- `ScalarFieldOutput.exported_fields_dense_grid` currently returns scalar values as gradients.
+- Existing integration assertions do not verify the projected surface residual or expected geometry.
+
+## Design Decisions
+
+- The finite-fault definition is declarative; serialized callables are not supported.
+- NumPy conversion happens at the numerical boundary, not in persisted fields.
+- Geometry belongs to the fault stack that defines it, not each destination stack affected by it.
+- Projection algorithm settings are evaluation concerns and are not part of geological input data.

gempy_engine/core/data/__init__.py

@@ -1,3 +1,4 @@
+from .finite_fault import FiniteFault, TaperType
 from .tensors_structure import TensorsStructure
 from .kernel_classes.orientations import Orientations, OrientationsInternals
 from .kernel_classes.surface_points import SurfacePoints, SurfacePointsInternals

gempy_engine/core/data/finite_fault.py

@@ -0,0 +1,86 @@
+from __future__ import annotations
+
+import math
+from enum import Enum
+from typing import Any
+
+import numpy as np
+from pydantic import ConfigDict, field_validator, model_validator
+from pydantic.dataclasses import dataclass
+
+
+class TaperType(str, Enum):
+    CUBIC = "cubic"
+    QUADRATIC = "quadratic"
+    SPLINE = "spline"
+
+
+Radius = float | tuple[float, float]
+SplineControlPoints = tuple[tuple[float, float], ...]
+
+
+@dataclass(frozen=True, config=ConfigDict(extra="forbid"))
+class FiniteFault:
+    """Serializable definition of an approximately planar finite fault."""
+
+    center: tuple[float, float, float]
+    strike_radius: Radius = 1.0
+    dip_radius: Radius = 1.0
+    taper: TaperType = TaperType.CUBIC
+    rotation_deg: float = 0.0
+    spline_control_points: SplineControlPoints | None = None
+
+    @field_validator("center", "strike_radius", "dip_radius", "spline_control_points", mode="before")
+    @classmethod
+    def _convert_numpy_arrays(cls, value: Any) -> Any:
+        return value.tolist() if isinstance(value, np.ndarray) else value
+
+    @field_validator("center")
+    @classmethod
+    def _validate_center(cls, value: tuple[float, float, float]) -> tuple[float, float, float]:
+        if not all(math.isfinite(component) for component in value):
+            raise ValueError("center coordinates must be finite")
+        return value
+
+    @field_validator("strike_radius", "dip_radius")
+    @classmethod
+    def _validate_radius(cls, value: Radius) -> Radius:
+        radii = value if isinstance(value, tuple) else (value,)
+        if not all(math.isfinite(radius) and radius > 0 for radius in radii):
+            raise ValueError("radii must be finite and greater than zero")
+        return value
+
+    @field_validator("rotation_deg")
+    @classmethod
+    def _validate_rotation(cls, value: float) -> float:
+        if not math.isfinite(value):
+            raise ValueError("rotation_deg must be finite")
+        return value
+
+    @model_validator(mode="after")
+    def _validate_spline_control_points(self) -> "FiniteFault":
+        points = self.spline_control_points
+        if points is None:
+            return self
+        if self.taper is not TaperType.SPLINE:
+            raise ValueError("spline_control_points require a spline taper")
+        if len(points) < 2:
+            raise ValueError("spline_control_points require at least two points")
+
+        distances = tuple(point[0] for point in points)
+        multipliers = tuple(point[1] for point in points)
+        if distances[0] != 0.0 or distances[-1] != 1.0:
+            raise ValueError("spline distances must start at 0 and end at 1")
+        if not all(math.isfinite(value) for point in points for value in point):
+            raise ValueError("spline control points must be finite")
+        if not all(left < right for left, right in zip(distances, distances[1:])):
+            raise ValueError("spline distances must be strictly increasing")
+        if not all(0.0 <= multiplier <= 1.0 for multiplier in multipliers):
+            raise ValueError("spline multipliers must be between 0 and 1")
+        return self
+
+    def calculate_slip(self, points: np.ndarray, normal: np.ndarray) -> np.ndarray:
+        """Calculate the prototype slip multiplier for the supplied points."""
+        from gempy_engine.modules.faults.finite_faults import calculate_slip
+
+        return calculate_slip(self, points=points, normal=normal)

gempy_engine/modules/faults/finite_faults.py

@@ -1,73 +1,38 @@
 import numpy as np
 
-from enum import Enum
-from dataclasses import dataclass
-from typing import Optional, Union, Tuple
 from scipy.interpolate import make_interp_spline
 
-
-class TaperType(Enum):
-    CUBIC = "cubic"
-    QUADRATIC = "quadratic"
-    SPLINE = "spline"
-
-
-@dataclass
-class FiniteFault:
-    """
-    Elegant API for defining finite faults.
-
-    Args:
-        center: Center of the fault in 3D space.
-        strike_radius: Radius along the strike direction (u). 
-            Can be a single float or a tuple (positive_u, negative_u) for anisotropy.
-        dip_radius: Radius along the dip direction (v).
-            Can be a single float or a tuple (positive_v, negative_v) for anisotropy.
-        normal_radius: Radius along the normal direction (w). Default 1.0.
-        taper: The tapering function to use.
-        spline_control_points: If taper is SPLINE, these points define the curve.
-    """
-    center: np.ndarray
-    strike_radius: Union[float, Tuple[float, float]] = 1.0
-    dip_radius: Union[float, Tuple[float, float]] = 1.0
-    normal_radius: Optional[Union[float, Tuple[float, float]]] = None
-    taper: TaperType = TaperType.CUBIC
-    rotation: float = 0.0
-    spline_control_points: Optional[np.ndarray] = None
-
-    def __post_init__(self):
-        if self.taper == TaperType.SPLINE and self.spline_control_points is None:
-            # Default bell-shaped spline if none provided
-            self.spline_control_points = np.array([
-                    [0.0, 1.0],
-                    [0.2, 0.95],
-                    [0.5, 0.5],
-                    [0.8, 0.05],
-                    [1.0, 0.0]
-            ])
-
-    def calculate_slip(self, points: np.ndarray, normal: np.ndarray) -> np.ndarray:
-        """
-        High-level method to calculate slip multiplier for given points.
-        """
-        u, v, w = get_local_frame(normal, angle_deg=self.rotation)
-        d = get_ellipsoid_distance(
-            points=points,
-            center=self.center,
-            u=u, v=v, w=w if self.normal_radius is not None else None,
-            a=self.strike_radius,
-            b=self.dip_radius,
-            c=self.normal_radius if self.normal_radius is not None else 1.0
-        )
-
-        if self.taper == TaperType.CUBIC:
-            return cubic_hermite_taper(d)
-        elif self.taper == TaperType.QUADRATIC:
-            return quadratic_taper(d)
-        elif self.taper == TaperType.SPLINE:
-            return spline_taper(d, self.spline_control_points)
-        else:
-            raise ValueError(f"Unknown taper type: {self.taper}")
+from gempy_engine.core.data.finite_fault import FiniteFault, TaperType
+
+_DEFAULT_SPLINE_CONTROL_POINTS = np.array([
+        [0.0, 1.0],
+        [0.2, 0.95],
+        [0.5, 0.5],
+        [0.8, 0.05],
+        [1.0, 0.0]
+])
+
+
+def calculate_slip(finite_fault: FiniteFault, points: np.ndarray, normal: np.ndarray) -> np.ndarray:
+    """Calculate a slip multiplier using the prototype local-plane model."""
+    u, v, _ = get_local_frame(normal, angle_deg=finite_fault.rotation_deg)
+    distance = get_ellipsoid_distance(
+        points=points,
+        center=np.asarray(finite_fault.center),
+        u=u,
+        v=v,
+        a=finite_fault.strike_radius,
+        b=finite_fault.dip_radius,
+    )
+
+    if finite_fault.taper is TaperType.CUBIC:
+        return cubic_hermite_taper(distance)
+    if finite_fault.taper is TaperType.QUADRATIC:
+        return quadratic_taper(distance)
+    if finite_fault.taper is TaperType.SPLINE:
+        control_points = finite_fault.spline_control_points or _DEFAULT_SPLINE_CONTROL_POINTS
+        return spline_taper(distance, np.asarray(control_points))
+    raise ValueError(f"Unknown taper type: {finite_fault.taper}")
 
 
 def get_local_frame(normal: np.ndarray, angle_deg: float = 0) -> tuple[np.ndarray, np.ndarray, np.ndarray]:

tests/test_common/test_api/test_faults/test_finite_fault_api.py

@@ -84,7 +84,7 @@ def test_finite_fault_rotation():
         center=center,
         strike_radius=2.0,
         dip_radius=1.0,
-        rotation=90 # Rotate 90 deg
+        rotation_deg=90 # Rotate 90 deg
     )
 
     # Normal along Z => default strike u=[1,0,0], dip v=[0,1,0]

tests/test_common/test_api/test_faults/test_finite_fault_integration.py

@@ -68,9 +68,9 @@ def test_finite_fault_integration_suite(finite_fault_setup):
             FiniteFault(center=center, strike_radius=0.8, dip_radius=0.4, taper=TaperType.CUBIC),
             FiniteFault(center=center, strike_radius=(1.0, 0.5), dip_radius=0.6, taper=TaperType.QUADRATIC),
             FiniteFault(center=center, strike_radius=0.7, dip_radius=0.7, taper=TaperType.SPLINE, spline_control_points=cp_bell),
-            FiniteFault(center=center, strike_radius=0.8, dip_radius=0.4, taper=TaperType.CUBIC, rotation=45),
-            FiniteFault(center=center, strike_radius=0.8, dip_radius=1, taper=TaperType.SPLINE, spline_control_points=cp_plateau, rotation=30),
-            FiniteFault(center=center, strike_radius=(1.2,0.8), dip_radius=(2,1), taper=TaperType.SPLINE, spline_control_points=cp_plateau, rotation=30)
+            FiniteFault(center=center, strike_radius=0.8, dip_radius=0.4, taper=TaperType.CUBIC, rotation_deg=45),
+            FiniteFault(center=center, strike_radius=0.8, dip_radius=1, taper=TaperType.SPLINE, spline_control_points=cp_plateau, rotation_deg=30),
+            FiniteFault(center=center, strike_radius=(1.2,0.8), dip_radius=(2,1), taper=TaperType.SPLINE, spline_control_points=cp_plateau, rotation_deg=30)
     ]
 
     if not plot_pyvista: