[ENH] Add finite-fault support to server schema, validation, and PyTorch backend

docs/implementation/finite_faults.md

@@ -41,11 +41,11 @@ for Phase 4.
 
 ### 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.
+- [x] Add a numerical integration test for a dependent stratigraphic stack.
+- [x] Verify that displacement reaches zero at the finite-fault tips.
+- [x] Document the approximately planar local-frame limitation.
+- [x] Define and test NumPy and PyTorch backend behavior.
+- [x] Add the finite-fault definition to the server payload.
 
 ## Input Contract
 
@@ -172,11 +172,74 @@ use the non-stacked symbolic evaluator. Ordinary stacks continue to use the
 optimized stacked evaluator. NumPy and PyKeOps reductions can differ slightly
 because their floating-point reduction orders are different.
 
+## Server Payload
+
+The server descriptor accepts the same finite-fault definition under an entry
+aligned with the fault stack. The fault dependency matrix is required when
+other stacks consume that fault:
+
+```json
+{
+  "input_data_descriptor": {
+    "number_of_points_per_surface": [9, 12],
+    "number_of_points_per_stack": [9, 12],
+    "number_of_orientations_per_stack": [1, 2],
+    "number_of_surfaces_per_stack": [1, 1],
+    "masking_descriptor": [3, 1],
+    "faults_relations": [
+      [false, true],
+      [false, false]
+    ],
+    "faults_input_data": [
+      {
+        "thickness": null,
+        "finite_fault": {
+          "center": [0.0, 0.0, 0.0],
+          "strike_radius": [2.0, 1.0],
+          "dip_radius": 0.75,
+          "taper": "quadratic",
+          "rotation_deg": 15.0,
+          "spline_control_points": null
+        }
+      },
+      null
+    ]
+  }
+}
+```
+
+`masking_descriptor` uses the existing numeric `StackRelationType` values;
+`3` denotes `FAULT` and `1` denotes `ERODE`.
+
+## Geometry Limitation
+
+The current implementation creates one strike/dip frame for the whole finite
+fault. Its normal is taken from the valid projected evaluation point nearest the
+configured center. This is appropriate for planar and approximately planar
+faults. On a strongly curved fault, one fixed frame distorts distances and does
+not constitute a true surface UV parameterization. Supporting those faults will
+require a spatially varying frame or a mesh-based UV map.
+
+## Backend Behavior
+
+Projection and taper evaluation use NumPy and SciPy. NumPy inputs remain NumPy.
+PyTorch CPU or GPU tensors are detached and transferred to CPU for geometric
+evaluation, and the resulting multiplier is converted back to the active
+backend, device, and dtype before it is applied to the fault drift. Flat
+NumPy/PyKeOps evaluation is also supported through an isolated, non-stacked
+finite-fault evaluation path.
+
+The CPU conversion means projection and taper geometry are not differentiable
+with respect to the fault scalar field, gradients, center, radii, rotation, or
+spline points. A differentiable finite-fault workflow would require backend-
+native projection and taper functions. NumPy and PyTorch CPU behavior are covered
+by tests; GPU placement follows the same conversion path and is exercised only
+when CUDA is available.
+
 ## Known Prototype Issues
 
 - Standalone `FiniteFault.calculate_slip` expects callers to project points separately; engine stack wiring performs that projection.
 - The local strike/dip frame is constant and therefore approximates curved faults.
-- Existing integration assertions do not verify the projected surface residual or expected geometry.
 
 ## Design Decisions
 

gempy_engine/core/data/input_data_descriptor.py

@@ -51,11 +51,24 @@ def from_schema(cls, schema: InputDataDescriptorSchema):
 
         # Convert list of ints into list of StackRelationType
         list_relations: list[StackRelationType] = [StackRelationType(x) for x in schema.masking_descriptor]
