(BUG) issues building folded fold frame

- trying to solve by being creative with how to build the fold frame

Author: lachlangrose

LoopStructural/interpolators/piecewiselinear_interpolator.py

@@ -72,7 +72,10 @@ def _setup_interpolator(self, **kwargs):
             self.interpolation_weights[key] = kwargs[key]
         if self.interpolation_weights['cgw'] > 0.:
             self.up_to_date = False
-            self.add_constant_gradient(self.interpolation_weights['cgw'])
+            self.add_constant_gradient(self.interpolation_weights['cgw'],
+            direction_feature=kwargs.get('direction_feature',None),
+            direction_vector=kwargs.get('direction_vector',None)
+            )
             logger.info("Using constant gradient regularisation w = %f"
                         %self.interpolation_weights['cgw'])
         logger.info("Added %i gradient constraints, %i normal constraints,"
@@ -87,7 +90,7 @@ def _setup_interpolator(self, **kwargs):
         if 'constant_norm' in kwargs:
             self.add_constant_norm(w=kwargs['constant_norm'])
 
-    def add_constant_gradient(self, w=0.1):
+    def add_constant_gradient(self, w= 0.1, direction_vector=None, direction_feature=None):
         """
         Add the constant gradient regularisation to the system
 
@@ -99,8 +102,59 @@ def add_constant_gradient(self, w=0.1):
         -------
 
         """
+        if direction_feature is not None:
+            print('dir fe')
+            direction_vector = direction_feature.evaluate_gradient(self.support.barycentre())
+        if direction_vector is not None:
+            if direction_vector.shape[0] == 1:
+                # if using a constant direction, tile array so it works for cg calc
+                direction_vector = np.tile(direction_vector,(self.support.barycentre().shape[0],1))
+
+            
+        # iterate over all elements
+        A, idc, B = self.support.get_constant_gradient(region=self.region,direction=direction_vector)
+        A = np.array(A)
+        B = np.array(B)
+        idc = np.array(idc)
+
+        gi = np.zeros(self.support.n_nodes)
+        gi[:] = -1
+        gi[self.region] = np.arange(0, self.nx)
+        idc = gi[idc]
+        outside = ~np.any(idc == -1, axis=1)
+
+        # w/=A.shape[0]
+        self.add_constraints_to_least_squares(A[outside, :] * w,
+                                              B[outside] * w, idc[outside, :],
+                                              name='regularisation')
+        return
+
+    def add_direction_constant_gradient(self, w= 0.1, direction_vector=None, direction_feature=None):
+        """
+        Add the constant gradient regularisation to the system where regularisation is projected
+        on a vector
+
+        Parameters
+        ----------
+        w (double) - weighting of the cg parameter
+        direction_vector
+        direction_feature
+
+        Returns
+        -------
+
+        """
+        if direction_feature:
+            print('dir fe')
+            direction_vector = direction_feature.evaluate_gradient(self.support.barycentre())
+        if direction_vector:
+            if direction_vector.shape[0] == 1:
+                # if using a constant direction, tile array so it works for cg calc
+                direction_vector = np.tile(direction_vector,(self.support.barycentre().shape[0],1))
+
+            
         # iterate over all elements
-        A, idc, B = self.support.get_constant_gradient(region=self.region)
+        A, idc, B = self.support.get_constant_gradient(region=self.region,direction=direction_vector)
         A = np.array(A)
         B = np.array(B)
         idc = np.array(idc)
@@ -116,6 +170,8 @@ def add_constant_gradient(self, w=0.1):
                                               B[outside] * w, idc[outside, :],
                                               name='regularisation')
         return
+
+
     def add_constant_norm(self, w=0.1):
         """
         Add the constant gradient regularisation to the system
@@ -128,7 +184,6 @@ def add_constant_norm(self, w=0.1):
         -------
 
         """
-        print('adding gradient norm')
         # iterate over all elements
         A, idc, B = self.support.get_constant_norm(region=self.region)
         A = np.array(A)
