diff --git a/SimPEG/TensorView.py b/SimPEG/TensorView.py
index ba4f1b63..687a520d 100644
--- a/SimPEG/TensorView.py
+++ b/SimPEG/TensorView.py
@@ -2,7 +2,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib
 from mpl_toolkits.mplot3d import Axes3D
-from utils import mkvc
+from SimPEG.utils import mkvc
 
 
 class TensorView(object):
diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py
index ea9faf33..e8946f88 100644
--- a/SimPEG/__init__.py
+++ b/SimPEG/__init__.py
@@ -1,5 +1,4 @@
-from TensorMesh import TensorMesh
-from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
+import mesh
 import utils
 import inverse
 from Solver import Solver
diff --git a/SimPEG/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem/DCProblem.py
index b8115598..1df8897b 100644
--- a/SimPEG/forward/DCProblem/DCProblem.py
+++ b/SimPEG/forward/DCProblem/DCProblem.py
@@ -1,4 +1,4 @@
-from SimPEG import TensorMesh
+from SimPEG.mesh import TensorMesh
 from SimPEG.forward import Problem, SyntheticProblem
 from SimPEG.tests import checkDerivative
 from SimPEG.utils import ModelBuilder, sdiag
diff --git a/SimPEG/BaseMesh.py b/SimPEG/mesh/BaseMesh.py
similarity index 99%
rename from SimPEG/BaseMesh.py
rename to SimPEG/mesh/BaseMesh.py
index c57828e9..29c43dd1 100644
--- a/SimPEG/BaseMesh.py
+++ b/SimPEG/mesh/BaseMesh.py
@@ -1,5 +1,5 @@
 import numpy as np
-from utils import mkvc
+from SimPEG.utils import mkvc
 
 
 class BaseMesh(object):
diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py
new file mode 100644
index 00000000..110fe9bd
--- /dev/null
+++ b/SimPEG/mesh/DiffOperators.py
@@ -0,0 +1,336 @@
+import numpy as np
+from scipy import sparse as sp
+from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros
+
+
+def ddx(n):
+    """Define 1D derivatives, inner, this means we go from n+1 to n+1"""
+    return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr")
+
+
+def checkBC(bc):
+    """ Checks if boundary condition 'bc' is valid. """
+    if(type(bc) is str):
+        bc = [bc, bc]
+    assert type(bc) is list, 'bc must be a list'
+    assert len(bc) == 2, 'bc must have two elements'
+
+    for bc_i in bc:
+        assert type(bc_i) is str, "each bc must be a string"
+        assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'"
+    return bc
+
+
+def ddxCellGrad(n, bc):
+    """Create 1D derivative operator from cell-centres to nodes this means we go from n to n+1"""
+    bc = checkBC(bc)
+
+    D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr")
+    # Set the first side
+    if(bc[0] == 'dirichlet'):
+        D[0, 0] = 2
+    elif(bc[0] == 'neumann'):
+        D[0, 0] = 0
+    # Set the second side
+    if(bc[1] == 'dirichlet'):
+        D[-1, -1] = -2
+    elif(bc[1] == 'neumann'):
+        D[-1, -1] = 0
+    return D
+
+
+def av(n):
+    """Define 1D averaging operator from cell-centres to nodes."""
+    return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
+
+
+class DiffOperators(object):
+    """
+        Class creates the differential operators that you need!
+    """
+    def __init__(self):
+        raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.')
+
+    def faceDiv():
+        doc = "Construct divergence operator (face-stg to cell-centres)."
+
+        def fget(self):
+            if(self._faceDiv is None):
+                # The number of cell centers in each direction
+                n = self.n
+                # Compute faceDivergence operator on faces
+                if(self.dim == 1):
+                    D = ddx(n[0])
+                elif(self.dim == 2):
+                    D1 = sp.kron(speye(n[1]), ddx(n[0]))
+                    D2 = sp.kron(ddx(n[1]), speye(n[0]))
+                    D = sp.hstack((D1, D2), format="csr")
+                elif(self.dim == 3):
+                    D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
+                    D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
+                    D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]))
+                    D = sp.hstack((D1, D2, D3), format="csr")
+                # Compute areas of cell faces & volumes
+                S = self.area
+                V = self.vol
+                self._faceDiv = sdiag(1/V)*D*sdiag(S)
+
+            return self._faceDiv
+        return locals()
+    _faceDiv = None
+    faceDiv = property(**faceDiv())
+
+    def nodalGrad():
+        doc = "Construct gradient operator (nodes to edges)."
+
+        def fget(self):
+            if(self._nodalGrad is None):
+                # The number of cell centers in each direction
+                n = self.n
+                # Compute divergence operator on faces
+                if(self.dim == 1):
+                    G = ddx(n[0])
+                elif(self.dim == 2):
+                    D1 = sp.kron(speye(n[1]+1), ddx(n[0]))
+                    D2 = sp.kron(ddx(n[1]), speye(n[0]+1))
+                    G = sp.vstack((D1, D2), format="csr")
+                elif(self.dim == 3):
+                    D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0]))
+                    D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1))
+                    D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1))
+                    G = sp.vstack((D1, D2, D3), format="csr")
+                # Compute lengths of cell edges
+                L = self.edge
+                self._nodalGrad = sdiag(1/L)*G
+            return self._nodalGrad
+        return locals()
+    _nodalGrad = None
+    nodalGrad = property(**nodalGrad())
+
+    def setCellGradBC(self, BC):
+        """
+        Function that sets the boundary conditions for cell-centred derivative operators.
+
+        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(type(BC) is str):
+            BC = [BC for _ in self.n]  # 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)
+
+        self._cellGrad = None  # ensure we create a new gradient next time we call it
+        self._cellGradBC = BC
+        return BC
+    _cellGradBC = 'neumann'
+
+    def cellGrad():
+        doc = "The cell centered Gradient, takes you to cell faces."
+
+        def fget(self):
+            if(self._cellGrad is None):
+                BC = self.setCellGradBC(self._cellGradBC)
+                n = self.n
+                if(self.dim == 1):
+                    G = ddxCellGrad(n[0], BC[0])
+                elif(self.dim == 2):
+                    G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0]))
+                    G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0]))
+                    G = sp.vstack((G1, G2), format="csr")
+                elif(self.dim == 3):
+                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0]))
+                    G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0]))
+                    G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0]))
+                    G = sp.vstack((G1, G2, G3), format="csr")
+                # Compute areas of cell faces & volumes
+                S = self.area
+                V = self.vol
+                self._cellGrad = sdiag(S)*G*sdiag(1/V)
+            return self._cellGrad
+        return locals()
+    _cellGrad = None
+    cellGrad = property(**cellGrad())
+
+    def edgeCurl():
+        doc = "Construct the 3D curl operator."
+
+        def fget(self):
+            if(self._edgeCurl is None):
+                # The number of cell centers in each direction
+                n1 = self.nCx
+                n2 = self.nCy
+                n3 = self.nCz
+
+                # Compute lengths of cell edges
+                L = self.edge
+
+                # Compute areas of cell faces
+                S = self.area
+
+                # Compute divergence operator on faces
+                d1 = ddx(n1)
+                d2 = ddx(n2)
+                d3 = ddx(n3)
+
+                D32 = kron3(d3, speye(n2), speye(n1+1))
+                D23 = kron3(speye(n3), d2, speye(n1+1))
+                D31 = kron3(d3, speye(n2+1), speye(n1))
+                D13 = kron3(speye(n3), speye(n2+1), d1)
+                D21 = kron3(speye(n3+1), d2, speye(n1))
+                D12 = kron3(speye(n3+1), speye(n2), d1)
+
+                O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1])
+                O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
+                O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])
+
+                C = sp.vstack((sp.hstack((O1, -D32, D23)),
+                               sp.hstack((D31, O2, -D13)),
+                               sp.hstack((-D21, D12, O3))), format="csr")
+
+                self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
+            return self._edgeCurl
+        return locals()
+    _edgeCurl = None
+    edgeCurl = property(**edgeCurl())
+
+    def aveF2CC():
+        doc = "Construct the averaging operator on cell faces to cell centers."
+
+        def fget(self):
+            if(self._aveF2CC is None):
+                n = self.n
+                if(self.dim == 1):
+                    self._aveF2CC = av(n[0])
+                elif(self.dim == 2):
+                    self._aveF2CC = sp.hstack((sp.kron(speye(n[1]), av(n[0])),
+                                               sp.kron(av(n[1]), speye(n[0]))), format="csr")
+                elif(self.dim == 3):
+                    self._aveF2CC = sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
+                                               kron3(speye(n[2]), av(n[1]), speye(n[0])),
+                                               kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
+            return self._aveF2CC
+        return locals()
+    _aveF2CC = None
+    aveF2CC = property(**aveF2CC())
+
+    def aveE2CC():
+        doc = "Construct the averaging operator on cell edges to cell centers."
+
+        def fget(self):
+            if(self._aveE2CC is None):
+                # The number of cell centers in each direction
+                n = self.n
+                if(self.dim == 1):
+                    raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
+                elif(self.dim == 2):
+                    self._aveE2CC = sp.hstack((sp.kron(av(n[1]), speye(n[0])),
+                                               sp.kron(speye(n[1]), av(n[0]))), format="csr")
+                elif(self.dim == 3):
+                    self._aveE2CC = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
+                                               kron3(av(n[2]), speye(n[1]), av(n[0])),
+                                               kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
+            return self._aveE2CC
+        return locals()
+    _aveE2CC = None
+    aveE2CC = property(**aveE2CC())
+
+    def aveN2CC():
+        doc = "Construct the averaging operator on cell nodes to cell centers."
+
+        def fget(self):
+            if(self._aveN2CC is None):
+                # The number of cell centers in each direction
+                n = self.n
+                if(self.dim == 1):
+                    self._aveN2CC = av(n[0])
+                elif(self.dim == 2):
+                    self._aveN2CC = sp.hstack((sp.kron(av(n[1]), av(n[0])),
+                                               sp.kron(av(n[1]), av(n[0]))), format="csr")
+                elif(self.dim == 3):
+                    self._aveN2CC = sp.hstack((kron3(av(n[2]), av(n[1]), av(n[0])),
+                                               kron3(av(n[2]), av(n[1]), av(n[0])),
+                                               kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr")
+            return self._aveN2CC
+        return locals()
+    _aveN2CC = None
+    aveN2CC = property(**aveN2CC())
+
+    def aveN2CCv():
+        doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension separate."
+
+        def fget(self):
+            if(self._aveN2CCv is None):
+                # The number of cell centers in each direction
+                n = self.n
+                if(self.dim == 1):
+                    self._aveN2CCv = av(n[0])
+                elif(self.dim == 2):
+                    self._aveN2CCv = sp.block_diag((sp.kron(av(n[1]), av(n[0])),
+                                                    sp.kron(av(n[1]), av(n[0]))), format="csr")
+                elif(self.dim == 3):
+                    self._aveN2CCv = sp.block_diag((kron3(av(n[2]), av(n[1]), av(n[0])),
+                                                    kron3(av(n[2]), av(n[1]), av(n[0])),
+                                                    kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr")
+            return self._aveN2CCv
+        return locals()
+    _aveN2CCv = None
+    aveN2CCv = property(**aveN2CCv())
+
+    def getMass(self, materialProp=None, loc='e'):
+        """ Produces mass matricies.
+
+        :param str loc: Average to location: 'e'-edges, 'f'-faces
+        :param None,float,numpy.ndarray materialProp: property to be averaged (see below)
+        :rtype: scipy.sparse.csr.csr_matrix
+        :return: M, the mass matrix
+
+        materialProp can be::
+
+            None            -> takes materialProp = 1 (default)
+            float           -> a constant value for entire domain
+            numpy.ndarray   -> if materialProp.size == self.nC
+                                    3D property model
+                               if materialProp.size = self.nCz
+                                    1D (layered eath) property model
+        """
+        if materialProp is None:
+            materialProp = np.ones(self.nC)
+        elif type(materialProp) is float:
+            materialProp = np.ones(self.nC)*materialProp
+        elif materialProp.shape == (self.nCz,):
+            materialProp = materialProp.repeat(self.nCx*self.nCy)
+        materialProp = mkvc(materialProp)
+        assert materialProp.shape == (self.nC,), "materialProp incorrect shape"
+
+        if loc=='e':
+            Av = self.aveE2CC
+        elif loc=='f':
+            Av = self.aveF2CC
+        else:
+            raise ValueError('Invalid loc')
+
+        diag = Av.T * (self.vol * mkvc(materialProp))
+
+        return sdiag(diag)
+
+    def getEdgeMass(self, materialProp=None):
+        """mass matrix for products of edge functions w'*M(materialProp)*e"""
+        return self.getMass(loc='e', materialProp=materialProp)
+
+    def getFaceMass(self, materialProp=None):
+        """mass matrix for products of face functions w'*M(materialProp)*f"""
+        return self.getMass(loc='f', materialProp=materialProp)
+
+    def getFaceMassDeriv(self):
+        Av = self.aveF2CC
+        return Av.T * sdiag(self.vol)
diff --git a/SimPEG/InnerProducts.py b/SimPEG/mesh/InnerProducts.py
similarity index 100%
rename from SimPEG/InnerProducts.py
rename to SimPEG/mesh/InnerProducts.py
diff --git a/SimPEG/LogicallyOrthogonalMesh.py b/SimPEG/mesh/LogicallyOrthogonalMesh.py
similarity index 99%
rename from SimPEG/LogicallyOrthogonalMesh.py
rename to SimPEG/mesh/LogicallyOrthogonalMesh.py
index 5d2b7ffa..7b1e1bca 100644
--- a/SimPEG/LogicallyOrthogonalMesh.py
+++ b/SimPEG/mesh/LogicallyOrthogonalMesh.py
@@ -3,7 +3,7 @@ from BaseMesh import BaseMesh
 from DiffOperators import DiffOperators
 from InnerProducts import InnerProducts
 from LomView import LomView
