Merge pull request #60 from simpeg/develop

Merge Develop and Get rid of develop.

SimPEG/Data.py

@@ -52,6 +52,7 @@ class BaseData(object):
             instead of recalculating the fields (which may be expensive!).
 
             .. math::
+
                 d_\\text{pred} = P(u(m))
 
             Where P is a projection of the fields onto the data space.
@@ -67,29 +68,22 @@ class BaseData(object):
 
 
             .. math::
+
                 d_\\text{pred} = \mathbf{P} u(m)
         """
         return u
 
-
     @Utils.count
-    def projectFieldsAdjoint(self, d):
+    def projectFieldsDeriv(self, u):
         """
-            This function is the adjoint of the projection.
-            **projectFieldsAdjoint** is used in the
-            calculation of the sensitivities.
+            This function projects the fields onto the data space.
+
 
             .. math::
-                u = \mathbf{P}^\\top d
 
-            :param numpy.array d: data
-            :param numpy.array u: fields (ish)
-            :rtype: fields like object
-            :return: data
+                \\frac{\partial d_\\text{pred}}{\partial u} = \mathbf{P}
         """
-        return d
-
-    #TODO: def projectFieldDeriv(self, u):  Does this need to be made??!
+        return sp.identity(u.size)
 
     @Utils.count
     def residual(self, m, u=None):

SimPEG/Mesh/BaseMesh.py

@@ -38,6 +38,371 @@ class BaseMesh(object):
         """
         return self._x0
 
+    @property
+    def dim(self):
+        """
+        The dimension of the mesh (1, 2, or 3).
+
+        :rtype: int
+        :return: dim
+        """
+        return len(self._n)
+
+    @property
+    def nC(self):
+        """
+        Total number of cells in the model.
+
+        :rtype: int
+        :return: nC
+
+        .. plot::
+            :include-source:
+
+            from SimPEG import Mesh, np
+            Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(centers=True,showIt=True)
+        """
+        return self._n.prod()
+
+    @property
+    def nN(self):
+        """
+        Total number of nodes
+
+        :rtype: int
+        :return: nN
+
+        .. plot::
+            :include-source:
+
+            from SimPEG import Mesh, np
+            Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(nodes=True,showIt=True)
+        """
+        return (self._n+1).prod()
+
+    @property
+    def nEx(self):
+        """
+        Number of x-edges
+
+        :rtype: int
+        :return: nEx
+        """
+        return (self._n + np.r_[0,1,1][:self.dim]).prod()
+
+    @property
+    def nEy(self):
+        """
+        Number of y-edges
+
+        :rtype: int
+        :return: nEy
+        """
+        if self.dim < 2: return None
+        return (self._n + np.r_[1,0,1][:self.dim]).prod()
+
+    @property
+    def nEz(self):
+        """
+        Number of z-edges
+
+        :rtype: int
+        :return: nEz
+        """
+        if self.dim < 3: return None
+        return (self._n + np.r_[1,1,0][:self.dim]).prod()
+
+    @property
+    def vnE(self):
+        """
+        Total number of edges in each direction
+
+        :rtype: numpy.array (dim, )
+        :return: [nEx, nEy, nEz]
+
+        .. plot::
+            :include-source:
+
+            from SimPEG import Mesh, np
+            Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(edges=True,showIt=True)
+        """
+        return np.array([x for x in [self.nEx, self.nEy, self.nEz] if not x is None])
+
+    @property
+    def nE(self):
+        """
+        Total number of edges.
+
+        :rtype: int
+        :return: sum([nEx, nEy, nEz])
+
+        """
+        return self.vnE.sum()
+
+    @property
+    def nFx(self):
+        """
+        Number of x-faces
+
+        :rtype: int
+        :return: nFx
+        """
+        return (self._n + np.r_[1,0,0][:self.dim]).prod()
+
+    @property
+    def nFy(self):
+        """
+        Number of y-faces
+
+        :rtype: int
+        :return: nFy
+        """
+        if self.dim < 2: return None
+        return (self._n + np.r_[0,1,0][:self.dim]).prod()
+
+    @property
+    def nFz(self):
+        """
+        Number of z-faces
+
+        :rtype: int
+        :return: nFz
+        """
+        if self.dim < 3: return None
+        return (self._n + np.r_[0,0,1][:self.dim]).prod()
+
+    @property
+    def vnF(self):
+        """
+        Total number of faces in each direction
+
+        :rtype: numpy.array (dim, )
+        :return: [nFx, nFy, nFz]
+
+        .. plot::
+            :include-source:
+
+            from SimPEG import Mesh, np
+            Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(faces=True,showIt=True)
+        """
+        return np.array([x for x in [self.nFx, self.nFy, self.nFz] if not x is None])
+
+    @property
+    def nF(self):
+        """
+        Total number of faces.
+
+        :rtype: int
+        :return: sum([nFx, nFy, nFz])
+
+        """
+        return self.vnF.sum()
+
+    @property
+    def normals(self):
+        """
+        Face Normals
+
+        :rtype: numpy.array (sum(nF), dim)
+        :return: normals
+        """
+        if self.dim == 2:
+            nX = np.c_[np.ones(self.nFx), np.zeros(self.nFx)]
+            nY = np.c_[np.zeros(self.nFy), np.ones(self.nFy)]
+            return np.r_[nX, nY]
+        elif self.dim == 3:
+            nX = np.c_[np.ones(self.nFx), np.zeros(self.nFx), np.zeros(self.nFx)]
+            nY = np.c_[np.zeros(self.nFy), np.ones(self.nFy), np.zeros(self.nFy)]
+            nZ = np.c_[np.zeros(self.nFz), np.zeros(self.nFz), np.ones(self.nFz)]
+            return np.r_[nX, nY, nZ]
+
+    @property
+    def tangents(self):
+        """
+        Edge Tangents
+
+        :rtype: numpy.array (sum(nE), dim)
+        :return: normals
+        """
+        if self.dim == 2:
+            tX = np.c_[np.ones(self.nEx), np.zeros(self.nEx)]
+            tY = np.c_[np.zeros(self.nEy), np.ones(self.nEy)]
+            return np.r_[tX, tY]
+        elif self.dim == 3:
+            tX = np.c_[np.ones(self.nEx), np.zeros(self.nEx), np.zeros(self.nEx)]
+            tY = np.c_[np.zeros(self.nEy), np.ones(self.nEy), np.zeros(self.nEy)]
+            tZ = np.c_[np.zeros(self.nEz), np.zeros(self.nEz), np.ones(self.nEz)]
+            return np.r_[tX, tY, tZ]
+
+    def projectFaceVector(self, fV):
+        """
+        Given a vector, fV, in cartesian coordinates, this will project it onto the mesh using the normals
+
+        :param numpy.array fV: face vector with shape (nF, dim)
+        :rtype: numpy.array with shape (nF, )
+        :return: projected face vector
+        """
+        assert type(fV) == np.ndarray, 'fV must be an ndarray'
+        assert len(fV.shape) == 2 and fV.shape[0] == self.nF and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)'
+        return np.sum(fV*self.normals, 1)
+
+    def projectEdgeVector(self, eV):
+        """
+        Given a vector, eV, in cartesian coordinates, this will project it onto the mesh using the tangents
+
+        :param numpy.array eV: edge vector with shape (nE, dim)
+        :rtype: numpy.array with shape (nE, )
+        :return: projected edge vector
+        """
+        assert type(eV) == np.ndarray, 'eV must be an ndarray'
+        assert len(eV.shape) == 2 and eV.shape[0] == self.nE and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)'
+        return np.sum(eV*self.tangents, 1)
+
+
+class BaseRectangularMesh(BaseMesh):
+    """BaseRectangularMesh"""
+    def __init__(self, n, x0=None):
+        BaseMesh.__init__(self, n, x0)
+
+    @property
+    def nCx(self):
+        """
+        Number of cells in the x direction
+
+        :rtype: int
+        :return: nCx
+        """
+        return self._n[0]
+
+    @property
+    def nCy(self):
+        """
+        Number of cells in the y direction
+
+        :rtype: int
+        :return: nCy or None if dim < 2
+        """
+        return None if self.dim < 2 else self._n[1]
+
+    @property
+    def nCz(self):
+        """Number of cells in the z direction
+
+        :rtype: int
+        :return: nCz or None if dim < 3
+        """
+        return None if self.dim < 3 else self._n[2]
+
+    @property
+    def vnC(self):
+        """
+        Total number of cells in each direction
+
+        :rtype: numpy.array (dim, )
+        :return: [nCx, nCy, nCz]
+        """
+        return np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None])
+
+    @property
+    def nNx(self):
+        """
+        Number of nodes in the x-direction
+
+        :rtype: int
+        :return: nNx
+        """
+        return self.nCx + 1
+
+    @property
+    def nNy(self):
+        """
+        Number of noes in the y-direction
+
+        :rtype: int
+        :return: nNy or None if dim < 2
+        """
+        return None if self.dim < 2 else self.nCy + 1
+
+    @property
+    def nNz(self):
+        """
+        Number of nodes in the z-direction
+
+        :rtype: int
+        :return: nNz or None if dim < 3
+        """
+        return None if self.dim < 3 else self.nCz + 1
+
+    @property
+    def vnN(self):
+        """
+        Total number of nodes in each direction
+
+        :rtype: numpy.array (dim, )
+        :return: [nNx, nNy, nNz]
+        """
+        return np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None])
+
+    @property
+    def vnEx(self):
+        """
+        Number of x-edges in each direction
+
+        :rtype: numpy.array (dim, )
+        :return: vnEx
+        """
+        return np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None])
+
+    @property
+    def vnEy(self):
+        """
+        Number of y-edges in each direction
+
+        :rtype: numpy.array (dim, )
+        :return: vnEy or None if dim < 2
+        """
+        return None if self.dim < 2 else np.array([x for x in [self.nNx, self.nCy, self.nNz] if not x is None])
+
+    @property
+    def vnEz(self):
+        """
+        Number of z-edges in each direction
+
+        :rtype: numpy.array (dim, )
+        :return: vnEz or None if dim < 3
+        """
+        return None if self.dim < 3 else np.array([x for x in [self.nNx, self.nNy, self.nCz] if not x is None])
+
+    @property
+    def vnFx(self):
+        """
+        Number of x-faces in each direction
+
+        :rtype: numpy.array (dim, )
+        :return: vnFx
+        """
+        return np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None])
+
+    @property
+    def vnFy(self):
+        """
+        Number of y-faces in each direction
+
+        :rtype: numpy.array (dim, )
+        :return: vnFy or None if dim < 2
+        """
+        return None if self.dim < 2 else np.array([x for x in [self.nCx, self.nNy, self.nCz] if not x is None])
+
+    @property
+    def vnFz(self):
+        """
+        Number of z-faces in each direction
+
+        :rtype: numpy.array (dim, )
+        :return: vnFz or None if dim < 3
+        """
+        return None if self.dim < 3 else np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None])
+
     def r(self, x, xType='CC', outType='CC', format='V'):
         """
         Mesh.r is a quick reshape command that will do the best it can at giving you what you want.
@@ -100,13 +465,13 @@ class BaseMesh(object):
             elif xType in ['F', 'E']:
                 # This will only deal with components of fields, not full 'F' or 'E'
                 xx = Utils.mkvc(xx)  # unwrap it in case it is a matrix
-                nn = self.nFv if xType == 'F' else self.nEv
+                nn = self.vnF if xType == 'F' else self.vnE
                 nn = np.r_[0, nn]
 
                 nx = [0, 0, 0]
-                nx[0] = self.nFx if xType == 'F' else self.nEx
-                nx[1] = self.nFy if xType == 'F' else self.nEy
-                nx[2] = self.nFz if xType == 'F' else self.nEz
+                nx[0] = self.vnFx if xType == 'F' else self.vnEx
+                nx[1] = self.vnFy if xType == 'F' else self.vnEy
+                nx[2] = self.vnFz if xType == 'F' else self.vnEz
 
                 for dim, dimName in enumerate(['x', 'y', 'z']):
                     if dimName in outType:
@@ -118,11 +483,11 @@ class BaseMesh(object):
             elif xTypeIsFExyz:
                 # This will deal with partial components (x, y or z) lying on edges or faces
                 if 'x' in xType:
-                    nn = self.nFx if 'F' in xType else self.nEx
+                    nn = self.vnFx if 'F' in xType else self.vnEx
                 elif 'y' in xType:
-                    nn = self.nFy if 'F' in xType else self.nEy
+                    nn = self.vnFy if 'F' in xType else self.vnEy
                 elif 'z' in xType:
-                    nn = self.nFz if 'F' in xType else self.nEz
+                    nn = self.vnFz if 'F' in xType else self.vnEz
                 assert xx.size == np.prod(nn), 'Vector is not the right size.'
                 return outKernal(xx, nn)
 
@@ -144,301 +509,3 @@ class BaseMesh(object):
             return out
         else:
             return switchKernal(x)
-
-    @property
-    def dim(self):
-        """
-        The dimension of the mesh (1, 2, or 3).
-
-        :rtype: int
-        :return: dim
-        """
-        return len(self._n)
-
-    @property
-    def nCx(self):
-        """
-        Number of cells in the x direction
-
-        :rtype: int
-        :return: nCx
-        """
-        return self._n[0]
-
-    @property
-    def nCy(self):
-        """
-        Number of cells in the y direction
-
-        :rtype: int
-        :return: nCy or None if dim < 2
-        """
-        return None if self.dim < 2 else self._n[1]
-
-    @property
-    def nCz(self):
-        """Number of cells in the z direction
-
-        :rtype: int
-        :return: nCz or None if dim < 3
-        """
-        return None if self.dim < 3 else self._n[2]
-
-    @property
-    def nCv(self):
-        """
-        Total number of cells in each direction
-
-        :rtype: numpy.array (dim, )
-        :return: [nCx, nCy, nCz]
-        """
-        return np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None])
-
-    @property
-    def nC(self):
-        doc = """
-        Total number of cells in the model.
-
-        :rtype: int
-        :return: nC
-
-        .. plot::
-            :include-source:
-
-            from SimPEG import Mesh, np
-            Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(centers=True,showIt=True)
-        """
-        return self.nCv.prod()
-
-    @property
-    def nNx(self):
-        """
-        Number of nodes in the x-direction
-
-        :rtype: int
-        :return: nNx
-        """
-        return self.nCx + 1
-
-    @property
-    def nNy(self):
-        """
-        Number of noes in the y-direction
-
-        :rtype: int
-        :return: nNy or None if dim < 2
-        """
-        return None if self.dim < 2 else self.nCy + 1
-
-    @property
-    def nNz(self):
-        """
-        Number of nodes in the z-direction
-
-        :rtype: int
-        :return: nNz or None if dim < 3
-        """
-        return None if self.dim < 3 else self.nCz + 1
-
-    @property
-    def nNv(self):
-        """
-        Total number of nodes in each direction
-
-        :rtype: numpy.array (dim, )
-        :return: [nNx, nNy, nNz]
-        """
-        return np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None])
-
-    @property
-    def nN(self):
-        doc = """
-        Total number of nodes
-
-        :rtype: int
-        :return: nN
-
-        .. plot::
-            :include-source:
-
-            from SimPEG import Mesh, np
-            Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(nodes=True,showIt=True)
-        """
-        return self.nNv.prod()
-
-    @property
-    def nEx(self):
-        """
-        Number of x-edges in each direction
-
-        :rtype: numpy.array (dim, )
-        :return: nEx
-        """
-        return np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None])
-
-    @property
-    def nEy(self):
-        """
-        Number of y-edges in each direction
-
-        :rtype: numpy.array (dim, )
-        :return: nEy or None if dim < 2
-        """
-        return None if self.dim < 2 else np.array([x for x in [self.nNx, self.nCy, self.nNz] if not x is None])
-
-    @property
-    def nEz(self):
-        """
-        Number of z-edges in each direction
-
-        :rtype: numpy.array (dim, )
-        :return: nEz or None if dim < 3
-        """
-        return None if self.dim < 3 else np.array([x for x in [self.nNx, self.nNy, self.nCz] if not x is None])
-
-    @property
-    def nEv(self):
-        """
-        Total number of edges in each direction
-
-        :rtype: numpy.array (dim, )
-        :return: [prod(nEx), prod(nEy), prod(nEz)]
-
-        .. plot::
-            :include-source:
-
-            from SimPEG import Mesh, np
-            Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(edges=True,showIt=True)
-        """
-        return np.array([np.prod(x) for x in [self.nEx, self.nEy, self.nEz] if not x is None])
-
-
-    @property
-    def nE(self):
-        """
-        Total number of edges.
-
-        :rtype: int
-        :return: sum([prod(nEx), prod(nEy), prod(nEz)])
-
-        """
-        return self.nEv.sum()
-
-    @property
-    def nFx(self):
-        """
-        Number of x-faces in each direction
-
-        :rtype: numpy.array (dim, )
-        :return: nFx
-        """
-        return np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None])
-
-    @property
-    def nFy(self):
-        """
-        Number of y-faces in each direction
-
-        :rtype: numpy.array (dim, )
-        :return: nFy or None if dim < 2
-        """
-        return None if self.dim < 2 else np.array([x for x in [self.nCx, self.nNy, self.nCz] if not x is None])
-
-    @property
-    def nFz(self):
-        """
-        Number of z-faces in each direction
-
-        :rtype: numpy.array (dim, )
-        :return: nFz or None if dim < 3
-        """
-        return None if self.dim < 3 else np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None])
-
-    @property
-    def nFv(self):
-        """
-        Total number of faces in each direction
-
-        :rtype: numpy.array (dim, )
-        :return: [prod(nFx), prod(nFy), prod(nFz)]
-
-        .. plot::
-            :include-source:
-
-            from SimPEG import Mesh, np
-            Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(faces=True,showIt=True)
-        """
-        return np.array([np.prod(x) for x in [self.nFx, self.nFy, self.nFz] if not x is None])
-
-
-    @property
-    def nF(self):
-        """
-        Total number of faces.
-
-        :rtype: int
-        :return: sum([nFx, nFy, nFz])
-
-        """
-        return self.nFv.sum()
-
-    @property
-    def normals(self):
-        """
-        Face Normals
-
-        :rtype: numpy.array (sum(nF), dim)
-        :return: normals
-        """
-        if self.dim == 2:
-            nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0])]
-            nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1])]
-            return np.r_[nX, nY]
-        elif self.dim == 3:
-            nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0]), np.zeros(self.nFv[0])]
-            nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1]), np.zeros(self.nFv[1])]
-            nZ = np.c_[np.zeros(self.nFv[2]), np.zeros(self.nFv[2]), np.ones(self.nFv[2])]
-            return np.r_[nX, nY, nZ]
-
-    @property
-    def tangents(self):
-        """
-        Edge Tangents
-
-        :rtype: numpy.array (sum(nE), dim)
-        :return: normals
-        """
-        if self.dim == 2:
-            tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0])]
-            tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1])]
-            return np.r_[tX, tY]
-        elif self.dim == 3:
-            tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0]), np.zeros(self.nEv[0])]
-            tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1]), np.zeros(self.nEv[1])]
-            tZ = np.c_[np.zeros(self.nEv[2]), np.zeros(self.nEv[2]), np.ones(self.nEv[2])]
-            return np.r_[tX, tY, tZ]
-
-
-    def projectFaceVector(self, fV):
-        """
-        Given a vector, fV, in cartesian coordinates, this will project it onto the mesh using the normals
-
-        :param numpy.array fV: face vector with shape (nF, dim)
-        :rtype: numpy.array with shape (nF, )
-        :return: projected face vector
-        """
-        assert type(fV) == np.ndarray, 'fV must be an ndarray'
-        assert len(fV.shape) == 2 and fV.shape[0] == np.sum(self.nF) and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)'
-        return np.sum(fV*self.normals, 1)
-
-    def projectEdgeVector(self, eV):
-        """
-        Given a vector, eV, in cartesian coordinates, this will project it onto the mesh using the tangents
-
-        :param numpy.array eV: edge vector with shape (nE, dim)
-        :rtype: numpy.array with shape (nE, )
-        :return: projected edge vector
-        """
-        assert type(eV) == np.ndarray, 'eV must be an ndarray'
-        assert len(eV.shape) == 2 and eV.shape[0] == np.sum(self.nE) and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)'
-        return np.sum(eV*self.tangents, 1)

SimPEG/Mesh/Cyl1DMesh.py

@@ -83,11 +83,11 @@ class Cyl1DMesh(object):
         return locals()
     nC = property(**nC())
 
-    def nCv():
+    def vnC():
         doc = "Total number of cells in each direction"
         fget = lambda self: np.array([self.nCx, self.nCz])
         return locals()
-    nCv = property(**nCv())
+    vnC = property(**vnC())
 
     def nNr():
         doc = "Number of nodes in the radial direction"
@@ -113,21 +113,21 @@ class Cyl1DMesh(object):
         return locals()
     nFr = property(**nFr())
 
-    def nFz():
+    def vnFz():
         doc = "Number of z faces"
         fget = lambda self: self.nNz * self.nCx
         return locals()
-    nFz = property(**nFz())
+    vnFz = property(**vnFz())
 
-    def nFv():
+    def vnF():
         doc = "Total number of faces in each direction"
-        fget = lambda self: np.array([self.nFr, self.nFz])
+        fget = lambda self: np.array([self.nFr, self.vnFz])
         return locals()
-    nFv = property(**nFv())
+    vnF = property(**vnF())
 
     def nF():
         doc = "Total number of faces"
-        fget = lambda self: self.nFr + self.nFz
+        fget = lambda self: self.nFr + self.vnFz
         return locals()
     nF = property(**nF())
 

SimPEG/Mesh/DiffOperators.py

@@ -129,7 +129,7 @@ class DiffOperators(object):
         def fget(self):
             if(self._faceDiv is None):
                 # The number of cell centers in each direction
-                n = self.nCv
+                n = self.vnC
                 # Compute faceDivergence operator on faces
                 if(self.dim == 1):
                     D = ddx(n[0])
@@ -158,7 +158,7 @@ class DiffOperators(object):
         def fget(self):
             if(self._faceDivx is None):
                 # The number of cell centers in each direction
-                n = self.nCv
+                n = self.vnC
                 # Compute faceDivergence operator on faces
                 if(self.dim == 1):
                     D1 = ddx(n[0])
@@ -183,7 +183,7 @@ class DiffOperators(object):
             if(self.dim < 2): return None
             if(self._faceDivy is None):
                 # The number of cell centers in each direction
-                n = self.nCv
+                n = self.vnC
                 # Compute faceDivergence operator on faces
                 if(self.dim == 2):
                     D2 = sp.kron(ddx(n[1]), speye(n[0]))
@@ -225,7 +225,7 @@ class DiffOperators(object):
         def fget(self):
             if(self._nodalGrad is None):
                 # The number of cell centers in each direction
-                n = self.nCv
+                n = self.vnC
                 # Compute divergence operator on faces
                 if(self.dim == 1):
                     G = ddx(n[0])
@@ -253,7 +253,7 @@ class DiffOperators(object):
             if(self._nodalLaplacian is None):
                 print 'Warning: Laplacian has not been tested rigorously.'
                 # The number of cell centers in each direction
-                n = self.nCv
+                n = self.vnC
                 # Compute divergence operator on faces
                 if(self.dim == 1):
                     D1 = sdiag(1./self.hx) * ddx(mesh.nCx)
@@ -291,8 +291,8 @@ class DiffOperators(object):
 
         """
         if(type(BC) is str):
