Face props mass matrices

discretize/base/base_tensor_mesh.py

@@ -797,18 +797,12 @@ def _fastInnerProduct(
 
         Parameters
         ----------
-        model : numpy.ndarray
-            material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
-
         projection_type : str
             'edges' or 'faces'
-
-        returnP : bool
-            returns the projection matrices
-
+        model : numpy.ndarray
+            material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
         invert_model : bool
             inverts the material property
-
         invert_matrix : bool
             inverts the matrix
 
@@ -875,9 +869,9 @@ def _fastInnerProductDeriv(
         Parameters
         ----------
         projection_type : str
-            'E' or 'F'
-        tensorType : TensorType
-            type of the tensor
+            'edges' or 'faces'
+        model : numpy.ndarray
+            material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
         invert_model : bool
             inverts the material property
         invert_matrix : bool

discretize/operators/inner_products.py

@@ -13,6 +13,8 @@
     inverse_3x3_block_diagonal,
     spzeros,
     sdinv,
+    mkvc,
+    is_scalar,
 )
 import numpy as np
 
@@ -95,6 +97,52 @@ def get_edge_inner_product(  # NOQA D102
             do_fast=do_fast,
         )
 
+    def get_edge_inner_product_surface(  # NOQA D102
+        self, model=None, invert_model=False, invert_matrix=False, **kwargs
+    ):
+        # Inherited documentation from discretize.base.BaseMesh
+        if model is None:
+            model = np.ones(self.nF)
+
+        if invert_model:
+            model = 1.0 / model
+
+        if is_scalar(model):
+            model = model * np.ones(self.nF)
+        # COULD ADD THIS CASE IF DESIRED
+        # elif len(model) == self.dim:
+        #     model = np.r_[
+        #         [model[ii]*self.vnF[ii] for ii in range(0, self.dim)]
+        #     ]
+
+        # Isotropic case only
+        if model.size != self.nF:
+            raise ValueError(
+                "Unexpected shape of tensor: {}".format(model.shape),
+                "Must be scalar or have length equal to total number of faces.",
+            )
+
+        # number of elements we are averaging (equals dim for regular
+        # meshes, but for cyl, where we use symmetry, it is 1 for edge
+        # variables and 2 for face variables)
+        if self._meshType == "CYL":
+            shape = self.vnE
+            if self.is_symmetric:
+                n_elements = 1
+            else:
+                n_elements = sum([1 if x != 0 else 0 for x in shape]) - 1
+        else:
+            n_elements = self.dim - 1
+
+        Aprop = self.face_areas * mkvc(model)
+        Av = self.average_edge_to_face
+        M = n_elements * sdiag(Av.T * Aprop)
+
+        if invert_matrix:
+            return sdinv(M)
+        else:
+            return M
+
     def _getInnerProduct(
         self,
         projection_type,
@@ -276,6 +324,87 @@ def get_edge_inner_product_deriv(  # NOQA D102
             invert_matrix=invert_matrix,
         )
 
+    def get_edge_inner_product_surface_deriv(  # NOQA D102
+        self, model, invert_model=False, invert_matrix=False, **kwargs
+    ):
+        # Inherited documentation from discretize.base.BaseMesh
+        if model is None:
+            tensorType = -1
+        elif is_scalar(model):
+            tensorType = 0
+        elif model.size == self.nF:
+            tensorType = 1
+        else:
+            raise ValueError(
+                "Unexpected shape of tensor: {}".format(model.shape),
+                "Must be scalar or have length equal to total number of faces.",
+            )
+
+        dMdprop = None
+
+        if invert_matrix or invert_model:
+            MI = self.get_edge_inner_product_surface(
+                model,
+                invert_model=invert_model,
+                invert_matrix=invert_matrix,
+            )
+
+        # number of elements we are averaging (equals dim for regular
+        # meshes, but for cyl, where we use symmetry, it is 1 for edge
+        # variables and 2 for face variables)
+        if self._meshType == "CYL":
+            shape = self.vnE
+            if self.is_symmetric:
+                n_elements = 1
+            else:
+                n_elements = sum([1 if x != 0 else 0 for x in shape]) - 1
+        else:
+            n_elements = self.dim - 1
+
+        A = sdiag(self.face_areas)
+        Av = self.average_edge_to_face
+
+        if tensorType == 0:  # isotropic, constant
+            ones = sp.csr_matrix(
+                (np.ones(self.nF), (range(self.nF), np.zeros(self.nF))),
+                shape=(self.nF, 1),
+            )
+            if not invert_matrix and not invert_model:
+                dMdprop = n_elements * Av.T * A * ones
+            elif invert_matrix and invert_model:
+                dMdprop = n_elements * (
+                    sdiag(MI.diagonal() ** 2)
+                    * Av.T
+                    * A
+                    * ones
+                    * sdiag(1.0 / model**2)
+                )
+            elif invert_model:
+                dMdprop = n_elements * Av.T * A * sdiag(-1.0 / model**2)
+            elif invert_matrix:
+                dMdprop = n_elements * (sdiag(-MI.diagonal() ** 2) * Av.T * A)
+
+        elif tensorType == 1:  # isotropic, variable in space
+            if not invert_matrix and not invert_model:
+                dMdprop = n_elements * Av.T * A
+            elif invert_matrix and invert_model:
+                dMdprop = n_elements * (
+                    sdiag(MI.diagonal() ** 2) * Av.T * A * sdiag(1.0 / model**2)
+                )
+            elif invert_model:
+                dMdprop = n_elements * Av.T * A * sdiag(-1.0 / model**2)
+            elif invert_matrix:
+                dMdprop = n_elements * (sdiag(-MI.diagonal() ** 2) * Av.T * A)
+
+        if dMdprop is not None:
+
+            def innerProductDeriv(v):
+                return sdiag(v) * dMdprop
+
+            return innerProductDeriv
+        else:
+            return None
+
     def _getInnerProductDeriv(
         self,
         model,

tests/base/test_basemesh.py

@@ -50,8 +50,14 @@ class TestBaseMesh(unittest.TestCase):
     not_implemented_functions = [
         "get_face_inner_product",
         "get_edge_inner_product",
+        "get_face_inner_product_surface",
+        "get_edge_inner_product_surface",
+        "get_edge_inner_product_line",
         "get_face_inner_product_deriv",
         "get_edge_inner_product_deriv",
+        "get_face_inner_product_surface_deriv",
+        "get_edge_inner_product_surface_deriv",
+        "get_edge_inner_product_line_deriv",
         "point2index",
         "get_interpolation_matrix",
     ]