From 9512b377c0cf8199052621a8a18b32aa639fb023 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Mon, 16 Jun 2014 12:17:01 -0600
Subject: [PATCH 01/14] InnerProducts working as an operator. Simplifications
 and generalizations in inner product code.

---
 SimPEG/Mesh/InnerProducts.py            | 257 +++++++++++-------------
 SimPEG/Mesh/TensorMesh.py               |  97 +++------
 SimPEG/Tests/test_massMatricesDerivs.py | 117 ++++-------
 SimPEG/Tests/test_utils.py              |  16 +-
 SimPEG/Utils/matutils.py                |  57 ++++--
 5 files changed, 230 insertions(+), 314 deletions(-)

diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py
index 8268b29a..026e406d 100644
--- a/SimPEG/Mesh/InnerProducts.py
+++ b/SimPEG/Mesh/InnerProducts.py
@@ -10,112 +10,43 @@ 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, prop=None, returnP=False,
-                            invProp=False, invMat=False, doFast=True):
+    def getFaceInnerProduct(self, prop=None, invProp=False, invMat=False, doFast=True):
         """
             :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
-            :param bool returnP: returns the projection matrices
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
             :param bool doFast: do a faster implementation if available.
             :rtype: scipy.csr_matrix
             :return: M, the inner product matrix (nF, nF)
         """
-        fast = None
+        return self._getInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat, doFast=True)
 
-        if returnP is False and hasattr(self, '_fastFaceInnerProduct') and doFast:
-            fast = self._fastFaceInnerProduct(prop=prop, invProp=invProp, invMat=invMat)
-
-        if fast is not None:
-            return fast
-
-        if invProp:
-            prop = invPropertyTensor(self, prop)
-
-        Mu = makePropertyTensor(self, prop)
-
-        d = self.dim
-        # We will multiply by sqrt on each side to keep symmetry
-        V = sp.kron(sp.identity(d), sdiag(np.sqrt((2**(-d))*self.vol)))
-
-        if d == 1:
-            fP = _getFacePx(self)
-            P000 = V*fP('fXm')
-            P100 = V*fP('fXp')
-        elif d == 2:
-            fP = _getFacePxx(self)
-            P000 = V*fP('fXm', 'fYm')
-            P100 = V*fP('fXp', 'fYm')
-            P010 = V*fP('fXm', 'fYp')
-            P110 = V*fP('fXp', 'fYp')
-        elif d == 3:
-            fP = _getFacePxxx(self)
-            P000 = V*fP('fXm', 'fYm', 'fZm')
-            P100 = V*fP('fXp', 'fYm', 'fZm')
-            P010 = V*fP('fXm', 'fYp', 'fZm')
-            P110 = V*fP('fXp', 'fYp', 'fZm')
-            P001 = V*fP('fXm', 'fYm', 'fZp')
-            P101 = V*fP('fXp', 'fYm', 'fZp')
-            P011 = V*fP('fXm', 'fYp', 'fZp')
-            P111 = V*fP('fXp', 'fYp', 'fZp')
-
-        A = P000.T*Mu*P000 + P100.T*Mu*P100
-        P = [P000, P100]
-
-        if d > 1:
-            A = A + P010.T*Mu*P010 + P110.T*Mu*P110
-            P += [P010, P110]
-        if d > 2:
-            A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111
-            P += [P001, P101, P011,  P111]
-
-        if invMat and tensorType(self, prop) < 3:
-            A = sdInv(A)
-        elif invMat and tensorType(self, prop) == 3:
-            raise Exception('Solver needed to invert A.')
-
-        if returnP:
-            return A, P
-        else:
-            return A
-
-    def getFaceInnerProductDeriv(self, prop=None, v=None, P=None, doFast=True):
+    def getEdgeInnerProduct(self, prop=None, invProp=False, invMat=False, doFast=True):
         """
             :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
-            :param numpy.array v: vector to multiply (required in the general implementation)
-            :param list P: list of projection matrices
-            :param bool doFast: do a faster implementation if available.
-            :rtype: scipy.csr_matrix
-            :return: dMdm, the derivative of the inner product matrix (nF, nC*nA)
-        """
-        fast = None
-
-        if hasattr(self, '_fastFaceInnerProductDeriv') and doFast:
-            fast = self._fastFaceInnerProductDeriv(prop=prop, v=v)
-
-        if fast is not None:
-            return fast
-
-        if P is None:
-            M, P = self.getFaceInnerProduct(prop=prop, returnP=True)
-
-        return self._getInnerProductDeriv(prop, v, P, self.nF)
-
-    def getEdgeInnerProduct(self, prop=None, returnP=False,
-                            invProp=False, invMat=False, doFast=True):
-        """
-            :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
-            :param bool returnP: returns the projection matrices
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
             :param bool doFast: do a faster implementation if available.
             :rtype: scipy.csr_matrix
             :return: M, the inner product matrix (nE, nE)
         """
+        return self._getInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat)
+
+    def _getInnerProduct(self, projType, prop=None, invProp=False, invMat=False, doFast=True):
+        """
+            :param str projType: 'F' for faces 'E' for edges
+            :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
+            :param bool invProp: inverts the material property
+            :param bool invMat: inverts the matrix
+            :param bool doFast: do a faster implementation if available.
+            :rtype: scipy.csr_matrix
+            :return: M, the inner product matrix (nE, nE)
+        """
+        assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
         fast = None
 
-        if returnP is False and hasattr(self, '_fastEdgeInnerProduct') and doFast:
-            fast = self._fastEdgeInnerProduct(prop=prop, invProp=invProp, invMat=invMat)
+        if hasattr(self, '_fastInnerProduct') and doFast:
+            fast = self._fastInnerProduct(projType, prop=prop, invProp=invProp, invMat=invMat)
 
         if fast is not None:
             return fast
@@ -123,72 +54,122 @@ class InnerProducts(object):
         if invProp:
             prop = invPropertyTensor(self, prop)
 
+        tensorType = TensorType(self, prop)
         Mu = makePropertyTensor(self, prop)
 
+        Ps = self._getInnerProductProjectionMatrices(projType, tensorType)
+
+        A = np.sum([P.T * Mu * P for P in Ps])
+
+        if invMat and tensorType < 3:
+            A = sdInv(A)
+        elif invMat and tensorType == 3:
+            raise Exception('Solver needed to invert A.')
+
+        return A
+
+    def _getInnerProductProjectionMatrices(self, projType, tensorType):
+        """
+            :param str projType: 'F' for faces 'E' for edges
+            :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
+        """
+        assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.'
+        assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
+
         d = self.dim
         # We will multiply by sqrt on each side to keep symmetry
         V = sp.kron(sp.identity(d), sdiag(np.sqrt((2**(-d))*self.vol)))
 
-        if d == 1:
-            raise NotImplementedError('getEdgeInnerProduct not implemented for 1D')
-        elif d == 2:
-            eP = _getEdgePxx(self)
-            P000 = V*eP('eX0', 'eY0')
-            P100 = V*eP('eX0', 'eY1')
-            P010 = V*eP('eX1', 'eY0')
-            P110 = V*eP('eX1', 'eY1')
-        elif d == 3:
-            eP = _getEdgePxxx(self)
-            P000 = V*eP('eX0', 'eY0', 'eZ0')
-            P100 = V*eP('eX0', 'eY1', 'eZ1')
-            P010 = V*eP('eX1', 'eY0', 'eZ2')
-            P110 = V*eP('eX1', 'eY1', 'eZ3')
-            P001 = V*eP('eX2', 'eY2', 'eZ0')
-            P101 = V*eP('eX2', 'eY3', 'eZ1')
-            P011 = V*eP('eX3', 'eY2', 'eZ2')
-            P111 = V*eP('eX3', 'eY3', 'eZ3')
+        nodes = ['000', '100', '010', '110', '001', '101', '011',  '111'][:2**d]
 
-        Mu = makePropertyTensor(self, prop)
-        A = P000.T*Mu*P000 + P100.T*Mu*P100 + P010.T*Mu*P010 + P110.T*Mu*P110
-        P = [P000, P100, P010, P110]
-        if d == 3:
-            A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111
-            P += [P001, P101, P011,  P111]
+        if projType == 'F':
+            locs = {
+                    '000': [('fXm',), ('fXm', 'fYm'), ('fXm', 'fYm', 'fZm')],
+                    '100': [('fXp',), ('fXp', 'fYm'), ('fXp', 'fYm', 'fZm')],
+                    '010': [  None  , ('fXm', 'fYp'), ('fXm', 'fYp', 'fZm')],
+                    '110': [  None  , ('fXp', 'fYp'), ('fXp', 'fYp', 'fZm')],
+                    '001': [  None  ,      None     , ('fXm', 'fYm', 'fZp')],
+                    '101': [  None  ,      None     , ('fXp', 'fYm', 'fZp')],
+                    '011': [  None  ,      None     , ('fXm', 'fYp', 'fZp')],
+                    '111': [  None  ,      None     , ('fXp', 'fYp', 'fZp')]
+                   }
+            if d == 1:
+                proj = _getFacePx(self)
+            elif d == 2:
+                proj = _getFacePxx(self)
+            elif d == 3:
+                proj = _getFacePxxx(self)
 
-        if invMat and tensorType(self, prop) < 3:
-            A = sdInv(A)
-        elif invMat and tensorType(self, prop) == 3:
-            raise Exception('Solver needed to invert A.')
+        elif projType == 'E':
+            locs = {
+                    '000': [  None  , ('eX0', 'eY0'), ('eX0', 'eY0', 'eZ0')],
+                    '100': [  None  , ('eX0', 'eY1'), ('eX0', 'eY1', 'eZ1')],
+                    '010': [  None  , ('eX1', 'eY0'), ('eX1', 'eY0', 'eZ2')],
+                    '110': [  None  , ('eX1', 'eY1'), ('eX1', 'eY1', 'eZ3')],
+                    '001': [  None  ,      None     , ('eX2', 'eY2', 'eZ0')],
+                    '101': [  None  ,      None     , ('eX2', 'eY3', 'eZ1')],
+                    '011': [  None  ,      None     , ('eX3', 'eY2', 'eZ2')],
+                    '111': [  None  ,      None     , ('eX3', 'eY3', 'eZ3')]
+                   }
+            if d == 1:
+                raise NotImplementedError('getEdgeInnerProduct not implemented for 1D')
+            elif d == 2:
+                proj = _getEdgePxx(self)
+            elif d == 3:
+                proj = _getEdgePxxx(self)
 
-        if returnP:
-            return A, P
-        else:
-            return A
+        return [V*proj(*locs[node][d-1]) for node in nodes]
 
-    def getEdgeInnerProductDeriv(self, prop=None, v=None, P=None, doFast=True):
+
+    def getFaceInnerProductDeriv(self, tensorType, P=None, doFast=True):
         """
-            :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
-            :param numpy.array v: vector to multiply (required in the general implementation)
+            :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
             :param list P: list of projection matrices
             :param bool doFast: do a faster implementation if available.
             :rtype: scipy.csr_matrix
-            :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
+            :return: dMdm, the derivative of the inner product matrix (nF, nC*nA)
         """
-
+        assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.'
         fast = None
 
-        if hasattr(self, '_fastEdgeInnerProductDeriv') and doFast:
-            fast = self._fastEdgeInnerProductDeriv(prop=prop, v=v)
+        if hasattr(self, '_fastInnerProductDeriv') and doFast:
+            fast = self._fastInnerProductDeriv('F', tensorType)
 
         if fast is not None:
             return fast
 
         if P is None:
-            M, P = self.getEdgeInnerProduct(prop=prop, returnP=True)
+            P = self._getInnerProductProjectionMatrices('F', tensorType=tensorType)
 
-        return self._getInnerProductDeriv(prop, v, P, self.nE)
+        def innerProductDeriv(v):
+            return self._getInnerProductDeriv(tensorType, P, self.nF, v)
+        return DerivOperator(innerProductDeriv)
 
-    def _getInnerProductDeriv(self, prop, v, P, n):
+
+    def getEdgeInnerProductDeriv(self, tensorType, P=None, doFast=True):
+        """
+            :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
+            :param list P: list of projection matrices
+            :param bool doFast: do a faster implementation if available.
+            :rtype: scipy.csr_matrix
+            :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
+        """
+        assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.'
+        fast = None
+
+        if hasattr(self, '_fastInnerProductDeriv') and doFast:
+            fast = self._fastInnerProductDeriv('E', tensorType)
+
+        if fast is not None:
+            return fast
+
+        if P is None:
+            P = self._getInnerProductProjectionMatrices('E', tensorType=tensorType)
+        def innerProductDeriv(v):
+            return self._getInnerProductDeriv(tensorType, P, self.nE, v)
+        return DerivOperator(innerProductDeriv)
+
+    def _getInnerProductDeriv(self, tensorType, P, n, v):
         """
             :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param numpy.array v: vector to multiply (required in the general implementation)
@@ -197,7 +178,7 @@ class InnerProducts(object):
             :rtype: scipy.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (n, nC*nA)
         """
-        if prop is None:
+        if tensorType == -1:
             return None
 
         if v is None:
@@ -206,24 +187,24 @@ class InnerProducts(object):
         d = self.dim
         Z = spzeros(self.nC, self.nC)
 
-        if isScalar(prop):
+        if tensorType == 0:
             dMdm = spzeros(n, 1)
             for i, p in enumerate(P):
                 dMdm = dMdm + sp.csr_matrix((p.T * (p * v), (range(n), np.zeros(n))), shape=(n,1))
         if d == 1:
-            if prop.size == self.nC:
+            if tensorType == 1:
                 dMdm = spzeros(n, self.nC)
                 for i, p in enumerate(P):
                     dMdm = dMdm + p.T * sdiag( p * v )
         elif d == 2:
-            if prop.size == self.nC:
+            if tensorType == 1:
                 dMdm = spzeros(n, self.nC)
                 for i, p in enumerate(P):
                     Y = p * v
                     y1 = Y[:self.nC]
                     y2 = Y[self.nC:]
                     dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 )))
-            elif prop.size == self.nC*2:
+            elif tensorType == 2:
                 dMdms = [spzeros(n, self.nC) for _ in range(2)]
                 for i, p in enumerate(P):
                     Y = p * v
@@ -232,7 +213,7 @@ class InnerProducts(object):
                     dMdms[0] = dMdms[0] + p.T * sp.vstack(( sdiag( y1 ), Z))
                     dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 )))
                 dMdm = sp.hstack(dMdms)
-            elif prop.size == self.nC*3:
+            elif tensorType == 3:
                 dMdms = [spzeros(n, self.nC) for _ in range(3)]
                 for i, p in enumerate(P):
                     Y = p * v
@@ -243,7 +224,7 @@ class InnerProducts(object):
                     dMdms[2] = dMdms[2] + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 )))
                 dMdm = sp.hstack(dMdms)
         elif d == 3:
-            if prop.size == self.nC:
+            if tensorType == 1:
                 dMdm = spzeros(n, self.nC)
                 for i, p in enumerate(P):
                     Y = p * v