-            BC = [BC for _ in self.nCv]  # Repeat the str self.dim times
-        elif(type(BC) is list):
+            BC = [BC]*self.dim
+        if(type(BC) is list):
             assert len(BC) == self.dim, 'BC list must be the size of your mesh'
         else:
             raise Exception("BC must be a str or a list.")
@@ -313,7 +313,7 @@ class DiffOperators(object):
         def fget(self):
             if(self._cellGrad is None):
                 BC = self.setCellGradBC(self._cellGradBC_list)
-                n = self.nCv
+                n = self.vnC
                 if(self.dim == 1):
                     G = ddxCellGrad(n[0], BC[0])
                 elif(self.dim == 2):
@@ -340,7 +340,7 @@ class DiffOperators(object):
         def fget(self):
             if(self._cellGradBC is None):
                 BC = self.setCellGradBC(self._cellGradBC_list)
-                n = self.nCv
+                n = self.vnC
                 if(self.dim == 1):
                     G = ddxCellGradBC(n[0], BC[0])
                 elif(self.dim == 2):
@@ -367,7 +367,7 @@ class DiffOperators(object):
         def fget(self):
             if getattr(self, '_cellGradx', None) is None:
                 BC = ['neumann', 'neumann']
-                n = self.nCv
+                n = self.vnC
                 if(self.dim == 1):
                     G1 = ddxCellGrad(n[0], BC)
                 elif(self.dim == 2):
@@ -388,7 +388,7 @@ class DiffOperators(object):
             if self.dim < 2: return None
             if getattr(self, '_cellGrady', None) is None:
                 BC = ['neumann', 'neumann']
-                n = self.nCv
+                n = self.vnC
                 if(self.dim == 2):
                     G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
                 elif(self.dim == 3):
@@ -460,13 +460,115 @@ class DiffOperators(object):
     _edgeCurl = None
     edgeCurl = property(**edgeCurl())
 
+    def getBCProjWF(self, BC, discretization='CC'):
+        """
+
+        The weak form boundary condition projection matrices.
+
+        Examples::
+
+            BC = 'neumann'                                            # Neumann in all directions
+            BC = ['neumann', 'dirichlet', 'neumann']                  # 3D, Dirichlet in y Neumann else
+            BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain,
+                                                                      #     Dirichlet else
+
+        """
+
+        if discretization is not 'CC':
+            raise NotImplementedError('Boundary conditions only implemented for CC discretization.')
+
+        if(type(BC) is str):
+            BC = [BC for _ in self.vnC]  # Repeat the str self.dim times
+        elif(type(BC) is list):
+            assert len(BC) == self.dim, 'BC list must be the size of your mesh'
+        else:
+            raise Exception("BC must be a str or a list.")
+
+        for i, bc_i in enumerate(BC):
+            BC[i] = checkBC(bc_i)
+
+
+        def projDirichlet(n, bc):
+            bc = checkBC(bc)
+            ij = ([0,n], [0,1])
+            vals = [0,0]
+            if(bc[0] == 'dirichlet'):
+                vals[0] = -1
+            if(bc[1] == 'dirichlet'):
+                vals[1] = 1
+            return sp.csr_matrix((vals, ij), shape=(n+1,2))
+
+        def projNeumannIn(n, bc):
+            bc = checkBC(bc)
+            P = sp.identity(n+1).tocsr()
+            if(bc[0] == 'neumann'):
+                P = P[1:,:]
+            if(bc[1] == 'neumann'):
+                P = P[:-1,:]
+            return P
+
+        def projNeumannOut(n, bc):
+            bc = checkBC(bc)
+            ij   = ([0, 1],[0, n])
+            vals = [0,0]
+            if(bc[0] == 'neumann'):
+                vals[0] = 1
+            if(bc[1] == 'neumann'):
+                vals[1] = 1
+            return sp.csr_matrix((vals, ij), shape=(2,n+1))
+
+        n = self.vnC
+        indF = self.faceBoundaryInd
+        if(self.dim == 1):
+            Pbc = projDirichlet(n[0], BC[0])
+            indF = indF[0] | indF[1]
+            Pbc = Pbc*sdiag(self.area[indF])
+
+            Pin = projNeumannIn(n[0], BC[0])
+
+            Pout = projNeumannOut(n[0], BC[0])
+        elif(self.dim == 2):
+            Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
+            Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
+            Pbc = sp.block_diag((Pbc1, Pbc2), format="csr")
+            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])]
+            Pbc = Pbc*sdiag(self.area[indF])
+
+            P1 = sp.kron(speye(n[1]), projNeumannIn(n[0], BC[0]))
+            P2 = sp.kron(projNeumannIn(n[1], BC[1]), speye(n[0]))
+            Pin = sp.block_diag((P1, P2), format="csr")
+
+            P1 = sp.kron(speye(n[1]), projNeumannOut(n[0], BC[0]))
+            P2 = sp.kron(projNeumannOut(n[1], BC[1]), speye(n[0]))
+            Pout = sp.block_diag((P1, P2), format="csr")
+        elif(self.dim == 3):
+            Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0]))
+            Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0]))
+            Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
+            Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
+            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])]
+            Pbc = Pbc*sdiag(self.area[indF])
+
+            P1 = kron3(speye(n[2]), speye(n[1]), projNeumannIn(n[0], BC[0]))
+            P2 = kron3(speye(n[2]), projNeumannIn(n[1], BC[1]), speye(n[0]))
+            P3 = kron3(projNeumannIn(n[2], BC[2]), speye(n[1]), speye(n[0]))
+            Pin = sp.block_diag((P1, P2, P3), format="csr")
+
+            P1 = kron3(speye(n[2]), speye(n[1]), projNeumannOut(n[0], BC[0]))
+            P2 = kron3(speye(n[2]), projNeumannOut(n[1], BC[1]), speye(n[0]))
+            P3 = kron3(projNeumannOut(n[2], BC[2]), speye(n[1]), speye(n[0]))
+            Pout = sp.block_diag((P1, P2, P3), format="csr")
+
+        return Pbc, Pin, Pout
+
+
     # --------------- Averaging ---------------------
 
     @property
     def aveF2CC(self):
         "Construct the averaging operator on cell faces to cell centers."
         if getattr(self, '_aveF2CC', None) is None:
-            n = self.nCv
+            n = self.vnC
             if(self.dim == 1):
                 self._aveF2CC = av(n[0])
             elif(self.dim == 2):
@@ -483,7 +585,7 @@ class DiffOperators(object):
     def aveF2CCV(self):
         "Construct the averaging operator on cell faces to cell centers."
         if getattr(self, '_aveF2CCV', None) is None:
-            n = self.nCv
+            n = self.vnC
             if(self.dim == 1):
                 self._aveF2CCV = av(n[0])
             elif(self.dim == 2):
@@ -499,7 +601,7 @@ class DiffOperators(object):
     def aveCC2F(self):
         "Construct the averaging operator on cell cell centers to faces."
         if getattr(self, '_aveCC2F', None) is None:
-            n = self.nCv
+            n = self.vnC
             if(self.dim == 1):
                 self._aveCC2F = avExtrap(n[0])
             elif(self.dim == 2):
@@ -516,7 +618,7 @@ class DiffOperators(object):
         "Construct the averaging operator on cell edges to cell centers."
         if getattr(self, '_aveE2CC', None) is None:
             # The number of cell centers in each direction
-            n = self.nCv
+            n = self.vnC
             if(self.dim == 1):
                 raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
             elif(self.dim == 2):
@@ -533,7 +635,7 @@ class DiffOperators(object):
         "Construct the averaging operator on cell edges to cell centers."
         if getattr(self, '_aveE2CCV', None) is None:
             # The number of cell centers in each direction
-            n = self.nCv
+            n = self.vnC
             if(self.dim == 1):
                 raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
             elif(self.dim == 2):
@@ -550,7 +652,7 @@ class DiffOperators(object):
         "Construct the averaging operator on cell nodes to cell centers."
         if getattr(self, '_aveN2CC', None) is None:
             # The number of cell centers in each direction
-            n = self.nCv
+            n = self.vnC
             if(self.dim == 1):
                 self._aveN2CC = av(n[0])
             elif(self.dim == 2):
@@ -565,7 +667,7 @@ class DiffOperators(object):
 
         if getattr(self, '_aveN2E', None) is None:
             # The number of cell centers in each direction
-            n = self.nCv
+            n = self.vnC
             if(self.dim == 1):
                 self._aveN2E = av(n[0])
             elif(self.dim == 2):
@@ -582,7 +684,7 @@ class DiffOperators(object):
         "Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate."
         if getattr(self, '_aveN2F', None) is None:
             # The number of cell centers in each direction
-            n = self.nCv
+            n = self.vnC
             if(self.dim == 1):
                 self._aveN2F = av(n[0])
             elif(self.dim == 2):

SimPEG/Mesh/InnerProducts.py

@@ -1,176 +1,55 @@
 from scipy import sparse as sp
-from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
+from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor
 import numpy as np
 
 
 class InnerProducts(object):
     """
-        Class creates the inner product matrices that you need!
-
-        InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class.
-
-
-        **Example problem for DC resistivity**
-
-        .. math::
-
-            \sigma^{-1}\mathbf{J} = \\nabla \phi
-
-        We can define in weak form by integrating with a general face function F:
-
-        .. math::
-
-            \int_{\\text{cell}}{\sigma^{-1}\mathbf{J} \cdot \mathbf{F}} = \int_{\\text{cell}}{\\nabla \phi  \cdot \mathbf{F}}
-
-            \int_{\\text{cell}}{\sigma^{-1}\mathbf{J} \cdot \mathbf{F}} = \int_{\\text{cell}}{(\\nabla \cdot \mathbf{F}) \phi   } + \int_{\partial \\text{cell}}{ \phi  \mathbf{F} \cdot \mathbf{n}}
-
-        We can then discretize for every cell:
-
-        .. math::
-
-            v_{\\text{cell}} \sigma^{-1} (\mathbf{J}_x \mathbf{F}_x +\mathbf{J}_y \mathbf{F}_y  + \mathbf{J}_z \mathbf{F}_z ) = -\phi^{\\top} v_{\\text{cell}} (\mathbf{D}_{\\text{cell}} \mathbf{F})  + \\text{BC}
-
-        We can represent this in vector form (again this is for every cell), and will generalize for the case of anisotropic (tensor) sigma.
-
-        .. math::
-
-            \mathbf{F}_c^{\\top} (\sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}})  \mathbf{J}_c = -\phi^{\\top} v_{\\text{cell}}( v_\\text{cell}^{-1} \mathbf{D}_{\\text{cell}} \mathbf{A} \mathbf{F})  + \\text{BC}
-
-        We multiply by volume on each side of the tensor conductivity to keep symmetry in the system. Here J_c is the Cartesian J (on the faces) and must be calculated differently depending on the mesh:
-
-        .. math::
-            \mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\\text{TENSOR} = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\\text{LOM}
-
-        Here the i index refers to where we choose to approximate this integral.
-        We will approximate this relation at every node of the cell, there are 8 in 3D, using a projection matrix Q_i to pick the appropriate fluxes.
-        We will then average to the cell center. For the TENSOR mesh, this looks like:
-
-        .. math::
-
-            \mathbf{F}^{\\top}
-                {1\over 8}
-                \left(\sum_{i=1}^8
-                \mathbf{Q}_{(i)}^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}}  \mathbf{Q}_{(i)}
-                \\right)
-                \mathbf{J}
-                =
-                -\mathbf{F}^{\\top} \mathbf{A} \mathbf{D}_{\\text{cell}}^{\\top} \phi   + \\text{BC}
-
-            \mathbf{M}(\Sigma^{-1}) \mathbf{J}
-                =
-                -\mathbf{A} \mathbf{D}_{\\text{cell}}^{\\top} \phi   + \\text{BC}
-
-            \mathbf{M}(\Sigma^{-1}) = {1\over 8}
-                \left(\sum_{i=1}^8
-                \mathbf{Q}_{(i)}^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma^{-1} \sqrt{v_{\\text{cell}}}  \mathbf{Q}_{(i)}
-                \\right)
-
-        The M is returned if mu is set equal  to \Sigma^{-1}.
-
-        If requested (returnP=True) the projection matricies are returned as well (ordered by nodes).
-        Here each P (3*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
-
-        .. math::
-            \mathbf{P}_{(i)} =  \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
-
-        Note that this is completed for each cell in the mesh at the same time.
+        This is a base for the SimPEG.Mesh classes. This mixIn creates the all the inner product matrices that you need!
     """
     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(M, mu=None, returnP=False):
+    def getFaceInnerProduct(self, materialProperty=None, returnP=False):
         """
-            :param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
+            :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param bool returnP: returns the projection matrices
             :rtype: scipy.csr_matrix
-            :return: M, the inner product matrix (sum(nF), sum(nF))
-
-            Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows:
-
-            .. math::
-                \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{1} & 0 \\\\ 0 & 0 & \mu_{1}  \end{matrix}\\right]
-
-                \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{2} & 0 \\\\ 0 & 0 & \mu_{3}  \end{matrix}\\right]
-
-                \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{4} & \mu_{5} \\\\ \mu_{4} & \mu_{2} & \mu_{6} \\\\ \mu_{5} & \mu_{6} & \mu_{3}  \end{matrix}\\right]
-
-                \mathbf{M}(\\vec{\mu}) = {1\over 8}
-                    \left(\sum_{i=1}^8
-                    \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \\vec{\mu} \sqrt{v_{\\text{cell}}}  \mathbf{J}_c
-                    \\right)
-
-            If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
-
-                P = [P000, P100, P010, P110, P001, P101, P011, P111]
-
-            Here each P (3*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
-
-            .. math::
-                \mathbf{P}_{(i)} =  \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
-
-            Note that this is completed for each cell in the mesh at the same time.
-
-            **For 2D:**
-
-             Depending on the number of columns (either 1, 2, or 3) of mu, the material property is interpreted as follows:
-
-            .. math::
-                \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{1} \end{matrix}\\right]
-
-                \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{2} \end{matrix}\\right]
-
-                \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{3} \\\\ \mu_{3} & \mu_{2} \end{matrix}\\right]
-
-
-            .. math::
-
-                \mathbf{M}(\\vec{\mu}) = {1\over 4}
-                    \left(\sum_{i=1}^4
-                    \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \\vec{\mu} \sqrt{v_{\\text{cell}}}  \mathbf{J}_c
-                    \\right)
-
-
-            If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
-
-                P = [P00, P10, P01, P11]
-
-            Here each P (2*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
-
-            .. math::
-                \mathbf{P}_{(i)} =  \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
-
-            Note that this is completed for each cell in the mesh at the same time.
-
+            :return: M, the inner product matrix (nF, nF)
         """
-        if M.dim == 2:
-            # Square root of cell volume multiplied by 1/4
-            v = np.sqrt(0.25*M.vol)
-            V2 = sdiag(np.r_[v, v])  # We will multiply on each side to keep symmetry
+        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)))
 
-            Pxx = _getFacePxx(M)
-            P000 = V2*Pxx('fXm', 'fYm')
-            P100 = V2*Pxx('fXp', 'fYm')
-            P010 = V2*Pxx('fXm', 'fYp')
-            P110 = V2*Pxx('fXp', 'fYp')
-        elif M.dim == 3:
-            # Square root of cell volume multiplied by 1/8
-            v = np.sqrt(0.125*M.vol)
-            V3 = sdiag(np.r_[v, v, v])  # We will multiply on each side to keep symmetry
+        if d == 1:
+            Px = _getFacePx(self)
+            P000 = V*Px('fXm')
+            P100 = V*Px('fXp')
+        elif d == 2:
+            Pxx = _getFacePxx(self)
+            P000 = V*Pxx('fXm', 'fYm')
+            P100 = V*Pxx('fXp', 'fYm')
+            P010 = V*Pxx('fXm', 'fYp')
+            P110 = V*Pxx('fXp', 'fYp')
+        elif d == 3:
+            Pxxx = _getFacePxxx(self)
+            P000 = V*Pxxx('fXm', 'fYm', 'fZm')
+            P100 = V*Pxxx('fXp', 'fYm', 'fZm')
+            P010 = V*Pxxx('fXm', 'fYp', 'fZm')
+            P110 = V*Pxxx('fXp', 'fYp', 'fZm')
+            P001 = V*Pxxx('fXm', 'fYm', 'fZp')
+            P101 = V*Pxxx('fXp', 'fYm', 'fZp')
+            P011 = V*Pxxx('fXm', 'fYp', 'fZp')
+            P111 = V*Pxxx('fXp', 'fYp', 'fZp')
 
-            Pxxx = _getFacePxxx(M)
-            P000 = V3*Pxxx('fXm', 'fYm', 'fZm')
-            P100 = V3*Pxxx('fXp', 'fYm', 'fZm')
-            P010 = V3*Pxxx('fXm', 'fYp', 'fZm')
-            P110 = V3*Pxxx('fXp', 'fYp', 'fZm')
-            P001 = V3*Pxxx('fXm', 'fYm', 'fZp')
-            P101 = V3*Pxxx('fXp', 'fYm', 'fZp')
-            P011 = V3*Pxxx('fXm', 'fYp', 'fZp')
-            P111 = V3*Pxxx('fXp', 'fYp', 'fZp')
+        Mu = makePropertyTensor(self, materialProperty)
+        A = P000.T*Mu*P000 + P100.T*Mu*P100
+        P = [P000, P100]
 
-        Mu = _makeTensor(M, mu)
-        A = P000.T*Mu*P000 + P100.T*Mu*P100 + P010.T*Mu*P010 + P110.T*Mu*P110
-        P = [P000, P100, P010, P110]
-        if M.dim == 3:
+        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 returnP:
@@ -178,89 +57,28 @@ class InnerProducts(object):
         else:
             return A
 
-    def getEdgeInnerProduct(M, sigma=None, returnP=False):
+    def getEdgeInnerProduct(self, materialProperty=None, returnP=False):
         """
-            :param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
+            :param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param bool returnP: returns the projection matrices
             :rtype: scipy.csr_matrix
-            :return: M, the inner product matrix (sum(nE), sum(nE))
-
-
-            Depending on the number of columns (either 1, 3, or 6) of sigma, the material property is interpreted as follows:
-
-            .. math::
-                \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{1} & 0 \\\\ 0 & 0 & \sigma_{1}  \end{matrix}\\right]
-
-                \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{2} & 0 \\\\ 0 & 0 & \sigma_{3}  \end{matrix}\\right]
-
-                \Sigma = \left[\\begin{matrix} \sigma_{1} & \sigma_{4} & \sigma_{5} \\\\ \sigma_{4} & \sigma_{2} & \sigma_{6} \\\\ \sigma_{5} & \sigma_{6} & \sigma_{3}  \end{matrix}\\right]
-
-            What is returned:
-
-            .. math::
-                \mathbf{M}(\Sigma) = {1\over 8}
-                    \left(\sum_{i=1}^8
-                    \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma \sqrt{v_{\\text{cell}}}  \mathbf{J}_c
-                    \\right)
-
-            If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
-
-                P = [P000, P100, P010, P110, P001, P101, P011, P111]
-
-            Here each P (3*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
-
-            .. math::
-                \mathbf{P}_{(i)} =  \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
-
-            Note that this is completed for each cell in the mesh at the same time.
-
-            **For 2D:**
-
-            Depending on the number of columns (either 1, 2, or 3) of sigma, the material property is interpreted as follows:
-
-            .. math::
-                \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{1} \end{matrix}\\right]
-
-                \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{2} \end{matrix}\\right]
-
-                \Sigma = \left[\\begin{matrix} \sigma_{1} & \sigma_{3} \\\\ \sigma_{3} & \sigma_{2} \end{matrix}\\right]
-
-
-            .. math::
-
-                \mathbf{M}(\Sigma) = {1\over 4}
-                    \left(\sum_{i=1}^4
-                    \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma \sqrt{v_{\\text{cell}}}  \mathbf{J}_c
-                    \\right)
-
-
-            If requested (returnP=True) the projection matricies are returned as well (ordered by nodes)::
-
-                P = [P00, P10, P01, P11]
-
-            Here each P (2*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates:
-
-            .. math::
-                \mathbf{P}_{(i)} =  \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)}
-
-            Note that this is completed for each cell in the mesh at the same time.
-
+            :return: M, the inner product matrix (nE, nE)
         """