-from utils import mkvc, ndgrid, volTetra, indexCube, faceInfo
+from SimPEG.utils import mkvc, ndgrid, volTetra, indexCube, faceInfo
 
 # Some helper functions.
 length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
diff --git a/SimPEG/LomView.py b/SimPEG/mesh/LomView.py
similarity index 99%
rename from SimPEG/LomView.py
rename to SimPEG/mesh/LomView.py
index 4b4f36a3..2a9b242b 100644
--- a/SimPEG/LomView.py
+++ b/SimPEG/mesh/LomView.py
@@ -2,7 +2,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib
 from mpl_toolkits.mplot3d import Axes3D
-from utils import mkvc
+from SimPEG.utils import mkvc
 
 
 class LomView(object):
diff --git a/SimPEG/TensorMesh.py b/SimPEG/mesh/TensorMesh.py
similarity index 96%
rename from SimPEG/TensorMesh.py
rename to SimPEG/mesh/TensorMesh.py
index 2015838c..a3329c14 100644
--- a/SimPEG/TensorMesh.py
+++ b/SimPEG/mesh/TensorMesh.py
@@ -3,7 +3,7 @@ from BaseMesh import BaseMesh
 from TensorView import TensorView
 from DiffOperators import DiffOperators
 from InnerProducts import InnerProducts