@@ -251,7 +232,7 @@ class InnerProducts(object):
                     y2 = Y[self.nC:self.nC*2]
                     y3 = Y[self.nC*2:]
                     dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 ), sdiag( y3 )))
-            elif prop.size == self.nC*3:
+            elif tensorType == 2:
                 dMdms = [spzeros(n, self.nC) for _ in range(3)]
                 for i, p in enumerate(P):
                     Y = p * v
@@ -262,7 +243,7 @@ class InnerProducts(object):
                     dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 ), Z))
                     dMdms[2] = dMdms[2] + p.T * sp.vstack(( Z, Z, sdiag( y3 )))
                 dMdm = sp.hstack(dMdms)
-            elif prop.size == self.nC*6:
+            elif tensorType == 3:
                 dMdms = [spzeros(n, self.nC) for _ in range(6)]
                 for i, p in enumerate(P):
                     Y = p * v
diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py
index ba5bd258..faf45cf2 100644
--- a/SimPEG/Mesh/TensorMesh.py
+++ b/SimPEG/Mesh/TensorMesh.py
@@ -260,49 +260,21 @@ class BaseTensorMesh(BaseRectangularMesh):
         return Q.tocsr()
 
 
-    def _fastFaceInnerProduct(self, prop=None, invProp=False, invMat=False):
+    def _fastInnerProduct(self, projType, prop=None, invProp=False, invMat=False):
         """
             Fast version of getFaceInnerProduct.
             This does not handle the case of a full tensor prop.
 
             :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
+            :param str projType: 'E' or 'F'
             :param bool returnP: returns the projection matrices
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
             :rtype: scipy.csr_matrix
             :return: M, the inner product matrix (nF, nF)
         """
-        return self._fastInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat)
+        assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
 
-
-    def _fastEdgeInnerProduct(self, prop=None, invProp=False, invMat=False):
-        """
-            Fast version of getEdgeInnerProduct.
-            This does not handle the case of a full tensor prop.
-
-            :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
-            :param bool returnP: returns the projection matrices
-            :param bool invProp: inverts the material property
-            :param bool invMat: inverts the matrix
-            :rtype: scipy.csr_matrix
-            :return: M, the inner product matrix (nE, nE)
-        """
-        return self._fastInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat)
-
-
-    def _fastInnerProduct(self, AvType, prop=None, invProp=False, invMat=False):
-        """
-            Fast version of getFaceInnerProduct.
-            This does not handle the case of a full tensor prop.
-
-            :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
-            :param str AvType: 'E' or 'F'
-            :param bool returnP: returns the projection matrices
-            :param bool invProp: inverts the material property
-            :param bool invMat: inverts the matrix
-            :rtype: scipy.csr_matrix
-            :return: M, the inner product matrix (nF, nF)
-        """
         if prop is None:
             prop = np.ones(self.nC)
 
@@ -313,11 +285,11 @@ class BaseTensorMesh(BaseRectangularMesh):
             prop = prop*np.ones(self.nC)
 
         if prop.size == self.nC:
-            Av = getattr(self, 'ave'+AvType+'2CC')
+            Av = getattr(self, 'ave'+projType+'2CC')
             Vprop = self.vol * Utils.mkvc(prop)
             M = self.dim * Utils.sdiag(Av.T * Vprop)
         elif prop.size == self.nC*self.dim:
-            Av = getattr(self, 'ave'+AvType+'2CCV')
+            Av = getattr(self, 'ave'+projType+'2CCV')
             V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
             M = Utils.sdiag(Av.T * V * Utils.mkvc(prop))
         else:
@@ -328,55 +300,40 @@ class BaseTensorMesh(BaseRectangularMesh):
         else:
             return M
 
-    def _fastFaceInnerProductDeriv(self, prop=None, v=None):
+    def _fastInnerProductDeriv(self, projType, tensorType):
         """
-            :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
+            :param str projType: 'E' or 'F'
+            :param TensorType tensorType: type of the tensor
             :rtype: scipy.csr_matrix
             :return: M, the inner product matrix (nF, nF)
         """
