fix: remove element volume scaling.

element volume scaling conflicts with row normalisation. Maybe need to look into this more.
There is an issue where the model's ability to solve is dependent on the

Author: lachlangrose

LoopStructural/interpolators/_finite_difference_interpolator.py

@@ -169,7 +169,7 @@ def add_value_constraints(self, w=1.0):
                 a,
                 points[inside, 3],
                 idc[inside, :],
-                w=w * points[inside, 4] * self.support.element_size,
+                w=w * points[inside, 4],
                 name="value",
             )
             if np.sum(inside) <= 0:
@@ -315,13 +315,9 @@ def add_gradient_constraints(self, w=1.0):
             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, B, idc[inside, :], w=w * self.vol, name="gradient"
-            )
+            self.add_constraints_to_least_squares(A, B, idc[inside, :], w=w, name="gradient")
             A = np.einsum("ij,ijk->ik", dip_vector.T, T)
-            self.add_constraints_to_least_squares(
-                A, B, idc[inside, :], w=w * self.vol, name="gradient"
-            )
+            self.add_constraints_to_least_squares(A, B, idc[inside, :], w=w, name="gradient")
             if np.sum(inside) <= 0:
                 logger.warning(
                     f" {np.sum(~inside)} \
@@ -370,21 +366,21 @@ def add_norm_constraints(self, w=1.0):
                 T[:, 0, :],
                 points[inside, 3],
                 idc[inside, :],
-                w=w * self.vol,
+                w=w,
                 name="norm",
             )
             self.add_constraints_to_least_squares(
                 T[:, 1, :],
                 points[inside, 4],
                 idc[inside, :],
-                w=w * self.vol,
+                w=w,
                 name="norm",
             )
             self.add_constraints_to_least_squares(
                 T[:, 2, :],
                 points[inside, 5],
                 idc[inside, :],
-                w=w * self.vol,
+                w=w,
                 name="norm",
             )
             if np.sum(inside) <= 0:

LoopStructural/interpolators/_p1interpolator.py

@@ -45,9 +45,9 @@ def add_norm_constraints(self, w=1.0):
         points = self.get_norm_constraints()
         if points.shape[0] > 0:
             grad, elements, inside = self.support.evaluate_shape_derivatives(points[:, :3])
-            size = self.support.element_size[elements[inside]]
+            size = self.support.element_scale[elements[inside]]
             wt = np.ones(size.shape[0])
-            wt *= w * size
+            wt *= w  # s* size
             elements = np.tile(self.support.elements[elements[inside]], (3, 1, 1))
 
             elements = elements.swapaxes(0, 1)
@@ -56,9 +56,10 @@ def add_norm_constraints(self, w=1.0):
             # elements = elements.swapaxes(1, 2)
 
             self.add_constraints_to_least_squares(
-                grad[inside, :, :] * wt[:, None, None],
-                points[inside, 3:6] * wt[:, None],
+                grad[inside, :, :],
+                points[inside, 3:6],
                 elements,
+                w=wt,
                 name="norm",
             )
             self.up_to_date = False
@@ -68,13 +69,14 @@ def add_value_constraints(self, w=1.0):
         points = self.get_value_constraints()
         if points.shape[0] > 1:
             N, elements, inside = self.support.evaluate_shape(points[:, :3])
-            size = self.support.element_size[elements[inside]]
+            size = self.support.element_scale[elements[inside]]
             wt = np.ones(size.shape[0])
-            wt *= w * size
+            wt *= w  # * size
             self.add_constraints_to_least_squares(
-                N[inside, :] * wt[:, None],
-                points[inside, 3] * wt,
+                N[inside, :],
+                points[inside, 3],
                 self.support.elements[elements[inside], :],
+                w=wt,
                 name="value",
             )
             self.up_to_date = False
@@ -88,7 +90,7 @@ def minimise_edge_jumps(self, w=0.1, vector_func=None, vector=None, name="edge j
         bc_t1 = self.support.barycentre[self.support.shared_element_relationships[:, 0]]
         bc_t2 = self.support.barycentre[self.support.shared_element_relationships[:, 1]]
         norm = self.support.shared_element_norm
-        shared_element_size = self.support.shared_element_size
+        shared_element_scale = self.support.shared_element_scale
 
         # evaluate normal if using vector func for cp2
         if vector_func:
@@ -116,9 +118,10 @@ def minimise_edge_jumps(self, w=0.1, vector_func=None, vector=None, name="edge j
         # tri_cp2 = np.hstack([self.support.elements[cp2_tri1],self.support.elements[tri2]])
         # add cp1 and cp2 to the least squares system
         self.add_constraints_to_least_squares(
-            const * shared_element_size[:, None] * w,
+            const,
             np.zeros(const.shape[0]),
             tri_cp1,
+            w=w,
             name=name,
         )
         self.up_to_date = False

LoopStructural/interpolators/supports/_3d_base_structured.py

@@ -446,6 +446,11 @@ def global_cell_indices(self, indexes) -> np.ndarray:
     def element_size(self):
         return np.prod(self.step_vector)
 
+    @property
+    def element_scale(self):
+        # all elements are the same size
+        return 1.0
+
     @property
     def vtk(self):
         try:

LoopStructural/interpolators/supports/_3d_structured_tetra.py

@@ -83,8 +83,13 @@ def element_size(self):
             self.nodes[self.elements[:, :4], :][:, 1:, :]
             - self.nodes[self.elements[:, :4], :][:, 0, None, :]
         )
+
         return np.abs(np.linalg.det(vecs)) / 6
 
+    @property
+    def element_scale(self):
+        return self.element_size / np.mean(self.element_size)
+
     @property
     def barycentre(self) -> np.ndarray:
         """
@@ -185,6 +190,10 @@ def shared_element_size(self):
         norm = self.shared_element_norm
         return 0.5 * np.linalg.norm(norm, axis=1)
 
+    @property
+    def shared_element_scale(self):
+        return self.shared_element_size / np.mean(self.shared_element_size)
+
     def evaluate_value(self, pos: np.ndarray, property_array: np.ndarray) -> np.ndarray:
         """
         Evaluate value of interpolant