Revert "IMprove smoothing at octree change"

This reverts commit 1c0c690.

Author: domfournier

simpeg_drivers/driver.py

@@ -560,7 +560,7 @@ def get_regularization(self):
                     f"aveCC2F{fun.orientation}",
                 )
                 fun.set_weights(**{comp: average_op @ weight})
-                fun.norm = np.round(average_op @ norm, decimals=3)
+                fun.norm = average_op @ norm
                 functions.append(fun)
 
                 if is_rotated:
@@ -583,7 +583,7 @@ def get_regularization(self):
                         f"aveCC2F{fun.orientation}",
                     )
                     backward_fun.set_weights(**{comp: average_op @ weight})
-                    backward_fun.norm = np.round(average_op @ norm, decimals=3)
+                    backward_fun.norm = average_op @ norm
                     functions.append(backward_fun)
 
             # Will avoid recomputing operators if the regularization mesh is the same

simpeg_drivers/utils/regularization.py

@@ -339,20 +339,20 @@ def gradient_operator(
     :param volumes: Partial volume array.
     :param n_cells: Number of cells in mesh.
     """
-    grad = ssp.csr_matrix(
+    Grad = ssp.csr_matrix(
         (volumes, (neighbors[:, 0], neighbors[:, 1])), shape=(n_cells, n_cells)
     )
 
     # Normalize rows
-    vol = mkvc(grad.sum(axis=1))
-    vol[vol > 0] = 1.0 / vol[vol > 0]
-    grad = -sdiag(vol) * grad
+    Vol = mkvc(Grad.sum(axis=1))
+    Vol[Vol > 0] = 1.0 / Vol[Vol > 0]
+    Grad = -sdiag(Vol) * Grad
 
     diag = np.ones(n_cells)
-    diag[vol == 0] = 0
-    grad = sdiag(diag) + grad
+    diag[Vol == 0] = 0
+    Grad = sdiag(diag) + Grad
 
-    return grad
+    return Grad
 
 
 def rotated_gradient(
@@ -399,8 +399,8 @@ def rotated_gradient(
     )
 
     unit_grad = gradient_operator(neighbors, volumes, n_cells)
-
-    return unit_grad
+    axes = "xyz" if dim == 3 else "xz"
+    return sdiag(1 / mesh.h_gridded[:, axes.find(axis)]) @ unit_grad
 
 
 def ensure_dip_direction_convention(
@@ -464,40 +464,26 @@ def set_rotated_operators(
     :param forward: Whether to use forward or backward difference for
         derivative approximations.
     """
-    mesh = function.regularization_mesh
-    axes = "xyz" if mesh.dim == 3 else "xz"
-
-    h_cell = mesh.mesh.h_gridded[:, axes.find(axis)]
-
-    unit_grad_op = rotated_gradient(mesh.mesh, neighbors, axis, dip, direction, forward)
-
-    count_op = unit_grad_op.copy()
-    count_op.data = np.ones_like(count_op.data)
-
-    grad_op_active = mesh.Pac.T @ (unit_grad_op @ mesh.Pac)
-    count_op = mesh.Pac.T @ (count_op @ mesh.Pac)
-    active_faces = np.isclose(grad_op_active @ np.ones(mesh.n_cells), 0)
-    active_faces &= grad_op_active.max(axis=1).toarray().ravel() != 0
-
-    count_op = count_op[active_faces, :]
-
-    h_op = sdiag(
-        (count_op @ (mesh.Pac.T @ h_cell) / np.asarray(count_op.sum(axis=1)).ravel())
-        ** -1.0
+    grad_op = rotated_gradient(
+        function.regularization_mesh.mesh, neighbors, axis, dip, direction, forward
+    )
+    grad_op_active = function.regularization_mesh.Pac.T @ (
+        grad_op @ function.regularization_mesh.Pac
+    )
+    active_faces = np.isclose(
+        grad_op_active @ np.ones(function.regularization_mesh.n_cells), 0
     )
+    active_faces &= grad_op_active.max(axis=1).toarray().ravel() != 0
 
-    grad_op = h_op @ grad_op_active[active_faces, :]
-    avg_op = abs(grad_op_active[active_faces, :])
-    avg_op = sdiag(np.asarray(avg_op.sum(axis=1)).ravel() ** -1.0) @ avg_op
     setattr(
-        mesh,
+        function.regularization_mesh,
         f"_cell_gradient_{function.orientation}",
-        grad_op,
+        grad_op_active[active_faces, :],
     )
     setattr(
-        mesh,
+        function.regularization_mesh,
         f"_aveCC2F{function.orientation}",
-        avg_op,
+        sdiag(np.ones(function.regularization_mesh.n_cells))[active_faces, :],
     )
 
     return function