-        return self._fastInnerProductDeriv('F', prop=prop, v=v)
-
-
-    def _fastEdgeInnerProductDeriv(self, prop=None, v=None):
-        """
-            :param numpy.array prop: 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 (nE, nE)
-        """
-        return self._fastInnerProductDeriv('E', prop=prop, v=v)
-
-
-    def _fastInnerProductDeriv(self, AvType, prop=None, v=None):
-        """
-            :param str AvType: 'E' or 'F'
-            :param numpy.array prop: 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 prop is None:
+        assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
+        if tensorType == -1:
             return None
 
-        if Utils.isScalar(prop):
-            Av = getattr(self, 'ave'+AvType+'2CC')
+        if tensorType == 0:
+            Av = getattr(self, 'ave'+projType+'2CC')
             V = Utils.sdiag(self.vol)
             ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1))
-            if v is None:
-                return self.dim * Av.T * V * ones
-            return Utils.sdiag(v) * self.dim * Av.T * V * ones
+            # if v is None:
+            #     return self.dim * Av.T * V * ones
+            def scalarInnerProductDeriv(v):
+                return Utils.sdiag(v) * self.dim * Av.T * V * ones
+            return Utils.DerivOperator(scalarInnerProductDeriv)
 
-        if prop.size == self.nC:
-            Av = getattr(self, 'ave'+AvType+'2CC')
+        if tensorType == 1:
+            Av = getattr(self, 'ave'+projType+'2CC')
             V = Utils.sdiag(self.vol)
-            if v is None:
-                return self.dim * Av.T * V
-            return Utils.sdiag(v) * self.dim * Av.T * V
+            def isotropicInnerProductDeriv(v):
+                return Utils.sdiag(v) * self.dim * Av.T * V
+            return Utils.DerivOperator(isotropicInnerProductDeriv)
 
-        if prop.size == self.nC*self.dim: # anisotropic
-            Av = getattr(self, 'ave'+AvType+'2CCV')
+        if tensorType == 2: # anisotropic
+            Av = getattr(self, 'ave'+projType+'2CCV')
             V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
-            if v is None:
-                return Av.T * V
-            return Utils.sdiag(v) * Av.T * V
+            def anisotropicInnerProductDeriv(v):
+                return Utils.sdiag(v) * Av.T * V
+            return Utils.DerivOperator(anisotropicInnerProductDeriv)
 
 
 
diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py
index 1da88173..b574c8f1 100644
--- a/SimPEG/Tests/test_massMatricesDerivs.py
+++ b/SimPEG/Tests/test_massMatricesDerivs.py
@@ -6,130 +6,93 @@ from TestUtils import checkDerivative
 
 class TestInnerProductsDerivs(unittest.TestCase):
 
-    def doTestFace(self, h, rep, vec, fast):
+    def doTestFace(self, h, rep, fast):
         mesh = Mesh.TensorMesh(h)
         v = np.random.rand(mesh.nF)
+        sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
+        Md = mesh.getFaceInnerProductDeriv(Utils.TensorType(mesh, sig), doFast=fast)
         def fun(sig):
             M = mesh.getFaceInnerProduct(sig)
-            if vec:
-                Md = mesh.getFaceInnerProductDeriv(sig, v=v, doFast=fast)
-                return M*v, Md
-            Md = mesh.getFaceInnerProductDeriv(sig, doFast=fast)
-            return M*v, Utils.sdiag(v)*Md
-        sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
+            return M*v, Md*v
         return checkDerivative(fun, sig, num=5, plotIt=False)
 
-    def doTestEdge(self, h, rep, vec, fast):
+    def doTestEdge(self, h, rep, fast):
         mesh = Mesh.TensorMesh(h)
         v = np.random.rand(mesh.nE)
+        sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
+        Md = mesh.getEdgeInnerProductDeriv(Utils.TensorType(mesh, sig), doFast=fast)
         def fun(sig):
             M = mesh.getEdgeInnerProduct(sig)
-            if vec:
-                Md = mesh.getEdgeInnerProductDeriv(sig, v=v, doFast=fast)
-                return M*v, Md
-            Md = mesh.getEdgeInnerProductDeriv(sig, doFast=fast)
-            return M*v, Utils.sdiag(v)*Md
-        sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
+            return M*v, Md*v
         return checkDerivative(fun, sig, num=5, plotIt=False)
 
     def test_FaceIP_1D_float(self):
-        self.assertTrue(self.doTestFace([10],0,True, False))
+        self.assertTrue(self.doTestFace([10],0, False))
     def test_FaceIP_2D_float(self):
-        self.assertTrue(self.doTestFace([10, 4],0,True, False))
+        self.assertTrue(self.doTestFace([10, 4],0, False))
     def test_FaceIP_3D_float(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],0,True, False))
+        self.assertTrue(self.doTestFace([10, 4, 5],0, False))
     def test_FaceIP_1D_isotropic(self):
-        self.assertTrue(self.doTestFace([10],1,True, False))
+        self.assertTrue(self.doTestFace([10],1, False))
     def test_FaceIP_2D_isotropic(self):
-        self.assertTrue(self.doTestFace([10, 4],1,True, False))
+        self.assertTrue(self.doTestFace([10, 4],1, False))
     def test_FaceIP_3D_isotropic(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],1,True, False))
+        self.assertTrue(self.doTestFace([10, 4, 5],1, False))
     def test_FaceIP_2D_anisotropic(self):
-        self.assertTrue(self.doTestFace([10, 4],2,True, False))
+        self.assertTrue(self.doTestFace([10, 4],2, False))
     def test_FaceIP_3D_anisotropic(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],3,True, False))
+        self.assertTrue(self.doTestFace([10, 4, 5],3, False))
     def test_FaceIP_2D_tensor(self):
-        self.assertTrue(self.doTestFace([10, 4],3,True, False))
+        self.assertTrue(self.doTestFace([10, 4],3, False))
     def test_FaceIP_3D_tensor(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],6,True, False))
+        self.assertTrue(self.doTestFace([10, 4, 5],6, False))
 
     def test_FaceIP_1D_float_fast(self):
-        self.assertTrue(self.doTestFace([10],0, False, True))
+        self.assertTrue(self.doTestFace([10],0, True))
     def test_FaceIP_2D_float_fast(self):
-        self.assertTrue(self.doTestFace([10, 4],0, False, True))
+        self.assertTrue(self.doTestFace([10, 4],0, True))
     def test_FaceIP_3D_float_fast(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],0, False, True))
+        self.assertTrue(self.doTestFace([10, 4, 5],0, True))
     def test_FaceIP_1D_isotropic_fast(self):
-        self.assertTrue(self.doTestFace([10],1, False, True))
+        self.assertTrue(self.doTestFace([10],1, True))
     def test_FaceIP_2D_isotropic_fast(self):
-        self.assertTrue(self.doTestFace([10, 4],1, False, True))
+        self.assertTrue(self.doTestFace([10, 4],1, True))
     def test_FaceIP_3D_isotropic_fast(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],1, False, True))
+        self.assertTrue(self.doTestFace([10, 4, 5],1, True))
     def test_FaceIP_2D_anisotropic_fast(self):
-        self.assertTrue(self.doTestFace([10, 4],2, False, True))
+        self.assertTrue(self.doTestFace([10, 4],2, True))
     def test_FaceIP_3D_anisotropic_fast(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],3, False, True))
-
-    def test_FaceIP_1D_float_fast_vec(self):
-        self.assertTrue(self.doTestFace([10],0, True, True))
-    def test_FaceIP_2D_float_fast_vec(self):
-        self.assertTrue(self.doTestFace([10, 4],0, True, True))
-    def test_FaceIP_3D_float_fast_vec(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],0, True, True))
-    def test_FaceIP_1D_isotropic_fast_vec(self):
-        self.assertTrue(self.doTestFace([10],1, True, True))
-    def test_FaceIP_2D_isotropic_fast_vec(self):
-        self.assertTrue(self.doTestFace([10, 4],1, True, True))
-    def test_FaceIP_3D_isotropic_fast_vec(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],1, True, True))
-    def test_FaceIP_2D_anisotropic_fast_vec(self):
-        self.assertTrue(self.doTestFace([10, 4],2, True, True))
-    def test_FaceIP_3D_anisotropic_fast_vec(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],3, True, True))
+        self.assertTrue(self.doTestFace([10, 4, 5],3, True))
 
     def test_EdgeIP_2D_float(self):
-        self.assertTrue(self.doTestEdge([10, 4],0,True, False))
+        self.assertTrue(self.doTestEdge([10, 4],0, False))
     def test_EdgeIP_3D_float(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],0,True, False))
+        self.assertTrue(self.doTestEdge([10, 4, 5],0, False))
     def test_EdgeIP_2D_isotropic(self):
-        self.assertTrue(self.doTestEdge([10, 4],1,True, False))
+        self.assertTrue(self.doTestEdge([10, 4],1, False))
     def test_EdgeIP_3D_isotropic(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],1,True, False))
+        self.assertTrue(self.doTestEdge([10, 4, 5],1, False))
     def test_EdgeIP_2D_anisotropic(self):
-        self.assertTrue(self.doTestEdge([10, 4],2,True, False))
+        self.assertTrue(self.doTestEdge([10, 4],2, False))
     def test_EdgeIP_3D_anisotropic(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],3,True, False))
+        self.assertTrue(self.doTestEdge([10, 4, 5],3, False))
     def test_EdgeIP_2D_tensor(self):
-        self.assertTrue(self.doTestEdge([10, 4],3,True, False))
+        self.assertTrue(self.doTestEdge([10, 4],3, False))
     def test_EdgeIP_3D_tensor(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],6,True, False))
+        self.assertTrue(self.doTestEdge([10, 4, 5],6, False))
 
     def test_EdgeIP_2D_float_fast(self):
-        self.assertTrue(self.doTestEdge([10, 4],0, False, True))
+        self.assertTrue(self.doTestEdge([10, 4],0, True))
     def test_EdgeIP_3D_float_fast(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],0, False, True))
+        self.assertTrue(self.doTestEdge([10, 4, 5],0, True))
     def test_EdgeIP_2D_isotropic_fast(self):
-        self.assertTrue(self.doTestEdge([10, 4],1, False, True))
+        self.assertTrue(self.doTestEdge([10, 4],1, True))
     def test_EdgeIP_3D_isotropic_fast(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],1, False, True))
+        self.assertTrue(self.doTestEdge([10, 4, 5],1, True))
     def test_EdgeIP_2D_anisotropic_fast(self):
-        self.assertTrue(self.doTestEdge([10, 4],2, False, True))
+        self.assertTrue(self.doTestEdge([10, 4],2, True))
     def test_EdgeIP_3D_anisotropic_fast(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],3, False, True))
-
-    def test_EdgeIP_2D_float_fast_vec(self):
-        self.assertTrue(self.doTestEdge([10, 4],0, True, True))
-    def test_EdgeIP_3D_float_fast_vec(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],0, True, True))
-    def test_EdgeIP_2D_isotropic_fast_vec(self):
-        self.assertTrue(self.doTestEdge([10, 4],1, True, True))
-    def test_EdgeIP_3D_isotropic_fast_vec(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],1, True, True))
-    def test_EdgeIP_2D_anisotropic_fast_vec(self):
-        self.assertTrue(self.doTestEdge([10, 4],2, True, True))
-    def test_EdgeIP_3D_anisotropic_fast_vec(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],3, True, True))
-
+        self.assertTrue(self.doTestEdge([10, 4, 5],3, True))
 
 
 if __name__ == '__main__':
diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py
index 6318527a..50cc0574 100644
--- a/SimPEG/Tests/test_utils.py
+++ b/SimPEG/Tests/test_utils.py
@@ -160,7 +160,7 @@ class TestSequenceFunctions(unittest.TestCase):
             Z = B2*A - sp.identity(M.nC*2)
             self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
 
-    def test_tensorType2D(self):
+    def test_TensorType2D(self):
         M = Mesh.TensorMesh([6, 6])
         a1 = np.random.rand(M.nC)
         a2 = np.random.rand(M.nC)
@@ -170,12 +170,12 @@ class TestSequenceFunctions(unittest.TestCase):
         prop3 = np.c_[a1, a2, a3]
 
         for ii, prop in enumerate([4, prop1, prop2, prop3]):
-            self.assertTrue(tensorType(M, prop) == ii)
+            self.assertTrue(TensorType(M, prop) == ii)
 
-        self.assertRaises(Exception, tensorType, M, np.c_[a1, a2, a3, a3])
-        self.assertTrue(tensorType(M, None) == -1)
+        self.assertRaises(Exception, TensorType, M, np.c_[a1, a2, a3, a3])
+        self.assertTrue(TensorType(M, None) == -1)
 
-    def test_tensorType3D(self):
+    def test_TensorType3D(self):
         M = Mesh.TensorMesh([6, 6, 7])
         a1 = np.random.rand(M.nC)
         a2 = np.random.rand(M.nC)
@@ -188,10 +188,10 @@ class TestSequenceFunctions(unittest.TestCase):
         prop3 = np.c_[a1, a2, a3, a4, a5, a6]
 
         for ii, prop in enumerate([4, prop1, prop2, prop3]):
-            self.assertTrue(tensorType(M, prop) == ii)
+            self.assertTrue(TensorType(M, prop) == ii)
 
-        self.assertRaises(Exception, tensorType, M, np.c_[a1, a2, a3, a3])
-        self.assertTrue(tensorType(M, None) == -1)
+        self.assertRaises(Exception, TensorType, M, np.c_[a1, a2, a3, a3])
+        self.assertTrue(TensorType(M, None) == -1)
 
 
     def test_invPropertyTensor3D(self):
diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py
index a93e6f07..1a4c4a0e 100644
--- a/SimPEG/Utils/matutils.py
+++ b/SimPEG/Utils/matutils.py
@@ -251,25 +251,34 @@ def inv2X2BlockDiagonal(a11, a12, a21, a22, returnMatrix=True):
     return sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))),
                       sp.hstack((sdiag(b21), sdiag(b22)))))
 
-def tensorType(M, tensor):
-    if tensor is None:  # default is ones
-        return -1
-
-    if isScalar(tensor):
-        return 0
-
-    if tensor.size == M.nC:
-        return 1
-
-    if ((M.dim == 2 and tensor.size == M.nC*2) or
-        (M.dim == 3 and tensor.size == M.nC*3)):
-        return 2
-
-    if ((M.dim == 2 and tensor.size == M.nC*3) or
-        (M.dim == 3 and tensor.size == M.nC*6)):
-        return 3
-
-    raise Exception('Unexpected shape of tensor')
+class TensorType(object):
+    def __init__(self, M, tensor):
+        if tensor is None:  # default is ones
+            self._tt  = -1
+            self._tts = 'none'
+        elif isScalar(tensor):
+            self._tt  = 0
+            self._tts = 'scalar'
+        elif tensor.size == M.nC:
+            self._tt  = 1
+            self._tts = 'isotropic'
+        elif ((M.dim == 2 and tensor.size == M.nC*2) or
+            (M.dim == 3 and tensor.size == M.nC*3)):
+            self._tt  = 2
+            self._tts = 'anisotropic'
+        elif ((M.dim == 2 and tensor.size == M.nC*3) or
+            (M.dim == 3 and tensor.size == M.nC*6)):
+            self._tt  = 3
+            self._tts = 'tensor'
+        else:
+            raise Exception('Unexpected shape of tensor')
+    def __str__(self):
+        return 'TensorType[%i]: %s' % (self._tt, self._tts)
+    def __eq__(self, v): return self._tt == v
+    def __le__(self, v): return self._tt <= v
+    def __ge__(self, v): return self._tt >= v
+    def __lt__(self, v): return self._tt < v
+    def __gt__(self, v): return self._tt > v
 
 def makePropertyTensor(M, tensor):
     if tensor is None:  # default is ones
@@ -278,7 +287,7 @@ def makePropertyTensor(M, tensor):
     if isScalar(tensor):
         tensor = tensor * np.ones(M.nC)
 
-    propType = tensorType(M, tensor)
+    propType = TensorType(M, tensor)
     if propType == 1: # Isotropic!
         Sigma = sp.kron(sp.identity(M.dim), sdiag(mkvc(tensor)))
     elif propType == 2: # Diagonal tensor
@@ -302,7 +311,7 @@ def makePropertyTensor(M, tensor):
 
 def invPropertyTensor(M, tensor, returnMatrix=False):
 
-    propType = tensorType(M, tensor)
+    propType = TensorType(M, tensor)
 
     if isScalar(tensor):
         T = 1./tensor
@@ -341,3 +350,9 @@ class SimPEGLinearOperator(LinearOperator):
     def T(self):
         return self.__class__((self.shape[1],self.shape[0]),self.rmatvec,rmatvec=self.matvec,matmat=self.matmat)
 
+
+class DerivOperator(object):
+    def __init__(self, f):
+        self.f   = f
+    def __mul__(self, v):
+        return self.f(v)

From 1d5ecbf303b32c66ef4426172034b87ee87fb890 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Mon, 16 Jun 2014 14:01:03 -0600
Subject: [PATCH 02/14] Edge inner products updates for 1D

---
 SimPEG/Mesh/DiffOperators.py            |  2 +-
 SimPEG/Mesh/InnerProducts.py            | 17 ++++++++++++++---
 SimPEG/Tests/test_massMatricesDerivs.py |  8 ++++++++
 3 files changed, 23 insertions(+), 4 deletions(-)

diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py
index f8d448d4..966397b3 100644
--- a/SimPEG/Mesh/DiffOperators.py
+++ b/SimPEG/Mesh/DiffOperators.py
@@ -620,7 +620,7 @@ class DiffOperators(object):
             # The number of cell centers in each direction
             n = self.vnC
             if(self.dim == 1):
-                raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
+                self._aveE2CC = speye(n[0])
             elif(self.dim == 2):
                 self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])),
                                                sp.kron(speye(n[1]), av(n[0]))), format="csr")
diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py
index 026e406d..1a1de3ed 100644
--- a/SimPEG/Mesh/InnerProducts.py
+++ b/SimPEG/Mesh/InnerProducts.py
@@ -102,8 +102,8 @@ class InnerProducts(object):
 
         elif projType == 'E':
             locs = {
-                    '000': [  None  , ('eX0', 'eY0'), ('eX0', 'eY0', 'eZ0')],
-                    '100': [  None  , ('eX0', 'eY1'), ('eX0', 'eY1', 'eZ1')],
+                    '000': [('eX0',), ('eX0', 'eY0'), ('eX0', 'eY0', 'eZ0')],
+                    '100': [('eX0',), ('eX0', 'eY1'), ('eX0', 'eY1', 'eZ1')],
                     '010': [  None  , ('eX1', 'eY0'), ('eX1', 'eY0', 'eZ2')],
                     '110': [  None  , ('eX1', 'eY1'), ('eX1', 'eY1', 'eZ3')],
                     '001': [  None  ,      None     , ('eX2', 'eY2', 'eZ0')],
@@ -112,7 +112,7 @@ class InnerProducts(object):
                     '111': [  None  ,      None     , ('eX3', 'eY3', 'eZ3')]
                    }
             if d == 1:
-                raise NotImplementedError('getEdgeInnerProduct not implemented for 1D')
+                proj = _getEdgePx(self)
             elif d == 2:
                 proj = _getEdgePxx(self)
             elif d == 3:
@@ -295,6 +295,10 @@ def _getFacePxxx(M):
 
     return _getFacePxxx_Rectangular(M)
 
+def _getEdgePx(M):
+    assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh'
+    return _getEdgePx_Rectangular(M)
+
 def _getEdgePxx(M):
     if M._meshType == 'TREE':
         return M._getEdgePxx
@@ -449,6 +453,13 @@ def _getFacePxxx_Rectangular(M):
         return PXXX
     return Pxxx
 
+def _getEdgePx_Rectangular(M):
+    """Returns a function for creating projection matrices"""
+    def Px(xEdge):
+        assert xEdge == 'eX0', 'xEdge = %s, not eX0' % xEdge
+        return sp.identity(M.nC)
+    return Px
+
 def _getEdgePxx_Rectangular(M):
     i, j = np.int64(range(M.nCx)), np.int64(range(M.nCy))
 
diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py
index b574c8f1..72e36f3e 100644
--- a/SimPEG/Tests/test_massMatricesDerivs.py
+++ b/SimPEG/Tests/test_massMatricesDerivs.py
@@ -64,10 +64,14 @@ class TestInnerProductsDerivs(unittest.TestCase):
     def test_FaceIP_3D_anisotropic_fast(self):
         self.assertTrue(self.doTestFace([10, 4, 5],3, True))
 
+    def test_EdgeIP_1D_float(self):
+        self.assertTrue(self.doTestEdge([10],0, False))
     def test_EdgeIP_2D_float(self):
         self.assertTrue(self.doTestEdge([10, 4],0, False))
     def test_EdgeIP_3D_float(self):
         self.assertTrue(self.doTestEdge([10, 4, 5],0, False))
+    def test_EdgeIP_1D_isotropic(self):
+        self.assertTrue(self.doTestEdge([10],1, False))
     def test_EdgeIP_2D_isotropic(self):
         self.assertTrue(self.doTestEdge([10, 4],1, False))
     def test_EdgeIP_3D_isotropic(self):
@@ -81,10 +85,14 @@ class TestInnerProductsDerivs(unittest.TestCase):
     def test_EdgeIP_3D_tensor(self):
         self.assertTrue(self.doTestEdge([10, 4, 5],6, False))
 
+    def test_EdgeIP_1D_float_fast(self):
+        self.assertTrue(self.doTestEdge([10],0, True))
     def test_EdgeIP_2D_float_fast(self):
         self.assertTrue(self.doTestEdge([10, 4],0, True))
     def test_EdgeIP_3D_float_fast(self):
         self.assertTrue(self.doTestEdge([10, 4, 5],0, True))
+    def test_EdgeIP_1D_isotropic_fast(self):
+        self.assertTrue(self.doTestEdge([10],1, True))
     def test_EdgeIP_2D_isotropic_fast(self):
         self.assertTrue(self.doTestEdge([10, 4],1, True))
     def test_EdgeIP_3D_isotropic_fast(self):

From fb90faae50fca9b622f46907c4dc14c70d5e0851 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Mon, 16 Jun 2014 14:16:19 -0600
Subject: [PATCH 03/14] more object oriented innerProducts

---
 SimPEG/Mesh/InnerProducts.py | 439 ++++++++++++++++-------------------
 1 file changed, 199 insertions(+), 240 deletions(-)

diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py
index 1a1de3ed..6ffb20b9 100644
--- a/SimPEG/Mesh/InnerProducts.py
+++ b/SimPEG/Mesh/InnerProducts.py
@@ -93,12 +93,7 @@ class InnerProducts(object):
                     '011': [  None  ,      None     , ('fXm', 'fYp', 'fZp')],
                     '111': [  None  ,      None     , ('fXp', 'fYp', 'fZp')]
                    }
-            if d == 1:
-                proj = _getFacePx(self)
-            elif d == 2:
-                proj = _getFacePxx(self)
-            elif d == 3:
-                proj = _getFacePxxx(self)
+            proj = getattr(self, '_getFaceP' + ('x'*d))()
 
         elif projType == 'E':
             locs = {
@@ -111,12 +106,7 @@ class InnerProducts(object):
                     '011': [  None  ,      None     , ('eX3', 'eY2', 'eZ2')],
                     '111': [  None  ,      None     , ('eX3', 'eY3', 'eZ3')]
                    }
-            if d == 1:
-                proj = _getEdgePx(self)
-            elif d == 2:
-                proj = _getEdgePxx(self)
-            elif d == 3:
-                proj = _getEdgePxxx(self)
+            proj = getattr(self, '_getEdgeP' + ('x'*d))()
 
         return [V*proj(*locs[node][d-1]) for node in nodes]
 
@@ -260,283 +250,252 @@ class InnerProducts(object):
 
         return dMdm
 
-# ------------------------ Geometries ------------------------------
-#
-#
-#         node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1)
-#              /                                    /
-#             /                                    / |
-#         edge3(i,j,k)     face1(i,j,k)        edge3(i,j+1,k)
-#           /                                    /   |
-#          /                                    /    |
-#    node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k)
-#         |                                     |    |
-#         |                                     |   node(i+1,j+1,k+1)
-#         |                                     |    /
-#    edge1(i,j,k)      face3(i,j,k)        edge1(i,j+1,k)
-#         |                                     |  /
-#         |                                     | /
-#         |                                     |/
-#    node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
+    # ------------------------ Geometries ------------------------------
+    #
+    #
+    #         node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1)
+    #              /                                    /
+    #             /                                    / |
+    #         edge3(i,j,k)     face1(i,j,k)        edge3(i,j+1,k)
+    #           /                                    /   |
+    #          /                                    /    |
+    #    node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k)
+    #         |                                     |    |
+    #         |                                     |   node(i+1,j+1,k+1)
+    #         |                                     |    /
+    #    edge1(i,j,k)      face3(i,j,k)        edge1(i,j+1,k)
+    #         |                                     |  /
+    #         |                                     | /
+    #         |                                     |/
+    #    node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
 
-def _getFacePx(M):
-    assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh'
-    return _getFacePx_Rectangular(M)
 
-def _getFacePxx(M):
-    if M._meshType == 'TREE':
-        return M._getFacePxx
-
-    return _getFacePxx_Rectangular(M)
-
-def _getFacePxxx(M):
-    if M._meshType == 'TREE':
-        return M._getFacePxxx
-
-    return _getFacePxxx_Rectangular(M)
-
-def _getEdgePx(M):
-    assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh'
-    return _getEdgePx_Rectangular(M)
-
-def _getEdgePxx(M):
-    if M._meshType == 'TREE':
-        return M._getEdgePxx
-
-    return _getEdgePxx_Rectangular(M)
-
-def _getEdgePxxx(M):
-    if M._meshType == 'TREE':
-        return M._getEdgePxxx
-
-    return _getEdgePxxx_Rectangular(M)
-
-def _getFacePx_Rectangular(M):
-    """Returns a function for creating projection matrices
-
-    """
-    ii = np.int64(range(M.nCx))
-
-    def Px(xFace):
-        """
-            xFace is 'fXp' or 'fXm'
-        """
-        posFx = 0 if xFace == 'fXm' else 1
-        IND = ii + posFx
-        PX = sp.csr_matrix((np.ones(M.nC), (range(M.nC), IND)), shape=(M.nC, M.nF))
-        return PX
-
-    return Px
-
-def _getFacePxx_Rectangular(M):
-    """returns a function for creating projection matrices
-
-        Mats takes you from faces a subset of all faces on only the
-        faces that you ask for.
-
-        These are centered around a single nodes.
-
-        For example, if this was your entire mesh:
-
-                        f3(Yp)
-                  2_______________3
-                  |               |
-                  |               |
-                  |               |
-          f0(Xm)  |       x       |  f1(Xp)
-                  |               |
-                  |               |
-                  |_______________|
-                  0               1
-                        f2(Ym)
-
-        Pxx('fXm','fYm') = | 1, 0, 0, 0 |
-                           | 0, 0, 1, 0 |
-
-        Pxx('fXp','fYm') = | 0, 1, 0, 0 |
-                           | 0, 0, 1, 0 |
+    def _getFacePx(M):
+        """Returns a function for creating projection matrices
 
         """