+        faults_relations = None if schema.faults_relations is None else np.array(schema.faults_relations, dtype=bool)
+        faults_input_data = None
+        if schema.faults_input_data is not None:
+            faults_input_data = [
+                None if fault_data is None else FaultsData.from_user_input(
+                    thickness=fault_data.thickness,
+                    finite_fault=fault_data.finite_fault,
+                )
+                for fault_data in schema.faults_input_data
+            ]
+
         stack_structure = StacksStructure(
             number_of_points_per_stack=np.array(schema.number_of_points_per_stack),
             number_of_orientations_per_stack=np.array(schema.number_of_orientations_per_stack),
             number_of_surfaces_per_stack=np.array(schema.number_of_surfaces_per_stack),
-            masking_descriptor=list_relations
+            masking_descriptor=list_relations,
+            faults_relations=faults_relations,
+            faults_input_data=faults_input_data,
         )
         return cls(tensors_structure=tensor_structure, stack_structure=stack_structure)
 

gempy_engine/core/data/kernel_classes/server/input_parser.py

@@ -1,6 +1,8 @@
 # Check if pydantic is installed and import it
 from typing import Optional
 
+from ...finite_fault import FiniteFault
+
 try:
     from pydantic import BaseModel, Field
 except ImportError:
@@ -27,12 +29,19 @@ class InterpolationInputSchema(BaseModel):
     grid: GridSchema
 
 
+class FaultsDataSchema(BaseModel):
+    thickness: Optional[float] = None
+    finite_fault: Optional[FiniteFault] = None
+
+
 class InputDataDescriptorSchema(BaseModel):
     number_of_points_per_surface: list[int]
     number_of_points_per_stack: list[int]
     number_of_orientations_per_stack: list[int]
     number_of_surfaces_per_stack: list[int]
     masking_descriptor: list[int]  # * StackRelationType
+    faults_relations: Optional[list[list[bool]]] = None
+    faults_input_data: Optional[list[Optional[FaultsDataSchema]]] = None
 
 
 class GemPyInput(BaseModel):

tests/test_common/test_api/test_faults/test_finite_fault_stack_wiring.py

@@ -9,7 +9,8 @@
     _options_with_finite_fault_gradients,
 )
 from gempy_engine.API.interp_single._multi_scalar_field_manager import _compute_independent_chunks
-from gempy_engine.core.data import FiniteFault, InterpolationOptions
+from gempy_engine.config import AvailableBackends
+from gempy_engine.core.data import FiniteFault, InterpolationOptions, TaperType
 from gempy_engine.core.backend_tensor import BackendTensor
 from gempy_engine.core.data.exported_fields import ExportedFields
 from gempy_engine.core.data.kernel_classes.faults import FaultsData
@@ -18,7 +19,7 @@
 from gempy_engine.core.data.stacks_structure import StacksStructure
 
 
