fix: isosurfacing of unconformity surfaces will now appear

* fix: disable unconformity with the same name as the feature when isosurfacing that feature.

for example when adding an unconformity to the bottom of a unit at value 0, we want to make sure this scalar field is never < 0 but if we do this we can't isosurface 0 because the cubes containing 0 isovalue have half the nodes being nan. To avoid this, when using the feature.surfaces() method, the region which is related to this feature can be disabled.

* style: style fixes by ruff and autoformatting by black

---------

Co-authored-by: lachlangrose <lachlangrose@users.noreply.github.com>

Author: lachlangrose

LoopStructural/modelling/features/_base_geological_feature.py

@@ -291,11 +291,24 @@ def surfaces(
             if self.model is None:
                 raise ValueError("Must specify bounding box")
             bounding_box = self.model.bounding_box
-        callable = lambda xyz: self.evaluate_value(self.model.scale(xyz))
-        isosurfacer = LoopIsosurfacer(bounding_box, callable=callable)
-        if name is None and self.name is not None:
-            name = self.name
-        return isosurfacer.fit(value, name)
+        regions = self.regions
+        try:
+            self.regions = [
+                r for r in self.regions if r.name != self.name and r.parent.name != self.name
+            ]
+            callable = lambda xyz: self.evaluate_value(self.model.scale(xyz))
+            isosurfacer = LoopIsosurfacer(bounding_box, callable=callable)
+            if name is None and self.name is not None:
+                name = self.name
+            surfaces = isosurfacer.fit(value, name)
+        except Exception as e:
+            logger.error(f"Failed to create surface for {self.name} at value {value}")
+            logger.error(e)
+            surfaces = []
+        finally:
+            self.regions = regions
+
+        return surfaces
 
     def scalar_field(self, bounding_box=None):
         """Create a scalar field for the feature

LoopStructural/modelling/features/_structural_frame.py

@@ -115,7 +115,7 @@ def get_feature(self, i):
         """
         return self.features[i]
 
-    def evaluate_value(self, evaluation_points, ignore_regions=False):
+    def evaluate_value(self, pos, ignore_regions=False):
         """
         Evaluate the value of the structural frame for the points.
         Can optionally only evaluate one coordinate
@@ -129,14 +129,14 @@ def evaluate_value(self, evaluation_points, ignore_regions=False):
         -------
 
         """
-        v = np.zeros(evaluation_points.shape)  # create new 3d array of correct length
+        v = np.zeros(pos.shape)  # create new 3d array of correct length
         v[:] = np.nan
-        v[:, 0] = self.features[0].evaluate_value(evaluation_points, ignore_regions=ignore_regions)
-        v[:, 1] = self.features[1].evaluate_value(evaluation_points, ignore_regions=ignore_regions)
-        v[:, 2] = self.features[2].evaluate_value(evaluation_points, ignore_regions=ignore_regions)
+        v[:, 0] = self.features[0].evaluate_value(pos, ignore_regions=ignore_regions)
+        v[:, 1] = self.features[1].evaluate_value(pos, ignore_regions=ignore_regions)
+        v[:, 2] = self.features[2].evaluate_value(pos, ignore_regions=ignore_regions)
         return v
 
-    def evaluate_gradient(self, evaluation_points, i=None, ignore_regions=False):
+    def evaluate_gradient(self, pos, i=None, ignore_regions=False):
         """
         Evaluate the gradient of the structural frame.
         Can optionally only evaluate the ith coordinate
@@ -151,19 +151,11 @@ def evaluate_gradient(self, evaluation_points, i=None, ignore_regions=False):
 
         """
         if i is not None:
-            return self.features[i].support.evaluate_gradient(
-                evaluation_points, ignore_regions=ignore_regions
-            )
+            return self.features[i].support.evaluate_gradient(pos, ignore_regions=ignore_regions)
         return (
-            self.features[0].support.evaluate_gradient(
-                evaluation_points, ignore_regions=ignore_regions
-            ),
-            self.features[1].support.evaluate_gradient(
-                evaluation_points, ignore_regions=ignore_regions
-            ),
-            self.features[2].support.evaluate_gradient(
-                evaluation_points, ignore_regions=ignore_regions
-            ),
+            self.features[0].support.evaluate_gradient(pos, ignore_regions=ignore_regions),
+            self.features[1].support.evaluate_gradient(pos, ignore_regions=ignore_regions),
+            self.features[2].support.evaluate_gradient(pos, ignore_regions=ignore_regions),
         )
 
     def get_data(self, value_map: Optional[dict] = None) -> List[Union[ValuePoints, VectorPoints]]:

LoopStructural/modelling/features/_unconformity_feature.py

@@ -19,7 +19,7 @@ def __init__(self, feature: GeologicalFeature, value: float, sign=True, onlap=Fa
         # just don't link the regions
         GeologicalFeature.__init__(
             self,
-            name=f"{feature.name}_unconformity",
+            name=f"__{feature.name}_unconformity",
             faults=feature.faults,
             regions=[],  # feature.regions.copy(),  # don't want to share regionsbetween unconformity and # feature.regions,
             builder=feature.builder,
@@ -75,9 +75,9 @@ def evaluate(self, pos: np.ndarray) -> np.ndarray:
             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
         if not self.sign:
-            return self.evaluate_value(pos) > self.value
+            return self.evaluate_value(pos) >= self.value
 
     def __call__(self, pos) -> np.ndarray:
         return self.evaluate(pos)

LoopStructural/modelling/features/fault/_fault_segment.py

@@ -197,7 +197,7 @@ def inside_volume(self, locations, threshold=0.001):
         return np.abs(v) > threshold
         # return np.all(np.logical_and(v > -1,v<1),axis=1)
 
-    def evaluate_value(self, locations):
+    def evaluate_value(self, locations, ignore_regions=False):
         """
         Return the value of the fault surface scalar field
 
@@ -222,7 +222,7 @@ def evaluate_value(self, locations):
         #     except:
         #         logger.error("nan slicing")
         # v[mask] = self.__getitem__(0).evaluate_value(locations[mask, :])
-        return self.__getitem__(0).evaluate_value(locations)
+        return self.__getitem__(0).evaluate_value(locations, ignore_regions=ignore_regions)
 
     def ellipsoid(self):
         pass