-    i, j = np.int64(range(M.nCx)), np.int64(range(M.nCy))
+        ii = np.arange(M.nCx)
 
-    iijj = ndgrid(i, j)
-    ii, jj = iijj[:, 0], iijj[:, 1]
+        def Px(xFace):
+            """
+                xFace is 'fXp' or 'fXm'
+            """
+            posFx = 0 if xFace == 'fXm' else 1
+            IND = ii + posFx
+            PX = sp.csr_matrix((np.ones(M.nC), (range(M.nC), IND)), shape=(M.nC, M.nF))
+            return PX
 
-    if M._meshType == 'LRM':
-        fN1 = M.r(M.normals, 'F', 'Fx', 'M')
-        fN2 = M.r(M.normals, 'F', 'Fy', 'M')
+        return Px
 
-    def Pxx(xFace, yFace):
-        """
-            xFace is 'fXp' or 'fXm'
-            yFace is 'fYp' or 'fYm'
-        """
-        # no | node      | f1     | f2
-        # 00 | i  ,j     | i  , j | i, j
-        # 10 | i+1,j     | i+1, j | i, j
-        # 01 | i  ,j+1   | i  , j | i, j+1
-        # 11 | i+1,j+1   | i+1, j | i, j+1
+    def _getFacePxx(M):
+        """returns a function for creating projection matrices
 
-        posFx = 0 if xFace == 'fXm' else 1
-        posFy = 0 if yFace == 'fYm' else 1
+            Mats takes you from faces a subset of all faces on only the
+            faces that you ask for.
 
-        ind1 = sub2ind(M.vnFx, np.c_[ii + posFx, jj])
-        ind2 = sub2ind(M.vnFy, np.c_[ii, jj + posFy]) + M.nFx
+            These are centered around a single nodes.
 
-        IND = np.r_[ind1, ind2].flatten()
+            For example, if this was your entire mesh:
 
-        PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF))
+                            f3(Yp)
+                      2_______________3
+                      |               |
+                      |               |
+                      |               |
+              f0(Xm)  |       x       |  f1(Xp)
+                      |               |
+                      |               |
+                      |_______________|
+                      0               1
+                            f2(Ym)
+
+            Pxx('fXm','fYm') = | 1, 0, 0, 0 |
+                               | 0, 0, 1, 0 |
+
+            Pxx('fXp','fYm') = | 0, 1, 0, 0 |
+                               | 0, 0, 1, 0 |
+
+            """
+        i, j = np.arange(M.nCx), np.arange(M.nCy)
+
+        iijj = ndgrid(i, j)
+        ii, jj = iijj[:, 0], iijj[:, 1]
 
         if M._meshType == 'LRM':
-            I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]),
-                                       getSubArray(fN2[0], [i, j + posFy]), getSubArray(fN2[1], [i, j + posFy]))
-            PXX = I2x2 * PXX
+            fN1 = M.r(M.normals, 'F', 'Fx', 'M')
+            fN2 = M.r(M.normals, 'F', 'Fy', 'M')
 
-        return PXX
+        def Pxx(xFace, yFace):
+            """
+                xFace is 'fXp' or 'fXm'
+                yFace is 'fYp' or 'fYm'
+            """
+            # no | node      | f1     | f2
+            # 00 | i  ,j     | i  , j | i, j
+            # 10 | i+1,j     | i+1, j | i, j
+            # 01 | i  ,j+1   | i  , j | i, j+1
+            # 11 | i+1,j+1   | i+1, j | i, j+1
 
-    return Pxx
+            posFx = 0 if xFace == 'fXm' else 1
+            posFy = 0 if yFace == 'fYm' else 1
 
-def _getFacePxxx_Rectangular(M):
-    """returns a function for creating projection matrices
+            ind1 = sub2ind(M.vnFx, np.c_[ii + posFx, jj])
+            ind2 = sub2ind(M.vnFy, np.c_[ii, jj + posFy]) + M.nFx
 
-        Mats takes you from faces a subset of all faces on only the
-        faces that you ask for.
+            IND = np.r_[ind1, ind2].flatten()
 
-        These are centered around a single nodes.
-    """
+            PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF))
 
-    i, j, k = np.int64(range(M.nCx)), np.int64(range(M.nCy)), np.int64(range(M.nCz))
+            if M._meshType == 'LRM':
+                I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]),
+                                           getSubArray(fN2[0], [i, j + posFy]), getSubArray(fN2[1], [i, j + posFy]))
+                PXX = I2x2 * PXX
 
-    iijjkk = ndgrid(i, j, k)
-    ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
+            return PXX
 
-    if M._meshType == 'LRM':
-        fN1 = M.r(M.normals, 'F', 'Fx', 'M')
-        fN2 = M.r(M.normals, 'F', 'Fy', 'M')
-        fN3 = M.r(M.normals, 'F', 'Fz', 'M')
+        return Pxx
 
-    def Pxxx(xFace, yFace, zFace):
-        """
-            xFace is 'fXp' or 'fXm'
-            yFace is 'fYp' or 'fYm'
-            zFace is 'fZp' or 'fZm'
+    def _getFacePxxx(M):
+        """returns a function for creating projection matrices
+
+            Mats takes you from faces a subset of all faces on only the
+            faces that you ask for.
+
+            These are centered around a single nodes.
         """
 
-        # no  | node        | f1        | f2        | f3
-        # 000 | i  ,j  ,k   | i  , j, k | i, j  , k | i, j, k
-        # 100 | i+1,j  ,k   | i+1, j, k | i, j  , k | i, j, k
-        # 010 | i  ,j+1,k   | i  , j, k | i, j+1, k | i, j, k
-        # 110 | i+1,j+1,k   | i+1, j, k | i, j+1, k | i, j, k
-        # 001 | i  ,j  ,k+1 | i  , j, k | i, j  , k | i, j, k+1
-        # 101 | i+1,j  ,k+1 | i+1, j, k | i, j  , k | i, j, k+1
-        # 011 | i  ,j+1,k+1 | i  , j, k | i, j+1, k | i, j, k+1
-        # 111 | i+1,j+1,k+1 | i+1, j, k | i, j+1, k | i, j, k+1
+        i, j, k = np.arange(M.nCx), np.arange(M.nCy), np.arange(M.nCz)
 
-        posX = 0 if xFace == 'fXm' else 1
-        posY = 0 if yFace == 'fYm' else 1
-        posZ = 0 if zFace == 'fZm' else 1
-
-        ind1 = sub2ind(M.vnFx, np.c_[ii + posX, jj, kk])
-        ind2 = sub2ind(M.vnFy, np.c_[ii, jj + posY, kk]) + M.nFx
-        ind3 = sub2ind(M.vnFz, np.c_[ii, jj, kk + posZ]) + M.nFx + M.nFy
-
-        IND = np.r_[ind1, ind2, ind3].flatten()
-
-        PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr()
+        iijjkk = ndgrid(i, j, k)
+        ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
 
         if M._meshType == 'LRM':
-            I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]),
-                                       getSubArray(fN2[0], [i, j + posY, k]), getSubArray(fN2[1], [i, j + posY, k]), getSubArray(fN2[2], [i, j + posY, k]),
-                                       getSubArray(fN3[0], [i, j, k + posZ]), getSubArray(fN3[1], [i, j, k + posZ]), getSubArray(fN3[2], [i, j, k + posZ]))
-            PXXX = I3x3 * PXXX
+            fN1 = M.r(M.normals, 'F', 'Fx', 'M')
+            fN2 = M.r(M.normals, 'F', 'Fy', 'M')
+            fN3 = M.r(M.normals, 'F', 'Fz', 'M')
 
-        return PXXX
-    return Pxxx
+        def Pxxx(xFace, yFace, zFace):
+            """
+                xFace is 'fXp' or 'fXm'
+                yFace is 'fYp' or 'fYm'
+                zFace is 'fZp' or 'fZm'
+            """
 
-def _getEdgePx_Rectangular(M):
-    """Returns a function for creating projection matrices"""
-    def Px(xEdge):
-        assert xEdge == 'eX0', 'xEdge = %s, not eX0' % xEdge
-        return sp.identity(M.nC)
-    return Px
+            # no  | node        | f1        | f2        | f3
+            # 000 | i  ,j  ,k   | i  , j, k | i, j  , k | i, j, k
+            # 100 | i+1,j  ,k   | i+1, j, k | i, j  , k | i, j, k
+            # 010 | i  ,j+1,k   | i  , j, k | i, j+1, k | i, j, k
+            # 110 | i+1,j+1,k   | i+1, j, k | i, j+1, k | i, j, k
+            # 001 | i  ,j  ,k+1 | i  , j, k | i, j  , k | i, j, k+1
+            # 101 | i+1,j  ,k+1 | i+1, j, k | i, j  , k | i, j, k+1
+            # 011 | i  ,j+1,k+1 | i  , j, k | i, j+1, k | i, j, k+1
+            # 111 | i+1,j+1,k+1 | i+1, j, k | i, j+1, k | i, j, k+1
 
-def _getEdgePxx_Rectangular(M):
-    i, j = np.int64(range(M.nCx)), np.int64(range(M.nCy))
+            posX = 0 if xFace == 'fXm' else 1
+            posY = 0 if yFace == 'fYm' else 1
+            posZ = 0 if zFace == 'fZm' else 1
 
-    iijj = ndgrid(i, j)
-    ii, jj = iijj[:, 0], iijj[:, 1]
+            ind1 = sub2ind(M.vnFx, np.c_[ii + posX, jj, kk])
+            ind2 = sub2ind(M.vnFy, np.c_[ii, jj + posY, kk]) + M.nFx
+            ind3 = sub2ind(M.vnFz, np.c_[ii, jj, kk + posZ]) + M.nFx + M.nFy
 
-    if M._meshType == 'LRM':
-        eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
-        eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
+            IND = np.r_[ind1, ind2, ind3].flatten()
 
-    def Pxx(xEdge, yEdge):
-        # no | node      | e1      | e2
-        # 00 | i  ,j     | i  ,j   | i  ,j
-        # 10 | i+1,j     | i  ,j   | i+1,j
-        # 01 | i  ,j+1   | i  ,j+1 | i  ,j
-        # 11 | i+1,j+1   | i  ,j+1 | i+1,j
-        posX = 0 if xEdge == 'eX0' else 1
-        posY = 0 if yEdge == 'eY0' else 1
+            PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr()
 
-        ind1 = sub2ind(M.vnEx, np.c_[ii, jj + posX])
-        ind2 = sub2ind(M.vnEy, np.c_[ii + posY, jj]) + M.nEx
+            if M._meshType == 'LRM':
+                I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]),
+                                           getSubArray(fN2[0], [i, j + posY, k]), getSubArray(fN2[1], [i, j + posY, k]), getSubArray(fN2[2], [i, j + posY, k]),
+                                           getSubArray(fN3[0], [i, j, k + posZ]), getSubArray(fN3[1], [i, j, k + posZ]), getSubArray(fN3[2], [i, j, k + posZ]))
+                PXXX = I3x3 * PXXX
 
-        IND = np.r_[ind1, ind2].flatten()
+            return PXXX
+        return Pxxx
 
-        PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr()
+    def _getEdgePx(M):
+        """Returns a function for creating projection matrices"""
+        def Px(xEdge):
+            assert xEdge == 'eX0', 'xEdge = %s, not eX0' % xEdge
+            return sp.identity(M.nC)
+        return Px
+
+    def _getEdgePxx(M):
+        i, j = np.arange(M.nCx), np.arange(M.nCy)
+
+        iijj = ndgrid(i, j)
+        ii, jj = iijj[:, 0], iijj[:, 1]
 
         if M._meshType == 'LRM':
-            I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]),
-                                       getSubArray(eT2[0], [i + posY, j]), getSubArray(eT2[1], [i + posY, j]))
-            PXX = I2x2 * PXX
+            eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
+            eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
 
-        return PXX
-    return Pxx
+        def Pxx(xEdge, yEdge):
+            # no | node      | e1      | e2
+            # 00 | i  ,j     | i  ,j   | i  ,j
+            # 10 | i+1,j     | i  ,j   | i+1,j
+            # 01 | i  ,j+1   | i  ,j+1 | i  ,j
+            # 11 | i+1,j+1   | i  ,j+1 | i+1,j
+            posX = 0 if xEdge == 'eX0' else 1
+            posY = 0 if yEdge == 'eY0' else 1
 
-def _getEdgePxxx_Rectangular(M):
-    i, j, k = np.int64(range(M.nCx)), np.int64(range(M.nCy)), np.int64(range(M.nCz))
+            ind1 = sub2ind(M.vnEx, np.c_[ii, jj + posX])
+            ind2 = sub2ind(M.vnEy, np.c_[ii + posY, jj]) + M.nEx
 
-    iijjkk = ndgrid(i, j, k)
-    ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
+            IND = np.r_[ind1, ind2].flatten()
 
-    if M._meshType == 'LRM':
-        eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
-        eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
-        eT3 = M.r(M.tangents, 'E', 'Ez', 'M')
+            PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr()
 
-    def Pxxx(xEdge, yEdge, zEdge):
+            if M._meshType == 'LRM':
+                I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]),
+                                           getSubArray(eT2[0], [i + posY, j]), getSubArray(eT2[1], [i + posY, j]))
+                PXX = I2x2 * PXX
 
-        # no  | node        | e1          | e2          | e3
-        # 000 | i  ,j  ,k   | i  ,j  ,k   | i  ,j  ,k   | i  ,j  ,k
-        # 100 | i+1,j  ,k   | i  ,j  ,k   | i+1,j  ,k   | i+1,j  ,k
-        # 010 | i  ,j+1,k   | i  ,j+1,k   | i  ,j  ,k   | i  ,j+1,k
-        # 110 | i+1,j+1,k   | i  ,j+1,k   | i+1,j  ,k   | i+1,j+1,k
-        # 001 | i  ,j  ,k+1 | i  ,j  ,k+1 | i  ,j  ,k+1 | i  ,j  ,k
-        # 101 | i+1,j  ,k+1 | i  ,j  ,k+1 | i+1,j  ,k+1 | i+1,j  ,k
-        # 011 | i  ,j+1,k+1 | i  ,j+1,k+1 | i  ,j  ,k+1 | i  ,j+1,k
-        # 111 | i+1,j+1,k+1 | i  ,j+1,k+1 | i+1,j  ,k+1 | i+1,j+1,k
+            return PXX
+        return Pxx
 
-        posX = [0,0] if xEdge == 'eX0' else [1, 0] if xEdge == 'eX1' else [0,1] if xEdge == 'eX2' else [1,1]
-        posY = [0,0] if yEdge == 'eY0' else [1, 0] if yEdge == 'eY1' else [0,1] if yEdge == 'eY2' else [1,1]
-        posZ = [0,0] if zEdge == 'eZ0' else [1, 0] if zEdge == 'eZ1' else [0,1] if zEdge == 'eZ2' else [1,1]
+    def _getEdgePxxx(M):
+        i, j, k = np.arange(M.nCx), np.arange(M.nCy), np.arange(M.nCz)
 
-        ind1 = sub2ind(M.vnEx, np.c_[ii, jj + posX[0], kk + posX[1]])
-        ind2 = sub2ind(M.vnEy, np.c_[ii + posY[0], jj, kk + posY[1]]) + M.nEx
-        ind3 = sub2ind(M.vnEz, np.c_[ii + posZ[0], jj + posZ[1], kk]) + M.nEx + M.nEy
-
-        IND = np.r_[ind1, ind2, ind3].flatten()
-
-        PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr()
+        iijjkk = ndgrid(i, j, k)
+        ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
 
         if M._meshType == 'LRM':
-            I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]),
-                                       getSubArray(eT2[0], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[1], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[2], [i + posY[0], j, k + posY[1]]),
-                                       getSubArray(eT3[0], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[1], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[2], [i + posZ[0], j + posZ[1], k]))
-            PXXX = I3x3 * PXXX
+            eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
+            eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
+            eT3 = M.r(M.tangents, 'E', 'Ez', 'M')
 
-        return PXXX
-    return Pxxx
+        def Pxxx(xEdge, yEdge, zEdge):
+
+            # no  | node        | e1          | e2          | e3
+            # 000 | i  ,j  ,k   | i  ,j  ,k   | i  ,j  ,k   | i  ,j  ,k
+            # 100 | i+1,j  ,k   | i  ,j  ,k   | i+1,j  ,k   | i+1,j  ,k
+            # 010 | i  ,j+1,k   | i  ,j+1,k   | i  ,j  ,k   | i  ,j+1,k
+            # 110 | i+1,j+1,k   | i  ,j+1,k   | i+1,j  ,k   | i+1,j+1,k
+            # 001 | i  ,j  ,k+1 | i  ,j  ,k+1 | i  ,j  ,k+1 | i  ,j  ,k
+            # 101 | i+1,j  ,k+1 | i  ,j  ,k+1 | i+1,j  ,k+1 | i+1,j  ,k
+            # 011 | i  ,j+1,k+1 | i  ,j+1,k+1 | i  ,j  ,k+1 | i  ,j+1,k
+            # 111 | i+1,j+1,k+1 | i  ,j+1,k+1 | i+1,j  ,k+1 | i+1,j+1,k
+
+            posX = [0,0] if xEdge == 'eX0' else [1, 0] if xEdge == 'eX1' else [0,1] if xEdge == 'eX2' else [1,1]
+            posY = [0,0] if yEdge == 'eY0' else [1, 0] if yEdge == 'eY1' else [0,1] if yEdge == 'eY2' else [1,1]
+            posZ = [0,0] if zEdge == 'eZ0' else [1, 0] if zEdge == 'eZ1' else [0,1] if zEdge == 'eZ2' else [1,1]
+
+            ind1 = sub2ind(M.vnEx, np.c_[ii, jj + posX[0], kk + posX[1]])
+            ind2 = sub2ind(M.vnEy, np.c_[ii + posY[0], jj, kk + posY[1]]) + M.nEx
+            ind3 = sub2ind(M.vnEz, np.c_[ii + posZ[0], jj + posZ[1], kk]) + M.nEx + M.nEy
+
+            IND = np.r_[ind1, ind2, ind3].flatten()
+
+            PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr()
+
+            if M._meshType == 'LRM':
+                I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]),
+                                           getSubArray(eT2[0], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[1], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[2], [i + posY[0], j, k + posY[1]]),
+                                           getSubArray(eT3[0], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[1], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[2], [i + posZ[0], j + posZ[1], k]))
+                PXXX = I3x3 * PXXX
+
+            return PXXX
+        return Pxxx
 
 if __name__ == '__main__':
     from TensorMesh import TensorMesh

From 6687a8324f5907741cfb28f1d5b483ca5d2bab92 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Mon, 16 Jun 2014 14:24:54 -0600
Subject: [PATCH 04/14] bug fixes in TreeMesh

---
 SimPEG/Mesh/InnerProducts.py |  7 +++--
 SimPEG/Mesh/TreeMesh.py      | 52 +++++++++++++++++++++---------------
 2 files changed, 33 insertions(+), 26 deletions(-)

diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py
index 6ffb20b9..79f617a5 100644
--- a/SimPEG/Mesh/InnerProducts.py
+++ b/SimPEG/Mesh/InnerProducts.py
@@ -19,7 +19,7 @@ class InnerProducts(object):
             :rtype: scipy.csr_matrix
             :return: M, the inner product matrix (nF, nF)
         """
