refactor: 🎨 black reformat

Author: Lachlan Grose

LoopStructural/__init__.py

@@ -9,6 +9,7 @@
 import tempfile
 from pathlib import Path
 from .version import __version__
+
 experimental = False
 ch = logging.StreamHandler()
 formatter = logging.Formatter(

LoopStructural/analysis/__init__.py

@@ -6,9 +6,12 @@
 """
 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")
+    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

@@ -22,9 +22,9 @@ def calculate_fault_topology_matrix(model, xyz=None, threshold=0.001, scale=True
     topology_matrix : np.array
         matrix containing nan (outside), 0 (footwall), 1 (hangingwall)
     """
-    if xyz is not None and scale==True:
+    if xyz is not None and scale == True:
         logger.warning("Scaling XYZ to model coordinate system")
-        xyz = model.scale(xyz,inplace=False)
+        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

@@ -6,7 +6,10 @@
 
 from LoopStructural.utils import getLogger
 import LoopStructural
+
 logger = getLogger(__name__)
+
+
 class InterpolatorType(IntEnum):
     """
     Enum for the different interpolator types
@@ -27,7 +30,15 @@ class InterpolatorType(IntEnum):
 
 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.supports import (
+    TetMesh,
+    StructuredGrid,
+    UnStructuredTetMesh,
+    P1Unstructured2d,
+    P2Unstructured2d,
+    StructuredGrid2D,
+    P2UnstructuredTetMesh,
+)
 
 
 from LoopStructural.interpolators._finite_difference_interpolator import (
@@ -39,7 +50,10 @@ class InterpolatorType(IntEnum):
 from LoopStructural.interpolators._discrete_fold_interpolator import (
     DiscreteFoldInterpolator,
 )
+
 if LoopStructural.experimental:
-    logger.warning("Using experimental interpolators: P1Interpolator and P2Interpolator")
+    logger.warning(
+        "Using experimental interpolators: P1Interpolator and P2Interpolator"
+    )
     from ._p1interpolator import P1Interpolator
     from ._p2interpolator import P2Interpolator

LoopStructural/interpolators/_discrete_interpolator.py

@@ -173,7 +173,7 @@ def add_constraints_to_least_squares(self, A, B, idc, w=1.0, name="undefined"):
 
         if len(A.shape) > 2:
             nr = A.shape[0] * A.shape[1]
-            if isinstance(w,np.ndarray):
+            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]))
@@ -280,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],
@@ -451,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,
@@ -502,7 +501,7 @@ def _solve_osqp(self, P, A, q, l, u, mkl=False):
         try:
             import osqp
         except ImportError:
-            raise LoopImportError("Missing osqp pip install osqp")      
+            raise LoopImportError("Missing osqp pip install osqp")
         prob = osqp.OSQP()
 
         # Setup workspace

LoopStructural/interpolators/_geological_interpolator.py

@@ -47,12 +47,12 @@ def __init__(self):
         self.valid = False
 
     def __str__(self):
-        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'
+        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):

LoopStructural/interpolators/_p1interpolator.py

@@ -46,7 +46,9 @@ def add_gradient_ctr_pts(self, w=1.0):
     def add_norm_ctr_pts(self, w=1.0):
         points = self.get_norm_constraints()
         if points.shape[0] > 0:
-            grad, elements, inside = self.support.evaluate_shape_derivatives(points[:, :3])
+            grad, elements, inside = self.support.evaluate_shape_derivatives(
+                points[:, :3]
+            )
             size = self.support.element_size[inside]
             wt = np.ones(size.shape[0])
             wt *= w * size
@@ -123,5 +125,3 @@ def minimize_edge_jumps(
             name="edge jump",
         )
         # p2.add_constraints_to_least_squares(const_cp2*e_len[:,None]*w,np.zeros(const_cp1.shape[0]),tri_cp2, name='edge jump cp2')
-
-   

LoopStructural/interpolators/_p2interpolator.py

@@ -40,6 +40,7 @@ def __init__(self, mesh):
             "tpw": 1.0,
             "ipw": 1.0,
         }
+
     def _setup_interpolator(self, **kwargs):
         """
         Searches through kwargs for any interpolation weights and updates
@@ -65,17 +66,19 @@ def _setup_interpolator(self, **kwargs):
             self.interpolation_weights[key] = kwargs[key]
         if self.interpolation_weights["cgw"] > 0.0:
             self.up_to_date = False
-            self.minimise_edge_jumps(
-                 self.interpolation_weights["cgw"])
+            self.minimise_edge_jumps(self.interpolation_weights["cgw"])
             #     direction_feature=kwargs.get("direction_feature", None),
             #     direction_vector=kwargs.get("direction_vector", None),
             # )
-            self.minimise_grad_steepness(w=self.interpolation_weights.get('steepness_weight',.01),wtfunc=self.interpolation_weights.get('steepness_wtfunc',None))
+            self.minimise_grad_steepness(
+                w=self.interpolation_weights.get("steepness_weight", 0.01),
+                wtfunc=self.interpolation_weights.get("steepness_wtfunc", None),
+            )
             logger.info(
                 "Using constant gradient regularisation w = %f"
                 % self.interpolation_weights["cgw"]
             )
-        
+
         logger.info(
             "Added %i gradient constraints, %i normal constraints,"
             "%i tangent constraints and %i value constraints"
@@ -87,7 +90,6 @@ def _setup_interpolator(self, **kwargs):
         self.add_tangent_constraints(self.interpolation_weights["tpw"])
         # self.add_interface_constraints(self.interpolation_weights["ipw"])
 
-
     def copy(self):
         return P2Interpolator(self.support)
 
@@ -103,8 +105,9 @@ def add_gradient_constraints(self, w=1.0):
             B = np.zeros(A.shape[0])
             elements = self.support[elements[inside]]
             self.add_constraints_to_least_squares(
-                A*wt[:,None], B, elements, name="gradient")
-         
+                A * wt[:, None], B, elements, name="gradient"
+            )
+
     def add_gradient_orthogonal_constraints(self, points, vector, w=1.0, B=0):
         """
         constraints scalar field to be orthogonal to a given vector
@@ -130,7 +133,8 @@ def add_gradient_orthogonal_constraints(self, points, vector, w=1.0, B=0):
             B = np.zeros(A.shape[0])
             elements = self.support.elements[elements[inside]]
             self.add_constraints_to_least_squares(
-                A*wt[:,None], B, elements, name="gradient orthogonal")
+                A * wt[:, None], B, elements, name="gradient orthogonal"
+            )
 
     def add_norm_constraints(self, w=1.0):
         points = self.get_norm_constraints()
@@ -156,16 +160,14 @@ def add_norm_constraints(self, w=1.0):
 
         pass
 
-     
-
     def add_value_constraints(self, w=1.0):
         points = self.get_value_constraints()
         if points.shape[0] > 1:
             N, elements, mask = self.support.evaluate_shape(points[:, :3])
             # mask = elements > 0
             size = self.support.element_size[elements[mask]]
             wt = np.ones(size.shape[0])
-            wt *= w 
+            wt *= w
             self.add_constraints_to_least_squares(
                 N[mask, :] * wt[:, None],
                 points[mask, 3] * wt,
@@ -197,10 +199,10 @@ def minimise_grad_steepness(self, w=0.1, maskall=True, wtfunc=None):
         d2 = self.support.evaluate_shape_d2(elements[mask])
         # d2 is [ele_idx, deriv, node]
         wt = np.ones(d2.shape[0])
-        
-        wt *= w #* self.support.element_size[mask]
+
+        wt *= w  # * self.support.element_size[mask]
         if callable(wtfunc):
-            logger.info('Using function to weight gradient steepness')
+            logger.info("Using function to weight gradient steepness")
             wt = wtfunc(self.support.barycentre) * self.support.element_size[mask]
         idc = self.support.elements[elements[mask]]
         for i in range(d2.shape[1]):
@@ -210,36 +212,36 @@ def minimise_grad_steepness(self, w=0.1, maskall=True, wtfunc=None):
                 idc[:, :],
                 name=f"grad_steepness_{i}",
             )
- 
+
     def minimise_edge_jumps(
         self, w=0.1, wtfunc=None, vector_func=None
     ):  # NOTE: imposes \phi_T1(xi)-\phi_T2(xi) dot n =0
         # iterate over all triangles
-        
+
         cp, weight = self.support.get_quadrature_points(3)
 
         norm = self.support.shared_element_norm
         shared_element_size = self.support.shared_element_size
-        
+
         # evaluate normal if using vector func for cp1
         for i in range(cp.shape[1]):
             if callable(vector_func):
-                norm = vector_func(cp[:,i,:])
+                norm = vector_func(cp[:, i, :])
             # evaluate the shape function for the edges for each neighbouring triangle
             cp_Dt, cp_tri1 = self.support.evaluate_shape_derivatives(
-                cp[:,i, :], elements=self.support.shared_element_relationships[:, 0]
+                cp[:, i, :], elements=self.support.shared_element_relationships[:, 0]
             )
             cp_Dn, cp_tri2 = self.support.evaluate_shape_derivatives(
-                cp[:,i, :], elements=self.support.shared_element_relationships[:, 1]
+                cp[:, i, :], elements=self.support.shared_element_relationships[:, 1]
             )
             # constraint for each cp is triangle - neighbour create a Nx12 matrix
             const_t_cp = np.einsum("ij,ijk->ik", norm, cp_Dt)
             const_n_cp = -np.einsum("ij,ijk->ik", norm, cp_Dn)
-           
+
             const_cp = np.hstack([const_t_cp, const_n_cp])
             tri_cp = np.hstack(
-                        [self.support.elements[cp_tri1], self.support.elements[cp_tri2]]
-                    )
+                [self.support.elements[cp_tri1], self.support.elements[cp_tri2]]
+            )
             wt = np.zeros(tri_cp.shape[0])
             wt[:] = w
             if wtfunc:
@@ -250,9 +252,6 @@ def minimise_edge_jumps(
                 tri_cp,
                 name=f"shared element jump cp{i}",
             )
-        
-
-    
 
     def evaluate_d2(self, evaluation_points):
         evaluation_points = np.array(evaluation_points)

LoopStructural/interpolators/supports/_2d_base_unstructured.py

@@ -23,7 +23,6 @@ def __init__(self, elements, vertices, neighbours):
         self.n_nodes = self.nx
         self.neighbours = neighbours
 
-        
         self.properties = {}
         # build an array of edges and edge relationships
         self.edges = np.zeros((self.nelements * 3, 2), dtype=int)
@@ -55,6 +54,7 @@ def ncps(self):
         Returns the number of nodes for an element in the mesh
         """
         return self.elements.shape[1]
+
     @property
     def nodes(self):
         """
@@ -66,6 +66,7 @@ def nodes(self):
             Fortran ordered
         """
         return self.vertices
+
     @property
     def barycentre(self):
         """