diff --git a/loopstructural/gui/dlg_settings.ui b/loopstructural/gui/dlg_settings.ui
index 998b597..477f844 100644
--- a/loopstructural/gui/dlg_settings.ui
+++ b/loopstructural/gui/dlg_settings.ui
@@ -274,7 +274,11 @@
-
-
+
+
+ 1000000
+
+
diff --git a/loopstructural/gui/modelling/fault_adjacency_tab.py b/loopstructural/gui/modelling/fault_adjacency_tab.py
index 6478f42..ffc910a 100644
--- a/loopstructural/gui/modelling/fault_adjacency_tab.py
+++ b/loopstructural/gui/modelling/fault_adjacency_tab.py
@@ -80,7 +80,7 @@ def _update(self, observable, event, *args, **kwargs):
self.update_fault_adjacency_table()
self.update_stratigraphic_units_table()
- def change_button_color(self, button, row, col):
+ def change_button_color(self, button, fault1, fault2):
"""Cycle the button color and update the fault relationship."""
current_color = button.styleSheet()
if "red" in current_color:
@@ -94,15 +94,23 @@ def change_button_color(self, button, row, col):
relationship = FaultRelationshipType.ABUTTING
button.setStyleSheet(f"background-color: {new_color};")
- f1 = self.data_manager._fault_topology.faults[row]
- f2 = self.data_manager._fault_topology.faults[col]
- self.data_manager._fault_topology.update_fault_relationship(f1, f2, relationship)
+ self.data_manager._fault_topology.update_fault_relationship(fault1, fault2, relationship)
+
+ def _displacement_fault_names(self):
+ """Fault names to show in these tables, excluding faults used as
+ stratigraphic-column domain boundaries (see `set_fault_boundary`).
+ Those are non-displacing splits, not faults that cut/abut other
+ faults or offset stratigraphic units, so FAULTED/ABUTTING and
+ fault-stratigraphy relationships don't apply to them.
+ """
+ domain_boundary_faults = self.data_manager.get_fault_boundary_fault_names()
+ return [
+ f for f in self.data_manager._fault_topology.faults if f not in domain_boundary_faults
+ ]
def update_fault_adjacency_table(self):
"""Update the fault adjacency table with QPushButtons."""
- faults = (
- self.data_manager._fault_topology.faults
- ) # Assuming faults is a list of fault names
+ faults = self._displacement_fault_names()
if not faults:
self.fault_table_group.hide()
return
@@ -145,15 +153,15 @@ def update_fault_adjacency_table(self):
else:
button.setStyleSheet("background-color: white;")
button.clicked.connect(
- lambda _, b=button, r=row, c=col: self.change_button_color(b, r, c)
+ lambda _, b=button, f1=faults[row], f2=faults[
+ col
+ ]: self.change_button_color(b, f1, f2)
)
self.table.setCellWidget(row, col, button)
def update_stratigraphic_units_table(self):
"""Update the stratigraphic units table with QPushButtons."""
- faults = (
- self.data_manager._fault_topology.faults
- ) # Assuming faults is a list of fault names
+ faults = self._displacement_fault_names()
group_units_pairs = self.data_manager._stratigraphic_column.get_group_unit_pairs()
if not faults or not group_units_pairs:
@@ -185,11 +193,13 @@ def update_stratigraphic_units_table(self):
# Default to white if no relationship or not faulted
button.setStyleSheet("background-color: white;")
button.clicked.connect(
- lambda _, b=button, r=row, c=col: self.change_button_colour_binary(b, r, c)
+ lambda _, b=button, u=units[row], f=faults[
+ col
+ ]: self.change_button_colour_binary(b, u, f)
)
self.stratigraphic_table.setCellWidget(row, col, button)
- def change_button_colour_binary(self, button, row, col):
+ def change_button_colour_binary(self, button, unit_name, fault_name):
"""Cycle the button color between red, green, and black."""
current_color = button.styleSheet()
@@ -199,8 +209,6 @@ def change_button_colour_binary(self, button, row, col):
else:
button.setStyleSheet("background-color: red;")
flag = True
- fault = self.data_manager._fault_topology.faults[col]
- unit = self.data_manager._stratigraphic_column.get_group_unit_pairs()[row]
self.data_manager._fault_topology.update_fault_stratigraphy_relationship(
- unit[1], fault, flag
+ unit_name, fault_name, flag
)
diff --git a/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py b/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py
index 13ac391..7dd559b 100644
--- a/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py
+++ b/loopstructural/gui/modelling/geological_model_tab/geological_model_tab.py
@@ -20,6 +20,7 @@
from .add_foliation_dialog import AddFoliationDialog
from .add_unconformity_dialog import AddUnconformityDialog
from .feature_details_panel import (
+ BaseFeatureDetailsPanel,
FaultFeatureDetailsPanel,
FoldedFeatureDetailsPanel,
FoliationFeatureDetailsPanel,
@@ -461,6 +462,17 @@ def on_feature_selected(self, item):
self.featureDetailsPanel = FoldedFeatureDetailsPanel(
feature=feature, model_manager=self.model_manager, data_manager=self.data_manager
)
+ elif feature.type == FeatureType.DOMAINFAULT:
+ # A domain fault is built by the same GeologicalFeatureBuilder
+ # as a foliation (see create_and_add_domain_fault), just with a
+ # different .type tag -- the generic base panel (interpolator
+ # settings, data layers, export/evaluate) already applies to it
+ # unchanged. Skip FoliationFeatureDetailsPanel's fold-frame
+ # attachment controls, which don't make sense for a domain
+ # boundary.
+ self.featureDetailsPanel = BaseFeatureDetailsPanel(
+ feature=feature, model_manager=self.model_manager, data_manager=self.data_manager
+ )
else:
self.featureDetailsPanel = QWidget() # Default empty panel
diff --git a/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_column.py b/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_column.py
index 0f1131b..27c66d6 100644
--- a/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_column.py
+++ b/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_column.py
@@ -232,7 +232,11 @@ def update_display(self):
widget, _ = self._widget_cache[unit.uuid]
# Update widget data without rebuilding
if hasattr(widget, 'setData'):
- widget.setData(unit.to_dict())
+ unit_data = unit.to_dict()
+ if isinstance(widget, UnconformityWidget):
+ unit_data = self._enrich_unconformity_data(unit_data)
+ widget.set_available_faults(self._get_available_fault_names())
+ widget.setData(unit_data)
return
# If order/content differs, do a full rebuild
@@ -263,7 +267,31 @@ def _full_rebuild_display(self, current_order):
if unit.element_type == StratigraphicColumnElementType.UNIT:
self.add_unit(unit_data=unit.to_dict(), create_new=False)
elif unit.element_type == StratigraphicColumnElementType.UNCONFORMITY:
- self.add_unconformity(unconformity_data=unit.to_dict(), create_new=False)
+ self.add_unconformity(
+ unconformity_data=self._enrich_unconformity_data(unit.to_dict()),
+ create_new=False,
+ )
+
+ def _enrich_unconformity_data(self, unconformity_data):
+ """Merge in the plugin-side fault-boundary link for an unconformity row.
+
+ `StratigraphicUnconformity.to_dict()` (core) only knows `erode`/
+ `onlap`; the fault link is tracked separately in the data manager
+ (see `ModellingDataManager.set_fault_boundary`), so it has to be
+ folded in here for display.
+ """
+ fault_name = self.data_manager.get_fault_boundary(unconformity_data.get('uuid'))
+ if fault_name:
+ unconformity_data = dict(unconformity_data)
+ unconformity_data['unconformity_type'] = 'fault'
+ unconformity_data['fault_name'] = fault_name
+ return unconformity_data
+
+ def _get_available_fault_names(self):
+ """Fault names offered when marking an unconformity as a domain boundary."""
+ if not self.data_manager:
+ return []
+ return list(self.data_manager._fault_topology.faults)
def init_stratigraphic_column_from_basal_contacts(self):
if self.data_manager:
@@ -482,11 +510,13 @@ def add_unconformity(self, *, unconformity_data=None, create_new=True):
widget, _ = self._widget_cache[unconformity.uuid]
# Just update the data, don't recreate the widget
if hasattr(widget, 'setData'):
+ widget.set_available_faults(self._get_available_fault_names())
widget.setData(unconformity_data)
return
unconformity_widget = UnconformityWidget(uuid=unconformity.uuid)
unconformity_widget.deleteRequested.connect(self.delete_unit)
+ unconformity_widget.dataChanged.connect(lambda: self.update_element(unconformity_widget))
unconformity_widget.dragHandlePressed.connect(
lambda: self._on_drag_start(unconformity_widget)
)
@@ -500,6 +530,8 @@ def add_unconformity(self, *, unconformity_data=None, create_new=True):
item.setSizeHint(unconformity_widget.sizeHint())
self.unitList.addItem(item)
self.unitList.setItemWidget(item, unconformity_widget)
+ unconformity_widget.set_available_faults(self._get_available_fault_names())
+ unconformity_widget.setData(unconformity_data)
# Cache the widget for efficient updates
self._widget_cache[unconformity.uuid] = (unconformity_widget, item)
@@ -605,6 +637,30 @@ def update_element(self, unit_widget):
"""
if self.data_manager:
unit_data = unit_widget.getData()
+ if isinstance(unit_widget, UnconformityWidget):
+ fault_name = unit_data.pop('fault_name', None)
+ is_fault_boundary = unit_data.get('unconformity_type') == 'fault'
+ if is_fault_boundary:
+ # The core stratigraphic column only knows erode/onlap --
+ # the fault link lives in the data manager's side table
+ # (see set_fault_boundary), so store it as a plain
+ # erosional boundary here.
+ unit_data['unconformity_type'] = 'erode'
+ if is_fault_boundary and fault_name:
+ self.data_manager.set_fault_boundary(unit_widget.uuid, fault_name)
+ if not self.data_manager.fault_spans_model_domain(fault_name):
+ QMessageBox.information(
+ self,
+ "Fault Domain Boundary",
+ f"Fault '{fault_name}' does not reach every edge of the model "
+ "bounding box.\n\nA fault used as a domain boundary crops the "
+ "whole model, so its digitised trace will automatically be "
+ "extended out to the domain edges along its overall trend when "
+ "the model is built. For best results the trace should still "
+ "roughly follow the fault's real direction across the gap.",
+ )
+ else:
+ self.data_manager.clear_fault_boundary(unit_widget.uuid)
self.data_manager._stratigraphic_column.update_element(unit_data)
# Trigger callback to notify all listeners of the change
if self.data_manager.stratigraphic_column_callback:
diff --git a/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_unit.ui b/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_unit.ui
index 3f9a3b9..306f738 100644
--- a/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_unit.ui
+++ b/loopstructural/gui/modelling/stratigraphic_column/stratigraphic_unit.ui
@@ -80,7 +80,7 @@
0.000000000000000
- 10000.000000000000000
+ 100000000.000000000000000
diff --git a/loopstructural/gui/modelling/stratigraphic_column/unconformity.py b/loopstructural/gui/modelling/stratigraphic_column/unconformity.py
index 115b9df..a828edd 100644
--- a/loopstructural/gui/modelling/stratigraphic_column/unconformity.py
+++ b/loopstructural/gui/modelling/stratigraphic_column/unconformity.py
@@ -10,6 +10,7 @@
class UnconformityWidget(QWidget):
deleteRequested = pyqtSignal(QWidget) # Signal to request deletion
+ dataChanged = pyqtSignal() # Type or fault-name changed
dragHandlePressed = pyqtSignal() # Drag handle mouse-down
dragHandleMoved = pyqtSignal(QPoint) # Drag handle mouse-move (global pos)
dragHandleReleased = pyqtSignal() # Drag handle mouse-up
@@ -25,9 +26,9 @@ def __init__(
self.buttonDelete.clicked.connect(self.request_delete)
self.uuid = uuid
self.unconformity_type = 'erode'
- # self.comboBoxUnconformityType.currentIndexChanged.connect(
- # lambda: setattr(self, 'unconformity_type', self.comboBoxUnconformityType.currentText())
- # )
+ self.fault_name = None
+ self.comboBoxUnconformityType.currentIndexChanged.connect(self._on_type_changed)
+ self.comboBoxFaultName.currentIndexChanged.connect(self._on_fault_name_changed)
# The row's combo box/buttons cover the whole widget, so a QListWidget's
# built-in drag-and-drop can never see a mouse press to start a
# reorder. Route presses on the dedicated grip label through here instead.
@@ -55,21 +56,73 @@ def request_delete(self):
self.deleteRequested.emit(self)
+ def _on_type_changed(self, _index):
+ self.unconformity_type = self.comboBoxUnconformityType.currentText()
+ self.comboBoxFaultName.setVisible(self.unconformity_type == 'fault')
+ if self.unconformity_type == 'fault':
+ self.fault_name = self.comboBoxFaultName.currentText() or None
+ else:
+ self.fault_name = None
+ self.dataChanged.emit()
+
+ def _on_fault_name_changed(self, _index):
+ if self.unconformity_type != 'fault':
+ return
+ self.fault_name = self.comboBoxFaultName.currentText() or None
+ self.dataChanged.emit()
+
+ def set_available_faults(self, fault_names):
+ """Populate the fault-name picker, keeping the current selection if
+ it is still available (e.g. after the fault trace layer changes).
+ """
+ fault_names = list(fault_names or [])
+ if [
+ self.comboBoxFaultName.itemText(i) for i in range(self.comboBoxFaultName.count())
+ ] == fault_names:
+ return
+ self.comboBoxFaultName.blockSignals(True)
+ try:
+ self.comboBoxFaultName.clear()
+ self.comboBoxFaultName.addItems(fault_names)
+ if self.fault_name and self.fault_name in fault_names:
+ self.comboBoxFaultName.setCurrentText(self.fault_name)
+ finally:
+ self.comboBoxFaultName.blockSignals(False)
+
def setData(self, data: Optional[dict] = None):
"""Set the data for the unconformity widget.
Parameters
----------
data : dict or None
- Dictionary containing 'unconformity_type' key. If None, defaults are used.
+ Dictionary with an 'unconformity_type' key ('erode', 'onlap' or
+ 'fault'), and a 'fault_name' key when the type is 'fault'. If
+ None, defaults are used.
"""
- if data:
- self.unconformity_type = data.get("unconformity_type", "")
- # self.unconformityTypeComboBox.setCurrentIndex(
- # self.unconformityTypeComboBox.findText(self.unconformity_type)
- # )
- else:
- self.unconformity_type = 'erode'
- # self.unconformityTypeComboBox.setCurrentIndex(
- # self.unconformityTypeComboBox.findText(self.unconformity_type)
- # )
+ self.unconformity_type = (data or {}).get("unconformity_type", "erode")
+ self.fault_name = (
+ (data or {}).get("fault_name") if self.unconformity_type == 'fault' else None
+ )
+
+ self.comboBoxUnconformityType.blockSignals(True)
+ self.comboBoxFaultName.blockSignals(True)
+ try:
+ index = self.comboBoxUnconformityType.findText(self.unconformity_type)
+ if index >= 0:
+ self.comboBoxUnconformityType.setCurrentIndex(index)
+ self.comboBoxFaultName.setVisible(self.unconformity_type == 'fault')
+ if self.fault_name:
+ self.comboBoxFaultName.setCurrentText(self.fault_name)
+ finally:
+ self.comboBoxUnconformityType.blockSignals(False)
+ self.comboBoxFaultName.blockSignals(False)
+
+ def getData(self):
+ """Return this row's data for the data manager: uuid, unconformity_type
+ and (when the boundary is fault-linked) fault_name.
+ """
+ return {
+ 'uuid': self.uuid,
+ 'unconformity_type': self.unconformity_type,
+ 'fault_name': self.fault_name,
+ }
diff --git a/loopstructural/gui/modelling/stratigraphic_column/unconformity.ui b/loopstructural/gui/modelling/stratigraphic_column/unconformity.ui
index cda9c0f..48239f4 100644
--- a/loopstructural/gui/modelling/stratigraphic_column/unconformity.ui
+++ b/loopstructural/gui/modelling/stratigraphic_column/unconformity.ui
@@ -61,6 +61,11 @@
onlap
+ -
+
+ fault
+
+
-
@@ -70,6 +75,16 @@
+ -
+
+
+ Fault whose surface realises this boundary as a domain split
+
+
+ false
+
+
+
diff --git a/loopstructural/gui/visualisation/feature_list_widget.py b/loopstructural/gui/visualisation/feature_list_widget.py
index 90c17a2..4c22961 100644
--- a/loopstructural/gui/visualisation/feature_list_widget.py
+++ b/loopstructural/gui/visualisation/feature_list_widget.py
@@ -456,8 +456,18 @@ def add_stratigraphic_surfaces(self):
stratigraphic_surfaces = self.model_manager.model.get_stratigraphic_surfaces()
for surface in stratigraphic_surfaces:
+ mesh = surface.vtk()
+ if mesh.n_points == 0:
+ # A unit with no digitised data of its own (e.g. an
+ # undigitised placeholder like "Top") can have no
+ # constrained geometry anywhere in the model, so its
+ # isovalue may not intersect the solved field at all --
+ # pyvista refuses to plot an empty mesh, so skip it rather
+ # than crashing every surface after it in this loop.
+ logger.info(f"Skipping '{surface.name}': isosurface has no geometry.")
+ continue
self.viewer.add_mesh_object(
- surface.vtk(),
+ mesh,
name=surface.name,
color=surface.colour,
source_feature=surface.name,
diff --git a/loopstructural/main/data_manager.py b/loopstructural/main/data_manager.py
index c91a4c8..6ef0f6a 100644
--- a/loopstructural/main/data_manager.py
+++ b/loopstructural/main/data_manager.py
@@ -89,6 +89,13 @@ def __init__(self, *, project=None, mapCanvas=None, logger=None):
self.logger = logger
self._stratigraphic_column = StratigraphicColumn()
self._fault_topology = FaultTopology(self._stratigraphic_column)
+ # Maps a stratigraphic-column unconformity's uuid to the name of an
+ # existing fault that should realise that boundary, instead of a
+ # flat isovalue surface -- see `set_fault_boundary`. Kept as a plain
+ # dict (not part of `StratigraphicColumn` itself) so this stays a
+ # plugin-side concept for now; passed by reference to the model
+ # manager the same way `_stratigraphic_column`/`_fault_topology` are.
+ self._fault_boundaries: dict[str, str] = {}
self._model_manager = None
self.bounding_box_callback = None
self.basal_contacts_callback = None
@@ -137,6 +144,7 @@ def set_model_manager(self, model_manager):
self._model_manager = model_manager
self._model_manager.set_stratigraphic_column(self._stratigraphic_column)
self._model_manager.set_fault_topology(self._fault_topology)
+ self._model_manager.set_fault_boundaries(self._fault_boundaries)
self._model_manager.update_bounding_box(self._bounding_box)
def set_bounding_box(
@@ -557,10 +565,103 @@ def add_to_stratigraphic_column(self, unit_data):
def remove_from_stratigraphic_column(self, unit_uuid):
"""Remove a unit or unconformity from the stratigraphic column."""
self._stratigraphic_column.remove_unit(uuid=unit_uuid)
+ self._fault_boundaries.pop(unit_uuid, None)
self.update_stratigraphy()
if self.stratigraphic_column_callback:
self.stratigraphic_column_callback()
+ def set_fault_boundary(self, unconformity_uuid, fault_name):
+ """Mark a stratigraphic-column unconformity as realised by an
+ existing fault instead of a flat isovalue surface.
+
+ Parameters
+ ----------
+ unconformity_uuid : str
+ uuid of the `StratigraphicUnconformity` element in the column.
+ fault_name : str
+ Name of an existing fault (as known to `_fault_topology`) whose
+ surface should be used as the domain boundary at this point in
+ the column.
+ """
+ self._fault_boundaries[unconformity_uuid] = fault_name
+ if self.stratigraphic_column_callback:
+ self.stratigraphic_column_callback()
+
+ def clear_fault_boundary(self, unconformity_uuid):
+ """Undo `set_fault_boundary`, reverting the unconformity to a plain isovalue boundary."""
+ if self._fault_boundaries.pop(unconformity_uuid, None) is not None:
+ if self.stratigraphic_column_callback:
+ self.stratigraphic_column_callback()
+
+ def get_fault_boundary(self, unconformity_uuid):
+ """Return the fault name linked to this unconformity, or None."""
+ return self._fault_boundaries.get(unconformity_uuid)
+
+ def get_fault_boundaries(self):
+ """Return the uuid -> fault_name mapping of all fault-linked boundaries."""
+ return dict(self._fault_boundaries)
+
+ def get_fault_boundary_fault_names(self):
+ """Return the set of fault names currently used as domain boundaries.
+
+ These faults are built as non-displacing domain splits (see
+ `GeologicalModelManager.update_foliation_features`), so they should
+ not also be offered in fault-fault (FAULTED/ABUTTING) or
+ fault-stratigraphy relationship editors, which assume a
+ displacement-modelled fault.
+ """
+ return set(self._fault_boundaries.values())
+
+ def fault_spans_model_domain(self, fault_name, *, tolerance=0.0):
+ """Check whether a fault's trace data reaches every edge of the
+ model's XY bounding box.
+
+ A fault used as a domain boundary crops the *entire* model on
+ either side of its interpolated surface (see
+ `LoopStructural.modelling.core._model_relationships`), so unlike an
+ ordinary local fault trace it needs to be constrained across the
+ whole domain -- otherwise the interpolator extrapolates the crop
+ surface into areas with no supporting data. Returns True if no
+ fault trace data is available yet (nothing to check against).
+
+ Parameters
+ ----------
+ fault_name : str
+ Name of the fault to check, as found in `get_fault_traces()`'s layer.
+ tolerance : float, optional
+ Allowed gap, in model units, between the trace's extent and the
+ bounding box edge before it is considered "not spanning".
+ """
+ if self._fault_traces is None or self._fault_traces['layer'] is None:
+ return True
+ layer = self._fault_traces['layer']
+ name_field = self._fault_traces['fault_name_field']
+ trace_extent = None
+ for feature in layer.getFeatures():
+ if name_field is not None and str(feature[name_field]) != str(fault_name):
+ continue
+ geom = feature.geometry()
+ if geom is None or geom.isEmpty():
+ continue
+ bbox = geom.boundingBox()
+ trace_extent = bbox if trace_extent is None else trace_extent.combineExtentWith(bbox)
+ if trace_extent is None:
+ return True
+ xmin, ymin = self._bounding_box.origin[0], self._bounding_box.origin[1]
+ xmax, ymax = self._bounding_box.maximum[0], self._bounding_box.maximum[1]
+ # A boundary only needs to fully cross the domain along one axis (an
+ # east-west or north-south cut) to split the whole model -- it does
+ # not need to cover the full bounding box in both directions.
+ spans_x = (
+ trace_extent.xMinimum() <= xmin + tolerance
+ and trace_extent.xMaximum() >= xmax - tolerance
+ )
+ spans_y = (
+ trace_extent.yMinimum() <= ymin + tolerance
+ and trace_extent.yMaximum() >= ymax - tolerance
+ )
+ return spans_x or spans_y
+
def update_stratigraphic_column_order(self, new_order):
"""Update the order of units in the stratigraphic column."""
if not isinstance(new_order, list):
@@ -956,6 +1057,7 @@ def to_dict(self):
'stratigraphic_column': (
self._stratigraphic_column.to_dict() if self._stratigraphic_column else None
),
+ 'fault_boundaries': dict(self._fault_boundaries),
'dem_layer': dem_layer_name if self.dem_layer else None,
'use_dem': self.use_dem,
'elevation': self.elevation,
@@ -997,7 +1099,12 @@ def from_dict(self, data):
self._structural_orientations = data['structural_orientations']
if 'stratigraphic_column' in data:
self._stratigraphic_column = StratigraphicColumn.from_dict(data['stratigraphic_column'])
+ # See the matching call in update_from_dict for why this is needed.
+ self._stratigraphic_column.update_unit_values()
self.stratigraphic_column_callback()
+ self._fault_boundaries.clear()
+ if data.get('fault_boundaries'):
+ self._fault_boundaries.update(data['fault_boundaries'])
if 'widget_settings' in data:
self.widget_settings = data['widget_settings']
@@ -1011,6 +1118,23 @@ def from_dict(self, data):
def update_from_dict(self, data):
"""Update the data manager from a dictionary."""
+ # Model CRS must be restored before anything below that reprojects
+ # a layer against it (basal_contacts/fault_traces/
+ # structural_orientations, all via get_model_crs()) -- restoring it
+ # last meant every reprojection during project load ran against
+ # whatever `_use_project_crs`/`_model_crs` still held from __init__
+ # (use_project_crs=True, i.e. the *project's* CRS) instead of the
+ # saved model CRS, silently skipping reprojection whenever a layer
+ # happened to already be in the project's CRS but not the model's.
+ if 'use_project_crs' in data:
+ self._use_project_crs = data['use_project_crs']
+ else:
+ self._use_project_crs = True
+ if 'model_crs' in data and data['model_crs'] is not None:
+ crs = QgsCoordinateReferenceSystem(data['model_crs'])
+ if crs.isValid():
+ self.set_model_crs(crs, use_project_crs=self._use_project_crs)
+
if 'bounding_box' in data:
self.set_bounding_box(
xmin=data['bounding_box']['origin'][0],
@@ -1077,25 +1201,30 @@ def update_from_dict(self, data):
)
if 'stratigraphic_column' in data:
self._stratigraphic_column.update_from_dict(data['stratigraphic_column'])
+ # update_from_dict restores elements via add_element, not
+ # add_unit -- only add_unit computes each unit's min/max
+ # scalar-field range as a side effect. Without this, every
+ # restored unit keeps the default (0, inf) range, so
+ # evaluate_model can't tell any unit in a group apart from any
+ # other and just labels every point with whichever unit was
+ # last in the group (see GeologicalModelManager.
+ # set_stratigraphic_column, which already does this for the
+ # very first load -- this covers every reload afterwards).
+ self._stratigraphic_column.update_unit_values()
else:
self._stratigraphic_column.clear()
+ # Mutate in place rather than reassign: `_fault_boundaries` is
+ # shared by reference with the model manager (see set_model_manager).
+ self._fault_boundaries.clear()
+ if data.get('fault_boundaries'):
+ self._fault_boundaries.update(data['fault_boundaries'])
+
if 'widget_settings' in data:
self.widget_settings = data['widget_settings']
else:
self.widget_settings = {}
- # Load model CRS settings
- if 'use_project_crs' in data:
- self._use_project_crs = data['use_project_crs']
- else:
- self._use_project_crs = True
-
- if 'model_crs' in data and data['model_crs'] is not None:
- crs = QgsCoordinateReferenceSystem(data['model_crs'])
- if crs.isValid():
- self.set_model_crs(crs, use_project_crs=self._use_project_crs)
-
if self.stratigraphic_column_callback:
self.stratigraphic_column_callback()
diff --git a/loopstructural/main/model_manager.py b/loopstructural/main/model_manager.py
index 4bea402..cbd6074 100644
--- a/loopstructural/main/model_manager.py
+++ b/loopstructural/main/model_manager.py
@@ -18,8 +18,11 @@
import pandas as pd
from LoopStructural.datatypes import BoundingBox
from LoopStructural.modelling.core.fault_topology import FaultRelationshipType
-from LoopStructural.modelling.core.stratigraphic_column import StratigraphicColumn
-from LoopStructural.modelling.features import FeatureType, StructuralFrame
+from LoopStructural.modelling.core.stratigraphic_column import (
+ StratigraphicColumn,
+ StratigraphicUnconformity,
+)
+from LoopStructural.modelling.features import FeatureType, StructuralFrame, UnconformityFeature
from LoopStructural.modelling.features.fold import FoldFrame
from LoopStructural.utils.observer import Observable
@@ -111,6 +114,10 @@ def __init__(self, debug_manager=None):
self.stratigraphy: Dict[str, StratigraphyEntry] = defaultdict(dict)
self.stratigraphic_column = None
self.fault_topology = None
+ # uuid (of a StratigraphicUnconformity in stratigraphic_column) ->
+ # fault name; see `set_fault_boundaries`. Shared by reference with
+ # ModellingDataManager, same as stratigraphic_column/fault_topology.
+ self.fault_boundaries: Dict[str, str] = {}
# Observers managed by Observable base class
self.dem_function = lambda x, y: 0
# internal flag to temporarily suppress notifications (used when
@@ -241,6 +248,12 @@ def set_fault_topology(self, fault_topology):
except Exception:
pass
+ def set_fault_boundaries(self, fault_boundaries: Dict[str, str]):
+ """Set the uuid -> fault_name mapping of fault-linked stratigraphic
+ column boundaries (see `ModellingDataManager.set_fault_boundary`).
+ """
+ self.fault_boundaries = fault_boundaries
+
# Topology events that change what actually feeds the interpolator (as
# opposed to just which side of an already-solved fault gets cropped
# away) and therefore need Initialize Model to be re-run before Solve
@@ -531,6 +544,230 @@ def update_stratigraphic_column(self, stratigraphic_column: StratigraphicColumn)
# def update_stratigraphic_unit(self, unit_data):
# self.data
+ def _closing_fault_boundary(self, group):
+ """Return the fault name that closes `group` from above, if any.
+
+ The boundary "closing" a group off from the next (younger) group is
+ the first `StratigraphicUnconformity` above the group's youngest
+ unit in `stratigraphic_column.order`. If that boundary has been
+ linked to a fault (`ModellingDataManager.set_fault_boundary`), the
+ group should be capped by that fault's surface -- a non-displacing
+ domain split, see `create_and_add_domain_fault` -- instead of the
+ flat isovalue-0 surface `add_unconformity` uses.
+ """
+ if not group.units or not self.fault_boundaries or self.stratigraphic_column is None:
+ return None
+ order = self.stratigraphic_column.order
+ youngest_uuid = group.units[0].uuid
+ start = next((i for i, e in enumerate(order) if e.uuid == youngest_uuid), None)
+ if start is None:
+ return None
+ for element in order[start + 1 :]:
+ if isinstance(element, StratigraphicUnconformity):
+ return self.fault_boundaries.get(element.uuid)
+ return None
+
+ def _clip_line_to_bounding_box(self, centroid, direction):
+ """Return the (t_min, t_max) range along `centroid + t*direction`
+ (XY only) that stays within the model's bounding box, or None if
+ the line never crosses it. Standard slab-method line/box clip.
+ """
+ origin_xy = np.array(self.model.bounding_box.origin[:2], dtype=float)
+ maximum_xy = np.array(self.model.bounding_box.maximum[:2], dtype=float)
+ t_min, t_max = -np.inf, np.inf
+ for axis in (0, 1):
+ d = direction[axis]
+ if abs(d) < 1e-12:
+ if centroid[axis] < origin_xy[axis] or centroid[axis] > maximum_xy[axis]:
+ return None
+ continue
+ t0 = (origin_xy[axis] - centroid[axis]) / d
+ t1 = (maximum_xy[axis] - centroid[axis]) / d
+ t0, t1 = min(t0, t1), max(t0, t1)
+ t_min = max(t_min, t0)
+ t_max = min(t_max, t1)
+ if t_min > t_max:
+ return None
+ return t_min, t_max
+
+ def _extend_fault_trace_to_domain(self, fault_data):
+ """Add two synthetic points that extend a fault's trace out to the
+ edges of the model's bounding box, and attach a `strike` column
+ derived from the trace's own *local* tangent at each point.
+
+ `create_and_add_domain_fault` interpolates a scalar field only from
+ the points it is given, over the model's exact bounding box (no
+ buffer, unlike a displacement fault's mesh) -- so a locally
+ digitised trace only reliably constrains the surface near itself,
+ and the domain crop can wander unpredictably further away.
+ Extending each end along its own local tangent, out to where it
+ meets the bounding box edge, keeps the interpolated surface
+ following the trace's actual trend all the way across the domain --
+ this is what makes the fault behave as an "infinite" domain
+ boundary rather than a locally-anchored patch.
+
+ Using each point's *local* tangent (rather than one global
+ best-fit line through the whole trace) matters for a genuinely
+ curved trace: fitting a single global line flattens that curvature
+ out, and extrapolating along it can land an extension point (or
+ bias the interpolated field generally) on the wrong side of the
+ real curve relative to data that's actually near the trace.
+ Confirmed on a live project: a global-line fit classified a
+ stratigraphic unit's own contact data as being on the opposite
+ side of the domain fault from a bounding-box corner that a proper
+ local (nearest-segment) classification put on the *same* side as
+ that data -- i.e. the global fit was extrapolating the wrong way.
+
+ Z at each synthetic point is extrapolated linearly against
+ distance along the local end segment, so a dipping trace keeps its
+ dip at the point it's extended from.
+ """
+ xy = fault_data[['X', 'Y']].to_numpy()
+ z = fault_data['Z'].to_numpy()
+ n = len(xy)
+ result = fault_data.copy()
+ if n < 2:
+ result['strike'] = np.nan
+ return result
+
+ # Local tangent per point: central difference for interior points,
+ # forward/backward difference at the ends. Assumes points follow
+ # the digitised line's vertex order (true for AllSampler-derived
+ # trace data, which walks each LineString's coords in order).
+ tangents = np.zeros((n, 2))
+ tangents[0] = xy[1] - xy[0]
+ tangents[-1] = xy[-1] - xy[-2]
+ if n > 2:
+ tangents[1:-1] = xy[2:] - xy[:-2]
+ # strikedip2vector's strike is degrees clockwise from North (+Y);
+ # atan2(dx, dy) matches that convention directly.
+ result['strike'] = np.degrees(np.arctan2(tangents[:, 0], tangents[:, 1])) % 360
+
+ new_rows = []
+ # Extend backward past the first point, continuing on its own
+ # local tangent (pointing away from the second point).
+ start_seg = xy[1] - xy[0]
+ start_len = np.linalg.norm(start_seg)
+ if start_len > 1e-9:
+ direction = -start_seg / start_len
+ clipped = self._clip_line_to_bounding_box(xy[0], direction)
+ if clipped is not None:
+ _, t_max = clipped
+ if t_max > 0:
+ point_xy = xy[0] + direction * t_max
+ z_slope = (z[1] - z[0]) / start_len
+ new_rows.append(
+ {
+ 'X': point_xy[0],
+ 'Y': point_xy[1],
+ 'Z': z[0] - z_slope * t_max,
+ 'strike': result['strike'].iloc[0],
+ }
+ )
+ # Extend forward past the last point, continuing on its own local
+ # tangent (pointing away from the second-to-last point).
+ end_seg = xy[-1] - xy[-2]
+ end_len = np.linalg.norm(end_seg)
+ if end_len > 1e-9:
+ direction = end_seg / end_len
+ clipped = self._clip_line_to_bounding_box(xy[-1], direction)
+ if clipped is not None:
+ _, t_max = clipped
+ if t_max > 0:
+ point_xy = xy[-1] + direction * t_max
+ z_slope = (z[-1] - z[-2]) / end_len
+ new_rows.append(
+ {
+ 'X': point_xy[0],
+ 'Y': point_xy[1],
+ 'Z': z[-1] + z_slope * t_max,
+ 'strike': result['strike'].iloc[-1],
+ }
+ )
+ if not new_rows:
+ return result
+ return pd.concat([result, pd.DataFrame(new_rows)], ignore_index=True)
+
+ def _domain_fault_dip(self, fault_entry):
+ """Dip (degrees from horizontal) for a domain-boundary fault.
+
+ Uses the ingested fault trace data's `dip` column if present
+ (matching how `update_fault_features` picks up dip for a
+ displacement fault), otherwise defaults to vertical (90 degrees).
+ """
+ raw_data = fault_entry.get('data') if fault_entry else None
+ if raw_data is not None and 'dip' in raw_data:
+ dip_values = raw_data['dip'].dropna()
+ if not dip_values.empty:
+ return float(dip_values.mean())
+ return 90.0
+
+ def _build_domain_fault_boundary(self, fault_name, groupname):
+ """Build `fault_name` as a domain-fault boundary in place of a flat unconformity.
+
+ `GeologicalModel.create_and_add_domain_fault` (unlike
+ `create_and_add_fault`/`create_and_add_foliation`) has no `data=`
+ parameter -- it always reads the fault's points from `model.data`
+ filtered by `feature_name`, so the fault's trace data is registered
+ there first. Registration replaces any rows already tagged with
+ this fault name so repeated Initialize Model runs stay idempotent
+ instead of accumulating duplicate points on every rebuild.
+
+ Returns True if the domain fault was built, False if it fell back
+ to a flat unconformity for lack of trace data (caller should then
+ call `self.model.add_unconformity` itself).
+ """
+ fault_entry = self.faults.get(fault_name)
+ fault_data = fault_entry.get('data') if fault_entry else None
+ if fault_data is None or fault_data.empty:
+ self._debug_manager and self._debug_manager.log(
+ f"Fault '{fault_name}' is linked as a domain boundary for group "
+ f"'{groupname}' but has no trace data; using a flat unconformity instead.",
+ log_level=2,
+ )
+ return False
+ extended = self._extend_fault_trace_to_domain(fault_data[['X', 'Y', 'Z']].copy())
+
+ value_rows = extended[['X', 'Y', 'Z']].copy()
+ value_rows['feature_name'] = fault_name
+ value_rows['val'] = 0
+
+ orientation_rows = extended.dropna(subset=['strike'])[['X', 'Y', 'Z', 'strike']].copy()
+ data_for_fault = value_rows
+ if not orientation_rows.empty:
+ orientation_rows['dip'] = self._domain_fault_dip(fault_entry)
+ orientation_rows['feature_name'] = fault_name
+ orientation_rows['val'] = np.nan
+ data_for_fault = pd.concat([value_rows, orientation_rows], ignore_index=True)
+ # Unlike create_and_add_foliation/create_and_add_fault (which
+ # normalise their own `data=` argument internally via
+ # model.prepare_data before building), create_and_add_domain_fault
+ # reads straight from model.data with no normalisation -- it
+ # expects every standard column (gx/gy/gz/nx/ny/nz/...) to already
+ # be present, or the interpolator crashes looking them up. Run it
+ # through prepare_data ourselves before writing it in.
+ data_for_fault = self.model.prepare_data(data_for_fault, include_feature_name=True)
+ existing_data = self.model.data
+ if existing_data is not None and not existing_data.empty:
+ existing_data = existing_data.loc[existing_data['feature_name'] != fault_name]
+ self.model.data = pd.concat([existing_data, data_for_fault], ignore_index=True)
+ else:
+ # An empty placeholder frame built with `pd.DataFrame(columns=...)`
+ # defaults every column to object dtype; concatenating that with
+ # `data_for_fault`'s float columns can leave the result as
+ # object dtype too, and `add_data_to_interpolator` then fails
+ # calling `np.isnan` on an object-dtype column. Assign directly
+ # instead of concatenating with a dtype-less placeholder.
+ self.model.data = data_for_fault
+ self.model.create_and_add_domain_fault(
+ fault_name,
+ nelements=PlgSettingsStructure.interpolator_nelements,
+ npw=PlgSettingsStructure.interpolator_npw,
+ cpw=PlgSettingsStructure.interpolator_cpw,
+ regularisation=PlgSettingsStructure.interpolator_regularisation,
+ )
+ return True
+
def update_foliation_features(self):
"""Builds the stratigraphic feature from the stratigraphic column data
and the basal contacts and structural orientations data.
@@ -550,6 +787,29 @@ def update_foliation_features(self):
groupname = group.name
stratigraphic_column[groupname] = {}
for u in group.units:
+ # A unit's own `val` is `u.min()` -- the cumulative
+ # thickness *before* this unit's own thickness is added --
+ # matching `StratigraphicColumn.update_unit_values` (a unit
+ # is only added to the column oldest-first via `where=
+ # 'top'`, so `min()` is the boundary shared with the
+ # next-*older* neighbour processed just before it, i.e.
+ # this unit's own base) and `get_isovalues()` (LoopStructural
+ # core). `group.units` (from `get_groups()`) is already in
+ # that oldest-after-youngest walk order, and `get_isovalues()`
+ # accumulates over it directly with no extra reversal, so
+ # this loop must not reverse it either -- doing so trains
+ # each unit with the wrong scalar value, see
+ # test_stratigraphic_value_consistency.py.
+ #
+ # `val` must accumulate every unit's thickness regardless of
+ # whether that unit has any digitised data -- get_isovalues()
+ # assigns each unit's isovalue purely from cumulative
+ # thickness, with no knowledge of which units were actually
+ # mapped. Skipping the increment for an unmapped unit (e.g.
+ # a "Top" placeholder with no contact points) would shift
+ # every val assigned to units after it in this loop, so
+ # extracted isosurfaces would get labelled with the wrong
+ # unit name even though the geometry itself is fine.
unit_data = self.stratigraphy.get(u.name, None)
if unit_data is not None:
if 'contact' in unit_data:
@@ -581,11 +841,51 @@ def update_foliation_features(self):
cpw=PlgSettingsStructure.interpolator_cpw,
regularisation=PlgSettingsStructure.interpolator_regularisation,
)
- self.model.add_unconformity(foliation, 0)
+ fault_name = self._closing_fault_boundary(group)
+ if fault_name is None or not self._build_domain_fault_boundary(fault_name, groupname):
+ self.model.add_unconformity(foliation, 0)
+ self._strip_spurious_regions_from_domain_faults()
self.model.stratigraphic_column = self.stratigraphic_column
# foliation features were rebuilt; let observers know
self._emit('foliation_features_updated')
+ def _strip_spurious_regions_from_domain_faults(self):
+ """Work around a LoopStructural core gap that corrupts a domain
+ fault's own scalar field.
+
+ `add_unconformity`'s backward crop walk (in LoopStructural's
+ `_model_relationships.FeatureRelationshipManager.add_unconformity`)
+ only recognises `FeatureType.FAULT`/`INACTIVEFAULT` as "already
+ handled, skip" -- it doesn't know about `FeatureType.DOMAINFAULT`.
+ So whenever a later group in the same Initialize Model run falls
+ back to a plain `add_unconformity` (no fault linked to its
+ boundary), that call walks straight through any domain-boundary
+ fault built earlier and incorrectly adds itself as a region on it.
+ Since `GeologicalFeature.evaluate_value` returns NaN wherever a
+ feature's regions don't hold, the domain fault's own field then
+ reads as NaN on whichever side that unrelated unconformity's
+ condition fails -- "NaN on one side" of an otherwise valid domain
+ boundary.
+
+ A domain fault is meant to crop other features, not be cropped
+ itself -- except by an *earlier* domain fault, which is a
+ legitimate, intentional cascade (`_add_domain_fault_above` adds
+ that as a plain lambda region, not an `UnconformityFeature`). So
+ strip only the `UnconformityFeature`-typed regions injected by the
+ bug, leaving any real domain-fault-vs-domain-fault crop intact.
+ """
+ for feature in self.model.features:
+ if getattr(feature, 'type', None) != FeatureType.DOMAINFAULT:
+ continue
+ kept = [r for r in feature.regions if not isinstance(r, UnconformityFeature)]
+ if len(kept) != len(feature.regions):
+ self._debug_manager and self._debug_manager.log(
+ f"Removing {len(feature.regions) - len(kept)} unconformity region(s) "
+ f"incorrectly applied to domain-boundary fault '{feature.name}'.",
+ log_level=2,
+ )
+ feature.regions = kept
+
def _report_progress(self, message: str):
"""Report progress on a long-running model update.
@@ -683,7 +983,13 @@ def _cutting_faults_for(self, fault_name):
def update_fault_features(self):
"""Update the fault features in the geological model."""
+ domain_boundary_faults = set(self.fault_boundaries.values())
for fault_name in self._fault_build_order():
+ if fault_name in domain_boundary_faults:
+ # Built as a non-displacing domain fault in
+ # update_foliation_features instead -- see
+ # `_build_domain_fault_boundary`.
+ continue
fault_data = self.faults[fault_name]
self._report_progress(f"Building fault '{fault_name}'")
if qgisAttributeIsNone(fault_name):
@@ -734,6 +1040,19 @@ def update_fault_features(self):
)
self.apply_fault_abutting_relationships()
+ def _get_feature_by_name_or_none(self, name):
+ """Non-raising counterpart to `GeologicalModel.get_feature_by_name`.
+
+ The wrapped `GeologicalModel` raises `ValueError` for a name that
+ hasn't been built yet rather than returning None, which several
+ call sites in this manager treat as "not built yet, skip it" --
+ `__contains__` (`name in self.model`) checks `feature_name_index`
+ directly, so this restores that non-raising lookup.
+ """
+ if name not in self.model:
+ return None
+ return self.model.get_feature_by_name(name)
+
def apply_fault_abutting_relationships(self):
"""Re-apply fault-fault ABUTTING relationships as region crops on the
already-built fault surfaces, for every pair currently in the topology.
@@ -754,20 +1073,30 @@ def apply_fault_abutting_relationships(self):
"""
if self.fault_topology is None:
return
+ # Domain-boundary faults (see `_build_domain_fault_boundary`) are
+ # non-displacing GeologicalFeatures, not FaultSegments, and are only
+ # built later in `update_foliation_features` -- calling this before
+ # that has run would look up a feature that doesn't exist yet.
+ # ABUTTING relationships don't apply to them either way (that's a
+ # FaultSegment-specific crop), matching their exclusion from the
+ # Fault Adjacency tab.
+ domain_boundary_faults = set(self.fault_boundaries.values())
for f in self.fault_topology.faults:
- fault_feature = self.model.get_feature_by_name(f)
+ if f in domain_boundary_faults:
+ continue
+ fault_feature = self._get_feature_by_name_or_none(f)
if fault_feature is None or not hasattr(fault_feature, 'abut'):
continue
coord0 = fault_feature.__getitem__(0)
for f2 in self.fault_topology.faults:
- if f == f2:
+ if f == f2 or f2 in domain_boundary_faults:
continue
relationship = self.fault_topology.get_fault_relationship(f, f2)
existing_region = fault_feature.abut.get(f2)
if relationship is FaultRelationshipType.ABUTTING:
if existing_region is not None:
continue # already cropped against f2
- f2_feature = self.model.get_feature_by_name(f2)
+ f2_feature = self._get_feature_by_name_or_none(f2)
if f2_feature is None:
continue
# Determine which side of f2 to keep ourselves, ignoring any
@@ -911,7 +1240,8 @@ def update_model(
len(self.stratigraphic_column.get_groups()) if self.stratigraphic_column else 0
)
self._progress_callback = progress_callback
- self._progress_total = len(self.faults) + group_count
+ displacement_fault_count = len(set(self.faults) - set(self.fault_boundaries.values()))
+ self._progress_total = displacement_fault_count + group_count
self._progress_current = 0
dbg = getattr(self, '_debug_manager', None)
if dbg is not None:
diff --git a/loopstructural/toolbelt/preferences.py b/loopstructural/toolbelt/preferences.py
index 99cfb6e..083ce8b 100644
--- a/loopstructural/toolbelt/preferences.py
+++ b/loopstructural/toolbelt/preferences.py
@@ -31,7 +31,7 @@ class PlgSettingsStructure:
debug_directory: str = ""
version: str = __version__
interpolator_type: str = 'FDI'
- interpolator_nelements: int = 10000
+ interpolator_nelements: int = 50000
interpolator_regularisation: float = 1.0
interpolator_cpw: float = 1.0
interpolator_npw: float = 1.0
diff --git a/tests/qgis/test_fault_domain_boundary.py b/tests/qgis/test_fault_domain_boundary.py
new file mode 100644
index 0000000..c4cfb35
--- /dev/null
+++ b/tests/qgis/test_fault_domain_boundary.py
@@ -0,0 +1,291 @@
+"""Regression tests for using a fault as a stratigraphic domain boundary.
+
+`GeologicalModelManager.update_foliation_features` normally closes each
+built stratigraphic group with a flat isovalue unconformity
+(`GeologicalModel.add_unconformity`). When the column boundary closing a
+group has been linked to a fault (via `fault_boundaries`, populated through
+`ModellingDataManager.set_fault_boundary`), the group should instead be
+capped by that fault's own surface -- a non-displacing domain split, built
+with `GeologicalModel.create_and_add_domain_fault` -- and that fault must be
+skipped by the ordinary displacement-fault build loop
+(`update_fault_features`), since a domain fault and a displacement fault are
+mutually exclusive roles for the same fault name.
+"""
+
+import numpy as np
+import pandas as pd
+import pytest
+from LoopStructural import FaultTopology, StratigraphicColumn
+from LoopStructural.datatypes import BoundingBox
+
+from loopstructural.main.model_manager import GeologicalModelManager
+from loopstructural.toolbelt.preferences import PlgSettingsStructure
+
+
+def _contact(unit_name):
+ return pd.DataFrame({'X': [0.0], 'Y': [0.0], 'Z': [0.0], 'source_unit': [unit_name]})
+
+
+def _fault_trace():
+ return pd.DataFrame({'X': [0.0, 10.0], 'Y': [0.0, 10.0], 'Z': [0.0, 0.0]})
+
+
+def _two_group_column():
+ """basin_floor/basin_fill, then a named unconformity, then cover_lower/cover_upper."""
+ column = StratigraphicColumn()
+ column.clear(basement=False)
+ column.add_unit(name='basin_floor', thickness=50.0, where='top')
+ column.add_unit(name='basin_fill', thickness=150.0, where='top')
+ boundary = column.add_unconformity(name='regional_unconformity', where='top')
+ column.add_unit(name='cover_lower', thickness=80.0, where='top')
+ column.add_unit(name='cover_upper', thickness=120.0, where='top')
+ return column, boundary
+
+
+@pytest.fixture
+def manager(monkeypatch):
+ manager = GeologicalModelManager()
+ calls = {'foliation': [], 'unconformity': [], 'domain_fault': [], 'fault': []}
+
+ def fake_create_and_add_foliation(name, data=None, **kwargs):
+ calls['foliation'].append((name, data))
+ return object()
+
+ def fake_add_unconformity(feature, value, **kwargs):
+ calls['unconformity'].append((feature, value))
+
+ def fake_create_and_add_domain_fault(fault_surface_data, **kwargs):
+ calls['domain_fault'].append(fault_surface_data)
+ return object()
+
+ def fake_create_and_add_fault(fault_name, displacement, **kwargs):
+ calls['fault'].append(fault_name)
+ return object()
+
+ monkeypatch.setattr(manager.model, 'create_and_add_foliation', fake_create_and_add_foliation)
+ monkeypatch.setattr(manager.model, 'add_unconformity', fake_add_unconformity)
+ monkeypatch.setattr(
+ manager.model, 'create_and_add_domain_fault', fake_create_and_add_domain_fault
+ )
+ monkeypatch.setattr(manager.model, 'create_and_add_fault', fake_create_and_add_fault)
+ manager._calls = calls
+ return manager
+
+
+class TestFaultDomainBoundary:
+ def test_group_boundary_linked_to_fault_builds_domain_fault(self, manager):
+ column, boundary = _two_group_column()
+ manager.stratigraphic_column = column
+ for name in ('basin_floor', 'basin_fill', 'cover_lower', 'cover_upper'):
+ manager.stratigraphy[name]['contact'] = _contact(name)
+ manager.faults['boundary_fault']['data'] = _fault_trace()
+ manager.fault_boundaries[boundary.uuid] = 'boundary_fault'
+
+ manager.update_foliation_features()
+
+ assert manager._calls['domain_fault'] == ['boundary_fault']
+ # Only the topmost group (nothing above it in the column) still
+ # falls back to a flat unconformity.
+ assert len(manager._calls['unconformity']) == 1
+
+ registered = manager.model.data
+ fault_rows = registered.loc[registered['feature_name'] == 'boundary_fault']
+ # The default bounding box here (never set explicitly) is smaller
+ # than the trace itself, so no synthetic edge-extension points get
+ # added (see TestExtendFaultTraceToDomain for that, with a
+ # realistic bounding box). What's registered is the 2 trace points
+ # as value (val=0) constraints, plus the same 2 points again as
+ # orientation (val=NaN, strike/dip) constraints -- a domain fault
+ # needs both, see _domain_fault_dip / _extend_fault_trace_to_domain.
+ assert len(fault_rows) == 4
+ assert len(fault_rows.loc[fault_rows['val'] == 0]) == 2
+ assert fault_rows['val'].isna().sum() == 2
+
+ def test_group_boundary_without_fault_link_uses_flat_unconformity(self, manager):
+ column, _boundary = _two_group_column()
+ manager.stratigraphic_column = column
+ for name in ('basin_floor', 'basin_fill', 'cover_lower', 'cover_upper'):
+ manager.stratigraphy[name]['contact'] = _contact(name)
+
+ manager.update_foliation_features()
+
+ assert manager._calls['domain_fault'] == []
+ assert len(manager._calls['unconformity']) == 2
+
+ def test_missing_fault_data_falls_back_to_flat_unconformity(self, manager):
+ column, boundary = _two_group_column()
+ manager.stratigraphic_column = column
+ for name in ('basin_floor', 'basin_fill', 'cover_lower', 'cover_upper'):
+ manager.stratigraphy[name]['contact'] = _contact(name)
+ # Linked to a fault, but no trace data was ever ingested for it.
+ manager.fault_boundaries[boundary.uuid] = 'boundary_fault'
+
+ manager.update_foliation_features()
+
+ assert manager._calls['domain_fault'] == []
+ assert len(manager._calls['unconformity']) == 2
+
+ def test_domain_boundary_fault_skipped_by_displacement_fault_build(self, manager):
+ _column, boundary = _two_group_column()
+ manager.faults['boundary_fault']['data'] = _fault_trace()
+ manager.fault_boundaries[boundary.uuid] = 'boundary_fault'
+ manager.faults['normal_fault']['data'] = _fault_trace()
+
+ manager.update_fault_features()
+
+ assert manager._calls['fault'] == ['normal_fault']
+
+ def test_abutting_relationships_skip_unbuilt_domain_boundary_fault(self, manager):
+ """`update_fault_features` calls `apply_fault_abutting_relationships`
+ before `update_foliation_features` has run, so a domain-boundary
+ fault's feature (built there, not here) doesn't exist in the model
+ yet. `GeologicalModel.get_feature_by_name` raises ValueError rather
+ than returning None for a name that isn't built -- this must not
+ propagate out of `apply_fault_abutting_relationships`.
+ """
+ column, boundary = _two_group_column()
+ manager.fault_topology = FaultTopology(column)
+ manager.fault_topology.add_fault('boundary_fault')
+ manager.fault_topology.add_fault('normal_fault')
+ manager.faults['boundary_fault']['data'] = _fault_trace()
+ manager.faults['normal_fault']['data'] = _fault_trace()
+ manager.fault_boundaries[boundary.uuid] = 'boundary_fault'
+
+ # Neither fault has been built into manager.model yet.
+ manager.apply_fault_abutting_relationships() # must not raise
+
+
+class TestDomainFaultBuildsAndSolves:
+ """End-to-end regression test against the real LoopStructural build/solve
+ path (nothing mocked here, unlike the other test classes in this module).
+
+ `create_and_add_domain_fault` reads its points straight from
+ `model.data` with no column normalisation, unlike
+ `create_and_add_foliation`/`create_and_add_fault` (which run their
+ `data=` argument through `GeologicalModel.prepare_data` internally) --
+ so registering raw X/Y/Z/feature_name/val rows for a domain-boundary
+ fault builds fine but crashes later, during Solve Model, with
+ `KeyError: "None of [Index(['gx', 'gy', 'gz']...` deep inside
+ `add_data_to_interpolator`. `_build_domain_fault_boundary` must run its
+ data through `model.prepare_data` itself first.
+
+ It also pins the fix for a second, more subtle bug: LoopStructural's
+ `add_unconformity` backward crop walk only recognises FAULT/
+ INACTIVEFAULT feature types as "already handled, skip" -- not
+ DOMAINFAULT. `_two_group_column()`'s topmost group (cover) has no fault
+ link, so it falls back to a plain `add_unconformity` call, which walks
+ straight through the already-built `boundary_fault` domain fault and
+ incorrectly adds itself as a region on it. Since
+ `GeologicalFeature.evaluate_value` returns NaN wherever a feature's
+ regions don't hold, this reads as the domain fault's own scalar field
+ being NaN on whichever side that unrelated unconformity's condition
+ fails. `_strip_spurious_regions_from_domain_faults` removes it.
+ """
+
+ def test_domain_fault_builds_and_solves(self, monkeypatch):
+ monkeypatch.setattr(PlgSettingsStructure, 'interpolator_nelements', 200)
+ manager = GeologicalModelManager()
+ manager.update_bounding_box(BoundingBox(origin=[0, 0, -50], maximum=[100, 100, 50]))
+ column, boundary = _two_group_column()
+ manager.stratigraphic_column = column
+ for name in ('basin_floor', 'basin_fill', 'cover_lower', 'cover_upper'):
+ manager.stratigraphy[name]['contact'] = _contact(name)
+ manager.faults['boundary_fault']['data'] = pd.DataFrame(
+ {'X': [10.0, 90.0], 'Y': [10.0, 90.0], 'Z': [0.0, 0.0]}
+ )
+ manager.fault_boundaries[boundary.uuid] = 'boundary_fault'
+
+ manager.update_model(notify_observers=False)
+ manager.update_all_features(notify_observers=False) # must not raise
+
+ domain_fault = manager.model.get_feature_by_name('boundary_fault')
+ assert domain_fault is not None
+ assert domain_fault.regions == []
+
+ xs, ys = np.meshgrid(np.linspace(0, 100, 11), np.linspace(0, 100, 11))
+ pts = np.column_stack([xs.ravel(), ys.ravel(), np.zeros(xs.size)])
+ values = domain_fault.evaluate_value(pts)
+ assert not np.any(np.isnan(values)), (
+ "domain fault scalar field is NaN somewhere in the model domain -- "
+ "an unconformity region was incorrectly left on it"
+ )
+
+
+class TestExtendFaultTraceToDomain:
+ """`_extend_fault_trace_to_domain` is what makes a domain-boundary fault
+ behave as an "infinite" cut: it adds two synthetic points extending the
+ trace's own trend out to the model's bounding box edges, so the
+ interpolated surface isn't only reliable near the digitised trace.
+ """
+
+ def test_extends_a_short_trace_to_the_domain_edges(self, manager):
+ manager.update_bounding_box(BoundingBox(origin=[0, 0, 0], maximum=[100, 100, 100]))
+ # A short trace entirely inside the domain, running along y=x.
+ trace = pd.DataFrame({'X': [40.0, 60.0], 'Y': [40.0, 60.0], 'Z': [0.0, 0.0]})
+
+ extended = manager._extend_fault_trace_to_domain(trace)
+
+ assert len(extended) == 4
+ xy = extended[['X', 'Y']].to_numpy()
+ # Two of the four points must sit on the bounding box boundary.
+ on_boundary = [
+ np.isclose(x, 0) or np.isclose(x, 100) or np.isclose(y, 0) or np.isclose(y, 100)
+ for x, y in xy
+ ]
+ assert sum(on_boundary) == 2
+
+ def test_extends_a_dipping_trace_keeping_its_trend(self, manager):
+ manager.update_bounding_box(BoundingBox(origin=[0, 0, 0], maximum=[100, 100, 100]))
+ # Trace along y=x with Z increasing 1:1 with distance along the line.
+ trace = pd.DataFrame({'X': [40.0, 60.0], 'Y': [40.0, 60.0], 'Z': [-10.0, 10.0]})
+
+ extended = manager._extend_fault_trace_to_domain(trace)
+
+ synthetic = extended.iloc[2:]
+ # The synthetic points should extrapolate the same Z-vs-along-line
+ # trend rather than reusing the nearest original Z value.
+ assert synthetic['Z'].max() > 10.0
+ assert synthetic['Z'].min() < -10.0
+
+ def test_single_point_trace_is_left_unchanged(self, manager):
+ manager.update_bounding_box(BoundingBox(origin=[0, 0, 0], maximum=[100, 100, 100]))
+ trace = pd.DataFrame({'X': [50.0], 'Y': [50.0], 'Z': [0.0]})
+
+ extended = manager._extend_fault_trace_to_domain(trace)
+
+ assert len(extended) == 1
+
+ def test_direction_parallel_to_axis_outside_domain_is_left_unchanged(self, manager):
+ manager.update_bounding_box(BoundingBox(origin=[0, 0, 0], maximum=[10, 10, 10]))
+ # A vertical (constant-X) trace entirely outside the domain in X.
+ trace = pd.DataFrame({'X': [500.0, 500.0], 'Y': [1.0, 2.0], 'Z': [0.0, 0.0]})
+
+ extended = manager._extend_fault_trace_to_domain(trace)
+
+ assert len(extended) == 2
+
+ def test_local_tangent_varies_along_a_curved_trace(self, manager):
+ """Regression test for a real bug: fitting one global best-fit line
+ through a curved trace (the old approach) flattens its curvature
+ out, and can extrapolate/orient the interpolated surface on the
+ wrong side of real nearby data. Confirmed on a live project where
+ a global-line fit put a stratigraphic unit's own contact data on
+ the opposite side of its domain-boundary fault from a bounding-box
+ corner that a correct local (nearest-segment) classification put
+ on the *same* side. Each point's strike must instead follow its
+ own local tangent.
+ """
+ manager.update_bounding_box(BoundingBox(origin=[0, 0, 0], maximum=[100, 100, 100]))
+ # An L-shaped trace: a horizontal leg then a vertical leg.
+ trace = pd.DataFrame(
+ {'X': [20.0, 50.0, 50.0], 'Y': [50.0, 50.0, 80.0], 'Z': [0.0, 0.0, 0.0]}
+ )
+
+ extended = manager._extend_fault_trace_to_domain(trace)
+
+ strikes = extended['strike'].to_numpy()
+ # Row 0's local tangent is horizontal (towards row 1); row 2's is
+ # vertical (away from row 1) -- these must differ substantially. A
+ # single global best-fit line would instead give every point close
+ # to the same strike.
+ assert abs(((strikes[0] - strikes[2] + 180) % 360) - 180) > 45