-        return self._getInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat, doFast=True)
+        return self._getInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat, doFast=doFast)
 
     def getEdgeInnerProduct(self, prop=None, invProp=False, invMat=False, doFast=True):
         """
@@ -30,7 +30,7 @@ class InnerProducts(object):
             :rtype: scipy.csr_matrix
             :return: M, the inner product matrix (nE, nE)
         """
-        return self._getInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat)
+        return self._getInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat, doFast=doFast)
 
     def _getInnerProduct(self, projType, prop=None, invProp=False, invMat=False, doFast=True):
         """
@@ -43,11 +43,10 @@ class InnerProducts(object):
             :return: M, the inner product matrix (nE, nE)
         """
         assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
-        fast = None
 
+        fast = None
         if hasattr(self, '_fastInnerProduct') and doFast:
             fast = self._fastInnerProduct(projType, prop=prop, invProp=invProp, invMat=invMat)
-
         if fast is not None:
             return fast
 
diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py
index dca39595..d2d18972 100644
--- a/SimPEG/Mesh/TreeMesh.py
+++ b/SimPEG/Mesh/TreeMesh.py
@@ -1022,31 +1022,39 @@ class TreeMesh(InnerProducts, BaseMesh):
             V += [1./lenj]*lenj
         return sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nE))
 
-    def _getFacePxx(self, xFace, yFace):
-        self.number()
-        xP = self._getFaceP(xFace, yFace, None)
-        yP = self._getFaceP(yFace, xFace, None)
-        return sp.vstack((xP, yP))
+    def _getFacePxx(self):
+        def Pxx(xFace, yFace):
+            self.number()
+            xP = self._getFaceP(xFace, yFace, None)
+            yP = self._getFaceP(yFace, xFace, None)
+            return sp.vstack((xP, yP))
+        return Pxx
 
-    def _getEdgePxx(self, xEdge, yEdge):
-        self.number()
-        xP = self._getEdgeP(xEdge, yEdge, None)
-        yP = self._getEdgeP(yEdge, xEdge, None)
-        return sp.vstack((xP, yP))
+    def _getEdgePxx(self):
+        def Pxx(xEdge, yEdge):
+            self.number()
+            xP = self._getEdgeP(xEdge, yEdge, None)
+            yP = self._getEdgeP(yEdge, xEdge, None)
+            return sp.vstack((xP, yP))
+        return Pxx
 
-    def _getFacePxxx(self, xFace, yFace, zFace):
-        self.number()
-        xP = self._getFaceP(xFace, yFace, zFace)
-        yP = self._getFaceP(yFace, xFace, zFace)
-        zP = self._getFaceP(zFace, xFace, yFace)
-        return sp.vstack((xP, yP, zP))
+    def _getFacePxxx(self):
+        def Pxxx(xFace, yFace, zFace):
+            self.number()
+            xP = self._getFaceP(xFace, yFace, zFace)
+            yP = self._getFaceP(yFace, xFace, zFace)
+            zP = self._getFaceP(zFace, xFace, yFace)
+            return sp.vstack((xP, yP, zP))
+        return Pxxx
 
-    def _getEdgePxxx(self, xEdge, yEdge, zEdge):
-        self.number()
-        xP = self._getEdgeP(xEdge, yEdge, zEdge)
-        yP = self._getEdgeP(yEdge, xEdge, zEdge)
-        zP = self._getEdgeP(zEdge, xEdge, yEdge)
-        return sp.vstack((xP, yP, zP))
+    def _getEdgePxxx(self):
+        def Pxxx(xEdge, yEdge, zEdge):
+            self.number()
+            xP = self._getEdgeP(xEdge, yEdge, zEdge)
+            yP = self._getEdgeP(yEdge, xEdge, zEdge)
+            zP = self._getEdgeP(zEdge, xEdge, yEdge)
+            return sp.vstack((xP, yP, zP))
+        return Pxxx
 
     def plotGrid(self, ax=None, text=False, centers=False, faces=False, edges=False, lines=True, nodes=False, showIt=False):
         self.number()

From d2bde11a2aad32ce192033026f6d00c11932524e Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Mon, 16 Jun 2014 17:09:06 -0600
Subject: [PATCH 05/14] Fixes #68

---
 SimPEG/Mesh/InnerProducts.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py
index 79f617a5..4ea9a3d2 100644
--- a/SimPEG/Mesh/InnerProducts.py
+++ b/SimPEG/Mesh/InnerProducts.py
@@ -54,10 +54,9 @@ class InnerProducts(object):
             prop = invPropertyTensor(self, prop)
 
         tensorType = TensorType(self, prop)
+
         Mu = makePropertyTensor(self, prop)
-
         Ps = self._getInnerProductProjectionMatrices(projType, tensorType)
-
         A = np.sum([P.T * Mu * P for P in Ps])
 
         if invMat and tensorType < 3:

From 7d34f596ff03078356570fe100dd5ef264d45da5 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Tue, 17 Jun 2014 09:47:03 -0600
Subject: [PATCH 06/14] More tests on more meshes, IPDeriv returns a function.

---
 SimPEG/Mesh/InnerProducts.py            |  12 +-
 SimPEG/Mesh/TensorMesh.py               |   6 +-
 SimPEG/Tests/test_massMatricesDerivs.py | 222 +++++++++++++++++++-----
 SimPEG/Utils/matutils.py                |   7 -
 4 files changed, 192 insertions(+), 55 deletions(-)

diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py
index 4ea9a3d2..f133410b 100644
--- a/SimPEG/Mesh/InnerProducts.py
+++ b/SimPEG/Mesh/InnerProducts.py
@@ -114,8 +114,14 @@ class InnerProducts(object):
             :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
             :param list P: list of projection matrices
             :param bool doFast: do a faster implementation if available.
+            :rtype: function
+            :return: dMdmu(u), the derivative of the inner product matrix (u)
+
+            Given u, dMdmu returns (nF, nC*nA)
+
+            :param np.ndarray u: vector that multiplies dMdmu
             :rtype: scipy.csr_matrix
-            :return: dMdm, the derivative of the inner product matrix (nF, nC*nA)
+            :return: dMdmu, the derivative of the inner product matrix for a certain u
         """
         assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.'
         fast = None
@@ -131,7 +137,7 @@ class InnerProducts(object):
 
         def innerProductDeriv(v):
             return self._getInnerProductDeriv(tensorType, P, self.nF, v)
-        return DerivOperator(innerProductDeriv)
+        return innerProductDeriv
 
 
     def getEdgeInnerProductDeriv(self, tensorType, P=None, doFast=True):
@@ -155,7 +161,7 @@ class InnerProducts(object):
             P = self._getInnerProductProjectionMatrices('E', tensorType=tensorType)
         def innerProductDeriv(v):
             return self._getInnerProductDeriv(tensorType, P, self.nE, v)
-        return DerivOperator(innerProductDeriv)
+        return innerProductDeriv
 
     def _getInnerProductDeriv(self, tensorType, P, n, v):
         """
diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py
index faf45cf2..96de4be5 100644
--- a/SimPEG/Mesh/TensorMesh.py
+++ b/SimPEG/Mesh/TensorMesh.py
@@ -319,21 +319,21 @@ class BaseTensorMesh(BaseRectangularMesh):
             #     return self.dim * Av.T * V * ones
             def scalarInnerProductDeriv(v):
                 return Utils.sdiag(v) * self.dim * Av.T * V * ones
-            return Utils.DerivOperator(scalarInnerProductDeriv)
+            return scalarInnerProductDeriv
 
         if tensorType == 1:
             Av = getattr(self, 'ave'+projType+'2CC')
             V = Utils.sdiag(self.vol)
             def isotropicInnerProductDeriv(v):
                 return Utils.sdiag(v) * self.dim * Av.T * V
-            return Utils.DerivOperator(isotropicInnerProductDeriv)
+            return isotropicInnerProductDeriv
 
         if tensorType == 2: # anisotropic
             Av = getattr(self, 'ave'+projType+'2CCV')
             V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
             def anisotropicInnerProductDeriv(v):
                 return Utils.sdiag(v) * Av.T * V
-            return Utils.DerivOperator(anisotropicInnerProductDeriv)
+            return anisotropicInnerProductDeriv
 
 
 
diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py
index 72e36f3e..dd208c68 100644
--- a/SimPEG/Tests/test_massMatricesDerivs.py
+++ b/SimPEG/Tests/test_massMatricesDerivs.py
@@ -6,102 +6,240 @@ from TestUtils import checkDerivative
 
 class TestInnerProductsDerivs(unittest.TestCase):
 
-    def doTestFace(self, h, rep, fast):
-        mesh = Mesh.TensorMesh(h)
+    def doTestFace(self, h, rep, fast, meshType):
+        if meshType == 'LRM':
+            hRect = Utils.exampleLrmGrid(h,'rotate')
+            mesh = Mesh.LogicallyRectMesh(hRect)
+        elif meshType == 'Tree':
+            mesh = Mesh.TreeMesh(h)
+        elif meshType == 'Tensor':
+            mesh = Mesh.TensorMesh(h)
         v = np.random.rand(mesh.nF)
         sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
         Md = mesh.getFaceInnerProductDeriv(Utils.TensorType(mesh, sig), doFast=fast)
         def fun(sig):
             M = mesh.getFaceInnerProduct(sig)