+        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')
         # We will multiply by V on each side to keep symmetry
-        if M.dim == 2:
-            # Square root of cell volume multiplied by 1/4
-            v = np.sqrt(0.25*M.vol)
-            V = sdiag(np.r_[v, v])
-            eP = _getEdgePxx(M)
+        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 M.dim == 3:
-            # Square root of cell volume multiplied by 1/8
-            v = np.sqrt(0.125*M.vol)
-            V = sdiag(np.r_[v, v, v])
-            eP = _getEdgePxxx(M)
+        elif d == 3:
+            eP = _getEdgePxxx(self)
             P000 = V*eP('eX0', 'eY0', 'eZ0')
             P100 = V*eP('eX0', 'eY1', 'eZ1')
             P010 = V*eP('eX1', 'eY0', 'eZ2')
@@ -270,11 +88,11 @@ class InnerProducts(object):
             P011 = V*eP('eX3', 'eY2', 'eZ2')
             P111 = V*eP('eX3', 'eY3', 'eZ3')
 
-        Sigma = _makeTensor(M, sigma)
-        A = P000.T*Sigma*P000 + P100.T*Sigma*P100 + P010.T*Sigma*P010 + P110.T*Sigma*P110
+        Mu = makePropertyTensor(self, materialProperty)
+        A = P000.T*Mu*P000 + P100.T*Mu*P100 + P010.T*Mu*P010 + P110.T*Mu*P110
         P = [P000, P100, P010, P110]
-        if M.dim == 3:
-            A = A + P001.T*Sigma*P001 + P101.T*Sigma*P101 + P011.T*Sigma*P011 + P111.T*Sigma*P111
+        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 returnP:
             return A, P
@@ -300,32 +118,10 @@ class InnerProducts(object):
 #         |                                     |/
 #    node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
 
-def _makeTensor(M, sigma):
-    if sigma is None:  # default is ones
-        sigma = np.ones((M.nC, 1))
 
-    if M.dim == 2:
-        if sigma.size == M.nC:  # Isotropic!
-            sigma = mkvc(sigma)  # ensure it is a vector.
-            Sigma = sdiag(np.r_[sigma, sigma])
-        elif sigma.shape[1] == 2:  # Diagonal tensor
-            Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]])
-        elif sigma.shape[1] == 3:  # Fully anisotropic
-            row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2])))
-            row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1])))
-            Sigma = sp.vstack((row1, row2))
-    elif M.dim == 3:
-        if sigma.size == M.nC:  # Isotropic!
-            sigma = mkvc(sigma)  # ensure it is a vector.
-            Sigma = sdiag(np.r_[sigma, sigma, sigma])
-        elif sigma.shape[1] == 3:  # Diagonal tensor
-            Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]])
-        elif sigma.shape[1] == 6:  # Fully anisotropic
-            row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4])))
-            row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5])))
-            row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2])))
-            Sigma = sp.vstack((row1, row2, row3))
-    return Sigma
+def _getFacePx(M):
+    assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh'
+    return _getFacePx_Rectangular(M)
 
 def _getFacePxx(M):
     if M._meshType == 'TREE':
@@ -351,6 +147,23 @@ def _getEdgePxxx(M):
 
     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
 
@@ -373,11 +186,11 @@ def _getFacePxx_Rectangular(M):
                   0               1
                         f2(Ym)
 
-        Pxx('m','m') = | 1, 0, 0, 0 |
-                       | 0, 0, 1, 0 |
+        Pxx('fXm','fYm') = | 1, 0, 0, 0 |
+                           | 0, 0, 1, 0 |
 
-        Pxx('p','m') = | 0, 1, 0, 0 |
-                       | 0, 0, 1, 0 |
+        Pxx('fXp','fYm') = | 0, 1, 0, 0 |
+                           | 0, 0, 1, 0 |
 
         """
     i, j = np.int64(range(M.nCx)), np.int64(range(M.nCy))
@@ -403,12 +216,12 @@ def _getFacePxx_Rectangular(M):
         posFx = 0 if xFace == 'fXm' else 1
         posFy = 0 if yFace == 'fYm' else 1
 
-        ind1 = sub2ind(M.nFx, np.c_[ii + posFx, jj])
-        ind2 = sub2ind(M.nFy, np.c_[ii, jj + posFy]) + M.nFv[0]
+        ind1 = sub2ind(M.vnFx, np.c_[ii + posFx, jj])
+        ind2 = sub2ind(M.vnFy, np.c_[ii, jj + posFy]) + M.nFx
 
         IND = np.r_[ind1, ind2].flatten()
 
-        PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, np.sum(M.nF)))
+        PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF))
 
         if M._meshType == 'LOM':
             I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]),
@@ -459,13 +272,13 @@ def _getFacePxxx_Rectangular(M):
         posY = 0 if yFace == 'fYm' else 1
         posZ = 0 if zFace == 'fZm' else 1
 
-        ind1 = sub2ind(M.nFx, np.c_[ii + posX, jj, kk])
-        ind2 = sub2ind(M.nFy, np.c_[ii, jj + posY, kk]) + M.nFv[0]
-        ind3 = sub2ind(M.nFz, np.c_[ii, jj, kk + posZ]) + M.nFv[0] + M.nFv[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, np.sum(M.nF))).tocsr()
+        PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr()
 
         if M._meshType == 'LOM':
             I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]),
@@ -495,12 +308,12 @@ def _getEdgePxx_Rectangular(M):
         posX = 0 if xEdge == 'eX0' else 1
         posY = 0 if yEdge == 'eY0' else 1
 
-        ind1 = sub2ind(M.nEx, np.c_[ii, jj + posX])
-        ind2 = sub2ind(M.nEy, np.c_[ii + posY, jj]) + M.nEv[0]
+        ind1 = sub2ind(M.vnEx, np.c_[ii, jj + posX])
+        ind2 = sub2ind(M.vnEy, np.c_[ii + posY, jj]) + M.nEx
 
         IND = np.r_[ind1, ind2].flatten()
 
-        PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, np.sum(M.nE))).tocsr()
+        PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr()
 
         if M._meshType == 'LOM':
             I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]),
@@ -537,13 +350,13 @@ def _getEdgePxxx_Rectangular(M):
         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.nEx, np.c_[ii, jj + posX[0], kk + posX[1]])
-        ind2 = sub2ind(M.nEy, np.c_[ii + posY[0], jj, kk + posY[1]]) + M.nEv[0]
-        ind3 = sub2ind(M.nEz, np.c_[ii + posZ[0], jj + posZ[1], kk]) + M.nEv[0] + M.nEv[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, np.sum(M.nE))).tocsr()
+        PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr()
 
         if M._meshType == 'LOM':
             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]]),

SimPEG/Mesh/LogicallyOrthogonalMesh.py

@@ -1,5 +1,5 @@
 from SimPEG import Utils, np
-from BaseMesh import BaseMesh
+from BaseMesh import BaseRectangularMesh
 from DiffOperators import DiffOperators
 from InnerProducts import InnerProducts
 from LomView import LomView
@@ -11,7 +11,7 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
 normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
 
 
-class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
+class LogicallyOrthogonalMesh(BaseRectangularMesh, DiffOperators, InnerProducts, LomView):
     """
     LogicallyOrthogonalMesh is a mesh class that deals with logically orthogonal meshes.
 
@@ -32,6 +32,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
 
     def __init__(self, nodes):
         assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray"
+        assert len(nodes) > 1, "len(node) must be greater than 1"
 
         for i, nodes_i in enumerate(nodes):
             assert type(nodes_i) == np.ndarray, ("nodes[%i] is not a numpy array." % i)
@@ -40,7 +41,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
         assert len(nodes[0].shape) == len(nodes), "Dimension mismatch"
         assert len(nodes[0].shape) > 1, "Not worth using LOM for a 1D mesh."
 
-        super(LogicallyOrthogonalMesh, self).__init__(np.array(nodes[0].shape)-1, None)
+        BaseRectangularMesh.__init__(self, np.array(nodes[0].shape)-1, None)
 
         # Save nodes to private variable _gridN as vectors
         self._gridN = np.ones((nodes[0].size, self.dim))
@@ -199,13 +200,13 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
         def fget(self):
             if(self._vol is None):
                 if self.dim == 2:
-                    A, B, C, D = Utils.indexCube('ABCD', self.nCv+1)
+                    A, B, C, D = Utils.indexCube('ABCD', self.vnC+1)
                     normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D)
                     self._vol = area
                 elif self.dim == 3:
                     # Each polyhedron can be decomposed into 5 tetrahedrons
                     # However, this presents a choice so we may as well divide in two ways and average.
-                    A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.nCv+1)
+                    A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.vnC+1)
 
                     vol1 = (Utils.volTetra(self.gridN, A, B, D, E) +  # cutted edge top
                             Utils.volTetra(self.gridN, B, E, F, G) +  # cutted edge top
@@ -233,11 +234,11 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
                 # Compute areas of cell faces
                 if(self.dim == 2):
                     xy = self.gridN
-                    A, B = Utils.indexCube('AB', self.nCv+1, np.array([self.nNx, self.nCy]))
+                    A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy]))
                     edge1 = xy[B, :] - xy[A, :]
                     normal1 = np.c_[edge1[:, 1], -edge1[:, 0]]
                     area1 = length2D(edge1)
-                    A, D = Utils.indexCube('AD', self.nCv+1, np.array([self.nCx, self.nNy]))
+                    A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy]))
                     # Note that we are doing A-D to make sure the normal points the right way.
                     # Think about it. Look at the picture. Normal points towards C iff you do this.
                     edge2 = xy[A, :] - xy[D, :]
@@ -247,13 +248,13 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
                     self._normals = [normalize2D(normal1), normalize2D(normal2)]
                 elif(self.dim == 3):
 
-                    A, E, F, B = Utils.indexCube('AEFB', self.nCv+1, np.array([self.nNx, self.nCy, self.nCz]))
+                    A, E, F, B = Utils.indexCube('AEFB', self.vnC+1, np.array([self.nNx, self.nCy, self.nCz]))
                     normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False)
 
-                    A, D, H, E = Utils.indexCube('ADHE', self.nCv+1, np.array([self.nCx, self.nNy, self.nCz]))
+                    A, D, H, E = Utils.indexCube('ADHE', self.vnC+1, np.array([self.nCx, self.nNy, self.nCz]))
                     normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False)
 
-                    A, B, C, D = Utils.indexCube('ABCD', self.nCv+1, np.array([self.nCx, self.nCy, self.nNz]))
+                    A, B, C, D = Utils.indexCube('ABCD', self.vnC+1, np.array([self.nCx, self.nCy, self.nNz]))
                     normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False)
 
                     self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)]
@@ -296,19 +297,19 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
             if(self._edge is None or self._tangents is None):
                 if(self.dim == 2):
                     xy = self.gridN
-                    A, D = Utils.indexCube('AD', self.nCv+1, np.array([self.nCx, self.nNy]))
+                    A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy]))
                     edge1 = xy[D, :] - xy[A, :]
-                    A, B = Utils.indexCube('AB', self.nCv+1, np.array([self.nNx, self.nCy]))
+                    A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy]))
                     edge2 = xy[B, :] - xy[A, :]
                     self._edge = np.r_[Utils.mkvc(length2D(edge1)), Utils.mkvc(length2D(edge2))]
                     self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge]
                 elif(self.dim == 3):
                     xyz = self.gridN
-                    A, D = Utils.indexCube('AD', self.nCv+1, np.array([self.nCx, self.nNy, self.nNz]))
+                    A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy, self.nNz]))
                     edge1 = xyz[D, :] - xyz[A, :]
-                    A, B = Utils.indexCube('AB', self.nCv+1, np.array([self.nNx, self.nCy, self.nNz]))
+                    A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy, self.nNz]))
                     edge2 = xyz[B, :] - xyz[A, :]
-                    A, E = Utils.indexCube('AE', self.nCv+1, np.array([self.nNx, self.nNy, self.nCz]))
+                    A, E = Utils.indexCube('AE', self.vnC+1, np.array([self.nNx, self.nNy, self.nCz]))
                     edge3 = xyz[E, :] - xyz[A, :]
                     self._edge = np.r_[Utils.mkvc(length3D(edge1)), Utils.mkvc(length3D(edge2)), Utils.mkvc(length3D(edge3))]
                     self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge]

SimPEG/Mesh/TensorMesh.py

@@ -1,10 +1,10 @@
 from SimPEG import Utils, np, sp
-from BaseMesh import BaseMesh
+from BaseMesh import BaseRectangularMesh
 from TensorView import TensorView
 from DiffOperators import DiffOperators
 from InnerProducts import InnerProducts
 
-class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
+class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
     """
     TensorMesh is a mesh class that deals with tensor product meshes.
 
@@ -48,7 +48,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
             assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i)
             h[i] = h_i[:] # make a copy.
 
-        BaseMesh.__init__(self, np.array([x.size for x in h]), x0)
+        BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0)
         assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)"
 
         # Ensure h contains 1D vectors
@@ -282,7 +282,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
                 # Ensure that we are working with column vectors
                 vh = self.h
                 # The number of cell centers in each direction
-                n = self.nCv
+                n = self.vnC
                 # Compute areas of cell faces
                 if(self.dim == 1):
                     self._area = np.ones(n[0]+1)
@@ -308,7 +308,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
                 # Ensure that we are working with column vectors
                 vh = self.h
                 # The number of cell centers in each direction
-                n = self.nCv
+                n = self.vnC
                 # Compute edge lengths
                 if(self.dim == 1):
                     self._edge = Utils.mkvc(vh[0])
@@ -367,7 +367,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
         return [t for t in ten if t is not None]
 
 
-    def isInside(self, pts):
+    def isInside(self, pts, locType='N'):
         """
         Determines if a set of points are inside a mesh.
 
@@ -376,15 +376,23 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
         :return inside, numpy array of booleans
         """
 
-        pts = np.atleast_2d(pts)
-        inside = (pts[:,0] >= self.vectorNx.min()) & (pts[:,0] <= self.vectorNx.max())
+        tensors = self.getTensor(locType)
+        if type(pts) == list:
+            pts = np.array(pts)
+        assert type(pts) == np.ndarray, "must be a numpy array"
         if self.dim > 1:
-            inside = inside & ((pts[:,1] >= self.vectorNy.min()) & (pts[:,1] <= self.vectorNy.max()))
-        if self.dim > 2:
-            inside = inside & ((pts[:,2] >= self.vectorNz.min()) & (pts[:,2] <= self.vectorNz.max()))
+            assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)"
+        elif len(pts.shape) == 1:
+            pts = pts[:,np.newaxis]
+        else:
+            assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)"
+
+        inside = np.ones(pts.shape[0],dtype=bool)
+        for i, tensor in enumerate(tensors):
+            inside = inside & (pts[:,i] >= tensor.min()) & (pts[:,i] <= tensor.max())
         return inside
 
-    def getInterpolationMat(self, loc, locType):
+    def getInterpolationMat(self, loc, locType, zerosOutside=False):
         """ Produces interpolation matrix
 
         :param numpy.ndarray loc: Location of points to interpolate to
