scale by element size

Author: lachlangrose

LoopStructural/interpolators/_discrete_interpolator.py

@@ -219,6 +219,7 @@ def add_constraints_to_least_squares(self, A, B, idc, w=1.0, name="undefined"):
         # normalise by rows of A
         # Should this be done? It should make the solution more stable
         length = np.linalg.norm(A, axis=1)
+        # length[length>0] = 1.
         B[length > 0] /= length[length > 0]
         # going to assume if any are nan they are all nan
         mask = np.any(np.isnan(A), axis=1)

LoopStructural/interpolators/_p1interpolator.py

@@ -45,7 +45,7 @@ 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  # s* size
             elements = np.tile(self.support.elements[elements[inside]], (3, 1, 1))
@@ -70,6 +70,7 @@ def add_value_constraints(self, w=1.0):
         if points.shape[0] > 1:
             N, elements, inside = self.support.evaluate_shape(points[:, :3])
             size = self.support.element_size[elements[inside]]
+
             wt = np.ones(size.shape[0])
             wt *= w  # * size
             self.add_constraints_to_least_squares(
@@ -110,13 +111,16 @@ def minimise_edge_jumps(self, w=0.1, vector_func=None, vector=None, name="edge j
         const_n = -np.einsum("ij,ijk->ik", norm, Dn)
         # const_t_cp2 = np.einsum('ij,ikj->ik',normal,cp2_Dt)
         # const_n_cp2 = -np.einsum('ij,ikj->ik',normal,cp2_Dn)
-
+        # shared_element_size = self.support.shared_element_size
+        # const_t /= shared_element_size[:, None]  # normalise by element size
+        # const_n /= shared_element_size[:, None]  # normalise by element size
         const = np.hstack([const_t, const_n])
 
         # get vertex indexes
         tri_cp1 = np.hstack([self.support.elements[tri1], self.support.elements[tri2]])
         # 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,
             np.zeros(const.shape[0]),

LoopStructural/interpolators/supports/_3d_structured_tetra.py

@@ -85,10 +85,14 @@ def element_size(self):
         )
 
         return np.abs(np.linalg.det(vecs)) / 6
-
+    
     @property
     def element_scale(self):
-        return self.element_size / np.mean(self.element_size)
+        size = self.element_size
+        size-= np.min(size)
+        size/= np.max(size)
+        size+=1.
+        return size
 
     @property
     def barycentre(self) -> np.ndarray:
@@ -189,7 +193,7 @@ 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)