Merge pull request #110 from Loop3D/2d_interpolation

2d interpolation

Author: lachlangrose

LoopStructural/__init__.py

@@ -1,6 +1,4 @@
 """
-
-
 LoopStructural API
 =======================
 
@@ -12,6 +10,7 @@
 from pathlib import Path
 from .version import __version__
 
+experimental = False
 ch = logging.StreamHandler()
 formatter = logging.Formatter(
     "%(levelname)s: %(asctime)s: %(filename)s:%(lineno)d -- %(message)s"

LoopStructural/analysis/__init__.py

@@ -4,6 +4,14 @@
 
 Various tools for analysing loopstructural models, including calculating fault intersections and fault toplogies
 """
-from ._fault_displacement import displacement_missfit
-from ._fault_intersection import calculate_fault_intersections
-from ._topology import calculate_fault_topology_matrix
+from LoopStructural.utils import getLogger
+import LoopStructural
+
+logger = getLogger(__name__)
+if LoopStructural.experimental:
+    logger.warning(
+        "LoopStructural.analysis is experimental and may not perform as expected"
+    )
+    from ._fault_displacement import displacement_missfit
+    from ._fault_intersection import calculate_fault_intersections
+    from ._topology import calculate_fault_topology_matrix

LoopStructural/analysis/_topology.py

@@ -4,7 +4,7 @@
 logger = getLogger(__name__)
 
 
-def calculate_fault_topology_matrix(model, xyz=None, threshold=0.001):
+def calculate_fault_topology_matrix(model, xyz=None, threshold=0.001, scale=True):
     """Calculate fault ellipsoid and hw/fw
 
     Parameters
@@ -13,12 +13,18 @@ def calculate_fault_topology_matrix(model, xyz=None, threshold=0.001):
         the model containing the faults
     xyz : np.array
         xyz locations in model coordinates
-
+    threshold : float
+        threshold for determining if point is inside fault volume
+    scale : bool
+        flag whether to rescale xyz to model coordinates
     Returns
     -------
     topology_matrix : np.array
         matrix containing nan (outside), 0 (footwall), 1 (hangingwall)
     """
+    if xyz is not None and scale == True:
+        logger.warning("Scaling XYZ to model coordinate system")
+        xyz = model.scale(xyz, inplace=False)
     if xyz is None:
         xyz = model.regular_grid(rescale=False, shuffle=False)
     topology_matrix = np.zeros((xyz.shape[0], len(model.faults)))

LoopStructural/interpolators/__init__.py

@@ -4,6 +4,11 @@
 """
 from enum import IntEnum
 
+from LoopStructural.utils import getLogger
+import LoopStructural
+
+logger = getLogger(__name__)
+
 
 class InterpolatorType(IntEnum):
     """
@@ -18,21 +23,37 @@ class InterpolatorType(IntEnum):
     FINITE_DIFFERENCE = 2
     DISCRETE_FOLD = 3
     PIECEWISE_LINEAR = 4
+    PIECEWISE_QUADRATIC = 5
     BASE_DATA_SUPPORTED = 10
     SURFE = 11
 
 