@@ -404,12 +412,25 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
             'CC'    -> scalar field defined on cell centers
         """
 
-        loc = np.atleast_2d(loc)
-        assert np.all(self.isInside(loc)), "Points outside of mesh"
+        if type(loc) == list:
+            loc = np.array(loc)
+        assert type(loc) == np.ndarray, "must be a numpy array"
+        if self.dim > 1:
+            assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)"
+        elif len(loc.shape) == 1:
+            loc = loc[:,np.newaxis]
+        else:
+            assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)"
+
+        if zerosOutside is False:
+            assert np.all(self.isInside(loc)), "Points outside of mesh"
+        else:
+            indZeros = np.logical_not(self.isInside(loc))
+            loc[indZeros, :] = np.array([v.mean() for v in self.getTensor('CC')])
 
         ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1
         if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind:
-            nF_nE = self.nFv if 'F' in locType else self.nEv
+            nF_nE = self.vnF if 'F' in locType else self.vnE
             components = [Utils.spzeros(loc.shape[0], n) for n in nF_nE]
             components[ind] = Utils.interpmat(loc, *self.getTensor(locType))
             Q = sp.hstack(components)
@@ -417,7 +438,59 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
             Q = Utils.interpmat(loc, *self.getTensor(locType))
         else:
             raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))
-        return Q
+        if zerosOutside:
+            Q[indZeros, :] = 0
+        return Q.tocsr()
+
+
+    @property
+    def faceBoundaryInd(self):
+        """
+            Find indices of boundary faces in each direction
+        """
+        if self.dim==1:
+            indxd = (self.gridFx==min(self.gridFx))
+            indxu = (self.gridFx==max(self.gridFx))
+            return indxd, indxu
+        elif self.dim==2:
+            indxd = (self.gridFx[:,0]==min(self.gridFx[:,0]))
+            indxu = (self.gridFx[:,0]==max(self.gridFx[:,0]))
+            indyd = (self.gridFy[:,1]==min(self.gridFy[:,1]))
+            indyu = (self.gridFy[:,1]==max(self.gridFy[:,1]))
+            return indxd, indxu, indyd, indyu
+        elif self.dim==3:
+            indxd = (self.gridFx[:,0]==min(self.gridFx[:,0]))
+            indxu = (self.gridFx[:,0]==max(self.gridFx[:,0]))
+            indyd = (self.gridFy[:,1]==min(self.gridFy[:,1]))
+            indyu = (self.gridFy[:,1]==max(self.gridFy[:,1]))
+            indzd = (self.gridFz[:,2]==min(self.gridFz[:,2]))
+            indzu = (self.gridFz[:,2]==max(self.gridFz[:,2]))
+            return indxd, indxu, indyd, indyu, indzd, indzu
+
+    @property
+    def cellBoundaryInd(self):
+        """
+            Find indices of boundary faces in each direction
+        """
+        if self.dim==1:
+            indxd = (self.gridCC==min(self.gridCC))
+            indxu = (self.gridCC==max(self.gridCC))
+            return indxd, indxu
+        elif self.dim==2:
+            indxd = (self.gridCC[:,0]==min(self.gridCC[:,0]))
+            indxu = (self.gridCC[:,0]==max(self.gridCC[:,0]))
+            indyd = (self.gridCC[:,1]==min(self.gridCC[:,1]))
+            indyu = (self.gridCC[:,1]==max(self.gridCC[:,1]))
+            return indxd, indxu, indyd, indyu
+        elif self.dim==3:
+            indxd = (self.gridCC[:,0]==min(self.gridCC[:,0]))
+            indxu = (self.gridCC[:,0]==max(self.gridCC[:,0]))
+            indyd = (self.gridCC[:,1]==min(self.gridCC[:,1]))
+            indyu = (self.gridCC[:,1]==max(self.gridCC[:,1]))
+            indzd = (self.gridCC[:,2]==min(self.gridCC[:,2]))
+            indzu = (self.gridCC[:,2]==max(self.gridCC[:,2]))
+            return indxd, indxu, indyd, indyu, indzd, indzu
+
 
 if __name__ == '__main__':
     print('Welcome to tensor mesh!')

SimPEG/Mesh/TensorView.py

@@ -61,20 +61,20 @@ class TensorView(object):
         elif imageType == 'N':
             assert I.size == self.nN, "Incorrect dimensions for N."
         elif imageType == 'Fx':
-            if I.size != np.prod(self.nFx): I, fy, fz = self.r(I,'F','F','M')
+            if I.size != np.prod(self.vnFx): I, fy, fz = self.r(I,'F','F','M')
         elif imageType == 'Fy':
-            if I.size != np.prod(self.nFy): fx, I, fz = self.r(I,'F','F','M')
+            if I.size != np.prod(self.vnFy): fx, I, fz = self.r(I,'F','F','M')
         elif imageType == 'Fz':
-            if I.size != np.prod(self.nFz): fx, fy, I = self.r(I,'F','F','M')
+            if I.size != np.prod(self.vnFz): fx, fy, I = self.r(I,'F','F','M')
         elif imageType == 'Ex':
-            if I.size != np.prod(self.nEx): I, ey, ez = self.r(I,'E','E','M')
+            if I.size != np.prod(self.vnEx): I, ey, ez = self.r(I,'E','E','M')
         elif imageType == 'Ey':
-            if I.size != np.prod(self.nEy): ex, I, ez = self.r(I,'E','E','M')
+            if I.size != np.prod(self.vnEy): ex, I, ez = self.r(I,'E','E','M')
         elif imageType == 'Ez':
-            if I.size != np.prod(self.nEz): ex, ey, I = self.r(I,'E','E','M')
+            if I.size != np.prod(self.vnEz): ex, ey, I = self.r(I,'E','E','M')
         elif imageType[0] == 'E':
             plotAll = len(imageType) == 1
-            options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt}
+            options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":False}
             fig = plt.figure(figNum)
             # Determine the subplot number: 131, 121
             numPlots = 130 if plotAll else len(imageType)/2*10+100
@@ -92,10 +92,11 @@ class TensorView(object):
                 ax_z = plt.subplot(numPlots+pltNum)
                 self.plotImage(ez, imageType='Ez', ax=ax_z, **options)
                 pltNum +=1
+            if showIt: plt.show()
             return
         elif imageType[0] == 'F':
             plotAll = len(imageType) == 1
-            options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt}
+            options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":False}
             fig = plt.figure(figNum)
             # Determine the subplot number: 131, 121
             numPlots = 130 if plotAll else len(imageType)/2*10+100
@@ -113,6 +114,7 @@ class TensorView(object):
                 ax_z = plt.subplot(numPlots+pltNum)
                 self.plotImage(fxyz[2], imageType='Fz', ax=ax_z, **options)
                 pltNum +=1
+            if showIt: plt.show()
             return
         else:
             raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'")
@@ -135,21 +137,21 @@ class TensorView(object):
             ax.axis('tight')
         elif self.dim == 2:
             if imageType == 'CC':
-                C = I[:].reshape(self.nCv, order='F')
+                C = I[:].reshape(self.vnC, order='F')
             elif imageType == 'N':
-                C = I[:].reshape(self.nNv, order='F')
+                C = I[:].reshape(self.vnN, order='F')
                 C = 0.25*(C[:-1, :-1] + C[1:, :-1] + C[:-1, 1:] + C[1:, 1:])
             elif imageType == 'Fx':
-                C = I[:].reshape(self.nFx, order='F')
+                C = I[:].reshape(self.vnFx, order='F')
                 C = 0.5*(C[:-1, :] + C[1:, :] )
             elif imageType == 'Fy':
-                C = I[:].reshape(self.nFy, order='F')
+                C = I[:].reshape(self.vnFy, order='F')
                 C = 0.5*(C[:, :-1] + C[:, 1:] )
             elif imageType == 'Ex':
-                C = I[:].reshape(self.nEx, order='F')
+                C = I[:].reshape(self.vnEx, order='F')
                 C = 0.5*(C[:,:-1] + C[:,1:] )
             elif imageType == 'Ey':
-                C = I[:].reshape(self.nEy, order='F')
+                C = I[:].reshape(self.vnEy, order='F')
                 C = 0.5*(C[:-1,:] + C[1:,:] )
 
             if clim is None:
@@ -164,27 +166,27 @@ class TensorView(object):
 
                 # get copy of image and average to cell-centres is necessary
                 if imageType == 'CC':
-                    Ic = I[:].reshape(self.nCv, order='F')
+                    Ic = I[:].reshape(self.vnC, order='F')
                 elif imageType == 'N':
-                    Ic = I[:].reshape(self.nNv, order='F')
+                    Ic = I[:].reshape(self.vnN, order='F')
                     Ic = .125*(Ic[:-1,:-1,:-1]+Ic[1:,:-1,:-1] + Ic[:-1,1:,:-1]+ Ic[1:,1:,:-1]+ Ic[:-1,:-1,1:]+Ic[1:,:-1,1:] + Ic[:-1,1:,1:]+ Ic[1:,1:,1:] )
                 elif imageType == 'Fx':
-                    Ic = I[:].reshape(self.nFx, order='F')
+                    Ic = I[:].reshape(self.vnFx, order='F')
                     Ic = .5*(Ic[:-1,:,:]+Ic[1:,:,:])
                 elif imageType == 'Fy':
-                    Ic = I[:].reshape(self.nFy, order='F')
+                    Ic = I[:].reshape(self.vnFy, order='F')
                     Ic = .5*(Ic[:,:-1,:]+Ic[:,1:,:])
                 elif imageType == 'Fz':
-                    Ic = I[:].reshape(self.nFz, order='F')
+                    Ic = I[:].reshape(self.vnFz, order='F')
                     Ic = .5*(Ic[:,:,:-1]+Ic[:,:,1:])
                 elif imageType == 'Ex':
-                    Ic = I[:].reshape(self.nEx, order='F')
+                    Ic = I[:].reshape(self.vnEx, order='F')
                     Ic = .25*(Ic[:,:-1,:-1]+Ic[:,1:,:-1]+Ic[:,:-1,1:]+Ic[:,1:,:1])
                 elif imageType == 'Ey':
-                    Ic = I[:].reshape(self.nEy, order='F')
+                    Ic = I[:].reshape(self.vnEy, order='F')
                     Ic = .25*(Ic[:-1,:,:-1]+Ic[1:,:,:-1]+Ic[:-1,:,1:]+Ic[1:,:,:1])
                 elif imageType == 'Ez':
-                    Ic = I[:].reshape(self.nEz, order='F')
+                    Ic = I[:].reshape(self.vnEz, order='F')
                     Ic = .25*(Ic[:-1,:-1,:]+Ic[1:,:-1,:]+Ic[:-1,1:,:]+Ic[1:,:1,:])
 
                 # determine number oE slices in x and y dimension
@@ -238,6 +240,135 @@ class TensorView(object):
         if showIt: plt.show()
         return ph
 
+    def plotSlice(self, v, vType='CC',
+                  normal='Z', ind=None, grid=False, view='real',
+                  ax=None, clim=None, showIt=False,
+                  pcolorOpts={},
+                  streamOpts={'color':'k'},
+                  gridOpts={'color':'k'}
+                  ):
+
+        """
+        Plots a slice of a 3D mesh.
+
+        .. plot::
+
+            from SimPEG import *
+            mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2)))
+            M = Mesh.TensorMesh(mT)
+            q = np.zeros(M.vnC)
+            q[[4,4],[4,4],[2,6]]=[-1,1]
+            q = Utils.mkvc(q)
+            A = M.faceDiv*M.cellGrad
+            b = Solver(A).solve(q)
+            M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
+
+        """
+        viewOpts = ['real','imag','abs','vec']
+        normalOpts = ['X', 'Y', 'Z']
+        vTypeOpts = ['CC', 'CCv','F','E']
+
+        # Some user error checking
+        assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
+        assert self.dim == 3, 'Must be a 3D mesh.'
+        assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts)
+        assert normal in normalOpts, "normal must be in ['%s']" % "','".join(normalOpts)
+        assert type(grid) is bool, 'grid must be a boolean'
+
+        szSliceDim = getattr(self, 'nC'+normal.lower()) #: Size of the sliced dimension
+        if ind is None: ind = int(szSliceDim/2)
+        assert type(ind) in [int, long], 'ind must be an integer'
+
+        if ax is None:
+            fig = plt.figure(1)
+            fig.clf()
+            ax = plt.subplot(111)
+        else:
+            assert isinstance(ax, matplotlib.axes.Axes), "ax must be an matplotlib.axes.Axes"
+            fig = ax.figure
+
+        # The slicing and plotting code!!
+
+        def getIndSlice(v):
+            if   normal == 'X': v = v[ind,:,:]
+            elif normal == 'Y': v = v[:,ind,:]
+            elif normal == 'Z': v = v[:,:,ind]
+            return v
+
+        def doSlice(v):
+            if vType == 'CC':
+                return getIndSlice(self.r(v,'CC','CC','M'))
+            elif vType == 'CCv':
+                v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M')
+                assert view == 'vec', 'Other types for CCv not yet supported'
+            else:
+                # Now just deal with 'F' and 'E'
+                aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC')
+                v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
+                v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M')
+            if view == 'vec':
+                outSlice = []
+                if 'X' not in normal: outSlice.append(getIndSlice(v[0]))
+                if 'Y' not in normal: outSlice.append(getIndSlice(v[1]))
+                if 'Z' not in normal: outSlice.append(getIndSlice(v[2]))
+                return outSlice
+            else:
+                return getIndSlice(self.r(v,'CC','CC','M'))
+
+        h2d = []
+        if 'X' not in normal: h2d.append(self.hx)
+        if 'Y' not in normal: h2d.append(self.hy)
+        if 'Z' not in normal: h2d.append(self.hz)
+        tM = self.__class__(h2d) #: Temp Mesh
+
+        out = ()
+        if view in ['real','imag','abs']:
+            v = getattr(np,view)(v) # e.g. np.real(v)
+            v = doSlice(v)
+            if clim is None:
+                clim = [v.min(),v.max()]
+            out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, v.T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
+        elif view in ['vec']:
+            U, V = doSlice(v)
+            if clim is None:
+                uv = np.r_[mkvc(U), mkvc(V)]
+                uv = np.sqrt(uv**2)
+                clim = [uv.min(),uv.max()]
+
+            # Matplotlib seems to not support irregular
+            # spaced vectors at the moment. So we will
+            # Interpolate down to a regular mesh at the
+            # smallest mesh size in this 2D slice.
+            nxi = int(tM.hx.sum()/tM.hx.min())
+            nyi = int(tM.hy.sum()/tM.hy.min())
+            tMi = self.__class__([np.ones(nxi)*tM.hx.sum()/nxi,
+                                  np.ones(nyi)*tM.hy.sum()/nyi])
+            P = tM.getInterpolationMat(tMi.gridCC,'CC',zerosOutside=True)
+            Ui = P*mkvc(U)
+            Vi = P*mkvc(V)
+            Ui = tMi.r(Ui, 'CC', 'CC', 'M')
+            Vi = tMi.r(Vi, 'CC', 'CC', 'M')
+            # End Interpolation
+
+            out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, np.sqrt(U**2+V**2).T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
+            out += (ax.streamplot(tMi.vectorCCx, tMi.vectorCCy, Ui.T, Vi.T, **streamOpts),)
+
+        if grid:
+            xXGrid = np.c_[tM.vectorNx,tM.vectorNx,np.nan*np.ones(tM.nNx)].flatten()
+            xYGrid = np.c_[tM.vectorNy[0]*np.ones(tM.nNx),tM.vectorNy[-1]*np.ones(tM.nNx),np.nan*np.ones(tM.nNx)].flatten()
+            yXGrid = np.c_[tM.vectorNx[0]*np.ones(tM.nNy),tM.vectorNx[-1]*np.ones(tM.nNy),np.nan*np.ones(tM.nNy)].flatten()
+            yYGrid = np.c_[tM.vectorNy,tM.vectorNy,np.nan*np.ones(tM.nNy)].flatten()
+            out += (ax.plot(np.r_[xXGrid,yXGrid],np.r_[xYGrid,yYGrid],**gridOpts)[0],)
+
+        ax.set_xlabel('y' if normal == 'X' else 'x')
+        ax.set_ylabel('y' if normal == 'Z' else 'z')
+        ax.set_title('Slice %d' % ind)
+        ax.set_xlim(*tM.vectorNx[[0,-1]])
+        ax.set_ylim(*tM.vectorNy[[0,-1]])
+
+        if showIt: plt.show()
+        return out
+
     def plotGrid(self, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False):
         """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
 
@@ -395,7 +526,7 @@ class TensorView(object):
             mesh.slicer(var, imageType=imageType, normal=normal, index=i, ax=ax, clim=clim)
             tlt.set_text(normal.upper()+('-Slice: %d, %4.4f' % (i,getattr(mesh,'vectorCC'+normal)[i])))
 
-        return animate(fig, animateFrame, frames=mesh.nCv['xyz'.index(normal)])
+        return animate(fig, animateFrame, frames=mesh.vnC['xyz'.index(normal)])
 
     def video(mesh, var, function, figsize=(10, 8), colorbar=True, skip=1):
         """
@@ -426,3 +557,13 @@ class TensorView(object):
         return animate(fig, animateFrame, frames=len(frames))
 
 
+if __name__ == '__main__':
+    from SimPEG import *
+    mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2)))
+    M = Mesh.TensorMesh(mT)
+    q = np.zeros(M.vnC)
+    q[[4,4],[4,4],[2,6]]=[-1,1]
+    q = Utils.mkvc(q)
+    A = M.faceDiv*M.cellGrad
+    b = Solver(A).solve(q)
+    M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})

SimPEG/Mesh/TreeMesh.py

@@ -679,6 +679,8 @@ class TreeMesh(InnerProducts, BaseMesh):
 
     def __init__(self, h_in, x0=None):
         assert type(h_in) is list, 'h_in must be a list'
+        assert len(h_in) > 1, "len(h_in) must be greater than 1"
+
         h = range(len(h_in))
         for i, h_i in enumerate(h_in):
             if type(h_i) in [int, long, float]:

SimPEG/Mesh/__init__.py

@@ -2,7 +2,7 @@ from Cyl1DMesh import Cyl1DMesh
 from TensorMesh import TensorMesh
 from TreeMesh import TreeMesh
 from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
-from BaseMesh import BaseMesh
+from BaseMesh import BaseMesh, BaseRectangularMesh
 from TensorView import TensorView
 from LomView import LomView
 from InnerProducts import InnerProducts

SimPEG/Model.py

@@ -64,6 +64,68 @@ class BaseModel(object):
             m = self.example()
         return checkDerivative(lambda m : [self.transform(m), self.transformDeriv(m)], m, plotIt=False)
 
+class BaseNonLinearModel(object):
+    """
+    SimPEG BaseNonLinearModel
+
+    """
+
+    __metaclass__ = Utils.SimPEGMetaClass
+
+    counter = None   #: A SimPEG.Utils.Counter object
+    mesh = None      #: A SimPEG Mesh
+
+    def __init__(self, mesh):
+        self.mesh = mesh
+
+    def transform(self, u, m):
+        """
+            :param numpy.array u: fields
+            :param numpy.array m: model
+            :rtype: numpy.array
+            :return: transformed model
+
+            The *transform* changes the model into the physical property.
+
+        """
+        return m
+
+    def transformDerivU(self, u, m):
+        """
+            :param numpy.array u: fields
+            :param numpy.array m: model
+            :rtype: scipy.csr_matrix
+            :return: derivative of transformed model
+
+            The *transform* changes the model into the physical property.
+            The *transformDerivU* provides the derivative of the *transform* with respect to the fields.
+        """
+        raise NotImplementedError('The transformDerivU is not implemented.')
+
+
+    def transformDerivM(self, u, m):
+        """
+            :param numpy.array u: fields
+            :param numpy.array m: model
+            :rtype: scipy.csr_matrix
+            :return: derivative of transformed model
+
+            The *transform* changes the model into the physical property.
+            The *transformDerivU* provides the derivative of the *transform* with respect to the model.
+        """
+        raise NotImplementedError('The transformDerivM is not implemented.')
+
+    @property
+    def nP(self):
+        """Number of parameters in the model."""
+        return self.mesh.nC
+
+    def example(self):
+        raise NotImplementedError('The example is not implemented.')
+
+    def test(self, m=None):
+        raise NotImplementedError('The test is not implemented.')
+
 
 class LogModel(BaseModel):
     """SimPEG LogModel"""
@@ -150,7 +212,7 @@ class Vertical1DModel(BaseModel):
 
            The number of cells in the
            last dimension of the mesh."""
-        return self.mesh.nCv[self.mesh.dim-1]
+        return self.mesh.vnC[self.mesh.dim-1]
 
     def transform(self, m):
         """
@@ -158,7 +220,7 @@ class Vertical1DModel(BaseModel):
             :rtype: numpy.array
             :return: transformed model
         """
-        repNum = self.mesh.nCv[:self.mesh.dim-1].prod()
+        repNum = self.mesh.vnC[:self.mesh.dim-1].prod()
         return Utils.mkvc(m).repeat(repNum)
 
     def transformDeriv(self, m):
@@ -167,19 +229,116 @@ class Vertical1DModel(BaseModel):
             :rtype: scipy.csr_matrix
             :return: derivative of transformed model
         """
-        repNum = self.mesh.nCv[:self.mesh.dim-1].prod()
+        repNum = self.mesh.vnC[:self.mesh.dim-1].prod()
         repVec = sp.csr_matrix(
                     (np.ones(repNum),
                     (range(repNum), np.zeros(repNum))
                     ), shape=(repNum, 1))
         return sp.kron(sp.identity(self.nP), repVec)
 
+class Mesh2Mesh(BaseModel):
+    """
+        Takes a model on one mesh are translates it to another mesh.
+
+        .. plot::
+
+            from SimPEG import *
+            M = Mesh.TensorMesh([100,100])
+            h1 = Utils.meshTensors(((7,6,1.5),(10,6),(7,6,1.5)))
+            h1 = h1/h1.sum()
+            M2 = Mesh.TensorMesh([h1,h1])
+            V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
+            v = Utils.mkvc(V)
+            modh = Model.Mesh2Mesh([M,M2])
+            modH = Model.Mesh2Mesh([M2,M])
+            H = modH.transform(v)
+            h = modh.transform(H)
+            ax = plt.subplot(131)
+            M.plotImage(v, ax=ax)
+            ax.set_title('Fine Mesh (Original)')
+            ax = plt.subplot(132)
+            M2.plotImage(H,clim=[0,1],ax=ax)
+            ax.set_title('Course Mesh')
+            ax = plt.subplot(133)
+            M.plotImage(h,clim=[0,1],ax=ax)
+            ax.set_title('Fine Mesh (Interpolated)')
+
+    """
+
+    def __init__(self, meshes, **kwargs):
+        Utils.setKwargs(self, **kwargs)
+
+        assert type(meshes) is list, "meshes must be a list of two meshes"
+        assert len(meshes) == 2, "meshes must be a list of two meshes"
+        assert meshes[0].dim == meshes[1].dim, """The two meshes must be the same dimension"""
+
+        self.mesh  = meshes[0]
+        self.mesh2 = meshes[1]
+
+        self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC,'CC',zerosOutside=True)
+
+    @property
+    def nP(self):
+        """Number of parameters in the model."""
+        return self.mesh2.nC
+    def transform(self, m):
+        return self.P*m
+    def transformDeriv(self, m):
+        return self.P
+
+
+class ActiveModel(BaseModel):
+    """
+        Active model parameters.
+
+    """
+
+    indActive   = None #: Active Cells
+    valInactive = None #: Values of inactive Cells
+    nC          = None #: Number of cells in the full model
+
+    def __init__(self, mesh, indActive, valInactive, nC=None):
+        self.mesh  = mesh
+
+        self.nC = nC or mesh.nC
+
+        if indActive.dtype is not bool:
+            z = np.zeros(self.nC,dtype=bool)
+            z[indActive] = True
+            indActive = z
+        self.indActive = indActive
+        self.indInactive = np.logical_not(indActive)
+        if type(valInactive) in [float, int, long]:
+            valInactive = np.ones(self.nC)*float(valInactive)
+
+        valInactive[self.indActive] = 0
+        self.valInactive = valInactive
+
+        inds = np.nonzero(self.indActive)[0]
+        self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
+
+    @property
+    def nP(self):
+        """Number of parameters in the model."""
+        return self.indActive.sum()
+
+    def transform(self, m):
+        return self.P*m + self.valInactive
+    def transformDeriv(self, m):
+        return self.P
+
 class ComboModel(BaseModel):
     """Combination of various models."""
 
     def __init__(self, mesh, models, **kwargs):
         BaseModel.__init__(self, mesh, **kwargs)
-        self.models = [m(mesh, **kwargs) for m in models]
+
+        self.models = []
+        for m in models:
+            if not isinstance(m, BaseModel):
+                self.models += [m(mesh, **kwargs)]
+            else:
+                self.models += [m]
 
     @property
     def nP(self):

SimPEG/ObjFunction.py

@@ -164,7 +164,7 @@ class BaseObjFunction(object):
 
         R = self.data.residualWeighted(m, u=u)
 
-        dmisfit = self.data.prob.Jt(m, self.data.Wd * R, u=u)
+        dmisfit = self.data.prob.Jtvec(m, self.data.Wd * R, u=u)
 
         return dmisfit
 
@@ -209,7 +209,7 @@ class BaseObjFunction(object):
         R = self.data.residualWeighted(m, u=u)
 
         # TODO: abstract to different norms a little cleaner.
-        #                                                 \/ it goes here. in l2 it is the identity.
-        dmisfit = self.data.prob.Jt_approx(m, self.data.Wd * self.data.Wd * self.data.prob.J_approx(m, v, u=u), u=u)
+        #                                                     \/ it goes here. in l2 it is the identity.
+        dmisfit = self.data.prob.Jtvec_approx(m, self.data.Wd * self.data.Wd * self.data.prob.Jvec_approx(m, v, u=u), u=u)
 
         return dmisfit

SimPEG/Problem.py

@@ -1,5 +1,5 @@
 import Utils, Data, numpy as np, scipy.sparse as sp
-
+import Model
 
 class BaseProblem(object):
     """
@@ -39,10 +39,12 @@ class BaseProblem(object):
     counter = None   #: A SimPEG.Utils.Counter object
 
     dataPair = Data.BaseData
+    modelPair = Model.BaseModel
 
-    def __init__(self, mesh, model, *args, **kwargs):
+    def __init__(self, mesh, model, **kwargs):
         Utils.setKwargs(self, **kwargs)
         self.mesh = mesh
+        assert isinstance(model, self.modelPair), "Model object must be an instance of a %s class."%(self.modelPair.__name__)
         self.model = model
 
     @property
@@ -70,7 +72,7 @@ class BaseProblem(object):
     def ispaired(self): return self.data is not None
 
     @Utils.timeIt
-    def J(self, m, v, u=None):
+    def Jvec(self, m, v, u=None):
         """
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
@@ -100,7 +102,7 @@ class BaseProblem(object):
         raise NotImplementedError('J is not yet implemented.')
 
     @Utils.timeIt
-    def Jt(self, m, v, u=None):
+    def Jtvec(self, m, v, u=None):
         """
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
@@ -114,7 +116,7 @@ class BaseProblem(object):
 
 
     @Utils.timeIt