-from utils import ndgrid, mkvc
+from SimPEG.utils import ndgrid, mkvc
 
 
 class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/mesh/TensorView.py
new file mode 100644
index 00000000..687a520d
--- /dev/null
+++ b/SimPEG/mesh/TensorView.py
@@ -0,0 +1,335 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib
+from mpl_toolkits.mplot3d import Axes3D
+from SimPEG.utils import mkvc
+
+
+class TensorView(object):
+    """
+    Provides viewing functions for TensorMesh
+
+    This class is inherited by TensorMesh
+    """
+    def __init__(self):
+        pass
+
+    def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False):
+        """
+        Mesh.plotImage(I)
+
+        Plots scalar fields on the given mesh.
+
+        Input:
+
+        :param numpy.array I: scalar field
+
+        Optional Input:
+
+        :param str imageType: type of image ('CC','N','F','Fx','Fy','Fz','E','Ex','Ey','Ez') or combinations, e.g. ExEy or FxFz
+        :param int figNum: number of figure to plot to
+        :param matplotlib.axes.Axes ax: axis to plot to
+        :param str direction: slice dimensions, 3D only ('x', 'y', 'z')
+        :param bool numbering: show numbering of slices, 3D only
+        :param str annotationColor: color of annotation, e.g. 'w', 'k', 'b'
+        :param bool showIt: call plt.show()
+
+        .. plot:: examples/mesh/plot_image_2D.py
+           :include-source:
+
+        .. plot:: examples/mesh/plot_image_3D.py
+           :include-source:
+        """
+        assert type(I) == np.ndarray, "I must be a numpy array"
+        assert type(numbering) == bool, "numbering must be a bool"
+        assert direction in ["x", "y","z"], "direction must be either x,y, or z"
+
+
+        if imageType == 'CC':
+            assert I.size == self.nC, "Incorrect dimensions for CC."
+        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')
+        elif imageType == 'Fy':
+            if I.size != np.prod(self.nFy): 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')
+        elif imageType == 'Ex':
+            if I.size != np.prod(self.nEx): 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')
+        elif imageType == 'Ez':
+            if I.size != np.prod(self.nEz): 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}
+            fig = plt.figure(figNum)
+            # Determine the subplot number: 131, 121
+            numPlots = 130 if plotAll else len(imageType)/2*10+100
+            pltNum = 1
+            ex, ey, ez = self.r(I,'E','E','M')
+            if plotAll or 'Ex' in imageType:
+                ax_x = plt.subplot(numPlots+pltNum)
+                self.plotImage(ex, imageType='Ex', ax=ax_x, **options)
+                pltNum +=1
+            if plotAll or 'Ey' in imageType:
+                ax_y = plt.subplot(numPlots+pltNum)
+                self.plotImage(ey, imageType='Ey', ax=ax_y, **options)
+                pltNum +=1
+            if plotAll or 'Ez' in imageType:
+                ax_z = plt.subplot(numPlots+pltNum)
+                self.plotImage(ez, imageType='Ez', ax=ax_z, **options)
+                pltNum +=1
+            return
+        elif imageType[0] == 'F':
+            plotAll = len(imageType) == 1
+            options = {"direction":direction,"numbering":numbering,"annotationColor":annotationColor,"showIt":showIt}
+            fig = plt.figure(figNum)
+            # Determine the subplot number: 131, 121
+            numPlots = 130 if plotAll else len(imageType)/2*10+100
+            pltNum = 1
+            fx, fy, fz = self.r(I,'F','F','M')
+            if plotAll or 'Fx' in imageType:
+                ax_x = plt.subplot(numPlots+pltNum)
+                self.plotImage(fx, imageType='Fx', ax=ax_x, **options)
+                pltNum +=1
+            if plotAll or 'Fy' in imageType:
+                ax_y = plt.subplot(numPlots+pltNum)
+                self.plotImage(fy, imageType='Fy', ax=ax_y, **options)
+                pltNum +=1
+            if plotAll or 'Fz' in imageType:
+                ax_z = plt.subplot(numPlots+pltNum)
+                self.plotImage(fz, imageType='Fz', ax=ax_z, **options)
+                pltNum +=1
+            return
+        else:
+            raise Exception("imageType must be 'CC', 'N','Fx','Fy','Fz','Ex','Ey','Ez'")
+
+
+        if ax is None:
+            fig = plt.figure(figNum)
+            fig.clf()
+            ax = plt.subplot(111)
+        else:
+            assert isinstance(ax,matplotlib.axes.Axes), "ax must be an Axes!"
+            fig = ax.figure
+
+        if self.dim == 1:
+            if imageType == 'CC':
+                ph = ax.plot(self.vectorCCx, I, '-ro')
+            elif imageType == 'N':
+                ph = ax.plot(self.vectorNx, I, '-bs')
+            ax.set_xlabel("x")
+            ax.axis('tight')
+        elif self.dim == 2:
+            if imageType == 'CC':
+                C = I[:].reshape(self.n, order='F')
+            elif imageType == 'N':
+                C = I[:].reshape(self.n+1, 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 = 0.5*(C[:-1, :] + C[1:, :] )
+            elif imageType == 'Fy':
+                C = I[:].reshape(self.nFy, order='F')
+                C = 0.5*(C[:, :-1] + C[:, 1:] )
+            elif imageType == 'Ex':
+                C = I[:].reshape(self.nEx, order='F')
+                C = 0.5*(C[:,:-1] + C[:,1:] )
+            elif imageType == 'Ey':
+                C = I[:].reshape(self.nEy, order='F')
+                C = 0.5*(C[:-1,:] + C[1:,:] )
+
+            ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T)
+            ax.axis('tight')
+            ax.set_xlabel("x")
+            ax.set_ylabel("y")
+
+        elif self.dim == 3:
+            if direction == 'z':
+
+                # get copy of image and average to cell-centres is necessary
+                if imageType == 'CC':
+                    Ic = I[:].reshape(self.n, order='F')
+                elif imageType == 'N':
+                    Ic = I[:].reshape(self.n+1, 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 = .5*(Ic[:-1,:,:]+Ic[1:,:,:])
+                elif imageType == 'Fy':
+                    Ic = I[:].reshape(self.nFy, order='F')
+                    Ic = .5*(Ic[:,:-1,:]+Ic[:,1:,:])
+                elif imageType == 'Fz':
+                    Ic = I[:].reshape(self.nFz, order='F')
+                    Ic = .5*(Ic[:,:,:-1]+Ic[:,:,1:])
+                elif imageType == 'Ex':
+                    Ic = I[:].reshape(self.nEx, 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 = .25*(Ic[:-1,:,:-1]+Ic[1:,:,:-1]+Ic[:-1,:,1:]+Ic[1:,:,:1])
+                elif imageType == 'Ez':
+                    Ic = I[:].reshape(self.nEz, 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
+                nX = np.ceil(np.sqrt(self.nCz))
+                nY = np.ceil(self.nCz/nX)
+
+                #  allocate space for montage
+                nCx = self.nCx
+                nCy = self.nCy
+
+                C = np.zeros((nX*nCx,nY*nCy))
+
+                for iy in range(int(nY)):
+                    for ix in range(int(nX)):
+                        iz = ix + iy*nX
+                        if iz < self.nCz:
+                            C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = Ic[:, :, iz]
+                        else:
+                            C[ix*nCx:(ix+1)*nCx, iy*nCy:(iy+1)*nCy] = np.nan
+
+                C = np.ma.masked_where(np.isnan(C), C)
+                xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())]
+                yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())]
+                # Plot the mesh
+                ph = ax.pcolormesh(xx, yy, C.T)
+                # Plot the lines
+                gx =  np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0])
+                gy =  np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1])
+                # Repeat and seperate with NaN
+                gxX = np.c_[gx, gx, gx+np.nan].ravel()
+                gxY = np.kron(np.ones((nX+1, 1)), np.array([0, sum(self.hy)*nY, np.nan])).ravel()
+                gyX = np.kron(np.ones((nY+1, 1)), np.array([0, sum(self.hx)*nX, np.nan])).ravel()
+                gyY = np.c_[gy, gy, gy+np.nan].ravel()
+                ax.plot(gxX, gxY, annotationColor+'-', linewidth=2)
+                ax.plot(gyX, gyY, annotationColor+'-', linewidth=2)
+                ax.axis('tight')
+
+                if numbering:
+                    pad = np.sum(self.hx)*0.04
+                    for iy in range(int(nY)):
+                        for ix in range(int(nX)):
+                            iz = ix + iy*nX
+                            if iz < self.nCz:
+                                ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad,
+                                         '#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large')
+
+        ax.set_title(imageType)
+        if showIt: plt.show()
+        return ph
+
+    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.
+
+        :param bool nodes: plot nodes
+        :param bool faces: plot faces
+        :param bool centers: plot centers
+        :param bool edges: plot edges
+        :param bool lines: plot lines connecting nodes
+        :param bool showIt: call plt.show()
+
+        .. plot:: examples/mesh/plot_grid_2D.py
+           :include-source:
+
+        .. plot:: examples/mesh/plot_grid_3D.py
+           :include-source:
+        """
+        if self.dim == 1:
+            fig = plt.figure(1)
+            fig.clf()
+            ax = plt.subplot(111)
+            xn = self.gridN
+            xc = self.gridCC
+            ax.hold(True)
+            ax.plot(xn, np.ones(np.shape(xn)), 'bs')
+            ax.plot(xc, np.ones(np.shape(xc)), 'ro')
+            ax.plot(xn, np.ones(np.shape(xn)), 'k--')
+            ax.grid(True)
+            ax.hold(False)
+            ax.set_xlabel('x1')
+            if showIt: plt.show()
+        elif self.dim == 2:
+            fig = plt.figure(2)
+            fig.clf()
+            ax = plt.subplot(111)
+            xn = self.gridN
+            xc = self.gridCC
+            xs1 = self.gridFx
+            xs2 = self.gridFy
+
+            ax.hold(True)
+            if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs')
+            if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro')
+            if faces:
+                ax.plot(xs1[:, 0], xs1[:, 1], 'g>')
+                ax.plot(xs2[:, 0], xs2[:, 1], 'g^')
+
+            # Plot the grid lines
+            if lines:
+                NN = self.r(self.gridN, 'N', 'N', 'M')
+                X1 = np.c_[mkvc(NN[0][0, :]), mkvc(NN[0][self.nCx, :]), mkvc(NN[0][0, :])*np.nan].flatten()
+                Y1 = np.c_[mkvc(NN[1][0, :]), mkvc(NN[1][self.nCx, :]), mkvc(NN[1][0, :])*np.nan].flatten()
+                X2 = np.c_[mkvc(NN[0][:, 0]), mkvc(NN[0][:, self.nCy]), mkvc(NN[0][:, 0])*np.nan].flatten()
+                Y2 = np.c_[mkvc(NN[1][:, 0]), mkvc(NN[1][:, self.nCy]), mkvc(NN[1][:, 0])*np.nan].flatten()
+                X = np.r_[X1, X2]
+                Y = np.r_[Y1, Y2]
+                plt.plot(X, Y)
+
+            ax.grid(True)
+            ax.hold(False)
+            ax.set_xlabel('x1')
+            ax.set_ylabel('x2')
+            if showIt: plt.show()
+        elif self.dim == 3:
+            fig = plt.figure(3)
+            fig.clf()
+            ax = fig.add_subplot(111, projection='3d')
+            xn = self.gridN
+            xc = self.gridCC
+            xfs1 = self.gridFx
+            xfs2 = self.gridFy
+            xfs3 = self.gridFz
+
+            xes1 = self.gridEx
+            xes2 = self.gridEy
+            xes3 = self.gridEz
+
+            ax.hold(True)
+            if nodes: ax.plot(xn[:, 0], xn[:, 1], 'bs', zs=xn[:, 2])
+            if centers: ax.plot(xc[:, 0], xc[:, 1], 'ro', zs=xc[:, 2])
+            if faces:
+                ax.plot(xfs1[:, 0], xfs1[:, 1], 'g>', zs=xfs1[:, 2])
+                ax.plot(xfs2[:, 0], xfs2[:, 1], 'g<', zs=xfs2[:, 2])
+                ax.plot(xfs3[:, 0], xfs3[:, 1], 'g^', zs=xfs3[:, 2])
+            if edges:
+                ax.plot(xes1[:, 0], xes1[:, 1], 'k>', zs=xes1[:, 2])
+                ax.plot(xes2[:, 0], xes2[:, 1], 'k<', zs=xes2[:, 2])
+                ax.plot(xes3[:, 0], xes3[:, 1], 'k^', zs=xes3[:, 2])
+
+            # Plot the grid lines
+            if lines:
+                NN = self.r(self.gridN, 'N', 'N', 'M')
+                X1 = np.c_[mkvc(NN[0][0, :, :]), mkvc(NN[0][self.nCx, :, :]), mkvc(NN[0][0, :, :])*np.nan].flatten()
+                Y1 = np.c_[mkvc(NN[1][0, :, :]), mkvc(NN[1][self.nCx, :, :]), mkvc(NN[1][0, :, :])*np.nan].flatten()
+                Z1 = np.c_[mkvc(NN[2][0, :, :]), mkvc(NN[2][self.nCx, :, :]), mkvc(NN[2][0, :, :])*np.nan].flatten()
+                X2 = np.c_[mkvc(NN[0][:, 0, :]), mkvc(NN[0][:, self.nCy, :]), mkvc(NN[0][:, 0, :])*np.nan].flatten()
+                Y2 = np.c_[mkvc(NN[1][:, 0, :]), mkvc(NN[1][:, self.nCy, :]), mkvc(NN[1][:, 0, :])*np.nan].flatten()
+                Z2 = np.c_[mkvc(NN[2][:, 0, :]), mkvc(NN[2][:, self.nCy, :]), mkvc(NN[2][:, 0, :])*np.nan].flatten()
+                X3 = np.c_[mkvc(NN[0][:, :, 0]), mkvc(NN[0][:, :, self.nCz]), mkvc(NN[0][:, :, 0])*np.nan].flatten()
+                Y3 = np.c_[mkvc(NN[1][:, :, 0]), mkvc(NN[1][:, :, self.nCz]), mkvc(NN[1][:, :, 0])*np.nan].flatten()
+                Z3 = np.c_[mkvc(NN[2][:, :, 0]), mkvc(NN[2][:, :, self.nCz]), mkvc(NN[2][:, :, 0])*np.nan].flatten()
+                X = np.r_[X1, X2, X3]
+                Y = np.r_[Y1, Y2, Y3]
+                Z = np.r_[Z1, Z2, Z3]
+                plt.plot(X, Y, 'b-', zs=Z)
+
+            ax.grid(True)
+            ax.hold(False)
+            ax.set_xlabel('x1')
+            ax.set_ylabel('x2')
+            ax.set_zlabel('x3')
+            if showIt: plt.show()
diff --git a/SimPEG/mesh/__init__.py b/SimPEG/mesh/__init__.py
new file mode 100644
index 00000000..9179efca
--- /dev/null
+++ b/SimPEG/mesh/__init__.py
@@ -0,0 +1,8 @@
+from TensorMesh import TensorMesh
+from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
+from BaseMesh import BaseMesh
+from TensorView import TensorView
+from LomView import LomView
+from InnerProducts import InnerProducts
+from DiffOperators import DiffOperators
+
diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py
index 79340da7..9b2158c4 100644
--- a/SimPEG/tests/TestUtils.py
+++ b/SimPEG/tests/TestUtils.py
@@ -2,7 +2,8 @@ import numpy as np
 import matplotlib.pyplot as plt
 from pylab import norm
 from SimPEG.utils import mkvc
-from SimPEG import TensorMesh, utils, LogicallyOrthogonalMesh
+from SimPEG import utils
+from SimPEG.mesh import TensorMesh, LogicallyOrthogonalMesh
 import numpy as np
 import unittest
 
diff --git a/SimPEG/tests/test_LogicallyOrthogonalMesh.py b/SimPEG/tests/test_LogicallyOrthogonalMesh.py
index 523b18d2..7f656aae 100644
--- a/SimPEG/tests/test_LogicallyOrthogonalMesh.py
+++ b/SimPEG/tests/test_LogicallyOrthogonalMesh.py
@@ -1,10 +1,7 @@
 import numpy as np
 import unittest
-import sys
-sys.path.append('../')
-from TensorMesh import TensorMesh
-from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
-from utils import ndgrid
+from SimPEG.mesh import TensorMesh, LogicallyOrthogonalMesh
+from SimPEG.utils import ndgrid
 
 
 class BasicLOMTests(unittest.TestCase):
diff --git a/SimPEG/tests/test_basemesh.py b/SimPEG/tests/test_basemesh.py
index 60373011..5762de57 100644
--- a/SimPEG/tests/test_basemesh.py
+++ b/SimPEG/tests/test_basemesh.py
@@ -1,7 +1,6 @@
 import unittest
 import sys
-sys.path.append('../')
-from BaseMesh import BaseMesh
+from SimPEG.mesh import BaseMesh
 import numpy as np
 
 
diff --git a/SimPEG/tests/test_forward_DCproblem.py b/SimPEG/tests/test_forward_DCproblem.py
index b5e6d844..6ffa4650 100644
--- a/SimPEG/tests/test_forward_DCproblem.py
+++ b/SimPEG/tests/test_forward_DCproblem.py
@@ -1,6 +1,6 @@
 import numpy as np
 import unittest
-from SimPEG import TensorMesh
+from SimPEG.mesh import TensorMesh
 from SimPEG.utils import ModelBuilder, sdiag
 from SimPEG.forward import Problem, SyntheticProblem
 from SimPEG.forward.DCProblem import DCProblem, DCutils
diff --git a/SimPEG/tests/test_forward_problem.py b/SimPEG/tests/test_forward_problem.py
index b7180dca..d913a697 100644
--- a/SimPEG/tests/test_forward_problem.py
+++ b/SimPEG/tests/test_forward_problem.py
@@ -1,6 +1,6 @@
 import numpy as np
 import unittest
-from SimPEG import TensorMesh
+from SimPEG.mesh import TensorMesh
 from SimPEG.forward import Problem
 from TestUtils import checkDerivative
 from scipy.sparse.linalg import dsolve
diff --git a/SimPEG/tests/test_tensorMesh.py b/SimPEG/tests/test_tensorMesh.py
index 300e6634..9544dab6 100644
--- a/SimPEG/tests/test_tensorMesh.py
+++ b/SimPEG/tests/test_tensorMesh.py
@@ -1,6 +1,6 @@
 import numpy as np
 import unittest
-from SimPEG import TensorMesh
+from SimPEG.mesh import TensorMesh
 from TestUtils import OrderTest
 from scipy.sparse.linalg import dsolve
 
diff --git a/SimPEG/utils/ModelBuilder.py b/SimPEG/utils/ModelBuilder.py
index 170da1bf..d0512680 100644
--- a/SimPEG/utils/ModelBuilder.py
+++ b/SimPEG/utils/ModelBuilder.py
@@ -15,6 +15,7 @@ def getIndecesBlock(p0,p1,ccMesh):
         ccMesh represents the cell-centered mesh
 
         The points p0 and p1 must live in the the same dimensional space as the mesh.
+
     """
 
     # Validation: p0 and p1 live in the same dimensional space
