updates to face innerProducts

Increase Speed
Add derivatives
Add tests

SimPEG/Mesh/InnerProducts.py

@@ -10,7 +10,8 @@ class InnerProducts(object):
     def __init__(self):
         raise Exception('InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class.')
 
-    def getFaceInnerProduct(self, materialProperty=None, returnP=False, invertProperty=False):
+    def getFaceInnerProduct(self, materialProperty=None, returnP=False,
+                            invertProperty=False):
         """
             :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param bool returnP: returns the projection matrices
@@ -18,6 +19,14 @@ class InnerProducts(object):
             :rtype: scipy.csr_matrix
             :return: M, the inner product matrix (nF, nF)
         """
+        fast = None
+
+        if returnP is False and hasattr(self, '_fastFaceInnerProduct'):
+            fast = self._fastFaceInnerProduct(materialProperty=materialProperty, invertProperty=invertProperty)
+
+        if fast is not None:
+            return fast
+
         if invertProperty:
             materialProperty = invPropertyTensor(self, materialProperty)
 
@@ -62,6 +71,34 @@ class InnerProducts(object):
         else:
             return A
 
+    def getFaceInnerProductDeriv(self, materialProperty=None, P=None):
+        """
+            :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
+            :rtype: scipy.csr_matrix
+            :return: M, the inner product matrix (nF, nF)
+        """
+
+        fast = None
+
+        if hasattr(self, '_fastFaceInnerProductDeriv'):
+            fast = self._fastFaceInnerProductDeriv(materialProperty=materialProperty)
+
+        if fast is not None:
+            return fast
+
+        raise NotImplementedError('Derivatives for the material property specified are not yet implemented.')
+
+        if P is None:
+            M, P = getFaceInnerProduct(self, materialProperty=materialProperty, returnP=True)
+
+        d = self.dim
+
+        if d == 1:
+            P[0].T * sp.identity(n) * P[0]
+
+
+
+
     def getEdgeInnerProduct(self, materialProperty=None, returnP=False, invertProperty=False):
         """
             :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))

SimPEG/Mesh/TensorMesh.py

@@ -39,6 +39,7 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
 
     def __init__(self, h_in, x0=None):
         assert type(h_in) is list, 'h_in must be a list'
+        assert len(h_in) in [1,2,3], 'h_in must be of dimension 1, 2, or 3'
         h = range(len(h_in))
         for i, h_i in enumerate(h_in):
             if type(h_i) in [int, long, float]:
@@ -491,6 +492,49 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
             indzu = (self.gridCC[:,2]==max(self.gridCC[:,2]))
             return indxd, indxu, indyd, indyu, indzd, indzu
 
+    def _fastFaceInnerProduct(self, materialProperty=None, invertProperty=False):
+        """
+            Fast version of getFaceInnerProduct.
+            This does not handle the case of a full tensor materialProperty.
+
+            :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
+            :param bool returnP: returns the projection matrices
+            :param bool invertProperty: inverts the material property
+            :rtype: scipy.csr_matrix
+            :return: M, the inner product matrix (nF, nF)
+        """
+        if materialProperty is None:
+            materialProperty = np.ones(self.nC)
+
+        if materialProperty.size == self.nC:
+            if invertProperty:
+                materialProperty = 1./materialProperty
+            Av = self.aveF2CC
+            V = Utils.sdiag(self.vol)
+            return self.dim * Utils.sdiag(Av.T * V * materialProperty)
+        if materialProperty.size == self.nC*self.dim:
+            if invertProperty:
+                materialProperty = 1./materialProperty
+            Av = self.aveF2CCV
+            V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
+            return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty))
+
+
+
+    def _fastFaceInnerProductDeriv(self, materialProperty=None):
+        """
+            :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
+            :rtype: scipy.csr_matrix
+            :return: M, the inner product matrix (nF, nF)
+        """
+        if materialProperty is None or materialProperty.size == self.nC:
+            Av = self.aveF2CC
+            return self.dim * Av.T * Utils.sdiag(self.vol)
+        if materialProperty.size == self.nC*self.dim: # anisotropic
+            Av = self.aveF2CCV
+            V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
+            return Av.T * V
+
 
 if __name__ == '__main__':
     print('Welcome to tensor mesh!')

SimPEG/Tests/test_massMatrices.py

@@ -30,7 +30,7 @@ class TestInnerProducts(OrderTest):
             sigma = np.c_[call(sigma1, Gc)]
             analytic = 647./360  # Found using sympy.
         elif self.sigmaTest == 3:
-            sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
+            sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
             analytic = 37./12  # Found using sympy.
         elif self.sigmaTest == 6:
             sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc),

SimPEG/Tests/test_massMatricesDerivs.py

@@ -0,0 +1,40 @@
+import numpy as np
+import unittest
+from SimPEG import *
+from TestUtils import checkDerivative
+
+
+class TestInnerProductsDerivs(unittest.TestCase):
+    def setUp(self):
+        pass
+
+    def test_FaceIP_derivs_isotropic(self):
+        for d in range(3):
+            mesh = Mesh.TensorMesh([10,5,4][d:])
+            M,Ps = mesh.getFaceInnerProduct(returnP=True)
+            v = np.random.rand(mesh.nF)
+            def fun(sig):
+                M = mesh.getFaceInnerProduct(sig)
+                Md = mesh.getFaceInnerProductDeriv(sig)
+                return M*v, Utils.sdiag(v)*Md
+            sig = np.random.rand(mesh.nC)
+            passed = checkDerivative(fun, sig, plotIt=False)
+            self.assertTrue(passed)
+
+    def test_FaceIP_derivs_anisotropic(self):
+        for d in range(3):
+            mesh = Mesh.TensorMesh([10,5,4][d:])
+            M,Ps = mesh.getFaceInnerProduct(returnP=True)
+            v = np.random.rand(mesh.nF)
+            def fun(sig):
+                M = mesh.getFaceInnerProduct(sig)
+                Md = mesh.getFaceInnerProductDeriv(sig)
+                return M*v, Utils.sdiag(v)*Md
+            sig = np.random.rand(mesh.nC*mesh.dim)
+            passed = checkDerivative(fun, sig, plotIt=False)
+            self.assertTrue(passed)
+
+
+
+if __name__ == '__main__':
+    unittest.main()

docs/api_InnerProducts.rst

@@ -131,6 +131,7 @@ If ``returnP=True`` is requested in any of these methods the projection matrices
 
 The derivation for ``edgeInnerProducts`` is exactly the same, however, when we approximate the integral using the fields around each node, the projection matrices look a bit different because we have 12 edges in 3D instead of just 6 faces. The interface to the code is exactly the same.
 
+
 Defining Tensor Properties
 --------------------------
 
@@ -158,6 +159,42 @@ Depending on the number of columns (either 1, 3, or 6) of mu, the material prope
     \vec{\mu} = \left[\begin{matrix} \mu_{1} & \mu_{3} \\ \mu_{3} & \mu_{2} \end{matrix}\right]
 
 
+Structure of Matrices
+---------------------
+
+Both the isotropic, and anisotropic material properties result in a diagonal mass matrix.
+Which is nice and easy to invert if necessary, however, in the fully anisotropic case which is not aligned with the grid, the matrix is not diagonal. This can be seen for a 3D mesh in the figure below.
+
+.. plot::
+
+    from SimPEG import *
+    mesh = Mesh.TensorMesh([10,50,3])
+    m1 = np.random.rand(mesh.nC)
+    m2 = np.random.rand(mesh.nC,3)
+    m3 = np.random.rand(mesh.nC,6)
+    M = range(3)
+    M[0] = mesh.getFaceInnerProduct(m1)
+    M[1] = mesh.getFaceInnerProduct(m2)
+    M[2] = mesh.getFaceInnerProduct(m3)
+    plt.figure(figsize=(13,5))
+    for i, lab in enumerate(['Isotropic','Anisotropic','Tensor']):
+        plt.subplot(131 + i)
+        plt.spy(M[i],ms=0.5,color='k')
+        plt.tick_params(axis='both',which='both',labeltop='off',labelleft='off')
+        plt.title(lab + ' Material Property')
+    plt.show()
+
+
+Taking Derivatives
+------------------
+
+TODO: Take the derivatives of the tensors.
+
+.. math::
+
+    \left[\begin{smallmatrix}0.5 \sigma_{1} & 0 & 0.25 \sigma_{3} & 0.25 \sigma_{3}\\0 & 0.5 \sigma_{1} & 0.25 \sigma_{3} & 0.25 \sigma_{3}\\0.25 \sigma_{3} & 0.25 \sigma_{3} & 0.5 \sigma_{2} & 0\\0.25 \sigma_{3} & 0.25 \sigma_{3} & 0 & 0.5 \sigma_{2}\end{smallmatrix}\right]
+
+
 The API
 -------