-    def J_approx(self, m, v, u=None):
+    def Jvec_approx(self, m, v, u=None):
         """
 
             :param numpy.array m: model
@@ -126,10 +128,10 @@ class BaseProblem(object):
             Approximate effect of J on a vector v
 
         """
-        return self.J(m, v, u)
+        return self.Jvec(m, v, u)
 
     @Utils.timeIt
-    def Jt_approx(self, m, v, u=None):
+    def Jtvec_approx(self, m, v, u=None):
         """
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
@@ -140,7 +142,7 @@ class BaseProblem(object):
             Approximate transpose of J*v
 
         """
-        return self.Jt(m, v, u)
+        return self.Jtvec(m, v, u)
 
     def fields(self, m):
         """

SimPEG/Regularization.py

@@ -59,7 +59,7 @@ class BaseRegularization(object):
 
     @Utils.timeIt
     def modelObj(self, m):
-        r = self.W * (m - self.mref)
+        r = self.W * self.model.transform(m - self.mref)
         return 0.5*r.dot(r)
 
     @Utils.timeIt
@@ -79,7 +79,7 @@ class BaseRegularization(object):
             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
         """
-        return self.W.T * ( self.W * (m - self.mref) )
+        return self.W.T * ( self.W * self.model.transform(m - self.mref) )
 
     @Utils.timeIt
     def modelObj2Deriv(self):
@@ -207,7 +207,7 @@ class Tikhonov(BaseRegularization):
     def Wx(self):
         """Regularization matrix Wx"""
         if getattr(self, '_Wx', None) is None:
-            Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFv[0]].T*self.mesh.vol
+            Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol
             self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx
         return self._Wx
 
@@ -215,7 +215,7 @@ class Tikhonov(BaseRegularization):
     def Wy(self):
         """Regularization matrix Wy"""
         if getattr(self, '_Wy', None) is None:
-            Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFv[0]:np.sum(self.mesh.nFv[:2])].T*self.mesh.vol
+            Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol
             self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady
         return self._Wy
 
@@ -223,7 +223,7 @@ class Tikhonov(BaseRegularization):
     def Wz(self):
         """Regularization matrix Wz"""
         if getattr(self, '_Wz', None) is None:
-            Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.nFv[:2]):].T*self.mesh.vol
+            Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol
             self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz
         return self._Wz
 

SimPEG/Solver.py

@@ -1,8 +1,7 @@
 import numpy as np
 import scipy.sparse as sp
 import scipy.sparse.linalg as linalg
-from Utils.matutils import mkvc
-from Utils.sputils import sdiag
+from Utils.matutils import mkvc, sdiag
 import warnings
 
 DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'triangular':'fortran', 'diagonal':'python'}

SimPEG/Tests/TestUtils.py

@@ -93,9 +93,9 @@ class OrderTest(unittest.TestCase):
                 h3 = np.ones(nc)/nc
                 h = [h1, h2, h3]
             elif 'random' in self._meshType:
-                h1 = np.random.rand(nc)
-                h2 = np.random.rand(nc)
-                h3 = np.random.rand(nc)
+                h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
+                h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
+                h3 = np.random.rand(nc)*nc*0.5 + nc*0.5
                 h = [hi/np.sum(hi) for hi in [h1, h2, h3]]  # normalize
             else:
                 raise Exception('Unexpected meshType')
@@ -111,10 +111,12 @@ class OrderTest(unittest.TestCase):
                 kwrd = 'rotate'
             else:
                 raise Exception('Unexpected meshType')
-            if self.meshDimension == 2:
+            if self.meshDimension == 1:
+                raise Exception('Lom not supported for 1D')
+            elif self.meshDimension == 2:
                 X, Y = Utils.exampleLomGird([nc, nc], kwrd)
                 self.M = LogicallyOrthogonalMesh([X, Y])
-            if self.meshDimension == 3:
+            elif self.meshDimension == 3:
                 X, Y, Z = Utils.exampleLomGird([nc, nc, nc], kwrd)
                 self.M = LogicallyOrthogonalMesh([X, Y, Z])
             return 1./nc

SimPEG/Tests/test_LogicallyOrthogonalMesh.py

@@ -44,43 +44,43 @@ class BasicLOMTests(unittest.TestCase):
 
     def test_tangents(self):
         T = self.LOM2.tangents
-        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM2.nEv[0])))
-        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM2.nEv[0])))
-        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM2.nEv[1])))
-        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM2.nEv[1])))
+        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM2.nEx)))
+        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM2.nEx)))
+        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM2.nEy)))
+        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM2.nEy)))
 
         T = self.LOM3.tangents
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM3.nEv[0])))
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM3.nEv[0])))
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LOM3.nEv[0])))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM3.nEx)))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM3.nEx)))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LOM3.nEx)))
 
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM3.nEv[1])))
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM3.nEv[1])))
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LOM3.nEv[1])))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM3.nEy)))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM3.nEy)))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LOM3.nEy)))
 
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LOM3.nEv[2])))
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LOM3.nEv[2])))
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LOM3.nEv[2])))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LOM3.nEz)))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LOM3.nEz)))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LOM3.nEz)))
 
     def test_normals(self):
         N = self.LOM2.normals
-        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM2.nFv[0])))
-        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM2.nFv[0])))
-        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM2.nFv[1])))
-        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM2.nFv[1])))
+        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM2.nFx)))
+        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM2.nFx)))
+        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM2.nFy)))
+        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM2.nFy)))
 
         N = self.LOM3.normals
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM3.nFv[0])))
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM3.nFv[0])))
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LOM3.nFv[0])))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM3.nFx)))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM3.nFx)))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LOM3.nFx)))
 
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM3.nFv[1])))
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM3.nFv[1])))
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LOM3.nFv[1])))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM3.nFy)))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM3.nFy)))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LOM3.nFy)))
 
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LOM3.nFv[2])))
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LOM3.nFv[2])))
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LOM3.nFv[2])))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LOM3.nFz)))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LOM3.nFz)))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LOM3.nFz)))
 
     def test_grid(self):
         self.assertTrue(np.all(self.LOM2.gridCC == self.TM2.gridCC))

SimPEG/Tests/test_basemesh.py

@@ -1,54 +1,46 @@
 import unittest
 import sys
-from SimPEG.Mesh import BaseMesh
+from SimPEG.Mesh import BaseRectangularMesh
 import numpy as np
 
 
 class TestBaseMesh(unittest.TestCase):
 
     def setUp(self):
-        self.mesh = BaseMesh([6, 2, 3])
+        self.mesh = BaseRectangularMesh([6, 2, 3])
 
     def test_meshDimensions(self):
         self.assertTrue(self.mesh.dim, 3)
 
     def test_mesh_nc(self):
-        self.assertTrue(np.all(self.mesh.nCv == [6, 2, 3]))
+        self.assertTrue(np.all(self.mesh.vnC == [6, 2, 3]))
 
     def test_mesh_nc_xyz(self):
-        x = np.all(self.mesh.nCx == 6)
-        y = np.all(self.mesh.nCy == 2)
-        z = np.all(self.mesh.nCz == 3)
-
-        self.assertTrue(np.all([x, y, z]))
+        self.assertTrue(np.all(self.mesh.nCx == 6))
+        self.assertTrue(np.all(self.mesh.nCy == 2))
+        self.assertTrue(np.all(self.mesh.nCz == 3))
 
     def test_mesh_nf(self):
-        x = np.all(self.mesh.nFx == [7, 2, 3])
-        y = np.all(self.mesh.nFy == [6, 3, 3])
-        z = np.all(self.mesh.nFz == [6, 2, 4])
-
-        self.assertTrue(np.all([x, y, z]))
+        self.assertTrue(np.all(self.mesh.vnFx == [7, 2, 3]))
+        self.assertTrue(np.all(self.mesh.vnFy == [6, 3, 3]))
+        self.assertTrue(np.all(self.mesh.vnFz == [6, 2, 4]))
 
     def test_mesh_ne(self):
-        x = np.all(self.mesh.nEx == [6, 3, 4])
-        y = np.all(self.mesh.nEy == [7, 2, 4])
-        z = np.all(self.mesh.nEz == [7, 3, 3])
-
-        self.assertTrue(np.all([x, y, z]))
+        self.assertTrue(np.all(self.mesh.vnEx == [6, 3, 4]))
+        self.assertTrue(np.all(self.mesh.vnEy == [7, 2, 4]))
+        self.assertTrue(np.all(self.mesh.vnEz == [7, 3, 3]))
 
     def test_mesh_numbers(self):
-        c = self.mesh.nC == 36
-        fv = np.all(self.mesh.nFv == [42, 54, 48])
-        ev = np.all(self.mesh.nEv == [72, 56, 63])
-        f = np.all(self.mesh.nF == np.sum([42, 54, 48]))
-        e = np.all(self.mesh.nE == np.sum([72, 56, 63]))
-
-        self.assertTrue(np.all([c, fv, ev, f, e]))
+        self.assertTrue(self.mesh.nC == 36)
+        self.assertTrue(np.all(self.mesh.vnF == [42, 54, 48]))
+        self.assertTrue(np.all(self.mesh.vnE == [72, 56, 63]))
+        self.assertTrue(np.all(self.mesh.nF == np.sum([42, 54, 48])))
+        self.assertTrue(np.all(self.mesh.nE == np.sum([72, 56, 63])))
 
     def test_mesh_r_E_V(self):
-        ex = np.ones(self.mesh.nEv[0])
-        ey = np.ones(self.mesh.nEv[1])*2
-        ez = np.ones(self.mesh.nEv[2])*3
+        ex = np.ones(self.mesh.nEx)
+        ey = np.ones(self.mesh.nEy)*2
+        ez = np.ones(self.mesh.nEz)*3
         e = np.r_[ex, ey, ez]
         tex = self.mesh.r(e, 'E', 'Ex', 'V')
         tey = self.mesh.r(e, 'E', 'Ey', 'V')
@@ -62,9 +54,9 @@ class TestBaseMesh(unittest.TestCase):
         self.assertTrue(np.all(tez == ez))
 
     def test_mesh_r_F_V(self):
-        fx = np.ones(self.mesh.nFv[0])
-        fy = np.ones(self.mesh.nFv[1])*2
-        fz = np.ones(self.mesh.nFv[2])*3
+        fx = np.ones(self.mesh.nFx)
+        fy = np.ones(self.mesh.nFy)*2
+        fz = np.ones(self.mesh.nFz)*3
         f = np.r_[fx, fy, fz]
         tfx = self.mesh.r(f, 'F', 'Fx', 'V')
         tfy = self.mesh.r(f, 'F', 'Fy', 'V')
@@ -78,25 +70,25 @@ class TestBaseMesh(unittest.TestCase):
         self.assertTrue(np.all(tfz == fz))
 
     def test_mesh_r_E_M(self):
-        g = np.ones((np.prod(self.mesh.nEx), 3))
+        g = np.ones((np.prod(self.mesh.vnEx), 3))
         g[:, 1] = 2
         g[:, 2] = 3
         Xex, Yex, Zex = self.mesh.r(g, 'Ex', 'Ex', 'M')
-        self.assertTrue(np.all(Xex.shape == self.mesh.nEx))
-        self.assertTrue(np.all(Yex.shape == self.mesh.nEx))
-        self.assertTrue(np.all(Zex.shape == self.mesh.nEx))
+        self.assertTrue(np.all(Xex.shape == self.mesh.vnEx))
+        self.assertTrue(np.all(Yex.shape == self.mesh.vnEx))
+        self.assertTrue(np.all(Zex.shape == self.mesh.vnEx))
         self.assertTrue(np.all(Xex == 1))
         self.assertTrue(np.all(Yex == 2))
         self.assertTrue(np.all(Zex == 3))
 
     def test_mesh_r_F_M(self):
-        g = np.ones((np.prod(self.mesh.nFx), 3))
+        g = np.ones((np.prod(self.mesh.vnFx), 3))
         g[:, 1] = 2
         g[:, 2] = 3
         Xfx, Yfx, Zfx = self.mesh.r(g, 'Fx', 'Fx', 'M')
-        self.assertTrue(np.all(Xfx.shape == self.mesh.nFx))
-        self.assertTrue(np.all(Yfx.shape == self.mesh.nFx))
-        self.assertTrue(np.all(Zfx.shape == self.mesh.nFx))
+        self.assertTrue(np.all(Xfx.shape == self.mesh.vnFx))
+        self.assertTrue(np.all(Yfx.shape == self.mesh.vnFx))
+        self.assertTrue(np.all(Zfx.shape == self.mesh.vnFx))
         self.assertTrue(np.all(Xfx == 1))
         self.assertTrue(np.all(Yfx == 2))
         self.assertTrue(np.all(Zfx == 3))
@@ -106,9 +98,9 @@ class TestBaseMesh(unittest.TestCase):
         g[:, 1] = 2
         g[:, 2] = 3
         Xc, Yc, Zc = self.mesh.r(g, 'CC', 'CC', 'M')
-        self.assertTrue(np.all(Xc.shape == self.mesh.nCv))
-        self.assertTrue(np.all(Yc.shape == self.mesh.nCv))
-        self.assertTrue(np.all(Zc.shape == self.mesh.nCv))
+        self.assertTrue(np.all(Xc.shape == self.mesh.vnC))
+        self.assertTrue(np.all(Yc.shape == self.mesh.vnC))
+        self.assertTrue(np.all(Zc.shape == self.mesh.vnC))
         self.assertTrue(np.all(Xc == 1))
         self.assertTrue(np.all(Yc == 2))
         self.assertTrue(np.all(Zc == 3))
@@ -117,47 +109,44 @@ class TestBaseMesh(unittest.TestCase):
 class TestMeshNumbers2D(unittest.TestCase):
 
     def setUp(self):
-        self.mesh = BaseMesh([6, 2])
+        self.mesh = BaseRectangularMesh([6, 2])
 
     def test_meshDimensions(self):
         self.assertTrue(self.mesh.dim, 2)
 
     def test_mesh_nc(self):
-        self.assertTrue(np.all(self.mesh.nCv == [6, 2]))
+        self.assertTrue(np.all(self.mesh.vnC == [6, 2]))
 
     def test_mesh_nc_xyz(self):
-        x = np.all(self.mesh.nCx == 6)
-        y = np.all(self.mesh.nCy == 2)
-        z = self.mesh.nCz is None
-
-        self.assertTrue(np.all([x, y, z]))
+        self.assertTrue(np.all(self.mesh.nCx == 6))
+        self.assertTrue(np.all(self.mesh.nCy == 2))
+        self.assertTrue(self.mesh.nCz is None)
 
     def test_mesh_nf(self):
-        x = np.all(self.mesh.nFx == [7, 2])
-        y = np.all(self.mesh.nFy == [6, 3])
-        z = self.mesh.nFz is None
-
-        self.assertTrue(np.all([x, y, z]))
+        self.assertTrue(np.all(self.mesh.vnFx == [7, 2]))
+        self.assertTrue(np.all(self.mesh.vnFy == [6, 3]))
+        self.assertTrue(self.mesh.vnFz is None)
 
     def test_mesh_ne(self):
-        x = np.all(self.mesh.nEx == [6, 3])
-        y = np.all(self.mesh.nEy == [7, 2])
-        z = self.mesh.nEz is None
-
-        self.assertTrue(np.all([x, y, z]))
+        self.assertTrue(np.all(self.mesh.vnEx == [6, 3]))
+        self.assertTrue(np.all(self.mesh.vnEy == [7, 2]))
+        self.assertTrue(self.mesh.vnEz is None)
 
     def test_mesh_numbers(self):
         c = self.mesh.nC == 12
-        fv = np.all(self.mesh.nFv == [14, 18])
-        ev = np.all(self.mesh.nEv == [18, 14])
-        f = np.all(self.mesh.nF == np.sum([14, 18]))
-        e = np.all(self.mesh.nE == np.sum([18, 14]))
-
-        self.assertTrue(np.all([c, fv, ev, f, e]))
+        self.assertTrue(np.all(self.mesh.vnF == [14, 18]))
+        self.assertTrue(np.all(self.mesh.nFx == 14))
+        self.assertTrue(np.all(self.mesh.nFy == 18))
+        self.assertTrue(np.all(self.mesh.nEx == 18))
+        self.assertTrue(np.all(self.mesh.nEy == 14))
+        self.assertTrue(np.all(self.mesh.vnE == [18, 14]))
+        self.assertTrue(np.all(self.mesh.vnE == [18, 14]))
+        self.assertTrue(np.all(self.mesh.nF == np.sum([14, 18])))
+        self.assertTrue(np.all(self.mesh.nE == np.sum([18, 14])))
 
     def test_mesh_r_E_V(self):
-        ex = np.ones(self.mesh.nEv[0])
-        ey = np.ones(self.mesh.nEv[1])*2
+        ex = np.ones(self.mesh.nEx)
+        ey = np.ones(self.mesh.nEy)*2
         e = np.r_[ex, ey]
         tex = self.mesh.r(e, 'E', 'Ex', 'V')
         tey = self.mesh.r(e, 'E', 'Ey', 'V')
@@ -169,8 +158,8 @@ class TestMeshNumbers2D(unittest.TestCase):
         self.assertRaises(AssertionError,   self.mesh.r, e, 'E', 'Ez', 'V')
 
     def test_mesh_r_F_V(self):
-        fx = np.ones(self.mesh.nFv[0])
-        fy = np.ones(self.mesh.nFv[1])*2
+        fx = np.ones(self.mesh.nFx)
+        fy = np.ones(self.mesh.nFy)*2
         f = np.r_[fx, fy]
         tfx = self.mesh.r(f, 'F', 'Fx', 'V')
         tfy = self.mesh.r(f, 'F', 'Fy', 'V')
@@ -182,20 +171,20 @@ class TestMeshNumbers2D(unittest.TestCase):
         self.assertRaises(AssertionError,   self.mesh.r, f, 'F', 'Fz', 'V')
 
     def test_mesh_r_E_M(self):
-        g = np.ones((np.prod(self.mesh.nEx), 2))
+        g = np.ones((np.prod(self.mesh.vnEx), 2))
         g[:, 1] = 2
         Xex, Yex = self.mesh.r(g, 'Ex', 'Ex', 'M')
-        self.assertTrue(np.all(Xex.shape == self.mesh.nEx))
-        self.assertTrue(np.all(Yex.shape == self.mesh.nEx))
+        self.assertTrue(np.all(Xex.shape == self.mesh.vnEx))
+        self.assertTrue(np.all(Yex.shape == self.mesh.vnEx))
         self.assertTrue(np.all(Xex == 1))
         self.assertTrue(np.all(Yex == 2))
 
     def test_mesh_r_F_M(self):
-        g = np.ones((np.prod(self.mesh.nFx), 2))
+        g = np.ones((np.prod(self.mesh.vnFx), 2))
         g[:, 1] = 2
         Xfx, Yfx = self.mesh.r(g, 'Fx', 'Fx', 'M')
-        self.assertTrue(np.all(Xfx.shape == self.mesh.nFx))
-        self.assertTrue(np.all(Yfx.shape == self.mesh.nFx))
+        self.assertTrue(np.all(Xfx.shape == self.mesh.vnFx))
+        self.assertTrue(np.all(Yfx.shape == self.mesh.vnFx))
         self.assertTrue(np.all(Xfx == 1))
         self.assertTrue(np.all(Yfx == 2))
 
@@ -203,8 +192,8 @@ class TestMeshNumbers2D(unittest.TestCase):
         g = np.ones((self.mesh.nC, 2))
         g[:, 1] = 2
         Xc, Yc = self.mesh.r(g, 'CC', 'CC', 'M')
-        self.assertTrue(np.all(Xc.shape == self.mesh.nCv))
-        self.assertTrue(np.all(Yc.shape == self.mesh.nCv))
+        self.assertTrue(np.all(Xc.shape == self.mesh.vnC))
+        self.assertTrue(np.all(Yc.shape == self.mesh.vnC))
         self.assertTrue(np.all(Xc == 1))
         self.assertTrue(np.all(Yc == 2))
 

SimPEG/Tests/test_boundaryPoisson.py