@@ -131,7 +132,7 @@ def scalarConductivity(ccMesh,pFunction):
 
 if __name__ == '__main__':
 
-    from SimPEG import TensorMesh
+    from SimPEG.mesh import TensorMesh
     from matplotlib import pyplot as plt
 
     # Define the mesh
diff --git a/docs/api_BaseMesh.rst b/docs/api_BaseMesh.rst
index 66249953..e5f81299 100644
--- a/docs/api_BaseMesh.rst
+++ b/docs/api_BaseMesh.rst
@@ -3,6 +3,6 @@
 Base Mesh
 *********
 
-.. automodule:: SimPEG.BaseMesh
+.. automodule:: SimPEG.mesh.BaseMesh
     :members:
     :undoc-members:
diff --git a/docs/api_DiffOperators.rst b/docs/api_DiffOperators.rst
index 57763a33..022fe564 100644
--- a/docs/api_DiffOperators.rst
+++ b/docs/api_DiffOperators.rst
@@ -3,6 +3,6 @@
 Differential Operators
 **********************
 
-.. automodule:: SimPEG.DiffOperators
+.. automodule:: SimPEG.mesh.DiffOperators
     :members:
     :undoc-members:
diff --git a/docs/api_InnerProducts.rst b/docs/api_InnerProducts.rst
index 6dfd2a77..90642807 100644
--- a/docs/api_InnerProducts.rst
+++ b/docs/api_InnerProducts.rst
@@ -3,6 +3,6 @@
 Inner Products
 **************
 
