[ENH] Refactor project_points_onto_surface for robustness and nonlinear projection support

docs/implementation/finite_faults.md

@@ -22,11 +22,11 @@ is not yet connected to the engine's older callable-based finite-fault hook.
 
 ### 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.
+- [x] Remove the incorrect fixed half-step from plane projection.
+- [x] Define explicit behavior for near-zero gradients.
+- [x] Test exact projection onto a plane and projection idempotency.
+- [x] Decide whether nonlinear fields need iterative re-evaluation.
+- [x] Fix the dense-grid gradient accessor before using it for projection.
 
 ### Phase 3: Stack wiring
 
@@ -117,14 +117,30 @@ Equivalent JSON:
 }
 ```
 
+## Projection Contract
+
+`project_points_onto_surface` performs one Newton step:
+
+```text
+P' = P - (F(P) - target) * grad(F(P)) / ||grad(F(P))||^2
+```
+
+This projection is exact for a linear scalar field. The function cannot iterate
+for a nonlinear field because its inputs contain scalar and gradient values only
+at the original points. Iterative projection, if required by integration tests,
+must re-evaluate the interpolator at each set of projected coordinates and will
+be implemented at the engine integration layer.
+
+A point whose gradient norm is below `gradient_tolerance` is left unchanged only
+when its scalar residual is within `surface_tolerance`. Otherwise projection is
+undefined and the function raises `ValueError`. Tolerances are keyword-only and
+must be non-negative.
+
 ## 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

gempy_engine/core/data/scalar_field_output.py

@@ -49,11 +49,16 @@ def scalar_field_at_sp(self):
     def exported_fields_dense_grid(self):
         slicer = self.grid.dense_grid_slice
         scalar_field = self.exported_fields.scalar_field[slicer]
-        gx_field = self.exported_fields.scalar_field[slicer]
-        gy_field = self.exported_fields.scalar_field[slicer]
-        gz_field = self.exported_fields.scalar_field[slicer]
-
-        return ExportedFields(scalar_field, gx_field, gy_field, gz_field)
+        gx_field = self.exported_fields.gx_field
+        gy_field = self.exported_fields.gy_field
+        gz_field = self.exported_fields.gz_field
+
+        return ExportedFields(
+            scalar_field,
+            None if gx_field is None else gx_field[slicer],
+            None if gy_field is None else gy_field[slicer],
+            None if gz_field is None else gz_field[slicer],
+        )
 
     @property
     def values_block_regular_grid(self):

gempy_engine/modules/faults/finite_faults.py

@@ -106,18 +106,47 @@ def get_radius_mask(val, radius):
     return d
 
 
-def project_points_onto_surface(points: np.ndarray, scalar_field_values: np.ndarray, gradient_fields: tuple[np.ndarray, np.ndarray, np.ndarray], target_scalar_value: float = 0.0) -> np.ndarray:
+def project_points_onto_surface(
+        points: np.ndarray,
+        scalar_field_values: np.ndarray,
+        gradient_fields: tuple[np.ndarray, np.ndarray, np.ndarray],
+        target_scalar_value: float = 0.0,
+        *,
+        gradient_tolerance: float = 1e-12,
+        surface_tolerance: float = 1e-12,
+) -> np.ndarray:
     """
-    Project points onto the surface F(x,y,z) = target_scalar_value.
-    Formula: P' = P - (F(P) - target) * grad(F) / ||grad(F)||^2
+    Apply one Newton projection step toward ``F(x, y, z) = target``.
+
+    The step is exact for a linear scalar field. Nonlinear fields require the
+    scalar field and gradient to be re-evaluated before applying another step.
+    A point away from the surface cannot be projected when its gradient is
+    effectively zero, so that case raises instead of silently leaving it in an
+    invalid position.
     """
+    points = np.asarray(points)
+    scalar_field_values = np.asarray(scalar_field_values)
     gx, gy, gz = gradient_fields
-    grad = np.stack([gx, gy, gz], axis=-1)
+    grad = np.stack([np.asarray(gx), np.asarray(gy), np.asarray(gz)], axis=-1)
+
+    if points.ndim != 2 or points.shape[1] != 3:
+        raise ValueError("points must have shape (n, 3)")
+    if scalar_field_values.shape != (len(points),) or grad.shape != points.shape:
+        raise ValueError("scalar and gradient fields must contain one value per point")
+    if gradient_tolerance < 0 or surface_tolerance < 0:
+        raise ValueError("projection tolerances must be non-negative")
+
     grad_norm_sq = np.sum(grad ** 2, axis=-1)
-    grad_norm_sq = np.where(grad_norm_sq < 1e-12, 1.0, grad_norm_sq)
-    f_p = scalar_field_values - target_scalar_value
-    projection = points - 0.5 * (f_p[:, np.newaxis] * grad) / grad_norm_sq[:, np.newaxis]
-    return projection
+    residual = scalar_field_values - target_scalar_value
+    near_zero_gradient = grad_norm_sq <= gradient_tolerance ** 2
+    unprojectable = near_zero_gradient & (np.abs(residual) > surface_tolerance)
+    if np.any(unprojectable):
+        count = int(np.count_nonzero(unprojectable))
+        raise ValueError(f"Cannot project {count} point(s) with near-zero scalar-field gradient")
+
+    safe_grad_norm_sq = np.where(near_zero_gradient, 1.0, grad_norm_sq)
+    correction = residual[:, np.newaxis] * grad / safe_grad_norm_sq[:, np.newaxis]
+    return points - correction
 
 
 def cubic_hermite_taper(d: np.ndarray) -> np.ndarray:

tests/test_common/test_api/test_faults/finite_fault_uv_design.md

@@ -25,8 +25,8 @@ To project the points we use the gradient of the scalar field of the fault surfa
 
 #### Phase 2: Point Projection ("Walking the Gradient")
 **Goal**: Project any 3D point $P$ onto the nearest point $P'$ on the fault surface $F(x,y,z)=c$.
-- **Task 2.1**: Implement the projection formula: $P' = P - 0.5 \cdot (F(P) - c) \frac{\nabla F(P)}{\|\nabla F(P)\|^2}$. 
-    - *Note*: We use a 0.5 factor because GemPy scalar fields typically behave quadratically ($F \approx d^2$) near the fault surface, so the gradient is twice as strong as a standard SDF gradient.
+- **Task 2.1**: Implement one Newton projection step: $P' = P - (F(P) - c) \frac{\nabla F(P)}{\|\nabla F(P)\|^2}$.
+    - *Note*: This is exact for a linear scalar field. Nonlinear fields require scalar and gradient re-evaluation before each additional step.
 - **Task 2.2**: Handle potential instabilities where $\|\nabla F\|$ is near zero.
 - **Tests**:
     - `test_projection_on_plane`: Verify projection works perfectly for a simple tilted plane.
@@ -84,6 +84,6 @@ $$Offset(d) = MaxSlip \times (1 - d^2)^2$$
 ### Step 3: Grid and other input points projection onto the fault surface
 1. The Projection Mechanism (Walking the Gradient)
 Your idea to interpolate the gradient of the fault's scalar field $F(x, y, z)$ is exactly the right path. The gradient $\nabla F$ acts as a vector field pointing perpendicularly toward/away from the fault surface.
-The Math: If your fault scalar field is a true (or approximate) Signed Distance Function (SDF), projecting a 3D point $P$ to its corresponding point on the fault surface $P'$ requires a single mathematical step:
-$$P' = P - F(P) \frac{\nabla F(P)}{\|\nabla F(P)\|}$$
-Note: If $F$ is not an SDF, a better approximation is $P' = P - 0.5 \cdot F(P) \frac{\nabla F(P)}{\|\nabla F(P)\|^2}$. The 0.5 factor accounts for quadratic behavior of GemPy's scalar field ($F \approx d^2$) where $\nabla F \approx 2d$.
+The Math: A Newton step projects a 3D point $P$ toward its corresponding point on the fault surface $P'$:
+$$P' = P - (F(P) - c) \frac{\nabla F(P)}{\|\nabla F(P)\|^2}$$
+For a signed distance field this reduces to a step of signed distance along the unit normal. For a general nonlinear field, evaluate $F$ and $\nabla F$ again at $P'$ before taking another step.

tests/test_common/test_api/test_faults/test_finite_fault_uv.py

@@ -10,6 +10,7 @@
 )
 from tests.conftest import plot_pyvista
 import numpy as np
+import pytest
 
 
 # --- Phase 1: Local Coordinate System & Analytical Ellipsoid UV ---
@@ -64,8 +65,6 @@ def test_ellipsoid_distance():
 def test_projection_on_plane():
     # Surface: F(x,y,z) = z - 5 = 0  => z = 5
     # grad(F) = [0, 0, 1]
-    # NOTE: Since we added a 0.5 factor to handle GemPy's quadratic scalar fields, 
-    # a single step on a linear field will only go halfway.
     points = np.array([
             [0, 0, 10.0],
             [1, 2, 0.0],
@@ -78,19 +77,50 @@ def test_projection_on_plane():
 
     projected = project_points_onto_surface(points, f_values, (gx, gy, gz), target_scalar_value=0.0)
 
-    # With 0.5 factor, it goes halfway:
-    # [0, 0, 10] -> [0, 0, 7.5]
-    # [1, 2, 0]  -> [1, 2, 2.5]
-    # [5, 5, 5]  -> [5, 5, 5.0]
     expected = np.array([
-            [0, 0, 7.5],
-            [1, 2, 2.5],
+            [0, 0, 5.0],
+            [1, 2, 5.0],
             [5, 5, 5.0]
     ])
 
     assert np.allclose(projected, expected)
 
 
+def test_projection_on_plane_is_idempotent():
+    points = np.array([[0.0, 0.0, 8.0], [1.0, 2.0, -2.0]])
+    gradients = (np.zeros(2), np.zeros(2), np.ones(2))
+
+    projected = project_points_onto_surface(points, points[:, 2] - 5.0, gradients)
+    projected_again = project_points_onto_surface(
+        projected,
+        projected[:, 2] - 5.0,
+        gradients,
+    )
+
+    assert np.allclose(projected_again, projected)
+
+
+def test_projection_rejects_near_zero_gradient_away_from_surface():
+    with pytest.raises(ValueError, match=r"Cannot project 1 point\(s\)"):
+        project_points_onto_surface(
+            points=np.array([[1.0, 2.0, 3.0]]),
+            scalar_field_values=np.array([1.0]),
+            gradient_fields=(np.zeros(1), np.zeros(1), np.zeros(1)),
+        )
+
+
+def test_projection_keeps_surface_point_with_near_zero_gradient():
+    points = np.array([[1.0, 2.0, 3.0]])
+
+    projected = project_points_onto_surface(
+        points=points,
+        scalar_field_values=np.array([0.0]),
+        gradient_fields=(np.zeros(1), np.zeros(1), np.zeros(1)),
+    )
+
+    assert np.array_equal(projected, points)
+
+
 # --- Phase 3: Slip Tapering Functions ---
 
 def test_taper_bounds():

tests/test_common/test_core/test_data_classes.py

@@ -3,7 +3,10 @@
 import numpy as np
 
 from gempy_engine.core.data import TensorsStructure
+from gempy_engine.core.data.exported_fields import ExportedFields
+from gempy_engine.core.data.scalar_field_output import ScalarFieldOutput
 from gempy_engine.core.data.solutions import Solutions
+from gempy_engine.core.data.stack_relation_type import StackRelationType
 
 
 def _make_stub_octree_level():
@@ -26,3 +29,44 @@ def test_solutions_repr_with_meshes():
 
     assert repr(solutions) == "Solutions(1 Octree Levels, 2 DualContouringMeshes)"
     assert solutions._repr_html_() == "<b>Solutions:</b> 1 Octree Levels, 2 DualContouringMeshes"
+
+
+def test_exported_fields_dense_grid_preserves_gradients():
+    exported_fields = ExportedFields(
+        _scalar_field=np.arange(5.0),
+        _gx_field=np.arange(10.0, 15.0),
+        _gy_field=np.arange(20.0, 25.0),
+        _gz_field=np.arange(30.0, 35.0),
+        _grid_size=5,
+    )
+    output = ScalarFieldOutput(
+        weights=None,
+        grid=SimpleNamespace(dense_grid_slice=slice(2, 5)),
+        exported_fields=exported_fields,
+        stack_relation=StackRelationType.ERODE,
+        values_block=None,
+    )
+
+    dense_fields = output.exported_fields_dense_grid
+
+    assert np.array_equal(dense_fields.scalar_field, np.arange(2.0, 5.0))
+    assert np.array_equal(dense_fields.gx_field, np.arange(12.0, 15.0))
+    assert np.array_equal(dense_fields.gy_field, np.arange(22.0, 25.0))
+    assert np.array_equal(dense_fields.gz_field, np.arange(32.0, 35.0))
+
+
+def test_exported_fields_dense_grid_preserves_missing_gradients():
+    exported_fields = ExportedFields(_scalar_field=np.arange(5.0), _grid_size=5)
+    output = ScalarFieldOutput(
+        weights=None,
+        grid=SimpleNamespace(dense_grid_slice=slice(2, 5)),
+        exported_fields=exported_fields,
+        stack_relation=StackRelationType.ERODE,
+        values_block=None,
+    )
+
+    dense_fields = output.exported_fields_dense_grid
+
+    assert dense_fields.gx_field is None
+    assert dense_fields.gy_field is None
+    assert dense_fields.gz_field is None