@@ -0,0 +1,501 @@
+import numpy as np
+import scipy.sparse as sp
+import unittest
+from TestUtils import OrderTest
+import matplotlib.pyplot as plt
+from SimPEG import Utils, Solver
+
+MESHTYPES = ['uniformTensorMesh']
+
+class Test1D_InhomogeneousDirichlet(OrderTest):
+    name = "1D - Dirichlet"
+    meshTypes = MESHTYPES
+    meshDimension = 1
+    expectedOrders = 2
+    meshSizes = [4, 8, 16, 32, 64, 128]
+
+    def getError(self):
+        #Test function
+        phi = lambda x: np.cos(np.pi*x)
+        j_fun = lambda x: -np.pi*np.sin(np.pi*x)
+        q_fun = lambda x: -(np.pi**2)*np.cos(np.pi*x)
+
+        xc_anal = phi(self.M.gridCC)
+        q_anal = q_fun(self.M.gridCC)
+        j_anal = j_fun(self.M.gridFx)
+
+        #TODO: Check where our boundary conditions are CCx or Nx
+        # vec = self.M.vectorNx
+        vec = self.M.vectorCCx
+
+        phi_bc = phi(vec[[0,-1]])
+        j_bc = j_fun(vec[[0,-1]])
+
+        P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'dirichlet']])
+
+        Mc = self.M.getFaceInnerProduct()
+        McI = Utils.sdInv(self.M.getFaceInnerProduct())
+        V = Utils.sdiag(self.M.vol)
+        G = -Pin.T*Pin*self.M.faceDiv.T * V
+        D = self.M.faceDiv
+        j = McI*(G*xc_anal + P*phi_bc)
+        q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
+
+        # Rearrange if we know q to solve for x
+        A = V*D*Pin.T*Pin*McI*G
+        rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
+        # A = D*McI*G
+        # rhs = q_anal - D*McI*P*phi_bc
+
+
+        if self.myTest == 'j':
+            err = np.linalg.norm((j-j_anal), np.inf)
+        elif self.myTest == 'q':
+            err = np.linalg.norm((q-V*q_anal), np.inf)
+        elif self.myTest == 'xc':
+            #TODO: fix the null space
+            solver = Solver(A, doDirect=False, options={'M':'J','iterSolver':'CG','backend':'scipy','maxIter':1000})
+            xc = solver.solve(rhs)
+            print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)
+            err = np.linalg.norm((xc-xc_anal), np.inf)
+        elif self.myTest == 'xcJ':
+            #TODO: fix the null space
+            xc = Solver(A).solve(rhs)
+            print np.linalg.norm(Utils.mkvc(A*xc) - rhs)
+            j = McI*(G*xc + P*phi_bc)
+            err = np.linalg.norm((j-j_anal), np.inf)
+
+        return err
+
+    def test_orderJ(self):
+        self.name = "1D - InhomogeneousDirichlet_Forward j"
+        self.myTest = 'j'
+        self.orderTest()
+
+    def test_orderQ(self):
+        self.name = "1D - InhomogeneousDirichlet_Forward q"
+        self.myTest = 'q'
+        self.orderTest()
+
+    def test_orderX(self):
+        self.name = "1D - InhomogeneousDirichlet_Inverse"
+        self.myTest = 'xc'
+        self.orderTest()
+
+    def test_orderXJ(self):
+        self.name = "1D - InhomogeneousDirichlet_Inverse J"
+        self.myTest = 'xcJ'
+        self.orderTest()
+
+
+class Test2D_InhomogeneousDirichlet(OrderTest):
+    name = "2D - Dirichlet"
+    meshTypes = MESHTYPES
+    meshDimension = 2
+    expectedOrders = 2
+    meshSizes = [4, 8, 16, 32]
+
+    def getError(self):
+        #Test function
+        phi = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
+        j_funX = lambda x: -np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
+        j_funY = lambda x: -np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
+        q_fun = lambda x: -2*(np.pi**2)*phi(x)
+
+        xc_anal = phi(self.M.gridCC)
+        q_anal = q_fun(self.M.gridCC)
+        jX_anal = j_funX(self.M.gridFx)
+        jY_anal = j_funY(self.M.gridFy)
+        j_anal = np.r_[jX_anal,jY_anal]
+
+        #TODO: Check where our boundary conditions are CCx or Nx
+        # fxm,fxp,fym,fyp = self.M.faceBoundaryInd
+        # gBFx = self.M.gridFx[(fxm|fxp),:]
+        # gBFy = self.M.gridFy[(fym|fyp),:]
+        fxm,fxp,fym,fyp = self.M.cellBoundaryInd
+        gBFx = self.M.gridCC[(fxm|fxp),:]
+        gBFy = self.M.gridCC[(fym|fyp),:]
+
+        bc = phi(np.r_[gBFx,gBFy])
+
+        # P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2))
+
+        P, Pin, Pout = self.M.getBCProjWF('dirichlet')
+
+        Mc = self.M.getFaceInnerProduct()
+        McI = Utils.sdInv(self.M.getFaceInnerProduct())
+        G = -self.M.faceDiv.T * Utils.sdiag(self.M.vol)
+        D = self.M.faceDiv
+        j = McI*(G*xc_anal + P*bc)
+        q = D*j
+
+        # self.M.plotImage(j, 'FxFy', showIt=True)
+
+        # Rearrange if we know q to solve for x
+        A = D*McI*G
+        rhs = q_anal - D*McI*P*bc
+
+        if self.myTest == 'j':
+            err = np.linalg.norm((j-j_anal), np.inf)
+        elif self.myTest == 'q':
+            err = np.linalg.norm((q-q_anal), np.inf)
+        elif self.myTest == 'xc':
+            xc = Solver(A).solve(rhs)
+            err = np.linalg.norm((xc-xc_anal), np.inf)
+        elif self.myTest == 'xcJ':
+            xc = Solver(A).solve(rhs)
+            j = McI*(G*xc + P*bc)
+            err = np.linalg.norm((j-j_anal), np.inf)
+
+        return err
+
+    def test_orderJ(self):
+        self.name = "2D - InhomogeneousDirichlet_Forward j"
+        self.myTest = 'j'
+        self.orderTest()
+
+    def test_orderQ(self):
+        self.name = "2D - InhomogeneousDirichlet_Forward q"
+        self.myTest = 'q'
+        self.orderTest()
+
+    def test_orderX(self):
+        self.name = "2D - InhomogeneousDirichlet_Inverse"
+        self.myTest = 'xc'
+        self.orderTest()
+
+    def test_orderXJ(self):
+        self.name = "2D - InhomogeneousDirichlet_Inverse J"
+        self.myTest = 'xcJ'
+        self.orderTest()
+
+class Test1D_InhomogeneousNeumann(OrderTest):
+    name = "1D - Neumann"
+    meshTypes = MESHTYPES
+    meshDimension = 1
+    expectedOrders = 2
+    meshSizes = [4, 8, 16, 32, 64, 128]
+
+    def getError(self):
+        #Test function
+        phi = lambda x: np.sin(np.pi*x)
+        j_fun = lambda x: np.pi*np.cos(np.pi*x)
+        q_fun = lambda x: -(np.pi**2)*np.sin(np.pi*x)
+
+        xc_anal = phi(self.M.gridCC)
+        q_anal = q_fun(self.M.gridCC)
+        j_anal = j_fun(self.M.gridFx)
+
+        #TODO: Check where our boundary conditions are CCx or Nx
+        vecN = self.M.vectorNx
+        vecC = self.M.vectorCCx
+
+        phi_bc = phi(vecC[[0,-1]])
+        j_bc = j_fun(vecN[[0,-1]])
+
+        P, Pin, Pout = self.M.getBCProjWF([['neumann', 'neumann']])
+
+        Mc = self.M.getFaceInnerProduct()
+        McI = Utils.sdInv(self.M.getFaceInnerProduct())
+        V = Utils.sdiag(self.M.vol)
+        G = -Pin.T*Pin*self.M.faceDiv.T * V
+        D = self.M.faceDiv
+        j = McI*(G*xc_anal + P*phi_bc)
+        q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
+
+        # Rearrange if we know q to solve for x
+        A = V*D*Pin.T*Pin*McI*G
+        rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
+        # A = D*McI*G
+        # rhs = q_anal - D*McI*P*phi_bc
+
+
+        if self.myTest == 'j':
+            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+        elif self.myTest == 'q':
+            err = np.linalg.norm((q-V*q_anal), np.inf)
+        elif self.myTest == 'xc':
+            #TODO: fix the null space
+            xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
+            err = np.linalg.norm((xc-xc_anal), np.inf)
+            if info > 0:
+                print 'Solve does not work well'
+                print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)            
+        elif self.myTest == 'xcJ':
+            #TODO: fix the null space
+            xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
+            j = McI*(G*xc + P*phi_bc)
+            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+            if info > 0:
+                print 'Solve does not work well'
+                print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)                        
+        return err
+
+    def test_orderJ(self):
+        self.name = "1D - InhomogeneousNeumann_Forward j"
+        self.myTest = 'j'
+        self.orderTest()
+
+    def test_orderQ(self):
+        self.name = "1D - InhomogeneousNeumann_Forward q"
+        self.myTest = 'q'
+        self.orderTest()
+
+    def test_orderXJ(self):
+        self.name = "1D - InhomogeneousNeumann_Inverse J"
+        self.myTest = 'xcJ'
+        self.orderTest()
+
+class Test2D_InhomogeneousNeumann(OrderTest):
+    name = "2D - Neumann"
+    meshTypes = MESHTYPES
+    meshDimension = 2
+    expectedOrders = 2
+    meshSizes = [4, 8, 16, 32]
+    # meshSizes = [4]
+
+    def getError(self):
+        #Test function
+        phi = lambda x: np.sin(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
+        j_funX = lambda x: np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
+        j_funY = lambda x: np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
+        q_fun = lambda x: -2*(np.pi**2)*phi(x)
+
+        xc_anal = phi(self.M.gridCC)
+        q_anal = q_fun(self.M.gridCC)
+        jX_anal = j_funX(self.M.gridFx)
+        jY_anal = j_funY(self.M.gridFy)
+        j_anal = np.r_[jX_anal,jY_anal]
+
+        #TODO: Check where our boundary conditions are CCx or Nx
+        
+        cxm,cxp,cym,cyp = self.M.cellBoundaryInd
+        fxm,fxp,fym,fyp = self.M.faceBoundaryInd
+
+        gBFx = self.M.gridFx[(fxm|fxp),:]
+        gBFy = self.M.gridFy[(fym|fyp),:]        
+
+        gBCx = self.M.gridCC[(cxm|cxp),:]
+        gBCy = self.M.gridCC[(cym|cyp),:]
+
+        phi_bc = phi(np.r_[gBFx,gBFy])
+        j_bc = np.r_[j_funX(gBFx), j_funY(gBFy)]
+
+        # P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2))
+
+        P, Pin, Pout = self.M.getBCProjWF('neumann')
+
+        Mc = self.M.getFaceInnerProduct()
+        McI = Utils.sdInv(self.M.getFaceInnerProduct())
+        V = Utils.sdiag(self.M.vol)
+        G = -Pin.T*Pin*self.M.faceDiv.T * V
+        D = self.M.faceDiv
+        j = McI*(G*xc_anal + P*phi_bc)
+        q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
+
+        # Rearrange if we know q to solve for x
+        A = V*D*Pin.T*Pin*McI*G
+        rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
+
+        if self.myTest == 'j':
+            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+        elif self.myTest == 'q':
+            err = np.linalg.norm((q-V*q_anal), np.inf)
+        elif self.myTest == 'xc':
+            #TODO: fix the null space
+            xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
+            err = np.linalg.norm((xc-xc_anal), np.inf)
+            if info > 0:
+                print 'Solve does not work well'
+                print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)            
+        elif self.myTest == 'xcJ':
+            #TODO: fix the null space
+            xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
+            j = McI*(G*xc + P*phi_bc)
+            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+            if info > 0:
+                print 'Solve does not work well'
+                print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)                        
+        return err
+
+    def test_orderJ(self):
+        self.name = "2D - InhomogeneousNeumann_Forward j"
+        self.myTest = 'j'
+        self.orderTest()
+
+    def test_orderQ(self):
+        self.name = "2D - InhomogeneousNeumann_Forward q"
+        self.myTest = 'q'
+        self.orderTest()
+
+    def test_orderXJ(self):
+        self.name = "2D - InhomogeneousNeumann_Inverse J"
+        self.myTest = 'xcJ'
+        self.orderTest()
+
+class Test1D_InhomogeneousMixed(OrderTest):
+    name = "1D - Mixed"
+    meshTypes = MESHTYPES
+    meshDimension = 1
+    expectedOrders = 2
+    meshSizes = [4, 8, 16, 32, 64, 128]
+
+    def getError(self):
+        #Test function
+        phi = lambda x: np.cos(0.5*np.pi*x)
+        j_fun = lambda x: -0.5*np.pi*np.sin(0.5*np.pi*x)
+        q_fun = lambda x: -0.25*(np.pi**2)*np.cos(0.5*np.pi*x)
+
+        xc_anal = phi(self.M.gridCC)
+        q_anal = q_fun(self.M.gridCC)
+        j_anal = j_fun(self.M.gridFx)
+
+        #TODO: Check where our boundary conditions are CCx or Nx
+        vecN = self.M.vectorNx
+        vecC = self.M.vectorCCx
+
+        phi_bc = phi(vecC[[0,-1]])
+        j_bc = j_fun(vecN[[0,-1]])
+
+        P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'neumann']])
+
+        Mc = self.M.getFaceInnerProduct()
+        McI = Utils.sdInv(self.M.getFaceInnerProduct())
+        V = Utils.sdiag(self.M.vol)
+        G = -Pin.T*Pin*self.M.faceDiv.T * V
+        D = self.M.faceDiv
+        j = McI*(G*xc_anal + P*phi_bc)
+        q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
+
+        # Rearrange if we know q to solve for x
+        A = V*D*Pin.T*Pin*McI*G
+        rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
+        # A = D*McI*G
+        # rhs = q_anal - D*McI*P*phi_bc
+
+
+        if self.myTest == 'j':
+            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+        elif self.myTest == 'q':
+            err = np.linalg.norm((q-V*q_anal), np.inf)
+        elif self.myTest == 'xc':
+            #TODO: fix the null space
+            xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
+            err = np.linalg.norm((xc-xc_anal), np.inf)
+            if info > 0:
+                print 'Solve does not work well'
+                print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)            
+        elif self.myTest == 'xcJ':
+            #TODO: fix the null space
+            xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
+            j = McI*(G*xc + P*phi_bc)
+            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+            if info > 0:
+                print 'Solve does not work well'
+                print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)                        
+        return err
+
+    def test_orderJ(self):
+        self.name = "1D - InhomogeneousMixed_Forward j"
+        self.myTest = 'j'
+        self.orderTest()
+
+    def test_orderQ(self):
+        self.name = "1D - InhomogeneousMixed_Forward q"
+        self.myTest = 'q'
+        self.orderTest()
+
+    def test_orderXJ(self):
+        self.name = "1D - InhomogeneousMixed_Inverse J"
+        self.myTest = 'xcJ'
+        self.orderTest()
+
+class Test2D_InhomogeneousMixed(OrderTest):
+    name = "2D - Mixed"
+    meshTypes = MESHTYPES
+    meshDimension = 2
+    expectedOrders = 2
+    meshSizes = [2, 4, 8, 16]
+    # meshSizes = [4]
+
+    def getError(self):
+        #Test function
+        phi = lambda x: np.cos(0.5*np.pi*x[:,0])*np.cos(0.5*np.pi*x[:,1])
+        j_funX = lambda x: -0.5*np.pi*np.sin(0.5*np.pi*x[:,0])*np.cos(0.5*np.pi*x[:,1])
+        j_funY = lambda x: -0.5*np.pi*np.cos(0.5*np.pi*x[:,0])*np.sin(0.5*np.pi*x[:,1])
+        q_fun = lambda x: -2*((0.5*np.pi)**2)*phi(x)
+
+        xc_anal = phi(self.M.gridCC)
+        q_anal = q_fun(self.M.gridCC)
+        jX_anal = j_funX(self.M.gridFx)
+        jY_anal = j_funY(self.M.gridFy)
+        j_anal = np.r_[jX_anal,jY_anal]
+
+        #TODO: Check where our boundary conditions are CCx or Nx
+        
+        cxm,cxp,cym,cyp = self.M.cellBoundaryInd
+        fxm,fxp,fym,fyp = self.M.faceBoundaryInd
+
+        gBFx = self.M.gridFx[(fxm|fxp),:]
+        gBFy = self.M.gridFy[(fym|fyp),:]        
+
+        gBCx = self.M.gridCC[(cxm|cxp),:]
+        gBCy = self.M.gridCC[(cym|cyp),:]
+
+        phi_bc = phi(np.r_[gBCx,gBCy])
+        j_bc = np.r_[j_funX(gBFx), j_funY(gBFy)]
+
+        # P = sp.csr_matrix(([-1,1],([0,self.M.nF-1],[0,1])), shape=(self.M.nF, 2))
+
+        P, Pin, Pout = self.M.getBCProjWF([['dirichlet', 'neumann'], ['dirichlet', 'neumann']])
+
+        Mc = self.M.getFaceInnerProduct()
+        McI = Utils.sdInv(self.M.getFaceInnerProduct())
+        V = Utils.sdiag(self.M.vol)
+        G = -Pin.T*Pin*self.M.faceDiv.T * V
+        D = self.M.faceDiv
+        j = McI*(G*xc_anal + P*phi_bc)
+        q = V*D*Pin.T*Pin*j + V*D*Pout.T*j_bc
+
+        # Rearrange if we know q to solve for x
+        A = V*D*Pin.T*Pin*McI*G
+        rhs = V*q_anal - V*D*Pin.T*Pin*McI*P*phi_bc - V*D*Pout.T*j_bc
+
+        if self.myTest == 'j':
+            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+        elif self.myTest == 'q':
+            err = np.linalg.norm((q-V*q_anal), np.inf)
+        elif self.myTest == 'xc':
+            #TODO: fix the null space
+            xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
+            err = np.linalg.norm((xc-xc_anal), np.inf)
+            if info > 0:
+                print 'Solve does not work well'
+                print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)            
+        elif self.myTest == 'xcJ':
+            #TODO: fix the null space
+            xc, info = sp.linalg.minres(A, rhs, tol = 1e-6)
+            j = McI*(G*xc + P*phi_bc)
+            err = np.linalg.norm((Pin*j-Pin*j_anal), np.inf)
+            if info > 0:
+                print 'Solve does not work well'
+                print 'ACCURACY', np.linalg.norm(Utils.mkvc(A*xc) - rhs)                        
+        return err
+
+    def test_orderJ(self):
+        self.name = "2D - InhomogeneousMixed_Forward j"
+        self.myTest = 'j'
+        self.orderTest()
+
+    def test_orderQ(self):
+        self.name = "2D - InhomogeneousMixed_Forward q"
+        self.myTest = 'q'
+        self.orderTest()
+
+    def test_orderXJ(self):
+        self.name = "2D - InhomogeneousMixed_Inverse J"
+        self.myTest = 'xcJ'
+        self.orderTest()
+
+if __name__ == '__main__':
+    unittest.main()

SimPEG/Tests/test_interpolation.py

@@ -2,6 +2,9 @@ import numpy as np
 import unittest
 from TestUtils import OrderTest
 from SimPEG.Utils import mkvc
+from SimPEG import Mesh
+import unittest
+
 
 MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh']
 TOLERANCES = [0.9, 0.5]
@@ -50,6 +53,19 @@ class TestInterpolation1D(OrderTest):
         self.name = 'Interpolation 1D: N'
         self.orderTest()
 
+class TestOutliersInterp1D(unittest.TestCase):
+
+    def setUp(self):
+        pass
+
+    def test_outliers(self):
+        M = Mesh.TensorMesh([4])
+        Q = M.getInterpolationMat(np.array([[0],[0.126],[0.127]]),'CC',zerosOutside=True)
+        x = np.arange(4)+1
+        self.assertTrue(np.all(Q*x == [1,1.004,1.008]))
+        Q = M.getInterpolationMat(np.array([[-1],[0.126],[0.127]]),'CC',zerosOutside=True)
+        self.assertTrue(np.all(Q*x == [0,1.004,1.008]))
+
 class TestInterpolation2d(OrderTest):
     name = "Interpolation 2D"
     LOCS = np.random.rand(50,2)*0.6+0.2

SimPEG/Tests/test_massMatrices.py

@@ -3,32 +3,6 @@ import unittest
 from TestUtils import OrderTest
 
 
-# MATLAB code:
-
-# syms x y z
-
-# ex = x.^2+y.*z;
-# ey = (z.^2).*x+y.*z;
-# ez = y.^2+x.*z;
-
-# e = [ex;ey;ez];
-
-# 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;
-
-# S1 = [sigma1,0,0;0,sigma1,0;0,0,sigma1];
-# S2 = [sigma1,0,0;0,sigma2,0;0,0,sigma3];
-# S3 = [sigma1,sigma4,sigma5;sigma4,sigma2,sigma6;sigma5,sigma6,sigma3];
-
-# i1 = int(int(int(e.'*S1*e,x,0,1),y,0,1),z,0,1);
-# i2 = int(int(int(e.'*S2*e,x,0,1),y,0,1),z,0,1);
-# i3 = int(int(int(e.'*S3*e,x,0,1),y,0,1),z,0,1);
-
-
 class TestInnerProducts(OrderTest):
     """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
 
@@ -54,14 +28,14 @@ class TestInnerProducts(OrderTest):
         Gc = self.M.gridCC
         if self.sigmaTest == 1:
             sigma = np.c_[call(sigma1, Gc)]
-            analytic = 647./360  # Found using matlab symbolic toolbox.
+            analytic = 647./360  # Found using sympy.
         elif self.sigmaTest == 3:
             sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
-            analytic = 37./12  # Found using matlab symbolic toolbox.
+            analytic = 37./12  # Found using sympy.
         elif self.sigmaTest == 6:
             sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc),
                           call(sigma4, Gc), call(sigma5, Gc), call(sigma6, Gc)]
-            analytic = 69881./21600  # Found using matlab symbolic toolbox.
+            analytic = 69881./21600  # Found using sympy.
 
         if self.location == 'edges':
             cart = lambda g: np.c_[call(ex, g), call(ey, g), call(ez, g)]
@@ -143,13 +117,13 @@ class TestInnerProducts2D(OrderTest):
         Gc = self.M.gridCC
         if self.sigmaTest == 1:
             sigma = np.c_[call(sigma1, Gc)]
-            analytic = 144877./360  # Found using matlab symbolic toolbox. z=5
+            analytic = 144877./360  # Found using sympy. z=5
         elif self.sigmaTest == 2:
             sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)]
-            analytic = 189959./120  # Found using matlab symbolic toolbox. z=5
+            analytic = 189959./120  # Found using sympy. z=5
         elif self.sigmaTest == 3:
             sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
-            analytic = 781427./360  # Found using matlab symbolic toolbox. z=5
+            analytic = 781427./360  # Found using sympy. z=5
 
         if self.location == 'edges':
             cart = lambda g: np.c_[call(ex, g), call(ey, g)]
@@ -206,5 +180,97 @@ class TestInnerProducts2D(OrderTest):
         self.orderTest()
 
 
+
+class TestInnerProducts1D(OrderTest):
+    """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
+
+    meshTypes = ['uniformTensorMesh']
+    meshDimension = 1
+    meshSizes = [4, 8, 16, 32, 64, 128]
+
+    def getError(self):
+
+        y = 12 # Because 12 is just such a great number.
+        z = 5  # Because 5 is just such a great number as well!
+
+        call = lambda fun, x: fun(x)
+
+        ex = lambda x: x**2+y*z
+
+        sigma1 = lambda x: x*y+1
+
+        Gc = self.M.gridCC
+        sigma = call(sigma1, Gc)
+        analytic = 128011./5  # Found using sympy. y=12, z=5
+
+        if self.location == 'faces':
+            F = call(ex, self.M.gridFx)
+            A = self.M.getFaceInnerProduct(sigma)
+            numeric = F.T.dot(A.dot(F))
+
+        err = np.abs(numeric - analytic)
+        return err
+
+    def test_order1_faces(self):
+        self.name = "1D Face Inner Product"
+        self.location = 'faces'
+        self.sigmaTest = 1
+        self.orderTest()
+
+
 if __name__ == '__main__':
     unittest.main()
+
+if __name__ == '__main__' and False:
+    import sympy
+
+    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])
+
+    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
+
+    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])
+
+    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]])
+
+    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
+
+    sigma1 = x*y+1
+
+    print '1D'
+    print sympy.integrate(e*sigma1*e, (x,0,1))

SimPEG/Tests/test_model.py

@@ -11,7 +11,8 @@ class ModelTests(unittest.TestCase):
 
         a = np.array([1, 1, 1])
         b = np.array([1, 2])
-        self.mesh2 = Mesh.TensorMesh([a, b], np.array([3, 5]))
+        self.mesh2 = Mesh.TensorMesh([a, b], x0=np.array([3, 5]))
+        self.mesh22 = Mesh.TensorMesh([b, a], x0=np.array([3, 5]))
 
     def test_modelTransforms(self):
         for M in dir(Model):
@@ -22,6 +23,10 @@ class ModelTests(unittest.TestCase):
                 continue
             self.assertTrue(model.test())
 
+    def test_Mesh2MeshModel(self):
+        model = Model.Mesh2Mesh([self.mesh22, self.mesh2])
+        self.assertTrue(model.test())
+
     def test_comboModels(self):
         combos = [(Model.LogModel, Model.Vertical1DModel)]
         for combo in combos:

SimPEG/Tests/test_utils.py

@@ -1,8 +1,9 @@
-import numpy as np
 import unittest
-from SimPEG.Utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
+from SimPEG.Utils import *
+from SimPEG import Mesh, np, sp
 from SimPEG.Tests import checkDerivative
 
+TOL = 1e-8
 
 class TestCheckDerivative(unittest.TestCase):
 
@@ -64,6 +65,19 @@ class TestSequenceFunctions(unittest.TestCase):
         self.assertTrue(np.all(XYZ[:, 1] == X2_test))
         self.assertTrue(np.all(XYZ[:, 2] == X3_test))
 
+    def test_sub2ind(self):
+        x = np.ones((5,2))
+        self.assertTrue(np.all(sub2ind(x.shape, [0,0]) == [0]))
+        self.assertTrue(np.all(sub2ind(x.shape, [4,0]) == [4]))
+        self.assertTrue(np.all(sub2ind(x.shape, [0,1]) == [5]))
+        self.assertTrue(np.all(sub2ind(x.shape, [4,1]) == [9]))
+        self.assertTrue(np.all(sub2ind(x.shape, [[0,0],[4,0],[0,1],[4,1]]) == [0,4,5,9]))
+
+    def test_ind2sub(self):
+        x = np.ones((5,2))
+        self.assertTrue(np.all(ind2sub(x.shape, [0,4,5,9])[0] == [0,4,0,4]))
+        self.assertTrue(np.all(ind2sub(x.shape, [0,4,5,9])[1] == [0,0,1,1]))
+
     def test_indexCube_2D(self):
         nN = np.array([3, 3])
         self.assertTrue(np.all(indexCube('A', nN) == np.array([0, 1, 3, 4])))
@@ -83,8 +97,6 @@ class TestSequenceFunctions(unittest.TestCase):
         self.assertTrue(np.all(indexCube('H', nN) == np.array([10, 11, 13, 14, 19, 20, 22, 23])))
 
     def test_invXXXBlockDiagonal(self):
-        import scipy.sparse as sp
-
         a = [np.random.rand(5, 1) for i in range(4)]
 
         B = inv2X2BlockDiagonal(*a)
@@ -93,7 +105,7 @@ class TestSequenceFunctions(unittest.TestCase):
                        sp.hstack((sdiag(a[2]), sdiag(a[3])))))
 
         Z2 = B*A - sp.eye(10, 10)
-        self.assertTrue(np.linalg.norm(Z2.todense().ravel(), 2) < 1e-12)
+        self.assertTrue(np.linalg.norm(Z2.todense().ravel(), 2) < TOL)
 
         a = [np.random.rand(5, 1) for i in range(9)]
         B = inv3X3BlockDiagonal(*a)
@@ -104,9 +116,53 @@ class TestSequenceFunctions(unittest.TestCase):
 
         Z3 = B*A - sp.eye(15, 15)
 
-        self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < 1e-12)
+        self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < TOL)
 
 
+    def test_invPropertyTensor2D(self):
+        M = Mesh.TensorMesh([6, 6])
+        a1 = np.random.rand(M.nC)
+        a2 = np.random.rand(M.nC)
+        a3 = np.random.rand(M.nC)
+        prop1 = a1
+        prop2 = np.c_[a1, a2]
+        prop3 = np.c_[a1, a2, a3]
+
+        for prop in [4, prop1, prop2, prop3]:
+            b = invPropertyTensor(M, prop)
+            A = makePropertyTensor(M, prop)
+            B1 = makePropertyTensor(M, b)
+            B2 = invPropertyTensor(M, prop, returnMatrix=True)
+
+            Z = B1*A - sp.identity(M.nC*2)
+            self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
+            Z = B2*A - sp.identity(M.nC*2)
+            self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
+
+
+    def test_invPropertyTensor3D(self):
+        M = Mesh.TensorMesh([6, 6, 6])
+        a1 = np.random.rand(M.nC)
+        a2 = np.random.rand(M.nC)
+        a3 = np.random.rand(M.nC)
+        a4 = np.random.rand(M.nC)
+        a5 = np.random.rand(M.nC)
+        a6 = np.random.rand(M.nC)
+        prop1 = a1
+        prop2 = np.c_[a1, a2, a3]
+        prop3 = np.c_[a1, a2, a3, a4, a5, a6]
+
+        for prop in [4, prop1, prop2, prop3]:
+            b = invPropertyTensor(M, prop)
+            A = makePropertyTensor(M, prop)
+            B1 = makePropertyTensor(M, b)
+            B2 = invPropertyTensor(M, prop, returnMatrix=True)
+
+            Z = B1*A - sp.identity(M.nC*3)
+            self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
+            Z = B2*A - sp.identity(M.nC*3)
+            self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
+
 
 if __name__ == '__main__':
     unittest.main()

SimPEG/Utils/__init__.py

@@ -1,7 +1,6 @@
-from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
-from sputils import spzeros, kron3, speye, sdiag, sdInv, ddx, av, avExtrap
+from matutils import *
 from meshutils import exampleLomGird, meshTensors
-from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube
+from lomutils import volTetra, faceInfo, indexCube
 from interputils import interpmat
 from ipythonutils import easyAnimate as animate
 import ModelBuilder

SimPEG/Utils/interputils.py

@@ -1,7 +1,6 @@
 import numpy as np
 import scipy.sparse as sp
-from sputils import spzeros
-from matutils import mkvc, sub2ind
+from matutils import mkvc, sub2ind, spzeros
 
 def _interp_point_1D(x, xr_i):
     """
@@ -20,9 +19,17 @@ def _interp_point_1D(x, xr_i):
     elif xr_i - x[im] < 0:  # Point on the right
         ind_x1 = im-1
         ind_x2 = im
+    ind_x1 = max(min(ind_x1, x.size-1), 0)
+    ind_x2 = max(min(ind_x2, x.size-1), 0)
     dx1 = xr_i - x[ind_x1]
     dx2 = x[ind_x2] - xr_i
-    return ind_x1, ind_x2, dx1, dx2
+
+    Dx =  x[ind_x2] - x[ind_x1]
+    if ind_x1 == ind_x2:
+        dx1 = 0.5
+        dx2 = 0.5
+        Dx = 1
+    return ind_x1, ind_x2, dx1, dx2, Dx
 
 
 def interpmat(locs, x, y=None, z=None):
@@ -70,14 +77,17 @@ def _interpmat1D(locs, x):
     Q = sp.lil_matrix((npts, nx))
 
     for i in range(npts):
-        ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i])
-        dv = (x[ind_x2] - x[ind_x1])
-        Dx =  x[ind_x2] - x[ind_x1]
+        ind_x1, ind_x2, dx1, dx2, Dx = _interp_point_1D(x, locs[i])
+        # dv = (x[ind_x2] - x[ind_x1])
         # Get the row in the matrix
         inds = [ind_x1, ind_x2]
         vals = [(1-dx1/Dx),(1-dx2/Dx)]
-        Q[i, inds] = vals
-    return Q.tocsr()
+
+        for I, V in zip(inds, vals):
+            Q[i, I] += V
+        # Q[i, mkvc(inds)] = vals
+
+    return Q
 
 
 
@@ -91,13 +101,10 @@ def _interpmat2D(locs, x, y):
 
 
     for i in range(npts):
-        ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i, 0])
-        ind_y1, ind_y2, dy1, dy2 = _interp_point_1D(y, locs[i, 1])
+        ind_x1, ind_x2, dx1, dx2, Dx = _interp_point_1D(x, locs[i, 0])
+        ind_y1, ind_y2, dy1, dy2, Dy = _interp_point_1D(y, locs[i, 1])
 
-        dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1])
-
-        Dx =  x[ind_x2] - x[ind_x1]
-        Dy =  y[ind_y2] - y[ind_y1]
+        # dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1])
 
         # Get the row in the matrix
 
@@ -112,9 +119,13 @@ def _interpmat2D(locs, x, y):
                 (1-dx2/Dx)*(1-dy1/Dy),
                 (1-dx2/Dx)*(1-dy2/Dy)]
 
-        Q[i, mkvc(inds)] = vals
 
-    return Q.tocsr()
+        for I, V in zip(mkvc(inds), vals):
+            Q[i, I] += V
+        # Q[i, mkvc(inds)] = vals
+
+
+    return Q
 
 
 
@@ -129,15 +140,11 @@ def _interpmat3D(locs, x, y, z):
 
 
     for i in range(npts):
-        ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i, 0])
-        ind_y1, ind_y2, dy1, dy2 = _interp_point_1D(y, locs[i, 1])
-        ind_z1, ind_z2, dz1, dz2 = _interp_point_1D(z, locs[i, 2])
+        ind_x1, ind_x2, dx1, dx2, Dx = _interp_point_1D(x, locs[i, 0])
+        ind_y1, ind_y2, dy1, dy2, Dy = _interp_point_1D(y, locs[i, 1])
+        ind_z1, ind_z2, dz1, dz2, Dz = _interp_point_1D(z, locs[i, 2])
 
-        dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1])
-
-        Dx =  x[ind_x2] - x[ind_x1]
-        Dy =  y[ind_y2] - y[ind_y1]
-        Dz =  z[ind_z2] - z[ind_z1]
+        # dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1])
 
         # Get the row in the matrix
 
@@ -160,20 +167,21 @@ def _interpmat3D(locs, x, y, z):
                 (1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz),
                 (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz)]
 
-        Q[i, mkvc(inds)] = vals
+        for I, V in zip(mkvc(inds), vals):
+            Q[i, I] += V
+        # Q[i, mkvc(inds)] = vals
 
-    return Q.tocsr()
+    return Q
 
 
 if __name__ == '__main__':
-    import SimPEG
-    import numpy as np
+    from SimPEG import *
     import matplotlib.pyplot as plt
     locs = np.random.rand(50)*0.8+0.1
     x = np.linspace(0,1,7)
     dense = np.linspace(0,1,200)
     fun = lambda x: np.cos(2*np.pi*x)
-    Q = SimPEG.Utils.interpmat(locs, x)
+    Q = Utils.interpmat(locs, x)
     plt.plot(x, fun(x), 'bs-')
     plt.plot(dense, fun(dense), 'y:')
     plt.plot(locs, Q*fun(x), 'mo')

SimPEG/Utils/lomutils.py

@@ -1,7 +1,6 @@
 import numpy as np
 from scipy import sparse as sp
-from matutils import mkvc, ndgrid, sub2ind
-from sputils import sdiag
+from matutils import mkvc, ndgrid, sub2ind, sdiag
 
 
 def volTetra(xyz, A, B, C, D):
@@ -188,78 +187,3 @@ def faceInfo(xyz, A, B, C, D, average=True, normalizeNormals=True):
 
     return N, area
 
-
-def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33):
-    """ B = inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33)
-
-    inverts a stack of 3x3 matrices
-
-    Input:
-     A   - a11, a12, a13, a21, a22, a23, a31, a32, a33
-
-    Output:
-     B   - inverse
-     """
-
-    a11 = mkvc(a11)
-    a12 = mkvc(a12)
-    a13 = mkvc(a13)
-    a21 = mkvc(a21)
-    a22 = mkvc(a22)
-    a23 = mkvc(a23)
-    a31 = mkvc(a31)
-    a32 = mkvc(a32)
-    a33 = mkvc(a33)
-
-    detA = a31*a12*a23 - a31*a13*a22 - a21*a12*a33 + a21*a13*a32 + a11*a22*a33 - a11*a23*a32
-
-    b11 = +(a22*a33 - a23*a32)/detA
-    b12 = -(a12*a33 - a13*a32)/detA
-    b13 = +(a12*a23 - a13*a22)/detA
-
-    b21 = +(a31*a23 - a21*a33)/detA
-    b22 = -(a31*a13 - a11*a33)/detA
-    b23 = +(a21*a13 - a11*a23)/detA
-
-    b31 = -(a31*a22 - a21*a32)/detA
-    b32 = +(a31*a12 - a11*a32)/detA
-    b33 = -(a21*a12 - a11*a22)/detA
-
-    B = sp.vstack((sp.hstack((sdiag(b11), sdiag(b12),  sdiag(b13))),
-                   sp.hstack((sdiag(b21), sdiag(b22),  sdiag(b23))),
-                   sp.hstack((sdiag(b31), sdiag(b32),  sdiag(b33)))))
-
-    return B
-
-
-def inv2X2BlockDiagonal(a11, a12, a21, a22):
-    """ B = inv2X2BlockDiagonal(a11, a12, a21, a22)
-
-    Inverts a stack of 2x2 matrices by using the inversion formula
-
-    inv(A) = (1/det(A)) * cof(A)^T
-
-    Input:
-    A   - a11, a12, a13, a21, a22, a23, a31, a32, a33
-
-    Output:
-    B   - inverse
-    """
-
-    a11 = mkvc(a11)
-    a12 = mkvc(a12)
-    a21 = mkvc(a21)
-    a22 = mkvc(a22)
-
-    # compute inverse of the determinant.
-    detAinv = 1./(a11*a22 - a21*a12)
-
-    b11 = +detAinv*a22
-    b12 = -detAinv*a12
-    b21 = -detAinv*a21
-    b22 = +detAinv*a11
-
-    B = sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))),
-                   sp.hstack((sdiag(b21), sdiag(b22)))))
-
-    return B

SimPEG/Utils/matutils.py

@@ -1,5 +1,5 @@
 import numpy as np
-
+import scipy.sparse as sp
 
 def mkvc(x, numDims=1):
     """Creates a vector with the number of dimension specified
@@ -30,6 +30,42 @@ def mkvc(x, numDims=1):
     elif numDims == 3:
         return x.flatten(order='F')[:, np.newaxis, np.newaxis]
 
+def sdiag(h):
+    """Sparse diagonal matrix"""
+    return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr")
+
+def sdInv(M):
+    "Inverse of a sparse diagonal matrix"
+    return sdiag(1/M.diagonal())
+
+def speye(n):
+    """Sparse identity"""
+    return sp.identity(n, format="csr")
+
+
+def kron3(A, B, C):
+    """Three kron prods"""
+    return sp.kron(sp.kron(A, B), C, format="csr")
+
+
+def spzeros(n1, n2):
+    """spzeros"""
+    return sp.coo_matrix((n1, n2)).tocsr()
+
+
+def ddx(n):
+    """Define 1D derivatives, inner, this means we go from n+1 to n"""
+    return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr")
+
+
+def av(n):
+    """Define 1D averaging operator from nodes to cell-centers."""
+    return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
+
+def avExtrap(n):
+    """Define 1D averaging operator from cell-centers to nodes."""
+    Av = sp.spdiags((0.5*np.ones((n, 1))*[1, 1]).T, [-1, 0], n+1, n, format="csr") + sp.csr_matrix(([0.5,0.5],([0,n],[0,n-1])),shape=(n+1,n))
+    return Av
 
 def ndgrid(*args, **kwargs):
     """
@@ -98,33 +134,23 @@ def ndgrid(*args, **kwargs):
             return XYZ[2], XYZ[1], XYZ[0]
 
 
-def ind2sub(shape, ind):
-    """From the given shape, returns the subscrips of the given index"""
-    revshp = []
-    revshp.extend(shape)
-    mult = [1]
-    for i in range(0, len(revshp)-1):
-        mult.extend([mult[i]*revshp[i]])
-    mult = np.array(mult).reshape(len(mult))
-
-    sub = []
-
-    for i in range(0, len(shape)):
-        sub.extend([np.math.floor(ind / mult[i])])
-        ind = ind - (np.math.floor(ind/mult[i]) * mult[i])
-    return sub
+def ind2sub(shape, inds):
+    """From the given shape, returns the subscripts of the given index"""
+    if type(inds) is not np.ndarray:
+        inds = np.array(inds)
+    assert len(inds.shape) == 1, 'Indexing must be done as a 1D row vector, e.g. [3,6,6,...]'
+    return np.unravel_index(inds, shape, order='F')
 
 
 def sub2ind(shape, subs):
     """From the given shape, returns the index of the given subscript"""
-    revshp = list(shape)
-    mult = [1]
-    for i in range(0, len(revshp)-1):
-        mult.extend([mult[i]*revshp[i]])
-    mult = np.array(mult).reshape(len(mult), 1)
-
-    idx = np.dot((subs), (mult))
-    return idx
+    if type(subs) is not np.ndarray:
+        subs = np.array(subs)
+    if subs.size == len(shape):
+        subs = subs[np.newaxis,:]
+    assert subs.shape[1] == len(shape), 'Indexing must be done as a column vectors. e.g. [[3,6],[6,2],...]'
+    inds = np.ravel_multi_index(subs.T, shape, order='F')
+    return mkvc(inds)
 
 
 def getSubArray(A, ind):
@@ -138,3 +164,159 @@ def getSubArray(A, ind):
         return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]]
     else:
         raise Exception("getSubArray does not support dimension asked.")
+
+
+def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33, returnMatrix=True):
+    """ B = inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33)
+
+    inverts a stack of 3x3 matrices
+
+    Input:
+     A   - a11, a12, a13, a21, a22, a23, a31, a32, a33
+
+    Output:
+     B   - inverse
+     """
+
+    a11 = mkvc(a11)
+    a12 = mkvc(a12)
+    a13 = mkvc(a13)
+    a21 = mkvc(a21)
+    a22 = mkvc(a22)
+    a23 = mkvc(a23)
+    a31 = mkvc(a31)
+    a32 = mkvc(a32)
+    a33 = mkvc(a33)
+
+    detA = a31*a12*a23 - a31*a13*a22 - a21*a12*a33 + a21*a13*a32 + a11*a22*a33 - a11*a23*a32
+
+    b11 = +(a22*a33 - a23*a32)/detA
+    b12 = -(a12*a33 - a13*a32)/detA
+    b13 = +(a12*a23 - a13*a22)/detA
+
+    b21 = +(a31*a23 - a21*a33)/detA
+    b22 = -(a31*a13 - a11*a33)/detA
+    b23 = +(a21*a13 - a11*a23)/detA
+
+    b31 = -(a31*a22 - a21*a32)/detA
+    b32 = +(a31*a12 - a11*a32)/detA
+    b33 = -(a21*a12 - a11*a22)/detA
+
+    if not returnMatrix:
+        return b11, b12, b13, b21, b22, b23, b31, b32, b33
+
+    return sp.vstack((sp.hstack((sdiag(b11), sdiag(b12),  sdiag(b13))),
+                      sp.hstack((sdiag(b21), sdiag(b22),  sdiag(b23))),
+                      sp.hstack((sdiag(b31), sdiag(b32),  sdiag(b33)))))
+
+
+
+def inv2X2BlockDiagonal(a11, a12, a21, a22, returnMatrix=True):
+    """ B = inv2X2BlockDiagonal(a11, a12, a21, a22)
+
+    Inverts a stack of 2x2 matrices by using the inversion formula
+
+    inv(A) = (1/det(A)) * cof(A)^T
+
+    Input:
+    A   - a11, a12, a21, a22
+
+    Output:
+    B   - inverse
+    """
+
+    a11 = mkvc(a11)
+    a12 = mkvc(a12)
+    a21 = mkvc(a21)
+    a22 = mkvc(a22)
+
+    # compute inverse of the determinant.
+    detAinv = 1./(a11*a22 - a21*a12)
+
+    b11 = +detAinv*a22
+    b12 = -detAinv*a12
+    b21 = -detAinv*a21
+    b22 = +detAinv*a11
+
+    if not returnMatrix:
+        return b11, b12, b21, b22
+
+    return sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))),
+                      sp.hstack((sdiag(b21), sdiag(b22)))))
+
+def makePropertyTensor(M, sigma):
+    if sigma is None:  # default is ones
+        sigma = np.ones(M.nC)
+
+    if type(sigma) in [float, int, long]:
+        sigma = sigma * np.ones(M.nC)
+
+    if M.dim == 1:
+        if sigma.size == M.nC:  # Isotropic!
+            sigma = mkvc(sigma)  # ensure it is a vector.
+            Sigma = sdiag(sigma)
+        else:
+            raise Exception('Unexpected shape of sigma')
+    elif M.dim == 2:
+        if sigma.size == M.nC:  # Isotropic!
+            sigma = mkvc(sigma)  # ensure it is a vector.
+            Sigma = sdiag(np.r_[sigma, sigma])
+        elif sigma.shape[1] == 2:  # Diagonal tensor
+            Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]])
+        elif sigma.shape[1] == 3:  # Fully anisotropic
+            row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2])))
+            row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1])))
+            Sigma = sp.vstack((row1, row2))
+        else:
+            raise Exception('Unexpected shape of sigma')
+    elif M.dim == 3:
+        if sigma.size == M.nC:  # Isotropic!
+            sigma = mkvc(sigma)  # ensure it is a vector.
+            Sigma = sdiag(np.r_[sigma, sigma, sigma])
+        elif sigma.shape[1] == 3:  # Diagonal tensor
+            Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]])
+        elif sigma.shape[1] == 6:  # Fully anisotropic
+            row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4])))
+            row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5])))
+            row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2])))
+            Sigma = sp.vstack((row1, row2, row3))
+        else:
+            raise Exception('Unexpected shape of sigma')
+    return Sigma
+
+
+def invPropertyTensor(M, tensor, returnMatrix=False):
+
+    T = None
+
+    if type(tensor) in [float, int, long]:
+        T = 1./tensor
+
+    elif tensor.size == M.nC:  # Isotropic!
+        T = 1./mkvc(tensor)  # ensure it is a vector.
+
+    elif M.dim == 2:
+        if tensor.shape[1] == 2:  # Diagonal tensor
+            T = 1./tensor
+        elif tensor.shape[1] == 3:  # Fully anisotropic
+            B = inv2X2BlockDiagonal(tensor[:,0], tensor[:,2],
+                                    tensor[:,2], tensor[:,1],
+                                    returnMatrix=False)
+            b11, b12, b21, b22 = B
+            T = np.c_[b11, b22, b12]
+    elif M.dim == 3:
+        if tensor.shape[1] == 3:  # Diagonal tensor
+            T = 1./tensor
+        elif tensor.shape[1] == 6:  # Fully anisotropic
+            B = inv3X3BlockDiagonal(tensor[:,0], tensor[:,3], tensor[:,4],
+                                    tensor[:,3], tensor[:,1], tensor[:,5],
+                                    tensor[:,4], tensor[:,5], tensor[:,2],
+                                    returnMatrix=False)
+            b11, b12, b13, b21, b22, b23, b31, b32, b33 = B
+            T = np.c_[b11, b22, b33, b12, b13, b23]
+
+    if T is None:
+        raise Exception('Unexpected shape of tensor')
+    if returnMatrix:
+        return makePropertyTensor(M, T)
+    return T

SimPEG/Utils/meshutils.py

@@ -1,7 +1,6 @@
 import numpy as np
 from scipy import sparse as sp
-from matutils import mkvc, ndgrid, sub2ind
-from sputils import sdiag
+from matutils import mkvc, ndgrid, sub2ind, sdiag
 
 def exampleLomGird(nC, exType):
     assert type(nC) == list, "nC must be a list containing the number of nodes"
@@ -30,7 +29,12 @@ def meshTensors(*args):
     """
         **meshTensors** takes any number of tuples that have the form::
 
-            h1 = ( (numPad, sizeStart [, increaseFactor]), (numCore, sizeCore), (numPad, sizeStart [, increaseFactor]) )
+            mT = ( (numPad, sizeStart [, increaseFactor]), (numCore, sizeCore), (numPad, sizeStart [, increaseFactor]) )
+
+        .. note::
+
+            The increaseFactor is an optional input.
+
 
         .. plot::
 

SimPEG/Utils/sputils.py

@@ -1,41 +0,0 @@
-from scipy import sparse as sp
-from matutils import mkvc
-import numpy as np
-
-
-def sdiag(h):
-    """Sparse diagonal matrix"""
-    return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr")
-
-def sdInv(M):
-    "Inverse of a sparse diagonal matrix"
-    return sdiag(1/M.diagonal())
-
-def speye(n):
-    """Sparse identity"""
-    return sp.identity(n, format="csr")
-
-
-def kron3(A, B, C):
-    """Three kron prods"""
-    return sp.kron(sp.kron(A, B), C, format="csr")
-
-
-def spzeros(n1, n2):
-    """spzeros"""
-    return sp.coo_matrix((n1, n2)).tocsr()
-
-
-def ddx(n):
-    """Define 1D derivatives, inner, this means we go from n+1 to n"""
-    return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr")
-
-
-def av(n):
-    """Define 1D averaging operator from nodes to cell-centers."""
-    return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
-
-def avExtrap(n):
-    """Define 1D averaging operator from cell-centers to nodes."""
-    Av = sp.spdiags((0.5*np.ones((n, 1))*[1, 1]).T, [-1, 0], n+1, n, format="csr") + sp.csr_matrix(([0.5,0.5],([0,n],[0,n-1])),shape=(n+1,n))
-    return Av

docs/api_BaseMesh.rst

@@ -1,8 +0,0 @@
-.. _api_BaseMesh:
-
-Base Mesh
-*********
-
-.. automodule:: SimPEG.Mesh.BaseMesh
-    :members:
-    :undoc-members:

docs/api_Cyl1DMesh.rst

@@ -1,8 +0,0 @@
-.. _api_Cyl1DMesh:
-
-Cylindrical 1D Mesh
-*******************
-
-.. automodule:: SimPEG.Mesh.Cyl1DMesh
-    :members:
-    :undoc-members:

docs/api_DiffOperators.rst

@@ -1,8 +0,0 @@
-.. _api_DiffOperators:
-
-Differential Operators
-**********************
-
-.. automodule:: SimPEG.Mesh.DiffOperators
-    :members:
-    :undoc-members:

docs/api_Forward.rst

@@ -2,7 +2,7 @@
 
 
 Model
-*****
+=====
 
 .. automodule:: SimPEG.Model
     :show-inheritance:
@@ -11,7 +11,7 @@ Model
     :inherited-members:
 
 Data
-****
+====
 
 .. automodule:: SimPEG.Data
     :show-inheritance:
@@ -20,7 +20,7 @@ Data
     :inherited-members:
 
 Problem
-*******
+=======
 
 .. automodule:: SimPEG.Problem
     :show-inheritance:

docs/api_InnerProducts.rst

@@ -1,8 +1,166 @@
 .. _api_InnerProducts:
 
+
 Inner Products
 **************
 
+By using the weak formulation of many of the PDEs in geophysical applications, we can rapidly develop discretizations. Much of this work, however, needs a good understanding of how to approximate inner products on our discretized meshes. We will define the inner product as:
+
+.. math::
+
+    \left(a,b\right) = \int_\Omega{a \cdot b}{\partial v}
+
+where a and b are either scalars or vectors.
+
+.. note::
+
+    The InnerProducts class is a base class providing inner product matrices for meshes and cannot run on its own.
+
+
+Example problem for DC resistivity
+----------------------------------
+We will start with the formulation of the Direct Current (DC) resistivity problem in geophysics.
+
+
+.. math::
+
+        \frac{1}{\sigma}\vec{j} = \nabla \phi \\
+
+        \nabla\cdot \vec{j} = q
+
+In the following discretization, \\\( \\sigma \\\) and \\\( \\phi \\\)
+will be discretized on the cell-centers and the flux, \\\(\\vec{j}\\\),
+will be on the faces. We will use the weak formulation to discretize
+the DC resistivity equation.
+
+We can define in weak form by integrating with a general face function \\\(\\vec{f}\\\):
+
+.. math::
+
+    \int_{\Omega}{\sigma^{-1}\vec{j} \cdot \vec{f}} = \int_{\Omega}{\nabla \phi  \cdot \vec{f}}
+
+Here we can integrate the right side by parts,
+
+.. math::
+
+    \nabla\cdot(\phi\vec{f})=\nabla\phi\cdot\vec{f} + \phi\nabla\cdot\vec{f}
+
+and rearrange it, and apply the Divergence theorem.
+
+.. math::
+
+    \int_{\Omega}{\sigma^{-1}\vec{j} \cdot \vec{f}} =
+    - \int_{\Omega}{(\phi \nabla \cdot \vec{f})}
+    + \int_{\partial \Omega}{ \phi  \vec{f} \cdot \mathbf{n}}
+
+We can then discretize for every cell:
+
+.. math::
+
+    v_{\text{cell}} \sigma^{-1} (\mathbf{J}_x \mathbf{F}_x +\mathbf{J}_y \mathbf{F}_y  + \mathbf{J}_z \mathbf{F}_z ) = -\phi^{\top} v_{\text{cell}} \mathbf{D}_{\text{cell}} \mathbf{F}  + \text{BC}
+
+.. note::
+
+    We have discretized the dot product above, but remember that we do not really have a single vector \\\(\\mathbf{J}\\\), but approximations of \\\(\\vec{j}\\\) on each face of our cell. In 2D that means 2 approximations of \\\(\\mathbf{J}_x\\\) and 2 approximations of \\\(\\mathbf{J}_y\\\). In 3D we also have 2 approximations of \\\(\\mathbf{J}_z\\\).
+
+Regardless of how we choose to approximate this dot product, we can represent this in vector form (again this is for every cell), and will generalize for the case of anisotropic (tensor) sigma.
+
+.. math::
+
+    \mathbf{F}_c^{\top} (\sqrt{v_{\text{cell}}} \Sigma^{-1} \sqrt{v_{\text{cell}}})  \mathbf{J}_c =
+    -\phi^{\top} v_{\text{cell}} \mathbf{D}_{\text{cell}} \mathbf{F})
+    + \text{BC}
+
+We multiply by  square-root of volume on each side of the tensor conductivity to keep symmetry in the system. Here \\\(\\mathbf{J}_c\\\) is the Cartesian \\\(\\mathbf{J}\\\) (on the faces that we choose to use in our approximation) and must be calculated differently depending on the mesh:
+
+.. math::
+    \mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\text{TENSOR} \\
+    \mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{LRM}
+
+Here the \\\(i\\\) index refers to where we choose to approximate this integral, as discussed in the note above.
+We will approximate this integral by taking the fluxes clustered around every node of the cell, there are 8 combinations in 3D, and 4 in 2D. We will use a projection matrix \\\( \\mathbf{Q}_{(i)} \\\) to pick the appropriate fluxes. So, now that we have 8 approximations of this integral, we will just take the average. For the TensorMesh, this looks like:
+
+.. math::
+
+    \mathbf{F}^{\top}
+        {1\over 8}
+        \left(\sum_{i=1}^8
+        \mathbf{Q}_{(i)}^{\top} \sqrt{v_{\text{cell}}} \Sigma^{-1} \sqrt{v_{\text{cell}}}  \mathbf{Q}_{(i)}
+        \right)
+        \mathbf{J}
+        =
+        -\mathbf{F}^{\top} \mathbf{D}_{\text{cell}}^{\top} v_{\text{cell}} \phi   + \text{BC}
+
+Or, when generalizing to the entire mesh and dropping our general face function:
+
+.. math::
+
+    \mathbf{M}^f_{\Sigma^{-1}} \mathbf{J}
+        =
+        - \mathbf{D}^{\top} \text{diag}(\mathbf{v}) \phi   + \text{BC}
+
+By defining the faceInnerProduct (8 combinations of fluxes in 3D, 4 in 2D, 2 in 1D) to be:
+
+.. math::
+
+    \mathbf{M}^f_{\Sigma^{-1}} =
+        \sum_{i=1}^{2^d}
+        \mathbf{P}_{(i)}^{\top} \Sigma^{-1} \mathbf{P}_{(i)}
+
+Where \\\(d\\\) is the dimension of the mesh.
+The \\\( \\mathbf{M}^f \\\) is returned when given the input of \\\( \\Sigma^{-1} \\\).
+
+Here each \\( \\mathbf{P} \\in \\mathbb{R}^{(d*nC, nF)} \\\) is a combination of the projection, volume, and any normalization to Cartesian coordinates (where the dot product is well defined):
+
+.. math::
+
+    \mathbf{P}_{(i)} =  \sqrt{ \frac{1}{2^d} \mathbf{I}^d \otimes \text{diag}(\mathbf{v})} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\text{LRM only}} \mathbf{Q}_{(i)}
+
+.. note::
+
+    This is actually completed for each cell in the mesh at the same time, and the full matrices are returned.
+
+If ``returnP=True`` is requested in any of these methods the projection matrices are returned as a list ordered by nodes around which the fluxes were approximated::
+
+    # In 3D
+    P = [P000, P100, P010, P110, P001, P101, P011, P111]
+    # In 2D
+    P = [P00, P10, P01, P11]
+    # In 1D
+    P = [P0, P1]
+
+The derivation for ``edgeInnerProducts`` is exactly the same, however, when we approximate the integral using the fields around each node, the projection matrices look a bit different because we have 12 edges in 3D instead of just 6 faces. The interface to the code is exactly the same.
+
+Defining Tensor Properties
+--------------------------
+
+**For 3D:**
+
+Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows:
+
+.. math::
+
+    \vec{\mu} = \left[\begin{matrix} \mu_{1} & 0 & 0 \\ 0 & \mu_{1} & 0 \\ 0 & 0 & \mu_{1}  \end{matrix}\right]
+
+    \vec{\mu} = \left[\begin{matrix} \mu_{1} & 0 & 0 \\ 0 & \mu_{2} & 0 \\ 0 & 0 & \mu_{3}  \end{matrix}\right]
+
+    \vec{\mu} = \left[\begin{matrix} \mu_{1} & \mu_{4} & \mu_{5} \\ \mu_{4} & \mu_{2} & \mu_{6} \\ \mu_{5} & \mu_{6} & \mu_{3}  \end{matrix}\right]
+
+**For 2D:**
+
+ Depending on the number of columns (either 1, 2, or 3) of mu, the material property is interpreted as follows:
+
+.. math::
+    \vec{\mu} = \left[\begin{matrix} \mu_{1} & 0 \\ 0 & \mu_{1} \end{matrix}\right]
+
+    \vec{\mu} = \left[\begin{matrix} \mu_{1} & 0 \\ 0 & \mu_{2} \end{matrix}\right]
+
+    \vec{\mu} = \left[\begin{matrix} \mu_{1} & \mu_{3} \\ \mu_{3} & \mu_{2} \end{matrix}\right]
+
+
+The API
+-------
+
 .. automodule:: SimPEG.Mesh.InnerProducts
     :members:
     :undoc-members:

docs/api_Inverse.rst

@@ -2,7 +2,7 @@
 
 
 Regularization
-**************
+==============
 
 .. automodule:: SimPEG.Regularization
     :show-inheritance:
@@ -11,7 +11,7 @@ Regularization
 
 
 Objective Function
-******************
+==================
 
 .. automodule:: SimPEG.ObjFunction
     :members:
@@ -19,7 +19,7 @@ Objective Function
 
 
 Optimize
-********
+========
 
 .. automodule:: SimPEG.Optimization
     :show-inheritance:
@@ -28,7 +28,7 @@ Optimize
     :undoc-members:
 
 Inversion
-*********
+=========
 
 .. automodule:: SimPEG.Inversion
     :show-inheritance:

docs/api_LogicallyOrthogonalMesh.rst

@@ -1,10 +0,0 @@
-.. _api_LogicallyOrthogonalMesh:
-
-Logically Orthogonal Mesh
-*************************
-
-.. automodule:: SimPEG.Mesh.LogicallyOrthogonalMesh
-    :show-inheritance:
-    :members:
-    :undoc-members:
-    :inherited-members:

docs/api_Mesh.rst

@@ -0,0 +1,45 @@
+.. _api_Mesh:
+
+
+Tensor Mesh
+===========
+
+.. automodule:: SimPEG.Mesh.TensorMesh
+    :show-inheritance:
+    :members:
+    :undoc-members:
+    :inherited-members:
+
+
+Cylindrical 1D Mesh
+===================
+
+.. automodule:: SimPEG.Mesh.Cyl1DMesh
+    :members:
+    :undoc-members:
+
+
+Logically Orthogonal Mesh
+=========================
+
+.. automodule:: SimPEG.Mesh.LogicallyOrthogonalMesh
+    :show-inheritance:
+    :members:
+    :undoc-members:
+    :inherited-members:
+
+
+Base Mesh
+=========
+
+.. automodule:: SimPEG.Mesh.BaseMesh
+    :members:
+    :undoc-members:
+
+
+Differential Operators
+======================
+
+.. automodule:: SimPEG.Mesh.DiffOperators
+    :members:
+    :undoc-members:

docs/api_Parameters.rst

@@ -2,7 +2,7 @@
 
 
 Parameters
-**********
+==========
 
 .. automodule:: SimPEG.Parameters
     :show-inheritance:

docs/api_TensorMesh.rst

@@ -1,10 +0,0 @@
-.. _api_TensorMesh:
-
-Tensor Mesh
-***********
-
-.. automodule:: SimPEG.Mesh.TensorMesh
-    :show-inheritance:
-    :members:
-    :undoc-members:
-    :inherited-members:

docs/api_Tests.rst

@@ -1,7 +1,7 @@
 .. _api_Tests:
 
 Testing SimPEG
-**************
+==============
 
 .. automodule:: SimPEG.Tests.TestUtils
     :members:

docs/api_Utils.rst

@@ -17,42 +17,35 @@ Utilities
     :undoc-members:
 
 Matrix Utilities
-****************
+================
 
 .. automodule:: SimPEG.Utils.matutils
     :members:
     :undoc-members:
 
-Sparse Utilities
-****************
-
-.. automodule:: SimPEG.Utils.sputils
-    :members:
-    :undoc-members:
-
 LOM Utilities
-*************
+=============
 
 .. automodule:: SimPEG.Utils.lomutils
     :members:
     :undoc-members:
 
 Mesh Utilities
-**************
+==============
 
 .. automodule:: SimPEG.Utils.meshutils
     :members:
     :undoc-members:
 
 Model Builder Utilities
-***********************
+=======================
 
 .. automodule:: SimPEG.Utils.ModelBuilder
     :members:
     :undoc-members:
 
 Interpolation Utilities
-***********************
+=======================
 
 .. automodule:: SimPEG.Utils.interputils
     :members:

docs/index.rst

@@ -11,10 +11,12 @@ The vision is to create a package for finite volume simulation and parameter est
 applications to geophysical imaging and subsurface flow. To enable
 these goals, this package has the following features:
 
-* is modular with respect to discretization, physics, optimization, and regularization
-* is built with the (large-scale) inverse problem in mind
-* provides a framework for geophysical and hydrogeologic problems
-* supports 1D, 2D and 3D problems
+
+   * is modular with respect to discretization, physics, optimization, and regularization
+   * is built with the (large-scale) inverse problem in mind
+   * provides a framework for geophysical and hydrogeologic problems
+   * supports 1D, 2D and 3D problems
+
 
 .. image:: simpeg-framework.png
    :width: 400 px
@@ -22,20 +24,16 @@ these goals, this package has the following features:
    :align: center
 
 Meshing & Operators
-===================
+*******************
 
 .. toctree::
    :maxdepth: 2
 
-   api_BaseMesh
-   api_TensorMesh
-   api_LogicallyOrthogonalMesh
-   api_Cyl1DMesh
-   api_DiffOperators
+   api_Mesh
    api_InnerProducts
 
 Forward Problems
-================
+****************
 
 .. toctree::
    :maxdepth: 2
@@ -43,7 +41,7 @@ Forward Problems
    api_Forward
 
 Inversion
-=========
+*********
 
 .. toctree::
    :maxdepth: 2
@@ -52,7 +50,7 @@ Inversion
    api_Parameters
 
 Testing SimPEG
-==============
+**************
 
 .. toctree::
    :maxdepth: 2
@@ -60,20 +58,20 @@ Testing SimPEG
    api_Tests
 
 * Master Branch
-   .. image:: https://travis-ci.org/simpeg/simpeg.png?branch=master
+   .. image:: https://travis-ci.org/simpeg/simpeg.png?branch*master
       :target: https://travis-ci.org/simpeg/simpeg
       :alt: Master Branch
       :align: center
 
 * Develop Branch
-   .. image:: https://travis-ci.org/simpeg/simpeg.png?branch=develop
+   .. image:: https://travis-ci.org/simpeg/simpeg.png?branch*develop
       :target: https://travis-ci.org/simpeg/simpeg
       :alt: Develop Branch
       :align: center
 
 
 Utility Codes
-=============
+*************
 
 .. toctree::
    :maxdepth: 2
@@ -82,7 +80,7 @@ Utility Codes
 
 
 Project Index & Search
-======================
+**********************
 
 * :ref:`genindex`
 * :ref:`modindex`