fix: bugfix for gradient constraints

Gradient constraint should be vector dot element matrix * x = 0
The vector should be a unit vector.
The vector was not being normalised and the element gradient matrix
was being scaled incorrectly. This problem was impacting both pli and
fdi for gradient orthogonal constraints and gradient constraints.

LoopStructural/interpolators/discrete_interpolator.py

@@ -294,7 +294,7 @@ def add_tangent_constraints(self, w=1.0):
         """
         points = self.get_tangent_constraints()
         if points.shape[0] > 1:
-            self.add_gradient_orthogonal_constraint(points[:,:3],points[:,3:6],w)
+            self.add_gradient_orthogonal_constraints(points[:,:3],points[:,3:6],w)
 
     def build_matrix(self, square=True, damp=True):
         """

LoopStructural/interpolators/finite_difference_interpolator.py

@@ -263,11 +263,13 @@ def add_gradient_constraints(self, w=1.):
             inside = np.logical_and(~np.any(idc == -1, axis=1), inside)
 
             T = self.support.calcul_T(points[inside, :3])
-            norm = np.linalg.norm(T,axis=2)
-            T/=norm[:,:,None]
+            # normalise constraint vector and scale element matrix by this
+            norm = np.linalg.norm(vector,axis=1)
+            vector/=norm[:,None]
+            T/=norm[:,None,None]
+            # calculate two orthogonal vectors to constraint (strike and dip vector)
             strike_vector, dip_vector = get_vectors(points[inside, 3:6])
             A = np.einsum('ij,ijk->ik', strike_vector.T, T)
-
             B = np.zeros(points[inside, :].shape[0])
             self.add_constraints_to_least_squares(A * w, B, idc[inside, :], name='gradient')
             A = np.einsum('ij,ijk->ik', dip_vector.T, T)
@@ -350,17 +352,15 @@ def add_gradient_orthogonal_constraints(self, points, vector, w=1.0,
 
             idc[inside, :] = gi[node_idx[inside, :]]
             inside = np.logical_and(~np.any(idc == -1, axis=1), inside)
-
+            # normalise vector and scale element gradient matrix by norm as well
+            norm = np.linalg.norm(vector,axis=1)
+            vector/=norm[:,None]
             #normalise element vector to unit vector for dot product
             T = self.support.calcul_T(points[inside, :3])
-            # norm = np.linalg.norm(T,axis=1)
-            # T/=norm[:,None,:]
-            # normalise vector to unit vector for dot product
-            # nro
-            # # vector[inside,:3] /= np.linalg.norm(vector[inside,:3],axis=1)[:,None]
-            # dot product of vector and element gradient 
-            A = np.einsum('ij,ijk->ik', vector[inside, :3], T)
+            T/=norm[:,None,None]
 
+            # dot product of vector and element gradient = 0
+            A = np.einsum('ij,ijk->ik', vector[inside, :3], T)
             B = np.zeros(points[inside, :].shape[0])
             self.add_constraints_to_least_squares(A * w, B, idc[inside, :], name='gradient orthogonal')
 

LoopStructural/interpolators/piecewiselinear_interpolator.py

@@ -200,9 +200,9 @@ def add_gradient_constraints(self, w=1.0):
             #nodes = self.support.nodes[self.support.elements[e]]
             vecs = vertices[:, 1:, :] - vertices[:, 0, None, :]
             vol = np.abs(np.linalg.det(vecs))  # / 6
-            # d_t = self.support.get_elements_gradients(e)
-            norm = np.linalg.norm(element_gradients, axis=2)
-            element_gradients /= norm[:, :, None]
+            norm = np.linalg.norm(vector,axis=1)
+            vector/=norm[:,None]
+            element_gradients /= norm[:, None, None]
             # d_t *= vol[:,None,None]
             strike_vector, dip_vector = get_vectors(points[:, 3:6])
             A = np.einsum('ji,ijk->ik', strike_vector, element_gradients)
@@ -221,7 +221,6 @@ def add_gradient_constraints(self, w=1.0):
                                                   name = 'gradient')
             A = np.einsum('ji,ijk->ik', dip_vector, element_gradients)
             A *= vol[:, None]
-            # A+=A2
             self.add_constraints_to_least_squares(A[outside, :] * w,
                                           B[outside], idc[outside, :],
                                                   name='gradient')
@@ -403,7 +402,7 @@ def add_interface_constraints(self, w=1.0):  # for now weight all value points t
         #                                               np.zeros(interface_A[outside,:].shape[0]),
         #                                               interface_idc[outside, :], name='interface_{}'.format(unique_id))
 
-    def add_gradient_orthogonal_constraint(self, points, vector, w=1.0,
+    def add_gradient_orthogonal_constraints(self, points, vector, w=1.0,
                                            B=0):
         """
         constraints scalar field to be orthogonal to a given vector
@@ -423,12 +422,11 @@ 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]
+            norm = np.linalg.norm(vector,axis=1)
+            vector /= norm[:,None]
             vecs = vertices[:, 1:, :] - vertices[:, 0, None, :]
             vol = np.abs(np.linalg.det(vecs))  # / 6
-            # d_t = self.support.get_elements_gradients(e)
-            norm = np.linalg.norm(element_gradients, axis=2)
-            element_gradients /= norm[:, :, None]
+            element_gradients /= norm[:,None,None]
 
             A = np.einsum('ij,ijk->ik', vector, element_gradients)
 

LoopStructural/modelling/features/geological_feature_builder.py

@@ -282,7 +282,7 @@ def install_gradient_constraint(self):
             vector[norm>0] /= norm[norm>0,None]
             element_idx = np.arange(self.interpolator.support.n_elements)
             np.random.shuffle(element_idx)
-            self.interpolator.add_gradient_orthogonal_constraint(
+            self.interpolator.add_gradient_orthogonal_constraints(
                 self.interpolator.support.barycentre()[element_idx[::step],:],
                 vector[element_idx[::step],:],
                 w=w,