fix: formatting

.gitignore

@@ -131,3 +131,4 @@ docs/source/test/*
 examples/*.png
 dev/scalar_field
 dev/unconf
+docs/source/sg_execution_times.rst

LoopStructural/interpolators/__init__.py

@@ -89,3 +89,4 @@ class InterpolatorType(IntEnum):
     InterpolatorType.BASE_DATA_SUPPORTED: GeologicalInterpolator,
     # InterpolatorType.SURFE: SurfeRBFInterpolator,
 }
+from ._interpolator_factory import InterpolatorFactory

LoopStructural/interpolators/_interpolator_factory.py

@@ -1,6 +1,6 @@
 from typing import Optional
 
-from . import interpolator_map, InterpolatorType, support_interpolator_map
+from . import interpolator_map, InterpolatorType  # , support_interpolator_map
 from LoopStructural.utils import BoundingBox
 from typing import Optional
 import numpy as np
@@ -32,9 +32,9 @@ def create_interpolator(
                 bounding_box=boundingbox,
                 nelements=nelements,
                 element_volume=element_volume,
-            )   
+            )
         return interpolator_map[interpolatortype](support)
- 
+
     @staticmethod
     def from_dict(d):
         d = d.copy()

LoopStructural/interpolators/_p2interpolator.py

@@ -2,6 +2,7 @@
 Piecewise linear interpolator
 """
 import logging
+from typing import Optional
 
 import numpy as np
 
@@ -169,7 +170,7 @@ def add_value_constraints(self, w=1.0):
             )
 
     def minimise_grad_steepness(
-        self, w: float = 0.1, maskall: bool = False, wtfunc: callable = None
+        self, w: float = 0.1, maskall: bool = False, wtfunc: Optional[callable] = None
     ):
         """This constraint minimises the second derivative of the gradient
         mimimising the 2nd derivative should prevent high curvature solutions

LoopStructural/interpolators/supports/__init__.py

@@ -39,6 +39,8 @@ class SupportType(IntEnum):
     SupportType.P2UnstructuredTetMesh: P2UnstructuredTetMesh,
 }
 
+from ._support_factory import SupportFactory
+
 __all__ = [
     "BaseUnstructured2d",
     "P1Unstructured2d",

LoopStructural/modelling/core/geological_model.py

@@ -696,202 +696,6 @@ def set_stratigraphic_column(self, stratigraphic_column, cmap="tab20"):
 
         self.stratigraphic_column = stratigraphic_column
 
-    # ## TODO CREATE SEPARATE API FOR THIS
-    # def get_interpolator(
-    #     self,
-    #     interpolatortype="FDI",
-    #     nelements=1e4,
-    #     buffer=0.2,
-    #     element_volume=None,
-    #     **kwargs,
-    # ):
-    #     """
-    #     Returns an interpolator given the arguments, also constructs a
-    #     support for a discrete interpolator
-
-    #     Parameters
-    #     ----------
-    #     interpolatortype : string
-    #         define the interpolator type
-    #     nelements : int
-    #         number of elements in the interpolator
-    #     buffer : double or numpy array 3x1
-    #         value(s) between 0,1 specifying the buffer around the bounding box
-    #     data_bb : bool
-    #         whether to use the model boundary or the boundary around
-    #     kwargs : no kwargs used, this just catches any additional arguments
-
-    #     Returns
-    #     -------
-    #     interpolator : GeologicalInterpolator
-    #         A geological interpolator
-
-    #     Notes
-    #     -----
-    #     This method will create a geological interpolator for the bounding box of the model. A
-    #     buffer area is added to the interpolation region to avoid boundaries and issues with faults.
-    #     This function wil create a :class:`LoopStructural.interpolators.GeologicalInterpolator` which can either be:
-    #     A discrete interpolator :class:`LoopStructural.interpolators.DiscreteInterpolator`
-
-    #     - 'FDI' :class:`LoopStructural.interpolators.FiniteDifferenceInterpolator`
-    #     - 'PLI' :class:`LoopStructural.interpolators.PiecewiseLinearInterpolator`
-    #     - 'P1'  :class:`LoopStructural.interpolators.P1Interpolator`
-    #     - 'DFI' :class:`LoopStructural.interpolators.DiscreteFoldInterpolator`
-    #     - 'P2'  :class:`LoopStructural.interpolators.P2Interpolator`
-    #     or
-
-    #     - 'surfe'  :class:`LoopStructural.interpolators.SurfeRBFInterpolator`
-
-    #     The discrete interpolators will require a support.
-
-    #     - 'PLI','DFI','P1Interpolator','P2Interpolator' :class:`LoopStructural.interpolators.supports.TetMesh` or you can provide another
-    #       mesh builder which returns :class:`LoopStructural.interpolators.support.UnStructuredTetMesh`
-
-    #     - 'FDI' :class:`LoopStructural.interpolators.supports.StructuredGrid`
-    #     """
-    #     bb = np.copy(self.bounding_box)
-    #     # add a buffer to the interpolation domain, this is necessary for
-    #     # faults but also generally a good
-    #     # idea to avoid boundary problems
-    #     # buffer = bb[1, :]
-    #     buffer = (np.min(bb[1, :] - bb[0, :])) * buffer
-    #     bb[0, :] -= buffer  # *(bb[1,:]-bb[0,:])
-    #     bb[1, :] += buffer  # *(bb[1,:]-bb[0,:])
-    #     box_vol = (bb[1, 0] - bb[0, 0]) * (bb[1, 1] - bb[0, 1]) * (bb[1, 2] - bb[0, 2])
-    #     if interpolatortype == "PLI" and pli:
-    #         if element_volume is None:
-    #             # nelements /= 5
-    #             element_volume = box_vol / nelements
-    #         # calculate the step vector of a regular cube
-    #         step_vector = np.zeros(3)
-    #         step_vector[:] = element_volume ** (1.0 / 3.0)
-    #         # step_vector /= np.array([1,1,2])
-    #         # number of steps is the length of the box / step vector
-    #         nsteps = np.ceil((bb[1, :] - bb[0, :]) / step_vector).astype(int)
-    #         if np.any(np.less(nsteps, 3)):
-    #             axis_labels = ["x", "y", "z"]
-    #             for i in range(3):
-    #                 if nsteps[i] < 3:
-    #                     logger.error(
-    #                         f"Number of steps in direction {axis_labels[i]} is too small, try increasing nelements"
-    #                     )
-    #             logger.error("Cannot create interpolator: number of steps is too small")
-    #             raise ValueError("Number of steps too small cannot create interpolator")
-    #         # create a structured grid using the origin and number of steps
-    #         if self.reuse_supports:
-    #             mesh_id = f"mesh_{nelements}"
-    #             mesh = self.support.get(
-    #                 mesh_id,
-    #                 TetMesh(origin=bb[0, :], nsteps=nsteps, step_vector=step_vector),
-    #             )
-    #             if mesh_id not in self.support:
-    #                 self.support[mesh_id] = mesh
-    #         else:
-    #             if "meshbuilder" in kwargs:
-    #                 mesh = kwargs["meshbuilder"](bb, nelements)
-    #             else:
-    #                 mesh = TetMesh(
-    #                     origin=bb[0, :], nsteps=nsteps, step_vector=step_vector
-    #                 )
-    #         logger.info(
-    #             "Creating regular tetrahedron mesh with %i elements \n"
-    #             "for modelling using PLI" % (mesh.ntetra)
-    #         )
-
-    #         return PLI(mesh)
-    #     if interpolatortype == "P2":
-    #         if element_volume is None:
-    #             # nelements /= 5
-    #             element_volume = box_vol / nelements
-    #         # calculate the step vector of a regular cube
-    #         step_vector = np.zeros(3)
-    #         step_vector[:] = element_volume ** (1.0 / 3.0)
-    #         # step_vector /= np.array([1,1,2])
-    #         # number of steps is the length of the box / step vector
-    #         nsteps = np.ceil((bb[1, :] - bb[0, :]) / step_vector).astype(int)
-    #         if "meshbuilder" in kwargs:
-    #             mesh = kwargs["meshbuilder"](bb, nelements)
-    #         else:
-    #             raise NotImplementedError(
-    #                 "Cannot use P2 interpolator without external mesh"
-    #             )
-    #         logger.info(
-    #             "Creating regular tetrahedron mesh with %i elements \n"
-    #             "for modelling using P2" % (mesh.ntetra)
-    #         )
-    #         return P2Interpolator(mesh)
-    #     if interpolatortype == "FDI":
-    #         # find the volume of one element
-    #         if element_volume is None:
-    #             element_volume = box_vol / nelements
-    #         # calculate the step vector of a regular cube
-    #         step_vector = np.zeros(3)
-    #         step_vector[:] = element_volume ** (1.0 / 3.0)
-    #         # number of steps is the length of the box / step vector
-    #         nsteps = np.ceil((bb[1, :] - bb[0, :]) / step_vector).astype(int)
-    #         if np.any(np.less(nsteps, 3)):
-    #             logger.error("Cannot create interpolator: number of steps is too small")
-    #             axis_labels = ["x", "y", "z"]
-    #             for i in range(3):
-    #                 if nsteps[i] < 3:
-    #                     logger.error(
-    #                         f"Number of steps in direction {axis_labels[i]} is too small, try increasing nelements"
-    #                     )
-    #             raise ValueError("Number of steps too small cannot create interpolator")
-    #         # create a structured grid using the origin and number of steps
-    #         if self.reuse_supports:
-    #             grid_id = "grid_{}".format(nelements)
-    #             grid = self.support.get(
-    #                 grid_id,
-    #                 StructuredGrid(
-    #                     origin=bb[0, :], nsteps=nsteps, step_vector=step_vector
-    #                 ),
-    #             )
-    #             if grid_id not in self.support:
-    #                 self.support[grid_id] = grid
-    #         else:
-    #             grid = StructuredGrid(
-    #                 origin=bb[0, :], nsteps=nsteps, step_vector=step_vector
-    #             )
-    #         logger.info(
-    #             f"Creating regular grid with {grid.n_elements} elements \n"
-    #             "for modelling using FDI"
-    #         )
-    #         return FDI(grid)
-
-    #     if interpolatortype == "DFI" and dfi is True:
-    #         if element_volume is None:
-    #             nelements /= 5
-    #             element_volume = box_vol / nelements
-    #         # calculate the step vector of a regular cube
-    #         step_vector = np.zeros(3)
-    #         step_vector[:] = element_volume ** (1.0 / 3.0)
-    #         # number of steps is the length of the box / step vector
-    #         nsteps = np.ceil((bb[1, :] - bb[0, :]) / step_vector).astype(int)
-    #         # create a structured grid using the origin and number of steps
-    #         if "meshbuilder" in kwargs:
-    #             mesh = kwargs["meshbuilder"].build(bb, nelements)
-    #         else:
-    #             mesh = kwargs.get(
-    #                 "mesh",
-    #                 TetMesh(origin=bb[0, :], nsteps=nsteps, step_vector=step_vector),
-    #             )
-    #         logger.info(
-    #             f"Creating regular tetrahedron mesh with {mesh.ntetra} elements \n"
-    #             "for modelling using DFI"
-    #         )
-    #         return DFI(mesh, kwargs["fold"])
-    #     if interpolatortype == "Surfe" or interpolatortype == "surfe":
-    #         # move import of surfe to where we actually try and use it
-    #         if not surfe:
-    #             logger.warning("Cannot import Surfe, try another interpolator")
-    #             raise ImportError("Cannot import surfepy, try pip install surfe")
-    #         method = kwargs.get("method", "single_surface")
-    #         logger.info("Using surfe interpolator")
-    #         return Surfe(method)
-    #     logger.warning("No interpolator")
-    #     raise InterpolatorError("Could not create interpolator")
-
     def create_and_add_foliation(
         self,
         series_surface_data: str,
@@ -1142,9 +946,7 @@ def create_and_add_folded_fold_frame(
             fold_frame = self.features[-1]
         assert type(fold_frame) == FoldFrame, "Please specify a FoldFrame"
         fold = FoldEvent(fold_frame, name=f"Fold_{fold_frame_data}")
-        # fold_interpolator = self.get_interpolator("DFI", fold=fold, **kwargs)
-        # gy_fold_interpolator = self.get_interpolator("DFI", fold=fold, **kwargs)
-        # frame_interpolator = self.get_interpolator(**kwargs)
+
         interpolatortypes = [
             "DFI",
             "FDI",
@@ -1461,7 +1263,9 @@ def add_onlap_unconformity(
 
         return uc_feature
 
-    def create_and_add_domain_fault(self, fault_surface_data, **kwargs):
+    def create_and_add_domain_fault(
+        self, fault_surface_data, nelements=10000, interpolatortype="FDI", **kwargs
+    ):
         """
         Parameters
         ----------
@@ -1478,10 +1282,14 @@ def create_and_add_domain_fault(self, fault_surface_data, **kwargs):
         * :meth:`LoopStructural.GeologicalModel.get_interpolator`
 
         """
-        interpolator = self.get_interpolator(**kwargs)
         domain_fault_feature_builder = GeologicalFeatureBuilder(
-            interpolator, name=fault_surface_data
+            bounding_box=self.bounding_box,
+            interpolatortype=interpolatortype,
+            nelements=nelements,
+            name=fault_surface_data,
+            **kwargs,
         )
+
         # add data
         unconformity_data = self.data[self.data["feature_name"] == fault_surface_data]
 
@@ -1495,7 +1303,7 @@ def create_and_add_domain_fault(self, fault_surface_data, **kwargs):
         # domain_fault = domain_fault_feature_builder.build(**kwargs)
         domain_fault = domain_fault_feature_builder.feature
         domain_fault_feature_builder.build_arguments = kwargs
-        domain_fault.type = "domain_fault"
+        domain_fault.type = FeatureType.DOMAINFAULT
         self._add_feature(domain_fault)
         self._add_domain_fault_below(domain_fault)
 
@@ -1800,7 +1608,6 @@ def regular_grid(self, nsteps=None, shuffle=True, rescale=False, order="C"):
     def evaluate_model(self, xyz, scale=True):
         """Evaluate the stratigraphic id at each location
 
-
         Parameters
         ----------
         xyz : np.array((N,3),dtype=float)

LoopStructural/modelling/features/__init__.py

@@ -20,6 +20,7 @@ class FeatureType(IntEnum):
     UNCONFORMITY = 7
     INTRUSION = 8
     FAULT = 9
+    DOMAINFAULT = 10
 
 
 from ._base_geological_feature import BaseFeature

setup.py

@@ -32,7 +32,7 @@
 setup(
     name="LoopStructural",
     install_requires=[
-        "numpy>=1.18",  # need to fix numpy to 1.18 because we build against it
+        "numpy>=1.18",
         "pandas",
         "scipy",
         "scikit-image",