-.. automodule:: SimPEG.InnerProducts
+.. automodule:: SimPEG.mesh.InnerProducts
     :members:
     :undoc-members:
diff --git a/docs/api_LOMView.rst b/docs/api_LOMView.rst
index 61630c26..cfbfa6fa 100644
--- a/docs/api_LOMView.rst
+++ b/docs/api_LOMView.rst
@@ -3,6 +3,6 @@
 LOM View
 ********
 
-.. automodule:: SimPEG.LomView
+.. automodule:: SimPEG.mesh.LomView
     :members:
     :undoc-members:
diff --git a/docs/api_LogicallyOrthogonalMesh.rst b/docs/api_LogicallyOrthogonalMesh.rst
index 3b8296d7..af9e73d0 100644
--- a/docs/api_LogicallyOrthogonalMesh.rst
+++ b/docs/api_LogicallyOrthogonalMesh.rst
@@ -3,6 +3,6 @@
 Logically Orthogonal Mesh
 *************************
 
-.. automodule:: SimPEG.LogicallyOrthogonalMesh
+.. automodule:: SimPEG.mesh.LogicallyOrthogonalMesh
     :members:
     :undoc-members:
diff --git a/docs/api_TensorMesh.rst b/docs/api_TensorMesh.rst
index 9eb346cb..36fb925d 100644
--- a/docs/api_TensorMesh.rst
+++ b/docs/api_TensorMesh.rst
@@ -3,6 +3,6 @@
 Tensor Mesh
 ***********
 
-.. automodule:: SimPEG.TensorMesh
+.. automodule:: SimPEG.mesh.TensorMesh
     :members:
     :undoc-members:
diff --git a/docs/api_TensorView.rst b/docs/api_TensorView.rst
index 6ae80276..c6b5a9b8 100644
--- a/docs/api_TensorView.rst
+++ b/docs/api_TensorView.rst
@@ -3,6 +3,6 @@
 Tensor View
 ***********
 
-.. automodule:: SimPEG.TensorView
+.. automodule:: SimPEG.mesh.TensorView
     :members:
     :undoc-members: