Hierarchical inner products for unstructured meshes

discretize/curvilinear_mesh.py

@@ -797,10 +797,13 @@ def _get_edge_surf_int_proj_mats(self, only_boundary=False, with_area=True):
             edge_dirs = self.edge_tangents[node_edges]
             t_for = np.concatenate((edge_dirs, face_normals[:, None, :]), axis=1)
             t_inv = np.linalg.inv(t_for)
-            t_inv = t_inv[:, :, :-1] / 4  # n_edges_per_thing
+            t_inv = t_inv[:, :, :-1]
 
             if with_area:
                 t_inv *= face_areas[:, None, None]
+                t_inv /= 4  # n_edges_per_thing
+            else:
+                t_inv /= 2  # sqrt n_edges_per_thing
 
             T = C2F @ sp.csr_matrix(
                 (t_inv.reshape(-1), T_col_inds, T_ind_ptr),

discretize/unstructured_mesh.py

@@ -519,7 +519,64 @@ def get_edge_inner_product(  # NOQA D102
                 "The inverse of the inner product matrix with a tetrahedral mesh is not supported."
             )
         return self.__get_inner_product("E", model, invert_model)
+
+    def get_edge_inner_product_surface(  # NOQA D102
+        self,
+        model=None,
+        invert_model=False,
+        invert_matrix=False
+    ):
+        # Documentation inherited from discretize.base.BaseMesh
+        if invert_matrix:
+            raise NotImplementedError(
+                "The inverse of the inner product matrix with a tetrahedral mesh is not supported."
+            )
+
+        # Edge inner product surface projection matrices
+        n_faces = self.n_faces
+        face_areas = self.face_areas
+        dim = self.dim
+
+        if dim == 2:
+
+            if model is None:
+                Mu = sp.diags(face_areas)
+            else:
+                if invert_model:
+                    model = 1.0 / model
+
+                if (model.size == 1) | (model.size == n_faces):
+                    Mu = sp.diags(model * face_areas)
+                else:
+                    raise ValueError(
+                        "Unrecognized size of model vector.",
+                        "Must be scalar or have length equal to total number of faces."
+                    )
+
+            return Mu        
+
+        else:
+            Ps = self._get_edge_surf_int_proj_mats(with_area=False)
+
+            if model is None:
+                Mu = sp.diags(np.tile(face_areas, dim))
+            else:
+                if invert_model:
+                    model = 1.0 / model
+
+                if (model.size == 1) | (model.size == n_faces):
+                    Mu = sp.diags(np.tile(model * face_areas, dim))
+                else:
+                    raise ValueError(
+                        "Unrecognized size of model vector.",
+                        "Must be scalar or have length equal to total number of faces."
+                )
 
+            A = np.sum([P.T @ Mu @ P for P in Ps])
+
+            return A
+
+
     def __get_inner_product_deriv_func(self, i_type, model):
         Ps, _ = self.__get_inner_product_projection_matrices(i_type)
         dim = self.dim
@@ -668,10 +725,13 @@ def _get_edge_surf_int_proj_mats(self, only_boundary=False, with_area=True):
             edge_dirs = self.edge_tangents[node_edges]
             t_for = np.concatenate((edge_dirs, face_normals[:, None, :]), axis=1)
             t_inv = invert_blocks(t_for)
-            t_inv = t_inv[:, :, :-1] / 3  # n_edges_per_thing
+            t_inv = t_inv[:, :, :-1]
 
             if with_area:
                 t_inv *= face_areas[:, None, None]
+                t_inv /= 3  # n_edges_per_thing
+            else:
+                t_inv /= np.sqrt(3)  # sqrt n_edges_per_thing
 
             T = C2F @ sp.csr_matrix(
                 (t_inv.reshape(-1), T_col_inds, T_ind_ptr),