From 5c47ef6020b2d2c1da53e771db5dac1308d870db Mon Sep 17 00:00:00 2001
From: Lachlan Grose <lachlan.grose@gmail.com>
Date: Mon, 1 Aug 2022 17:15:39 +1000
Subject: [PATCH] fix: basal unconformities cropping lowest surface

---
 .../modelling/core/geological_model.py         | 15 +++++++++------
 .../features/_unconformity_feature.py          | 18 ++++++++++++++----
 2 files changed, 23 insertions(+), 10 deletions(-)

diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py
index b5f645400..a01621cfb 100644
--- a/LoopStructural/modelling/core/geological_model.py
+++ b/LoopStructural/modelling/core/geological_model.py
@@ -115,6 +115,7 @@ def __init__(
         reuse_supports=False,
         logfile=None,
         loglevel="info",
+        epsilon=0.04,
     ):
         """
         Parameters
@@ -125,7 +126,8 @@ def __init__(
             specifying the maximum extent of the model
         rescale : bool
             whether to rescale the model to between 0/1
-
+        epsion : float
+            a fudge factor for isosurfacing, used to make sure surfaces appear
         Examples
         --------
         Demo data
@@ -184,6 +186,7 @@ def __init__(
             )
 
         self.bounding_box /= self.scale_factor
+        self.epsilon = epsilon * float(np.max(lengths)) / self.scale_factor
         self.support = {}
         self.reuse_supports = reuse_supports
         if self.reuse_supports:
@@ -1252,7 +1255,7 @@ def _add_unconformity_above(self, feature):
 
         for f in reversed(self.features):
             if f.type == FeatureType.UNCONFORMITY and f.name != feature.name:
-                feature.add_region(lambda pos: f.evaluate(pos))
+                feature.add_region(f)
                 break
 
     def create_and_add_unconformity(self, unconformity_surface_data, **kwargs):
@@ -1327,7 +1330,7 @@ def add_unconformity(
             logger.debug(f"Adding {uc_feature.name} as unconformity to {f.name}")
             if f.type == FeatureType.FAULT:
                 continue
-            f.add_region(lambda pos: ~uc_feature.evaluate(pos))
+            f.add_region(uc_feature.inverse())
         # now add the unconformity to the feature list
         self._add_feature(uc_feature)
         return uc_feature
@@ -1353,12 +1356,12 @@ def add_onlap_unconformity(
         """
 
         uc_feature = UnconformityFeature(feature, value, True)
-        feature.add_region(lambda pos: ~uc_feature.evaluate(pos))
+        feature.add_region(uc_feature.inverse())
         for f in reversed(self.features):
             if f.type == FeatureType.UNCONFORMITY:
                 break
             if f != feature:
-                f.add_region(lambda pos: uc_feature.evaluate(pos))
+                f.add_region(uc_feature)
         self._add_feature(uc_feature)
 
         return uc_feature
@@ -1593,7 +1596,7 @@ def create_and_add_fault(
 
         for f in reversed(self.features):
             if f.type == FeatureType.UNCONFORMITY:
-                fault.add_region(lambda pos: f.evaluate_value(pos) <= 0)
+                fault.add_region(f)
                 break
         if displacement == 0:
             fault.type = "fault_inactive"
diff --git a/LoopStructural/modelling/features/_unconformity_feature.py b/LoopStructural/modelling/features/_unconformity_feature.py
index 40b0ac9a8..33c30967d 100644
--- a/LoopStructural/modelling/features/_unconformity_feature.py
+++ b/LoopStructural/modelling/features/_unconformity_feature.py
@@ -1,6 +1,8 @@
 from LoopStructural.modelling.features import GeologicalFeature
 from LoopStructural.modelling.features import FeatureType
 
+import numpy as np
+
 
 class UnconformityFeature(GeologicalFeature):
     """ """
@@ -28,7 +30,12 @@ def __init__(self, feature: GeologicalFeature, value: float, sign=True):
         self.type = FeatureType.UNCONFORMITY
         self.sign = sign
 
-    def evaluate(self, pos):
+    def inverse(self):
+        uc = UnconformityFeature(self, self.value, sign=not self.sign)
+        uc.name = self.name + "_inverse"
+        return uc
+
+    def evaluate(self, pos: np.ndarray) -> np.ndarray:
         """
 
         Parameters
@@ -38,10 +45,13 @@ def evaluate(self, pos):
 
         Returns
         -------
-        boolean
+        np.ndarray.dtype(bool)
             true if above the unconformity, false if below
         """
         if self.sign:
-            return self.evaluate_value(pos) < self.value
+            return self.evaluate_value(pos) < self.value + self.model.epsilon
         if not self.sign:
-            return self.evaluate_value(pos) > self.value
+            return self.evaluate_value(pos) > self.value - self.model.epsilon
+
+    def __call__(self, pos) -> np.ndarray:
+        return self.evaluate(pos)