-def test_modify_fault_values_applies_projected_taper():
+def test_modify_fault_values_reaches_zero_at_fault_tip():
     finite_fault = FiniteFault(center=(0.0, 0.0, 0.0), strike_radius=1.0, dip_radius=1.0)
     fault_input = FaultsData.from_user_input(thickness=None, finite_fault=finite_fault)
     points = np.array([
@@ -47,6 +48,51 @@ def test_modify_fault_values_applies_projected_taper():
     assert np.allclose(tapered_values, [[0.0, 2.0, 0.0]])
 
 
+@pytest.mark.parametrize("use_gpu", [False, True])
+def test_modify_fault_values_preserves_pytorch_backend(use_gpu):
+    torch = pytest.importorskip("torch")
+    if use_gpu and not torch.cuda.is_available():
+        pytest.skip("CUDA is not available")
+    BackendTensor._change_backend(AvailableBackends.PYTORCH, use_gpu=use_gpu)
+    dtype = BackendTensor.dtype_obj
+    device = BackendTensor.device
+
+    finite_fault = FiniteFault(
+        center=(0.0, 0.0, 0.0),
+        strike_radius=1.0,
+        dip_radius=1.0,
+        taper=TaperType.SPLINE,
+    )
+    fault_input = FaultsData.from_user_input(thickness=None, finite_fault=finite_fault)
+    points = torch.tensor([
+            [2.0, 0.0, 2.0],
+            [0.0, 0.0, 2.0],
+            [1.0, 0.0, 2.0],
+    ], dtype=dtype, device=device)
+    exported_fields = ExportedFields(
+        _scalar_field=torch.full((3,), 2.0, dtype=dtype, device=device),
+        _gx_field=torch.zeros(3, dtype=dtype, device=device),
+        _gy_field=torch.zeros(3, dtype=dtype, device=device),
+        _gz_field=torch.ones(3, dtype=dtype, device=device),
+        _grid_size=3,
+        _scalar_field_at_surface_points=torch.tensor([0.0], dtype=dtype, device=device),
+    )
+    output = ScalarFieldOutput(
+        weights=None,
+        grid=None,
+        exported_fields=exported_fields,
+        stack_relation=StackRelationType.FAULT,
+        values_block=torch.tensor([[1.0, 3.0, 3.0]], dtype=dtype, device=device),
+    )
+
+    tapered_values = _modify_faults_values_output(fault_input, output, points)
+
+    assert isinstance(tapered_values, torch.Tensor)
+    assert tapered_values.dtype == dtype
+    assert tapered_values.device == points.device
+    assert torch.allclose(tapered_values, torch.tensor([[0.0, 2.0, 0.0]], dtype=dtype, device=device))
+
+
 def test_finite_fault_gradient_options_do_not_mutate_input():
     options = InterpolationOptions.from_args(range=1.0, c_o=1.0, mesh_extraction=False)
     fault_input = FaultsData.from_user_input(

tests/test_server/test_finite_fault_input.py

@@ -0,0 +1,50 @@
+from gempy_engine.core.data.finite_fault import TaperType
+from gempy_engine.core.data.input_data_descriptor import InputDataDescriptor
+from gempy_engine.core.data.kernel_classes.server.input_parser import InputDataDescriptorSchema
+
+
+def test_finite_fault_server_schema_converts_to_runtime_data():
+    schema = InputDataDescriptorSchema.model_validate({
+            "number_of_points_per_surface"      : [3, 3],
+            "number_of_points_per_stack"        : [3, 3],
+            "number_of_orientations_per_stack"  : [1, 1],
+            "number_of_surfaces_per_stack"      : [1, 1],
+            "masking_descriptor"                : [3, 1],
+            "faults_relations"                  : [[False, True], [False, False]],
+            "faults_input_data"                 : [
+                    {
+                        "finite_fault": {
+                            "center"       : [0.0, 0.0, 0.0],
+                            "strike_radius": [2.0, 1.0],
+                            "dip_radius"   : 0.75,
+                            "taper"        : "quadratic",
+                            "rotation_deg" : 15.0,
+                        }
+                    },
+                    None,
+            ],
+    })
+
+    descriptor = InputDataDescriptor.from_schema(schema)
+    fault_data = descriptor.stack_structure.faults_input_data[0]
+
+    assert descriptor.stack_structure.faults_relations.tolist() == [[False, True], [False, False]]
+    assert fault_data.finite_fault.center == (0.0, 0.0, 0.0)
+    assert fault_data.finite_fault.strike_radius == (2.0, 1.0)
+    assert fault_data.finite_fault.taper is TaperType.QUADRATIC
+
+
+def test_finite_fault_server_schema_serializes_as_json_values():
+    schema = InputDataDescriptorSchema.model_validate({
+            "number_of_points_per_surface"      : [3],
+            "number_of_points_per_stack"        : [3],
+            "number_of_orientations_per_stack"  : [1],
+            "number_of_surfaces_per_stack"      : [1],
+            "masking_descriptor"                : [3],
+            "faults_input_data"                 : [{"finite_fault": {"center": [1.0, 2.0, 3.0]}}],
+    })
+
+    payload = schema.model_dump(mode="json")
+
+    assert payload["faults_input_data"][0]["finite_fault"]["center"] == [1.0, 2.0, 3.0]
+    assert payload["faults_input_data"][0]["finite_fault"]["taper"] == "cubic"