-            return M*v, Md*v
+            return M*v, Md(v)
+        print meshType, 'Face', h, rep, fast
         return checkDerivative(fun, sig, num=5, plotIt=False)
 
-    def doTestEdge(self, h, rep, fast):
-        mesh = Mesh.TensorMesh(h)
+    def doTestEdge(self, h, rep, fast, meshType):
+        if meshType == 'LRM':
+            hRect = Utils.exampleLrmGrid(h,'rotate')
+            mesh = Mesh.LogicallyRectMesh(hRect)
+        elif meshType == 'Tree':
+            mesh = Mesh.TreeMesh(h)
+        elif meshType == 'Tensor':
+            mesh = Mesh.TensorMesh(h)
         v = np.random.rand(mesh.nE)
         sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
         Md = mesh.getEdgeInnerProductDeriv(Utils.TensorType(mesh, sig), doFast=fast)
         def fun(sig):
             M = mesh.getEdgeInnerProduct(sig)
-            return M*v, Md*v
+            return M*v, Md(v)
+        print meshType, 'Edge', h, rep, fast
         return checkDerivative(fun, sig, num=5, plotIt=False)
 
     def test_FaceIP_1D_float(self):
-        self.assertTrue(self.doTestFace([10],0, False))
+        self.assertTrue(self.doTestFace([10],0, False, 'Tensor'))
     def test_FaceIP_2D_float(self):
-        self.assertTrue(self.doTestFace([10, 4],0, False))
+        self.assertTrue(self.doTestFace([10, 4],0, False, 'Tensor'))
     def test_FaceIP_3D_float(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],0, False))
+        self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Tensor'))
     def test_FaceIP_1D_isotropic(self):
-        self.assertTrue(self.doTestFace([10],1, False))
+        self.assertTrue(self.doTestFace([10],1, False, 'Tensor'))
     def test_FaceIP_2D_isotropic(self):
-        self.assertTrue(self.doTestFace([10, 4],1, False))
+        self.assertTrue(self.doTestFace([10, 4],1, False, 'Tensor'))
     def test_FaceIP_3D_isotropic(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],1, False))
+        self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Tensor'))
     def test_FaceIP_2D_anisotropic(self):
-        self.assertTrue(self.doTestFace([10, 4],2, False))
+        self.assertTrue(self.doTestFace([10, 4],2, False, 'Tensor'))
     def test_FaceIP_3D_anisotropic(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],3, False))
+        self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Tensor'))
     def test_FaceIP_2D_tensor(self):
-        self.assertTrue(self.doTestFace([10, 4],3, False))
+        self.assertTrue(self.doTestFace([10, 4],3, False, 'Tensor'))
     def test_FaceIP_3D_tensor(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],6, False))
+        self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Tensor'))
 
     def test_FaceIP_1D_float_fast(self):
-        self.assertTrue(self.doTestFace([10],0, True))
+        self.assertTrue(self.doTestFace([10],0, True, 'Tensor'))
     def test_FaceIP_2D_float_fast(self):
-        self.assertTrue(self.doTestFace([10, 4],0, True))
+        self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor'))
     def test_FaceIP_3D_float_fast(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],0, True))
+        self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tensor'))
     def test_FaceIP_1D_isotropic_fast(self):
-        self.assertTrue(self.doTestFace([10],1, True))
+        self.assertTrue(self.doTestFace([10],1, True, 'Tensor'))
     def test_FaceIP_2D_isotropic_fast(self):
-        self.assertTrue(self.doTestFace([10, 4],1, True))
+        self.assertTrue(self.doTestFace([10, 4],1, True, 'Tensor'))
     def test_FaceIP_3D_isotropic_fast(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],1, True))
+        self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tensor'))
     def test_FaceIP_2D_anisotropic_fast(self):
-        self.assertTrue(self.doTestFace([10, 4],2, True))
+        self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor'))
     def test_FaceIP_3D_anisotropic_fast(self):
-        self.assertTrue(self.doTestFace([10, 4, 5],3, True))
+        self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tensor'))
 
     def test_EdgeIP_1D_float(self):
-        self.assertTrue(self.doTestEdge([10],0, False))
+        self.assertTrue(self.doTestEdge([10],0, False, 'Tensor'))
     def test_EdgeIP_2D_float(self):
-        self.assertTrue(self.doTestEdge([10, 4],0, False))
+        self.assertTrue(self.doTestEdge([10, 4],0, False, 'Tensor'))
     def test_EdgeIP_3D_float(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],0, False))
+        self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'Tensor'))
     def test_EdgeIP_1D_isotropic(self):
-        self.assertTrue(self.doTestEdge([10],1, False))
+        self.assertTrue(self.doTestEdge([10],1, False, 'Tensor'))
     def test_EdgeIP_2D_isotropic(self):
-        self.assertTrue(self.doTestEdge([10, 4],1, False))
+        self.assertTrue(self.doTestEdge([10, 4],1, False, 'Tensor'))
     def test_EdgeIP_3D_isotropic(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],1, False))
+        self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'Tensor'))
     def test_EdgeIP_2D_anisotropic(self):
-        self.assertTrue(self.doTestEdge([10, 4],2, False))
+        self.assertTrue(self.doTestEdge([10, 4],2, False, 'Tensor'))
     def test_EdgeIP_3D_anisotropic(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],3, False))
+        self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'Tensor'))
     def test_EdgeIP_2D_tensor(self):
-        self.assertTrue(self.doTestEdge([10, 4],3, False))
+        self.assertTrue(self.doTestEdge([10, 4],3, False, 'Tensor'))
     def test_EdgeIP_3D_tensor(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],6, False))
+        self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Tensor'))
 
     def test_EdgeIP_1D_float_fast(self):
-        self.assertTrue(self.doTestEdge([10],0, True))
+        self.assertTrue(self.doTestEdge([10],0, True, 'Tensor'))
     def test_EdgeIP_2D_float_fast(self):
-        self.assertTrue(self.doTestEdge([10, 4],0, True))
+        self.assertTrue(self.doTestEdge([10, 4],0, True, 'Tensor'))
     def test_EdgeIP_3D_float_fast(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],0, True))
+        self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'Tensor'))
     def test_EdgeIP_1D_isotropic_fast(self):
-        self.assertTrue(self.doTestEdge([10],1, True))
+        self.assertTrue(self.doTestEdge([10],1, True, 'Tensor'))
     def test_EdgeIP_2D_isotropic_fast(self):
-        self.assertTrue(self.doTestEdge([10, 4],1, True))
+        self.assertTrue(self.doTestEdge([10, 4],1, True, 'Tensor'))
     def test_EdgeIP_3D_isotropic_fast(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],1, True))
+        self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'Tensor'))
     def test_EdgeIP_2D_anisotropic_fast(self):
-        self.assertTrue(self.doTestEdge([10, 4],2, True))
+        self.assertTrue(self.doTestEdge([10, 4],2, True, 'Tensor'))
     def test_EdgeIP_3D_anisotropic_fast(self):
-        self.assertTrue(self.doTestEdge([10, 4, 5],3, True))
+        self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Tensor'))
 
 
+
+    def test_FaceIP_2D_float_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4],0, False, 'LRM'))
+    def test_FaceIP_3D_float_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'LRM'))
+    def test_FaceIP_2D_isotropic_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4],1, False, 'LRM'))
+    def test_FaceIP_3D_isotropic_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'LRM'))
+    def test_FaceIP_2D_anisotropic_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4],2, False, 'LRM'))
+    def test_FaceIP_3D_anisotropic_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'LRM'))
+    def test_FaceIP_2D_tensor_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4],3, False, 'LRM'))
+    def test_FaceIP_3D_tensor_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'LRM'))
+
+    def test_FaceIP_2D_float_fast_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4],0, True, 'LRM'))
+    def test_FaceIP_3D_float_fast_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'LRM'))
+    def test_FaceIP_2D_isotropic_fast_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4],1, True, 'LRM'))
+    def test_FaceIP_3D_isotropic_fast_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'LRM'))
+    def test_FaceIP_2D_anisotropic_fast_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4],2, True, 'LRM'))
+    def test_FaceIP_3D_anisotropic_fast_LRM(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'LRM'))
+
+    def test_EdgeIP_2D_float_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4],0, False, 'LRM'))
+    def test_EdgeIP_3D_float_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'LRM'))
+    def test_EdgeIP_2D_isotropic_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4],1, False, 'LRM'))
+    def test_EdgeIP_3D_isotropic_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'LRM'))
+    def test_EdgeIP_2D_anisotropic_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4],2, False, 'LRM'))
+    def test_EdgeIP_3D_anisotropic_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'LRM'))
+    def test_EdgeIP_2D_tensor_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4],3, False, 'LRM'))
+    def test_EdgeIP_3D_tensor_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'LRM'))
+
+    def test_EdgeIP_2D_float_fast_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4],0, True, 'LRM'))
+    def test_EdgeIP_3D_float_fast_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'LRM'))
+    def test_EdgeIP_2D_isotropic_fast_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4],1, True, 'LRM'))
+    def test_EdgeIP_3D_isotropic_fast_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'LRM'))
+    def test_EdgeIP_2D_anisotropic_fast_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4],2, True, 'LRM'))
+    def test_EdgeIP_3D_anisotropic_fast_LRM(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'LRM'))
+
+
+
+
+    def test_FaceIP_2D_float_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4],0, False, 'Tree'))
+    def test_FaceIP_3D_float_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Tree'))
+    def test_FaceIP_2D_isotropic_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4],1, False, 'Tree'))
+    def test_FaceIP_3D_isotropic_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Tree'))
+    def test_FaceIP_2D_anisotropic_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4],2, False, 'Tree'))
+    def test_FaceIP_3D_anisotropic_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Tree'))
+    def test_FaceIP_2D_tensor_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4],3, False, 'Tree'))
+    def test_FaceIP_3D_tensor_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Tree'))
+
+    def test_FaceIP_2D_float_fast_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4],0, True, 'Tree'))
+    def test_FaceIP_3D_float_fast_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tree'))
+    def test_FaceIP_2D_isotropic_fast_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4],1, True, 'Tree'))
+    def test_FaceIP_3D_isotropic_fast_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tree'))
+    def test_FaceIP_2D_anisotropic_fast_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4],2, True, 'Tree'))
+    def test_FaceIP_3D_anisotropic_fast_Tree(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tree'))
+
+    def test_EdgeIP_2D_float_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4],0, False, 'Tree'))
+    def test_EdgeIP_3D_float_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'Tree'))
+    def test_EdgeIP_2D_isotropic_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4],1, False, 'Tree'))
+    def test_EdgeIP_3D_isotropic_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'Tree'))
+    def test_EdgeIP_2D_anisotropic_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4],2, False, 'Tree'))
+    def test_EdgeIP_3D_anisotropic_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'Tree'))
+    def test_EdgeIP_2D_tensor_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4],3, False, 'Tree'))
+    def test_EdgeIP_3D_tensor_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Tree'))
+
+    def test_EdgeIP_2D_float_fast_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4],0, True, 'Tree'))
+    def test_EdgeIP_3D_float_fast_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'Tree'))
+    def test_EdgeIP_2D_isotropic_fast_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4],1, True, 'Tree'))
+    def test_EdgeIP_3D_isotropic_fast_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'Tree'))
+    def test_EdgeIP_2D_anisotropic_fast_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4],2, True, 'Tree'))
+    def test_EdgeIP_3D_anisotropic_fast_Tree(self):
+        self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Tree'))
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py
index 1a4c4a0e..691aaabb 100644
--- a/SimPEG/Utils/matutils.py
+++ b/SimPEG/Utils/matutils.py
@@ -349,10 +349,3 @@ class SimPEGLinearOperator(LinearOperator):
     @property
     def T(self):
         return self.__class__((self.shape[1],self.shape[0]),self.rmatvec,rmatvec=self.matvec,matmat=self.matmat)
-
-
-class DerivOperator(object):
-    def __init__(self, f):
-        self.f   = f
-    def __mul__(self, v):
-        return self.f(v)

From a93b25b4881e9cd3fe6962bc17491891ec969b44 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Wed, 18 Jun 2014 10:21:10 -0700
Subject: [PATCH 07/14] pre multiplication in innerproduct code.

---
 SimPEG/Mesh/TensorMesh.py | 18 ++++++++++++------
 1 file changed, 12 insertions(+), 6 deletions(-)

diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py
index 96de4be5..854a111a 100644
--- a/SimPEG/Mesh/TensorMesh.py
+++ b/SimPEG/Mesh/TensorMesh.py
@@ -304,8 +304,8 @@ class BaseTensorMesh(BaseRectangularMesh):
         """
             :param str projType: 'E' or 'F'
             :param TensorType tensorType: type of the tensor
-            :rtype: scipy.csr_matrix
-            :return: M, the inner product matrix (nF, nF)
+            :rtype: function
+            :return: dMdmu, the derivative of the inner product matrix
         """
         assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
         if tensorType == -1:
@@ -317,22 +317,28 @@ class BaseTensorMesh(BaseRectangularMesh):
             ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1))
             # if v is None:
             #     return self.dim * Av.T * V * ones
+            dim_x_AvT_x_V_x_ones = self.dim * Av.T * V * ones
             def scalarInnerProductDeriv(v):
-                return Utils.sdiag(v) * self.dim * Av.T * V * ones
+                return Utils.sdiag(v) * dim_x_AvT_x_V_x_ones
             return scalarInnerProductDeriv
 
         if tensorType == 1:
             Av = getattr(self, 'ave'+projType+'2CC')
             V = Utils.sdiag(self.vol)
-            def isotropicInnerProductDeriv(v):
-                return Utils.sdiag(v) * self.dim * Av.T * V
+            dim_x_AvT_x_V = self.dim * Av.T * V
+            def isotropicInnerProductDeriv(v=None):
+                if v is None:
+                    print 'Depreciation Warning: TensorMesh.isotropicInnerProductDeriv. You should be supplying a vector. Returning: d * Av.T * V'
+                    return dim_x_AvT_x_V
+                return Utils.sdiag(v) * dim_x_AvT_x_V
             return isotropicInnerProductDeriv
 
         if tensorType == 2: # anisotropic
             Av = getattr(self, 'ave'+projType+'2CCV')
             V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
+            AvT_x_V = Av.T * V
             def anisotropicInnerProductDeriv(v):
-                return Utils.sdiag(v) * Av.T * V
+                return Utils.sdiag(v) * AvT_x_V
             return anisotropicInnerProductDeriv
 
 

From 32fa4c95dfe61439e12e22d727471ebee64e3bd5 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Wed, 18 Jun 2014 11:08:42 -0700
Subject: [PATCH 08/14] Add tensorType to problem.

---
 SimPEG/Problem.py                 |  4 ++
 SimPEG/Tests/test_massMatrices.py | 86 ++++++++++++++++---------------
 2 files changed, 49 insertions(+), 41 deletions(-)

diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py
index dea83b1a..f60d027b 100644
--- a/SimPEG/Problem.py
+++ b/SimPEG/Problem.py
@@ -338,6 +338,10 @@ class BaseProblem(object):
             if hasattr(self, prop):
                 delattr(self, prop)
 
+    @property
+    def tensorType(self):
+        return Utils.TensorType(self.mesh, self.curModel.transform)
+
     @property
     def ispaired(self):
         """True if the problem is paired to a survey."""
diff --git a/SimPEG/Tests/test_massMatrices.py b/SimPEG/Tests/test_massMatrices.py
index b167ce9e..c0bd85a0 100644
--- a/SimPEG/Tests/test_massMatrices.py
+++ b/SimPEG/Tests/test_massMatrices.py
@@ -260,56 +260,60 @@ class TestInnerProducts1D(OrderTest):
 if __name__ == '__main__':
     unittest.main()
 
-if __name__ == '__main__' and False:
-    import sympy
+###################################################
+#### Uncomment to Reevaluate the InnerProducts ####
+###################################################
 
-    x,y,z = sympy.symbols(['x','y','z'])
-    ex = x**2+y*z
-    ey = (z**2)*x+y*z
-    ez = y**2+x*z
-    e = sympy.Matrix([ex,ey,ez])
+# if __name__ == '__main__':
+    # import sympy
 
-    sigma1 = x*y+1
-    sigma2 = x*z+2
-    sigma3 = 3+z*y
-    sigma4 = 0.1*x*y*z
-    sigma5 = 0.2*x*y
-    sigma6 = 0.1*z
+    # x,y,z = sympy.symbols(['x','y','z'])
+    # ex = x**2+y*z
+    # ey = (z**2)*x+y*z
+    # ez = y**2+x*z
+    # e = sympy.Matrix([ex,ey,ez])
 
-    S1 = sympy.Matrix([[sigma1,0,0],[0,sigma1,0],[0,0,sigma1]])
-    S2 = sympy.Matrix([[sigma1,0,0],[0,sigma2,0],[0,0,sigma3]])
-    S3 = sympy.Matrix([[sigma1,sigma4,sigma5],[sigma4,sigma2,sigma6],[sigma5,sigma6,sigma3]])
+    # sigma1 = x*y+1
+    # sigma2 = x*z+2
+    # sigma3 = 3+z*y
+    # sigma4 = 0.1*x*y*z
+    # sigma5 = 0.2*x*y
+    # sigma6 = 0.1*z
 
-    print '3D'
-    print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1)), (z,0,1))
-    print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1)), (z,0,1))
-    print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1)), (z,0,1))
+    # S1 = sympy.Matrix([[sigma1,0,0],[0,sigma1,0],[0,0,sigma1]])
+    # S2 = sympy.Matrix([[sigma1,0,0],[0,sigma2,0],[0,0,sigma3]])
+    # S3 = sympy.Matrix([[sigma1,sigma4,sigma5],[sigma4,sigma2,sigma6],[sigma5,sigma6,sigma3]])
+
+    # print '3D'
+    # print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1)), (z,0,1))
+    # print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1)), (z,0,1))
+    # print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1)), (z,0,1))
 
 
-    z = 5
-    ex = x**2+y*z
-    ey = (z**2)*x+y*z
-    e = sympy.Matrix([ex,ey])
+    # z = 5
+    # ex = x**2+y*z
+    # ey = (z**2)*x+y*z
+    # e = sympy.Matrix([ex,ey])
 
-    sigma1 = x*y+1
-    sigma2 = x*z+2
-    sigma3 = 3+z*y
+    # sigma1 = x*y+1
+    # sigma2 = x*z+2
+    # sigma3 = 3+z*y
 
-    S1 = sympy.Matrix([[sigma1,0],[0,sigma1]])
-    S2 = sympy.Matrix([[sigma1,0],[0,sigma2]])
-    S3 = sympy.Matrix([[sigma1,sigma3],[sigma3,sigma2]])
+    # S1 = sympy.Matrix([[sigma1,0],[0,sigma1]])
+    # S2 = sympy.Matrix([[sigma1,0],[0,sigma2]])
+    # S3 = sympy.Matrix([[sigma1,sigma3],[sigma3,sigma2]])
 
-    print '2D'
-    print sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1))
-    print sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1))
-    print sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1))
+    # print '2D'
+    # print sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1))
+    # print sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1))
+    # print sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1))
 
-    y = 12
-    z = 5
-    ex = x**2+y*z
-    e = ex
+    # y = 12
+    # z = 5
+    # ex = x**2+y*z
+    # e = ex
 
-    sigma1 = x*y+1
+    # sigma1 = x*y+1
 
-    print '1D'
-    print sympy.integrate(e*sigma1*e, (x,0,1))
+    # print '1D'
+    # print sympy.integrate(e*sigma1*e, (x,0,1))

From e77f43137092e0e1a89e26d44a66fd9f4a9122f4 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Wed, 2 Jul 2014 16:22:33 -0700
Subject: [PATCH 09/14] closest inds for 1d mesh

---
 SimPEG/Utils/meshutils.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py
index 9f7890ee..7eeac690 100644
--- a/SimPEG/Utils/meshutils.py
+++ b/SimPEG/Utils/meshutils.py
@@ -94,7 +94,10 @@ def closestPoints(mesh, pts, gridLoc='CC'):
     nodeInds = np.empty(pts.shape[0], dtype=int)
 
     for i, pt in enumerate(pts):
-        nodeInds[i] = ((np.tile(pt, (grid.shape[0],1)) - grid)**2).sum(axis=1).argmin()
+        if mesh.dim == 1:
+            nodeInds[i] = ((pt - grid)**2).argmin()
+        else:
+            nodeInds[i] = ((np.tile(pt, (grid.shape[0],1)) - grid)**2).sum(axis=1).argmin()
 
     return nodeInds
 
@@ -107,7 +110,7 @@ def readUBCTensorMesh(fileName):
 
         Output:
         :param SimPEG TensorMesh object
-        :return 
+        :return
     """
 
     # Interal function to read cell size lines for the UBC mesh files.
@@ -119,8 +122,8 @@ def readUBCTensorMesh(fileName):
                 re = np.array(sp[0],dtype=int)*(' ' + sp[1])
                 line = line.replace(st,re.strip())
         return np.array(line.split(),dtype=float)
-    
-    # Read the file as line strings, remove lines with comment = ! 
+
+    # Read the file as line strings, remove lines with comment = !
     msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!')
 
     # Fist line is the size of the model
@@ -134,7 +137,7 @@ def readUBCTensorMesh(fileName):
     h3 = h3temp[::-1] # Invert the indexing of the vector to start from the bottom.
     # Adjust the reference point to the bottom south west corner
     x0[2] = x0[2] - np.sum(h3)
-    # Make the mesh 
+    # Make the mesh
     from SimPEG import Mesh
     tensMsh = Mesh.TensorMesh([h1,h2,h3],x0)
     return tensMsh

From 9d048f28843dea197f48cd697cd86a489efd7253 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Wed, 2 Jul 2014 16:52:46 -0700
Subject: [PATCH 10/14] generalize the derivs for innerProducts

---
 SimPEG/Mesh/InnerProducts.py | 46 ++++++++++++++++--------------------
 1 file changed, 21 insertions(+), 25 deletions(-)

diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py
index f133410b..11f5d496 100644
--- a/SimPEG/Mesh/InnerProducts.py
+++ b/SimPEG/Mesh/InnerProducts.py
@@ -109,10 +109,9 @@ class InnerProducts(object):
         return [V*proj(*locs[node][d-1]) for node in nodes]
 
 
-    def getFaceInnerProductDeriv(self, tensorType, P=None, doFast=True):
+    def getFaceInnerProductDeriv(self, tensorType, doFast=True):
         """
             :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
-            :param list P: list of projection matrices
             :param bool doFast: do a faster implementation if available.
             :rtype: function
             :return: dMdmu(u), the derivative of the inner product matrix (u)
@@ -123,27 +122,22 @@ class InnerProducts(object):
             :rtype: scipy.csr_matrix
             :return: dMdmu, the derivative of the inner product matrix for a certain u
         """
-        assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.'
-        fast = None
-
-        if hasattr(self, '_fastInnerProductDeriv') and doFast:
-            fast = self._fastInnerProductDeriv('F', tensorType)
-
-        if fast is not None:
-            return fast
-
-        if P is None:
-            P = self._getInnerProductProjectionMatrices('F', tensorType=tensorType)
-
-        def innerProductDeriv(v):
-            return self._getInnerProductDeriv(tensorType, P, self.nF, v)
-        return innerProductDeriv
+        return self._getInnerProductDeriv(tensorType, 'F', doFast=doFast)
 
 
-    def getEdgeInnerProductDeriv(self, tensorType, P=None, doFast=True):
+    def getEdgeInnerProductDeriv(self, tensorType, doFast=True):
         """
             :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
-            :param list P: list of projection matrices
+            :param bool doFast: do a faster implementation if available.
+            :rtype: scipy.csr_matrix
+            :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
+        """
+        return self._getInnerProductDeriv(tensorType, 'E', doFast=doFast)
+
+    def _getInnerProductDeriv(self, tensorType, projType, doFast=True):
+        """
+            :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
+            :param str projType: 'F' for faces 'E' for edges
             :param bool doFast: do a faster implementation if available.
             :rtype: scipy.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
@@ -152,26 +146,28 @@ class InnerProducts(object):
         fast = None
 
         if hasattr(self, '_fastInnerProductDeriv') and doFast:
-            fast = self._fastInnerProductDeriv('E', tensorType)
+            fast = self._fastInnerProductDeriv(projType, tensorType)
 
         if fast is not None:
             return fast
 
-        if P is None:
-            P = self._getInnerProductProjectionMatrices('E', tensorType=tensorType)
+        P = self._getInnerProductProjectionMatrices(projType, tensorType=tensorType)
         def innerProductDeriv(v):
-            return self._getInnerProductDeriv(tensorType, P, self.nE, v)
+            return self._getInnerProductDerivFunction(tensorType, P, projType, v)
         return innerProductDeriv
 
-    def _getInnerProductDeriv(self, tensorType, P, n, v):
+    def _getInnerProductDerivFunction(self, tensorType, P, projType, v):
         """
             :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param numpy.array v: vector to multiply (required in the general implementation)
             :param list P: list of projection matrices
-            :param int n: nF or nE
+            :param str projType: 'F' for faces 'E' for edges
             :rtype: scipy.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (n, nC*nA)
         """
+        assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
+        n = getattr(self,'n'+projType)
+
         if tensorType == -1:
             return None
 

From 8c65ddacae8edc8990f7cdebd286bb60f220e0cd Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Thu, 3 Jul 2014 10:08:39 -0700
Subject: [PATCH 11/14] Change innerProduct to take sigma not tensorType.

---
 SimPEG/Mesh/InnerProducts.py            | 18 +++++++++---------
 SimPEG/Tests/test_massMatricesDerivs.py |  4 ++--
 2 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py
index 11f5d496..f8bcee46 100644
--- a/SimPEG/Mesh/InnerProducts.py
+++ b/SimPEG/Mesh/InnerProducts.py
@@ -109,9 +109,9 @@ class InnerProducts(object):
         return [V*proj(*locs[node][d-1]) for node in nodes]
 
 
-    def getFaceInnerProductDeriv(self, tensorType, doFast=True):
+    def getFaceInnerProductDeriv(self, prop, doFast=True):
         """
-            :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
+            :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param bool doFast: do a faster implementation if available.
             :rtype: function
             :return: dMdmu(u), the derivative of the inner product matrix (u)
@@ -122,27 +122,27 @@ class InnerProducts(object):
             :rtype: scipy.csr_matrix
             :return: dMdmu, the derivative of the inner product matrix for a certain u
         """
-        return self._getInnerProductDeriv(tensorType, 'F', doFast=doFast)
+        return self._getInnerProductDeriv(prop, 'F', doFast=doFast)
 
 
-    def getEdgeInnerProductDeriv(self, tensorType, doFast=True):
+    def getEdgeInnerProductDeriv(self, prop, doFast=True):
         """
-            :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
+            :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param bool doFast: do a faster implementation if available.
             :rtype: scipy.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
         """
-        return self._getInnerProductDeriv(tensorType, 'E', doFast=doFast)
+        return self._getInnerProductDeriv(prop, 'E', doFast=doFast)
 
-    def _getInnerProductDeriv(self, tensorType, projType, doFast=True):
+    def _getInnerProductDeriv(self, prop, projType, doFast=True):
         """
-            :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
+            :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param str projType: 'F' for faces 'E' for edges
             :param bool doFast: do a faster implementation if available.
             :rtype: scipy.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
         """
-        assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.'
+        tensorType = TensorType(self, prop)
         fast = None
 
         if hasattr(self, '_fastInnerProductDeriv') and doFast:
diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py
index dd208c68..21ac8c41 100644
--- a/SimPEG/Tests/test_massMatricesDerivs.py
+++ b/SimPEG/Tests/test_massMatricesDerivs.py
@@ -16,7 +16,7 @@ class TestInnerProductsDerivs(unittest.TestCase):
             mesh = Mesh.TensorMesh(h)
         v = np.random.rand(mesh.nF)
         sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
-        Md = mesh.getFaceInnerProductDeriv(Utils.TensorType(mesh, sig), doFast=fast)
+        Md = mesh.getFaceInnerProductDeriv(sig, doFast=fast)
         def fun(sig):
             M = mesh.getFaceInnerProduct(sig)
             return M*v, Md(v)
@@ -33,7 +33,7 @@ class TestInnerProductsDerivs(unittest.TestCase):
             mesh = Mesh.TensorMesh(h)
         v = np.random.rand(mesh.nE)
         sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
-        Md = mesh.getEdgeInnerProductDeriv(Utils.TensorType(mesh, sig), doFast=fast)
+        Md = mesh.getEdgeInnerProductDeriv(sig, doFast=fast)
         def fun(sig):
             M = mesh.getEdgeInnerProduct(sig)
             return M*v, Md(v)

From bb0b9c1de4da6de5aee25edc6c00890f94ca68a3 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Thu, 3 Jul 2014 10:31:04 -0700
Subject: [PATCH 12/14] write test for invProp and invMat (failing)

---
 SimPEG/Mesh/InnerProducts.py            | 21 ++++++++++++-----
 SimPEG/Tests/test_massMatricesDerivs.py | 31 ++++++++++++++++++++-----
 2 files changed, 40 insertions(+), 12 deletions(-)

diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py
index f8bcee46..9bc092d3 100644
--- a/SimPEG/Mesh/InnerProducts.py
+++ b/SimPEG/Mesh/InnerProducts.py
@@ -109,10 +109,12 @@ class InnerProducts(object):
         return [V*proj(*locs[node][d-1]) for node in nodes]
 
 
