fix: interpolator support is rescaled for fault displacement.

Author: lachlangrose

LoopStructural/datatypes/_bounding_box.py

@@ -36,7 +36,8 @@ def __init__(
             origin = np.zeros(3)
         if maximum is None and nsteps is not None and step_vector is not None:
             maximum = origin + nsteps * step_vector
-
+        if origin is not None and global_origin is None:
+            global_origin = origin
         self._origin = np.array(origin)
         self._maximum = np.array(maximum)
         if global_origin is None:
@@ -248,7 +249,9 @@ def with_buffer(self, buffer: float = 0.2) -> BoundingBox:
         # local coordinates, rescale into the original bounding boxes global coordinates
         origin = self.origin - buffer * (self.maximum - self.origin)
         maximum = self.maximum + buffer * (self.maximum - self.origin)
-        return BoundingBox(origin=origin, maximum=maximum, global_origin=self.global_origin)
+        return BoundingBox(
+            origin=origin, maximum=maximum, global_origin=self.global_origin + origin
+        )
 
     def get_value(self, name):
         ix, iy = self.name_map.get(name, (-1, -1))

LoopStructural/modelling/core/geological_model.py

@@ -692,6 +692,7 @@ def create_and_add_foliation(
             interpolatortype=interpolatortype,
             nelements=nelements,
             name=series_surface_data,
+            model=self,
             **kwargs,
         )
         # add data
@@ -706,6 +707,9 @@ def create_and_add_foliation(
         # series_feature = series_builder.build(**kwargs)
         series_feature = series_builder.feature
         series_builder.build_arguments = kwargs
+        # this support is built for the entire model domain? Possibly would
+        # could just pass a regular grid of points - mask by any above unconformities??
+        series_builder.build_arguments['domain'] = True
         series_builder.build_arguments["tol"] = tol
         series_feature.type = FeatureType.INTERPOLATED
         self._add_feature(series_feature)
@@ -746,6 +750,7 @@ def create_and_add_fold_frame(
             name=foldframe_data,
             frame=FoldFrame,
             nelements=nelements,
+            model=self,
             **kwargs,
         )
         # add data
@@ -827,6 +832,7 @@ def create_and_add_folded_foliation(
             fold=fold,
             name=foliation_data,
             svario=svario,
+            model=self,
             **kwargs,
         )
 
@@ -912,6 +918,7 @@ def create_and_add_folded_fold_frame(
             name=fold_frame_data,
             fold=fold,
             frame=FoldFrame,
+            model=self,
             **kwargs,
         )
         fold_frame_builder.add_data_from_data_frame(
@@ -1237,6 +1244,7 @@ def create_and_add_domain_fault(
             interpolatortype=interpolatortype,
             nelements=nelements,
             name=fault_surface_data,
+            model=self,
             **kwargs,
         )
 

LoopStructural/modelling/features/builders/_base_builder.py

@@ -4,7 +4,7 @@
 
 
 class BaseBuilder:
-    def __init__(self, name="Feature"):
+    def __init__(self, model, name: str = "Feature"):
         """Base builder that provides a template for
         implementing different builders.
 
@@ -23,11 +23,16 @@ def __init__(self, name="Feature"):
         If the build arguments are changed, this will flag that the feature needs to be rebuilt
         """
         self._name = name
+        self._model = model
         self._feature = None
         self._up_to_date = False
         self._build_arguments = {}
         self.faults = []
 
+    @property
+    def model(self):
+        return self._model
+
     @property
     def feature(self):
         return self._feature

LoopStructural/modelling/features/builders/_geological_feature_builder.py

@@ -39,6 +39,7 @@ def __init__(
         nelements: int = 1000,
         name="Feature",
         interpolation_region=None,
+        model=None,
         **kwargs,
     ):
         """
@@ -52,7 +53,7 @@ def __init__(
             defining whether the location (xyz) should be included in the
         kwargs - name of the feature, region to interpolate the feature
         """
-        BaseBuilder.__init__(self, name)
+        BaseBuilder.__init__(self, model, name)
         interpolator = InterpolatorFactory.create_interpolator(
             interpolatortype=interpolatortype,
             boundingbox=bounding_box,
@@ -451,16 +452,21 @@ def set_interpolation_geometry(self, origin, maximum, rotation=None):
         while self.interpolator.nx < 100:
             self.interpolator.support.step_vector = self.interpolator.support.step_vector * 0.9
 
-    def check_interpolation_geometry(self, data):
+    def check_interpolation_geometry(self, data, buffer=0.3):
         """Check the interpolation support geometry
-        to data to make sure everything fits"""
+        to data to make sure everything fits
+        Apply the fault to the model grid to ensure that the support
+        is big enough to capture the faulted feature.
+        """
+
         origin = self.interpolator.support.origin
         maximum = self.interpolator.support.maximum
-        origin[origin < np.min(data, axis=0)] = np.min(data, axis=0)[origin < np.min(data, axis=0)]
-        maximum[maximum < np.max(data, axis=0)] = np.max(data, axis=0)[
-            maximum < np.max(data, axis=0)
-        ]
+        pts = self.model.bounding_box.with_buffer(buffer).regular_grid(local=True)
+        for f in self.faults:
+            pts = f.apply_to_points(pts)
 
+        origin[origin > np.min(pts, axis=0)] = np.min(pts, axis=0)[origin > np.min(pts, axis=0)]
+        maximum[maximum < np.max(pts, axis=0)] = np.max(pts, axis=0)[maximum < np.max(pts, axis=0)]
         self.interpolator.support.origin = origin
         self.interpolator.support.maximum = maximum
 
@@ -485,6 +491,9 @@ def build(self, fold=None, fold_weights={}, data_region=None, **kwargs):
         # self.get_interpolator(**kwargs)
         for f in self.faults:
             f.builder.update()
+        domain = kwargs.get("domain", None)
+        if domain:
+            self.check_interpolation_geometry(None)
         self.add_data_to_interpolator(**kwargs)
         if data_region is not None:
             xyz = self.interpolator.get_data_locations()

LoopStructural/modelling/intrusions/intrusion_builder.py

@@ -38,11 +38,10 @@ def __init__(
             a geological interpolator, by default None
         """
 
-        BaseBuilder.__init__(self, name=name)
+        BaseBuilder.__init__(self, model, name=name)
 
         self.intrusion_frame = frame
         self._up_to_date = False
-        self.model = model
         self._feature = IntrusionFeature(
             frame=frame,
             builder=self,