@@ -358,6 +413,7 @@ def add_gradient_orthogonal_constraint(self, points, vector, w=1.0,
             vertices, element_gradients, tetras, inside = self.support.get_tetra_gradient_for_location(points[:,:3])
             #e, inside = self.support.elements_for_array(points[:, :3])
             #nodes = self.support.nodes[self.support.elements[e]]
+            vector /= np.linalg.norm(vector,axis=1)[:,None]
             vecs = vertices[:, 1:, :] - vertices[:, 0, None, :]
             vol = np.abs(np.linalg.det(vecs))  # / 6
             # d_t = self.support.get_elements_gradients(e)
@@ -373,7 +429,7 @@ def add_gradient_orthogonal_constraint(self, points, vector, w=1.0,
             gi[self.region] = np.arange(0, self.nx).astype(int)
             w /= 3
             idc = gi[tetras]
-            B = np.zeros(idc.shape[0])
+            B = np.zeros(idc.shape[0])+B
             outside = ~np.any(idc == -1, axis=1)
             self.add_constraints_to_least_squares(A[outside, :] * w,
                                                   B[outside], idc[outside, :])

LoopStructural/interpolators/structured_tetra.py

@@ -4,7 +4,7 @@
 import logging
 
 import numpy as np
-from LoopStructural.interpolators.cython.dsi_helper import cg, constant_norm
+from LoopStructural.interpolators.cython.dsi_helper import cg, constant_norm, fold_cg
 
 logger = logging.getLogger(__name__)
 
@@ -233,7 +233,7 @@ def get_tetra_for_location(self, pos):
         tetra_return[inside,:] = tetras[mask,:]
         return vertices_return, c_return, tetra_return, inside
 
-    def get_constant_gradient(self, region):
+    def get_constant_gradient(self, region, direction=None):
         """
         Get the constant gradient for the specified nodes
 
@@ -246,6 +246,34 @@ def get_constant_gradient(self, region):
         -------
 
         """
+        """
+        Add the constant gradient regularisation to the system
+
+        Parameters
+        ----------
+        w (double) - weighting of the cg parameter
+
+        Returns
+        -------
+
+        """
+        if direction is not None:
+            print('using cg direction')
+            logger.info("Running constant gradient")
+            elements_gradients = self.get_element_gradients(np.arange(self.ntetra))
+            if elements_gradients.shape[0] != direction.shape[0]:
+                logger.error('Cannot add directional CG, vector field is not the correct length')
+                return
+            region = region.astype('int64')
+
+            neighbours = self.get_neighbours()
+            elements = self.get_elements()
+            idc, c, ncons = fold_cg(elements_gradients, direction, neighbours.astype('int64'), elements.astype('int64'), self.nodes)
+
+            idc = np.array(idc[:ncons, :])
+            c = np.array(c[:ncons, :])
+            B = np.zeros(c.shape[0])
+            return c, idc, B
         if self.cg is None:
             logger.info("Running constant gradient")
             elements_gradients = self.get_element_gradients(np.arange(self.ntetra))

LoopStructural/modelling/core/geological_model.py

@@ -718,8 +718,10 @@ def create_and_add_folded_fold_frame(self, fold_frame_data,
         assert type(fold_frame) == FoldFrame, "Please specify a FoldFrame"
         fold = FoldEvent(fold_frame,name='Fold_{}'.format(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)
-        interpolators = [fold_interpolator, frame_interpolator,
+        interpolators = [fold_interpolator, gy_fold_interpolator,
                          frame_interpolator.copy()]
         fold_frame_builder = StructuralFrameBuilder(
             interpolators=interpolators, name=fold_frame_data, **kwargs)

LoopStructural/modelling/features/cross_product_geological_feature.py

@@ -34,7 +34,7 @@ def __init__(self, name, geological_feature_a, geological_feature_b):
         super().__init__(name, None)
         self.geological_feature_a = geological_feature_a
         self.geological_feature_b = geological_feature_b
-
+        self.value_feature = None
     def evaluate_gradient(self, locations):
         """
         Calculate the gradient of the geological feature by using numpy to
@@ -49,8 +49,11 @@ def evaluate_gradient(self, locations):
         -------
 
         """
-        return np.cross(self.geological_feature_a.evaluate_gradient(locations),
-                        self.geological_feature_b.evaluate_gradient(locations),
+        v1 = self.geological_feature_a.evaluate_gradient(locations)
+        # v1 /= np.linalg.norm(v1,axis=1)[:,None]
+        v2 = self.geological_feature_b.evaluate_gradient(locations)
+        # v2 /= np.linalg.norm(v2,axis=1)[:,None]
+        return np.cross(v1,v2,
                         axisa=1,
                         axisb=1)
 
@@ -65,13 +68,22 @@ def evaluate_value(self, evaluation_points):
         -------
 
         """
-        return np.zeros(evaluation_points.shape[0])
+        values = np.zeros(evaluation_points.shape[0])
+        if self.value_feature is not None:
+            values[:] = self.value_feature.evaluate_value(evaluation_points)
+        return values
 
     def mean(self):
+        if self.value_feature:
+            return self.value_feature.mean()
         return 0.
 
     def min(self):
+        if self.value_feature:
+            return self.value_feature.min()
         return 0.
 
     def max(self):
+        if self.value_feature:
+            return self.value_feature.max()
         return 0.

LoopStructural/modelling/features/structural_frame_builder.py

@@ -128,7 +128,7 @@ def build(self, w1=1., w2=1., w3=1., frame=StructuralFrame, **kwargs):
                                                     0], **kwargs)
             # remove fold from kwargs
 
-            fold = kwargs.pop('fold', False)
+            fold = kwargs.pop('fold', None)
         if gx_feature is None:
             logger.warning(
                 "Not enough constraints for structural frame coordinate 0, \n"
@@ -159,7 +159,22 @@ def build(self, w1=1., w2=1., w3=1., frame=StructuralFrame, **kwargs):
                 self.builders[1].add_orthogonal_feature(gx_feature, gxxgy,step=step)
             if gz_feature is not None and gyxgz>0:
                 self.builders[1].add_orthogonal_feature(gz_feature, gyxgz,step=step)
-            gy_feature = self.builders[1].build(regularisation=regularisation[1], **kwargs)
+            gy_const_norm = kwargs.get('gy_const_norm',0.)
+
+            ## bit of an ugly hack, adding in norm constraints for the norm we are forcing 
+            tmp = CrossProductGeologicalFeature('tmp',gx_feature,gz_feature)
+            self.builders[1].add_orthogonal_feature(tmp,10.,step=step,B=1)
+            # gy_feature.value_feature = fold.fold frame[0]
+            # vector = feature.evaluate_gradient(self.builders[1].interpolator.support.barycentre())
+            # vector /= np.linalg.norm(vector,axis=1)[:,None]
+            # element_idx = np.arange(self.builders[1].interpolator.support.n_elements)
+            # np.random.shuffle(element_idx)
+            # norm_pts = np.hstack([self.builders[1].interpolator.support.barycentre()[element_idx[::step],:],vector[element_idx[::step],:],np.ones((vector[element_idx[::step],:].shape[0],1))])
+
+            # self.builders[1].interpolator.set_normal_constraints(norm_pts)
+            self.builders[1].data_added=True
+            
+            gy_feature = self.builders[1].build(regularisation=regularisation[1],**kwargs)
 
         if gy_feature is None:
             logger.warning(

LoopStructural/modelling/fold/fold.py

@@ -26,7 +26,7 @@ def __init__(self, foldframe, fold_axis_rotation=None, fold_limb_rotation=None,
         self.fold_limb_rotation = fold_limb_rotation
         self.fold_axis = fold_axis
         self.name = name
-
+        
     def get_fold_axis_orientation(self, points):
         """
         gets the fold axis orientation for evaluation points