-from LoopStructural.interpolators.geological_interpolator import GeologicalInterpolator
-from LoopStructural.interpolators.discrete_interpolator import DiscreteInterpolator
-from LoopStructural.interpolators.structured_tetra import TetMesh
-from LoopStructural.interpolators.unstructured_tetra import UnStructuredTetMesh
-from LoopStructural.interpolators.structured_grid import StructuredGrid
-from LoopStructural.interpolators.finite_difference_interpolator import (
+from LoopStructural.interpolators._geological_interpolator import GeologicalInterpolator
+from LoopStructural.interpolators._discrete_interpolator import DiscreteInterpolator
+from LoopStructural.interpolators.supports import (
+    TetMesh,
+    StructuredGrid,
+    UnStructuredTetMesh,
+    P1Unstructured2d,
+    P2Unstructured2d,
+    StructuredGrid2D,
+    P2UnstructuredTetMesh,
+)
+
+
+from LoopStructural.interpolators._finite_difference_interpolator import (
     FiniteDifferenceInterpolator,
 )
 from LoopStructural.interpolators.piecewiselinear_interpolator import (
     PiecewiseLinearInterpolator,
 )
-from LoopStructural.interpolators.discrete_fold_interpolator import (
+from LoopStructural.interpolators._discrete_fold_interpolator import (
     DiscreteFoldInterpolator,
 )
+
+if LoopStructural.experimental:
+    logger.warning(
+        "Using experimental interpolators: P1Interpolator and P2Interpolator"
+    )
+    from ._p1interpolator import P1Interpolator
+    from ._p2interpolator import P2Interpolator

LoopStructural/interpolators/discrete_fold_interpolator.py renamed to LoopStructural/interpolators/_discrete_fold_interpolator.py

@@ -4,14 +4,15 @@
 import logging
 
 import numpy as np
-from LoopStructural.interpolators.cython.dsi_helper import fold_cg
 
 from LoopStructural.interpolators import PiecewiseLinearInterpolator, InterpolatorType
 
 from LoopStructural.utils import getLogger
 
 logger = getLogger(__name__)
 
+from .cython.dsi_helper import fold_cg
+
 
 class DiscreteFoldInterpolator(PiecewiseLinearInterpolator):
     """ """
@@ -29,7 +30,7 @@ def __init__(self, support, fold):
         """
 
         PiecewiseLinearInterpolator.__init__(self, support)
-        self.type = InterpolatorType.FINITE_DIFFERENCE
+        self.type = InterpolatorType.DISCRETE_FOLD
         self.fold = fold
 
     @classmethod
@@ -126,7 +127,7 @@ def add_fold_constraints(
         ]
         # calculate the fold geometry for the elements barycentre
         deformed_orientation, fold_axis, dgz = self.fold.get_deformed_orientation(
-            self.support.barycentre()
+            self.support.barycentre
         )
         element_idx = np.arange(self.support.n_elements)
         np.random.shuffle(element_idx)
@@ -249,5 +250,5 @@ def add_fold_constraints(
             B = np.zeros(A.shape[0])
             idc = np.array(idc[:ncons, :])
             self.add_constraints_to_least_squares(
-                A, B, fold_regularisation[1], idc, name="fold regularisation 3"
+                A, B, fold_regularisation[2], idc, name="fold regularisation 3"
             )

LoopStructural/interpolators/discrete_interpolator.py renamed to LoopStructural/interpolators/_discrete_interpolator.py

@@ -16,6 +16,8 @@
 
 logger = getLogger(__name__)
 
+from ._geological_interpolator import GeologicalInterpolator
+
 
 class DiscreteInterpolator(GeologicalInterpolator):
     """ """
@@ -70,7 +72,7 @@ def nx(self):
 
     @property
     def region(self):
-        return self.region_function(self.support.nodes)
+        return self.region_function(self.support.nodes).astype(bool)
 
     @property
     def region_map(self):
@@ -166,11 +168,13 @@ def add_constraints_to_least_squares(self, A, B, idc, w=1.0, name="undefined"):
         # logger.debug('Adding constraints to interpolator: {} {} {}'.format(A.shape[0]))
         # print(A.shape,B.shape,idc.shape)
         if A.shape != idc.shape:
+            logger.error("Cannot add constraints: A and indexes have different shape")
             return
 
         if len(A.shape) > 2:
             nr = A.shape[0] * A.shape[1]
-            w = np.tile(w, (A.shape[1]))
+            if isinstance(w, np.ndarray):
+                w = np.tile(w, (A.shape[1]))
             A = A.reshape((A.shape[0] * A.shape[1], A.shape[2]))
             idc = idc.reshape((idc.shape[0] * idc.shape[1], idc.shape[2]))
             B = B.reshape((A.shape[0]))
@@ -184,7 +188,6 @@ def add_constraints_to_least_squares(self, A, B, idc, w=1.0, name="undefined"):
         A[length > 0, :] /= length[length > 0, None]
         if isinstance(w, (float, int)):
             w = np.ones(A.shape[0]) * w
-
         if isinstance(w, np.ndarray) == False:
             raise BaseException("w must be a numpy array")
 
@@ -277,7 +280,7 @@ def add_equality_constraints(self, node_idx, values, name="undefined"):
         gi[self.region] = np.arange(0, self.nx)
         idc = gi[node_idx]
         outside = ~(idc == -1)
-        
+
         self.equal_constraints[name] = {
             "A": np.ones(idc[outside].shape[0]),
             "B": values[outside],
@@ -423,7 +426,7 @@ def build_matrix(self, square=True, damp=0.0, ie=False):
         # can help speed up solving, but might also introduce some errors
 
         if len(self.equal_constraints) > 0:
-            logger.info("Equality block is %i x %i" % (self.eq_const_c, self.nx))
+            logger.info(f"Equality block is {self.eq_const_c} x {self.nx}")
             # solving constrained least squares using
             # | ATA CT | |c| = b
             # | C   0  | |y|   d
@@ -448,7 +451,6 @@ def build_matrix(self, square=True, damp=0.0, ie=False):
                 cols.extend(c["col"].flatten()[~mask].tolist())
 
             C = coo_matrix(
-
                 (np.array(a), (np.array(rows), cols)),
                 shape=(self.eq_const_c, self.nx),
                 dtype=float,
@@ -500,30 +502,6 @@ def _solve_osqp(self, P, A, q, l, u, mkl=False):
             import osqp
         except ImportError:
             raise LoopImportError("Missing osqp pip install osqp")
-        # m = A.shape[0]
-        # n = A.shape[1]
-        # Ad = sparse.random(m, n, density=0.7, format='csc')
-        # b = np.random.randn(m)
-
-        # # OSQP data
-        # P = sparse.block_diag([sparse.csc_matrix((n, n)), sparse.eye(m)], format='csc')
-        # q = np.zeros(n+m)
-        # A = sparse.vstack([
-        #         sparse.hstack([Ad, -sparse.eye(m)]),
-        #         sparse.hstack([sparse.eye(n), sparse.csc_matrix((n, m))])], format='csc')
-        # l = np.hstack([b, np.zeros(n)])
-        # u = np.hstack([b, np.ones(n)])
-
-        # # Create an OSQP object
-        # prob = osqp.OSQP()
-
-        # # Setup workspace
-        # prob.setup(P, q, A, l, u)
-
-        # # Solve problem
-        # res = prob.solve()
-
-        # Create an OSQP object
         prob = osqp.OSQP()
 
         # Setup workspace
@@ -727,6 +705,7 @@ def _solve(self, solver="cg", **kwargs):
             self.c[self.region] = self._solve_chol(A, B)
         if solver == "lu":
             logger.info("Solving using scipy LU")
+            print(self.region)
             self.c[self.region] = self._solve_lu(A, B)
         if solver == "pyamg":
             try:

LoopStructural/interpolators/finite_difference_interpolator.py renamed to LoopStructural/interpolators/_finite_difference_interpolator.py

@@ -6,10 +6,10 @@
 import numpy as np
 
 from LoopStructural.utils.helper import get_vectors
-from .discrete_interpolator import DiscreteInterpolator
+from ._discrete_interpolator import DiscreteInterpolator
 from LoopStructural.interpolators import InterpolatorType
 
-from .operator import Operator
+from ._operator import Operator
 
 from LoopStructural.utils import getLogger
 

LoopStructural/interpolators/geological_interpolator.py renamed to LoopStructural/interpolators/_geological_interpolator.py

@@ -47,8 +47,13 @@ def __init__(self):
         self.valid = False
 
     def __str__(self):
-
-        return self.__str
+        name = f"{self.type} \n"
+        name += f"{self.n_g} gradient points\n"
+        name += f"{self.n_i} interface points\n"
+        name += f"{self.n_n} normal points\n"
+        name += f"{self.n_t} tangent points\n"
+        name += f"{self.n_g + self.n_i + self.n_n + self.n_t} total points\n"
+        return name
 
     def set_region(self, **kwargs):
         pass