fix: actually fixing divide by zero

LoopStructural/interpolators/base_structured_3d_support.py

@@ -27,6 +27,8 @@ def __init__(
         # we use property decorators to update these when different parts of
         # the geometry need to change
         # inisialise the private attributes
+        if np.any(step_vector==0):
+            logger.warning(f"Step vector {step_vector} has zero values")
         self._nsteps = np.array(nsteps, dtype=int)
         self._step_vector = np.array(step_vector)
         self._origin = np.array(origin)
@@ -72,7 +74,7 @@ def step_vector(self, step_vector):
         newsteps = self._nsteps / change_factor
         self._nsteps = np.ceil(newsteps).astype(int)
         self._step_vector = step_vector
-
+        
     @property
     def origin(self):
         return self._origin

LoopStructural/interpolators/piecewiselinear_interpolator.py

@@ -435,11 +435,12 @@ def add_gradient_orthogonal_constraints(self, points, vector, w=1.0, B=0):
             # e, inside = self.support.elements_for_array(points[:, :3])
             # nodes = self.support.nodes[self.support.elements[e]]
             norm = np.linalg.norm(vector, axis=1)
-            mask = ~np.isnan(norm)
-            vector /= norm[mask, None]
+            mask = np.isnan(norm)
+            mask = ~np.logical_or(mask, norm==0)
+            vector[mask,:] /= norm[mask, None]
             vecs = vertices[:, 1:, :] - vertices[:, 0, None, :]
-            vol = np.abs(np.linalg.det(vecs)) / 6
-            element_gradients /= norm[mask, None, None]
+            # vol = np.abs(np.linalg.det(vecs)) / 6
+            element_gradients[mask,:,:] /= norm[mask, None, None]
 
             A = np.einsum("ij,ijk->ik", vector, element_gradients)