-    def getFaceInnerProductDeriv(self, prop, doFast=True):
+    def getFaceInnerProductDeriv(self, prop, doFast=True, invProp=False, invMat=False):
         """
             :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param bool doFast: do a faster implementation if available.
+            :param bool invProp: inverts the material property
+            :param bool invMat: inverts the matrix
             :rtype: function
             :return: dMdmu(u), the derivative of the inner product matrix (u)
 
@@ -122,23 +124,27 @@ class InnerProducts(object):
             :rtype: scipy.csr_matrix
             :return: dMdmu, the derivative of the inner product matrix for a certain u
         """
-        return self._getInnerProductDeriv(prop, 'F', doFast=doFast)
+        return self._getInnerProductDeriv(prop, 'F', doFast=doFast, invProp=invProp, invMat=invMat)
 
 
-    def getEdgeInnerProductDeriv(self, prop, doFast=True):
+    def getEdgeInnerProductDeriv(self, prop, doFast=True, invProp=False, invMat=False):
         """
             :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param bool doFast: do a faster implementation if available.
+            :param bool invProp: inverts the material property
+            :param bool invMat: inverts the matrix
             :rtype: scipy.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
         """
-        return self._getInnerProductDeriv(prop, 'E', doFast=doFast)
+        return self._getInnerProductDeriv(prop, 'E', doFast=doFast, invProp=invProp, invMat=invMat)
 
-    def _getInnerProductDeriv(self, prop, projType, doFast=True):
+    def _getInnerProductDeriv(self, prop, projType, doFast=True, invProp=False, invMat=False):
         """
             :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param str projType: 'F' for faces 'E' for edges
             :param bool doFast: do a faster implementation if available.
+            :param bool invProp: inverts the material property
+            :param bool invMat: inverts the matrix
             :rtype: scipy.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
         """
@@ -146,7 +152,10 @@ class InnerProducts(object):
         fast = None
 
         if hasattr(self, '_fastInnerProductDeriv') and doFast:
-            fast = self._fastInnerProductDeriv(projType, tensorType)
+            fast = self._fastInnerProductDeriv(projType, tensorType, invProp=invProp, invMat=invMat)
+
+        if invProp or invMat:
+            raise NotImplementedError('inverting the property or the matrix is not yet implemented for this mesh/tensorType. You should write it!')
 
         if fast is not None:
             return fast
diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py
index 21ac8c41..5eac4ca8 100644
--- a/SimPEG/Tests/test_massMatricesDerivs.py
+++ b/SimPEG/Tests/test_massMatricesDerivs.py
@@ -6,7 +6,7 @@ from TestUtils import checkDerivative
 
 class TestInnerProductsDerivs(unittest.TestCase):
 
-    def doTestFace(self, h, rep, fast, meshType):
+    def doTestFace(self, h, rep, fast, meshType, invProp=False, invMat=False):
         if meshType == 'LRM':
             hRect = Utils.exampleLrmGrid(h,'rotate')
             mesh = Mesh.LogicallyRectMesh(hRect)
@@ -16,14 +16,14 @@ class TestInnerProductsDerivs(unittest.TestCase):
             mesh = Mesh.TensorMesh(h)
         v = np.random.rand(mesh.nF)
         sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
-        Md = mesh.getFaceInnerProductDeriv(sig, doFast=fast)
         def fun(sig):
-            M = mesh.getFaceInnerProduct(sig)
+            M  = mesh.getFaceInnerProduct(sig, invProp=invProp, invMat=invMat)
+            Md = mesh.getFaceInnerProductDeriv(sig, invProp=invProp, invMat=invMat, doFast=fast)
             return M*v, Md(v)
         print meshType, 'Face', h, rep, fast
         return checkDerivative(fun, sig, num=5, plotIt=False)
 
-    def doTestEdge(self, h, rep, fast, meshType):
+    def doTestEdge(self, h, rep, fast, meshType, invProp=False, invMat=False):
         if meshType == 'LRM':
             hRect = Utils.exampleLrmGrid(h,'rotate')
             mesh = Mesh.LogicallyRectMesh(hRect)
@@ -33,9 +33,9 @@ class TestInnerProductsDerivs(unittest.TestCase):
             mesh = Mesh.TensorMesh(h)
         v = np.random.rand(mesh.nE)
         sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
-        Md = mesh.getEdgeInnerProductDeriv(sig, doFast=fast)
         def fun(sig):
-            M = mesh.getEdgeInnerProduct(sig)
+            M  = mesh.getEdgeInnerProduct(sig, invProp=invProp, invMat=invMat)
+            Md = mesh.getEdgeInnerProductDeriv(sig, invProp=invProp, invMat=invMat, doFast=fast)
             return M*v, Md(v)
         print meshType, 'Edge', h, rep, fast
         return checkDerivative(fun, sig, num=5, plotIt=False)
@@ -118,6 +118,25 @@ class TestInnerProductsDerivs(unittest.TestCase):
 
 
 
+    # def test_FaceIP_1D_float_fast_harmonic(self):
+    #     self.assertTrue(self.doTestFace([10],0, True, 'Tensor', invProp=True, invMat=True))
+    # def test_FaceIP_2D_float_fast_harmonic(self):
+    #     self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor', invProp=True, invMat=True))
+    # def test_FaceIP_3D_float_fast_harmonic(self):
+    #     self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tensor', invProp=True, invMat=True))
+    def test_FaceIP_1D_isotropic_fast_harmonic(self):
+        self.assertTrue(self.doTestFace([10],1, True, 'Tensor', invProp=True, invMat=True))
+    def test_FaceIP_2D_isotropic_fast_harmonic(self):
+        self.assertTrue(self.doTestFace([10, 4],1, True, 'Tensor', invProp=True, invMat=True))
+    def test_FaceIP_3D_isotropic_fast_harmonic(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tensor', invProp=True, invMat=True))
+    # def test_FaceIP_2D_anisotropic_fast_harmonic(self):
+    #     self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor', invProp=True, invMat=True))
+    # def test_FaceIP_3D_anisotropic_fast_harmonic(self):
+    #     self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tensor', invProp=True, invMat=True))
+
+
+
     def test_FaceIP_2D_float_LRM(self):
         self.assertTrue(self.doTestFace([10, 4],0, False, 'LRM'))
     def test_FaceIP_3D_float_LRM(self):

From 42278d35fe3e81f1848b8a7c56f97d05da263053 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Thu, 3 Jul 2014 11:00:25 -0700
Subject: [PATCH 13/14] implement and test harmonic averaged derivatives

---
 SimPEG/Mesh/InnerProducts.py            | 10 ++---
 SimPEG/Mesh/TensorMesh.py               | 51 +++++++++++++++----------
 SimPEG/Tests/test_massMatricesDerivs.py | 24 ++++++------
 3 files changed, 47 insertions(+), 38 deletions(-)

diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py
index 9bc092d3..f0e45665 100644
--- a/SimPEG/Mesh/InnerProducts.py
+++ b/SimPEG/Mesh/InnerProducts.py
@@ -148,18 +148,16 @@ class InnerProducts(object):
             :rtype: scipy.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
         """
-        tensorType = TensorType(self, prop)
         fast = None
-
         if hasattr(self, '_fastInnerProductDeriv') and doFast:
-            fast = self._fastInnerProductDeriv(projType, tensorType, invProp=invProp, invMat=invMat)
+            fast = self._fastInnerProductDeriv(projType, prop, invProp=invProp, invMat=invMat)
+        if fast is not None:
+            return fast
 
         if invProp or invMat:
             raise NotImplementedError('inverting the property or the matrix is not yet implemented for this mesh/tensorType. You should write it!')
 
-        if fast is not None:
-            return fast
-
+        tensorType = TensorType(self, prop)
         P = self._getInnerProductProjectionMatrices(projType, tensorType=tensorType)
         def innerProductDeriv(v):
             return self._getInnerProductDerivFunction(tensorType, P, projType, v)
diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py
index 854a111a..d9a14aaa 100644
--- a/SimPEG/Mesh/TensorMesh.py
+++ b/SimPEG/Mesh/TensorMesh.py
@@ -300,46 +300,57 @@ class BaseTensorMesh(BaseRectangularMesh):
         else:
             return M
 
-    def _fastInnerProductDeriv(self, projType, tensorType):
+    def _fastInnerProductDeriv(self, projType, prop, invProp=False, invMat=False):
         """
             :param str projType: 'E' or 'F'
             :param TensorType tensorType: type of the tensor
+            :param bool invProp: inverts the material property
+            :param bool invMat: inverts the matrix
             :rtype: function
             :return: dMdmu, the derivative of the inner product matrix
         """
         assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
-        if tensorType == -1:
-            return None
+        tensorType = Utils.TensorType(self, prop)
+
+        dMdprop = None
+
+        if invMat:
+            MI = self._fastInnerProduct(projType, prop, invProp=invProp, invMat=invMat)
 
         if tensorType == 0:
             Av = getattr(self, 'ave'+projType+'2CC')
             V = Utils.sdiag(self.vol)
             ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1))
-            # if v is None:
-            #     return self.dim * Av.T * V * ones
-            dim_x_AvT_x_V_x_ones = self.dim * Av.T * V * ones
-            def scalarInnerProductDeriv(v):
-                return Utils.sdiag(v) * dim_x_AvT_x_V_x_ones
-            return scalarInnerProductDeriv
+            if not invMat and not invProp:
+                dMdprop = self.dim * Av.T * V * ones
+            elif invMat and invProp:
+                dMdprop =  self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * ones * Utils.sdiag(1./prop**2)
 
         if tensorType == 1:
             Av = getattr(self, 'ave'+projType+'2CC')
             V = Utils.sdiag(self.vol)
-            dim_x_AvT_x_V = self.dim * Av.T * V
-            def isotropicInnerProductDeriv(v=None):
-                if v is None:
-                    print 'Depreciation Warning: TensorMesh.isotropicInnerProductDeriv. You should be supplying a vector. Returning: d * Av.T * V'
-                    return dim_x_AvT_x_V
-                return Utils.sdiag(v) * dim_x_AvT_x_V
-            return isotropicInnerProductDeriv
+            if not invMat and not invProp:
+                dMdprop = self.dim * Av.T * V
+            elif invMat and invProp:
+                dMdprop =  self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2)
 
         if tensorType == 2: # anisotropic
             Av = getattr(self, 'ave'+projType+'2CCV')
             V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
-            AvT_x_V = Av.T * V
-            def anisotropicInnerProductDeriv(v):
-                return Utils.sdiag(v) * AvT_x_V
-            return anisotropicInnerProductDeriv
+            if not invMat and not invProp:
+                dMdprop = Av.T * V
+            elif invMat and invProp:
+                dMdprop =  Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2)
+
+        if dMdprop is not None:
+            def innerProductDeriv(v=None):
+                if v is None:
+                    print 'Depreciation Warning: TensorMesh.innerProductDeriv. You should be supplying a vector. Use: sdiag(u)*dMdprop'
+                    return dMdprop
+                return Utils.sdiag(v) * dMdprop
+            return innerProductDeriv
+        else:
+            return None
 
 
 
diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py
index 5eac4ca8..4b7d63cd 100644
--- a/SimPEG/Tests/test_massMatricesDerivs.py
+++ b/SimPEG/Tests/test_massMatricesDerivs.py
@@ -20,7 +20,7 @@ class TestInnerProductsDerivs(unittest.TestCase):
             M  = mesh.getFaceInnerProduct(sig, invProp=invProp, invMat=invMat)
             Md = mesh.getFaceInnerProductDeriv(sig, invProp=invProp, invMat=invMat, doFast=fast)
             return M*v, Md(v)
-        print meshType, 'Face', h, rep, fast
+        print meshType, 'Face', h, rep, fast, ('harmonic' if invProp and invMat else 'standard')
         return checkDerivative(fun, sig, num=5, plotIt=False)
 
     def doTestEdge(self, h, rep, fast, meshType, invProp=False, invMat=False):
@@ -37,7 +37,7 @@ class TestInnerProductsDerivs(unittest.TestCase):
             M  = mesh.getEdgeInnerProduct(sig, invProp=invProp, invMat=invMat)
             Md = mesh.getEdgeInnerProductDeriv(sig, invProp=invProp, invMat=invMat, doFast=fast)
             return M*v, Md(v)
-        print meshType, 'Edge', h, rep, fast
+        print meshType, 'Edge', h, rep, fast, ('harmonic' if invProp and invMat else 'standard')
         return checkDerivative(fun, sig, num=5, plotIt=False)
 
     def test_FaceIP_1D_float(self):
@@ -118,22 +118,22 @@ class TestInnerProductsDerivs(unittest.TestCase):
 
 
 
-    # def test_FaceIP_1D_float_fast_harmonic(self):
-    #     self.assertTrue(self.doTestFace([10],0, True, 'Tensor', invProp=True, invMat=True))
-    # def test_FaceIP_2D_float_fast_harmonic(self):
-    #     self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor', invProp=True, invMat=True))
-    # def test_FaceIP_3D_float_fast_harmonic(self):
-    #     self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tensor', invProp=True, invMat=True))
+    def test_FaceIP_1D_float_fast_harmonic(self):
+        self.assertTrue(self.doTestFace([10],0, True, 'Tensor', invProp=True, invMat=True))
+    def test_FaceIP_2D_float_fast_harmonic(self):
+        self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor', invProp=True, invMat=True))
+    def test_FaceIP_3D_float_fast_harmonic(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tensor', invProp=True, invMat=True))
     def test_FaceIP_1D_isotropic_fast_harmonic(self):
         self.assertTrue(self.doTestFace([10],1, True, 'Tensor', invProp=True, invMat=True))
     def test_FaceIP_2D_isotropic_fast_harmonic(self):
         self.assertTrue(self.doTestFace([10, 4],1, True, 'Tensor', invProp=True, invMat=True))
     def test_FaceIP_3D_isotropic_fast_harmonic(self):
         self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tensor', invProp=True, invMat=True))
-    # def test_FaceIP_2D_anisotropic_fast_harmonic(self):
-    #     self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor', invProp=True, invMat=True))
-    # def test_FaceIP_3D_anisotropic_fast_harmonic(self):
-    #     self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tensor', invProp=True, invMat=True))
+    def test_FaceIP_2D_anisotropic_fast_harmonic(self):
+        self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor', invProp=True, invMat=True))
+    def test_FaceIP_3D_anisotropic_fast_harmonic(self):
+        self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tensor', invProp=True, invMat=True))
 
 
 

From 0de54dc2ae5152e098c643e2f1bb350255528ab3 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Thu, 3 Jul 2014 11:14:29 -0700
Subject: [PATCH 14/14] remove tensor type property from problem.

---
 SimPEG/Problem.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py
index f60d027b..dea83b1a 100644
--- a/SimPEG/Problem.py
+++ b/SimPEG/Problem.py
@@ -338,10 +338,6 @@ class BaseProblem(object):
             if hasattr(self, prop):
                 delattr(self, prop)
 
-    @property
-    def tensorType(self):
-        return Utils.TensorType(self.mesh, self.curModel.transform)
-
     @property
     def ispaired(self):
         """True if the problem is paired to a survey."""