diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py
index 9a8a9db9..d50cc5ba 100644
--- a/SimPEG/DiffOperators.py
+++ b/SimPEG/DiffOperators.py
@@ -1,6 +1,6 @@
 import numpy as np
 from scipy import sparse as sp
-from utils import mkvc, sdiag, speye, kron3, spzeros
+from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros
 
 
 def ddx(n):
@@ -287,15 +287,19 @@ class DiffOperators(object):
     nodalVectorAve = property(**nodalVectorAve())
 
     def getEdgeMass(self, materialProp=None):
-        """mass matix for products of edge functions w'*M(materialProp)*e"""
+        """mass matrix for products of edge functions w'*M(materialProp)*e"""
         if(materialProp is None):
             materialProp = np.ones(self.nC)
         Av = self.edgeAve
         return sdiag(Av.T * (self.vol * mkvc(materialProp)))
 
     def getFaceMass(self, materialProp=None):
-        """mass matix for products of edge functions w'*M(materialProp)*e"""
+        """mass matrix for products of face functions w'*M(materialProp)*f"""
         if(materialProp is None):
             materialProp = np.ones(self.nC)
         Av = self.faceAve
-        return sdiag(Av.T*(self.vol*mkvc(materialProp)))
+        return sdiag(Av.T * (self.vol * mkvc(materialProp)))
+
+    def getFaceMassDeriv(self):
+        Av = self.faceAve
+        return Av.T * sdiag(self.vol)
diff --git a/SimPEG/InnerProducts.py b/SimPEG/InnerProducts.py
index fca632f4..9fe84ac3 100644
--- a/SimPEG/InnerProducts.py
+++ b/SimPEG/InnerProducts.py
@@ -1,5 +1,5 @@
 from scipy import sparse as sp
-from utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
+from SimPEG.utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
 import numpy as np
 
 
diff --git a/SimPEG/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem/DCProblem.py
new file mode 100644
index 00000000..b94b35b5
--- /dev/null
+++ b/SimPEG/forward/DCProblem/DCProblem.py
@@ -0,0 +1,93 @@
+from SimPEG import TensorMesh
+from SimPEG.forward import Problem, SyntheticProblem
+from SimPEG.utils import ModelBuilder
+import numpy as np
+import scipy.sparse.linalg as linalg
+import DCutils
+
+class DCProblem(Problem):
+    """docstring for DCProblem"""
+    def __init__(self, mesh):
+        super(DCProblem, self).__init__(mesh)
+        self.mesh.setCellGradBC('neumann')
+
+    def createMatrix(self, m):
+        D = self.mesh.faceDiv
+        G = self.mesh.cellGrad
+        sigma = self.modelTransform(m)
+        Msig = self.mesh.getFaceMass(sigma)
+        A = D*Msig*G
+        return A.tocsc()
+
+    def field(self, m):
+        A = self.createMatrix(m)
+        solve = linalg.factorized(A)
+
+        nRHSs = self.RHS.shape[1]  # Number of RHSs
+        phi = np.zeros((self.mesh.nC, nRHSs)) + np.nan
+        for ii in range(nRHSs):
+            phi[:,ii] = solve(self.RHS[:,ii])
+
+        return phi
+
+    def J(self, m, v, u=None, RHSii=0, solve=None):
+        P = self.P
+        D = self.mesh.faceDiv
+        G = self.mesh.cellGrad
+        A = self.createMatrix(m)
+        Av_dm = self.mesh.getFaceMassDeriv()
+        mT_dm = self.modelTransform(m)
+
+        dCdu = A
+        dCdm = - D * ( sdiag( G * u[:, RHSii] ) * ( Av_dm * ( mT_dm * v ) ) )
+
+        if solve is None:
+            solve = linalg.factorized(dCdu)
+
+        return - P * solve(dCdm)
+
+
+
+if __name__ == '__main__':
+    # Create the mesh
+    h1 = np.ones(100)
+    h2 = np.ones(100)
+    mesh = TensorMesh([h1,h2])
+
+    # Create some parameters for the model
+    sig1 = 1
+    sig2 = 0.01
+
+    # Create a synthetic model from a block in a half-space
+    p0 = [20, 20]
+    p1 = [50, 50]
+    condVals = [sig1, sig2]
+    mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals)
+    mesh.plotImage(mSynth, showIt=False)
+
+
+    # Set up the projection
+    nelec = 50
+    spacelec = 2
+    surfloc = 0.5
+    elecini = 0.5
+    elecend = 0.5+spacelec*(nelec-1)
+    elecLocR = np.linspace(elecini, elecend, nelec)
+    rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
+    q, Q, rxmidloc = DCutils.genTxRxmat(nelec, spacelec, surfloc, elecini, mesh)
+
+
+    # Create some data
+    class syntheticDCProblem(DCProblem, SyntheticProblem):
+        pass
+
+    synthetic = syntheticDCProblem(mesh);
+    synthetic.P = Q.T
+    synthetic.RHS = q
+    dobs, Wd = synthetic.createData(mSynth)
+
+    # Now set up the problem to do some minimization
+    problem = DCProblem(mesh)
+
+
+
diff --git a/SimPEG/forward/DCProblem/DCutils.py b/SimPEG/forward/DCProblem/DCutils.py
new file mode 100644
index 00000000..f3445096
--- /dev/null
+++ b/SimPEG/forward/DCProblem/DCutils.py
@@ -0,0 +1,29 @@
+import numpy as np
+import scipy.sparse as sp
+
+def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh):
+    """ Generate projection matrix (Q) and """
+    elecend = 0.5+spacelec*(nelec-1)
+    elecLocR = np.linspace(elecini, elecend, nelec)
+    elecLocT = elecLocR+1
+    nrx = nelec-1
+    ntx = nelec-1
+    q = np.zeros((mesh.nC, ntx))
+    Q = np.zeros((mesh.nC, nrx))
+
+    for i in range(nrx):
+
+        rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i]))
+        rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1]))
+
+        txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i]))
+        txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1]))
+
+        q[txind1,i] = 1
+        q[txind2,i] = -1
+        Q[rxind1,i] = 1
+        Q[rxind2,i] = -1
+
+    Q = sp.csr_matrix(Q)
+    rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
+    return q, Q, rxmidLoc
diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py
index 48848a07..ff8bdb00 100644
--- a/SimPEG/forward/Problem.py
+++ b/SimPEG/forward/Problem.py
@@ -84,8 +84,15 @@ class Problem(object):
         self._P = value
 
 
-    def J(self, u):
+    def J(self, m, v, u=None, RHSii=0):
         """
+            :param numpy.array m: model
+            :param numpy.array v: vector to multiply
+            :param numpy.array u: fields
+            :param int RHSii: which RHS to calculate sensitivity too
+            :rtype: numpy.array
+            :return: Jv
+
             Working with the general PDE, c(m, u) = 0, where m is the model and u is the field,
             the sensitivity is defined as:
 
@@ -107,15 +114,26 @@ class Problem(object):
         """
         pass
 
-    def Jt(self, v):
+    def Jt(self, m, v, u=None, RHSii=0):
         """
+            :param numpy.array m: model
+            :param numpy.array v: vector to multiply
+            :param numpy.array u: fields
+            :param int RHSii: which RHS to calculate sensitivity too
+            :rtype: numpy.array
+            :return: JTv
+
             Transpose of J
         """
         pass
 
     def field(self, m):
         """
-            The fields.
+            The field given the model.
+
+            .. math::
+                u(m)
+
         """
         pass
 
@@ -179,10 +197,10 @@ class Problem(object):
         m = np.random.rand(5)
         return checkDerivative(lambda m : [self.modelTransform(m), self.modelTransformDeriv(m)], m)
 
-    def misfit(self, m, R=None):
+    def misfit(self, m, u=None):
         """
             :param numpy.array m: geophysical model
-            :param numpy.array R: residual, R = W o (dpred - dobs)
+            :param numpy.array u: fields
             :rtype: float
             :return: data misfit
 
@@ -195,15 +213,15 @@ class Problem(object):
             Where P is a projection matrix that brings the field on the full domain to the data measurement locations;
             u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
         """
-        if R is None:
-            R = self.W*(self.dpred(m) - self.dobs)
 
+        R = self.W*(self.dpred(m, u=u) - self.dobs)
         R = mkvc(R)
         return 0.5*R.inner(R)
 
-    def misfitDeriv(self, m, R=None, u=None):
+    def misfitDeriv(self, m, u=None):
         """
             :param numpy.array m: geophysical model
+            :param numpy.array u: fields
             :rtype: numpy.array
             :return: data misfit derivative
 
@@ -213,6 +231,12 @@ class Problem(object):
 
                 \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2
 
+            If the field, u, is provided, the calculation of the data is fast:
+
+            .. math::
+
+                \mathbf{d}_\\text{pred} = \mathbf{Pu(m)}
+
                 \mathbf{R} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}
 
                 \mu_\\text{data} = {1\over 2}\left| \mathbf{W \circ R} \\right|_2^2
@@ -230,8 +254,7 @@ class Problem(object):
         if u is None:
             u = self.field(m)
 
-        if R is None:
-            R = self.W*(self.dpred(m, u=u) - self.dobs)
+        R = self.W*(self.dpred(m, u=u) - self.dobs)
 
         dmisfit = 0
         for i in range(self.RHS.shape[1]): # Loop over each right hand side
@@ -240,9 +263,40 @@ class Problem(object):
         return dmisfit
 
 
+class SyntheticProblem(object):
+    """
+        Has helpful functions when dealing with synthetic problems
+
+        To use this class, inherit to your problem::
+
+            class mySyntheticExample(Problem, SyntheticProblem):
+                pass
+    """
+    def createData(self, m, std=0.05):
+        """
+            :param numpy.array m: geophysical model
+            :param numpy.array std: standard deviation
+            :rtype: numpy.array, numpy.array
+            :return: dobs, Wd
+
+            Create synthetic data given a model, and a standard deviation.
+
+            Returns the observed data with random Gaussian noise
+            and Wd which is the same size as data, and can be used to weight the inversion.
+        """
+        dobs = self.dpred(m)
+        dobs = dobs
+        noise = std*abs(dobs)*np.random.randn(*dobs.shape)
+        dobs = dobs+noise
+        eps = np.linalg.norm(mkvc(dobs),2)*1e-5
+        Wd = 1/(abs(dobs)*std+eps)
+        return dobs, Wd
+
+
+
 if __name__ == '__main__':
     from SimPEG.inverse import checkDerivative
 
     p = Problem(None)
     m = np.random.rand(5)
-    checkDerivative(lambda m : [p.modelTransform(m), p.modelTransformDeriv(m)], m)
+    checkDerivative(lambda m : [p.modelTransform(m), p.modelTransformDeriv(m)], m, plotIt=False)
diff --git a/SimPEG/utils/ModelBuilder.py b/SimPEG/utils/ModelBuilder.py
index 3b1f977f..527d5eef 100644
--- a/SimPEG/utils/ModelBuilder.py
+++ b/SimPEG/utils/ModelBuilder.py
@@ -15,10 +15,6 @@ def getIndecesBlock(p0,p1,ccMesh):
     The points p0 and p1 must live in the the same dimensional space as the mesh.
     """
 
-    # Validation of the input
-    assert type(p0) == np.ndarray, "Vector must be a numpy array"
-    assert type(p1) == np.ndarray, "Vector must be a numpy array"
-
     # Validation: p0 and p1 live in the same dimensional space
     assert len(p0) == len(p1), "Dimension mismatch. len(p0) != len(p1)"
 
@@ -47,7 +43,7 @@ def getIndecesBlock(p0,p1,ccMesh):
 
         ind  = np.where(indX & indY)
 
-    else:
+    elif dimMesh == 3:
         # Define the points
         x1 = p0[0]
         y1 = p0[1]
@@ -98,13 +94,16 @@ def defineTwoLayeredConductivity(depth,ccMesh,condVals):
 
     # Identify 1st cell centered reference point
     p0[0] = ccMesh[0,0]
-    p0[1] = ccMesh[0,1]
-    p0[2] = ccMesh[0,2]
+    if dim>1: p0[1] = ccMesh[0,1]
+    if dim>2: p0[2] = ccMesh[0,2]
 
     # Identify the last cell-centered reference point
     p1[0] = ccMesh[-1,0]
-    p1[1] = ccMesh[-1,1]
-    p1[2] = ccMesh[-1,2] - depth;
+    if dim>1: p1[1] = ccMesh[-1,1]
+    if dim>2: p1[2] = ccMesh[-1,2]
+
+    # The depth is always defined on the last one.
+    p1[len(p1)-1] -= depth
 
     ind   = getIndecesBlock(p0,p1,ccMesh)
 
@@ -117,23 +116,24 @@ def scalarConductivity(ccMesh,pFunction):
     Define the distribution conductivity in the mesh according to the
     analytical expression given in pFunction
     """
-    xCC = ccMesh[:,0]
-    yCC = ccMesh[:,1]
-    zCC = ccMesh[:,2]
+    dim = np.size(ccMesh[0,:])
+    CC = [ccMesh[:,0]]
+    if dim>1: CC.append(ccMesh[:,1])
+    if dim>2: CC.append(ccMesh[:,2])
 
-    sigma = pFunction(xCC,yCC,zCC)
+
+    sigma = pFunction(*CC)
 
     return sigma
 
 if __name__ == '__main__':
 
-    import sys
-    sys.path.append('../')
-    from TensorMesh import TensorMesh
+    from SimPEG import TensorMesh
+    from matplotlib import pyplot as plt
 
     # Define the mesh
 
-    testDim = 3
+    testDim = 2
     h1 = 0.3*np.ones(7)
     h1[0] = 0.5
     h1[-1] = 0.6
@@ -157,8 +157,8 @@ if __name__ == '__main__':
     # ------------------- Test conductivities! --------------------------
     print('Testing 1 block conductivity')
 
-    p0 = np.array([0.5,0.5,0.5])
-    p1 = np.array([1.0,1.0,1.0])
+    p0 = np.array([0.5,0.5,0.5])[:testDim]
+    p1 = np.array([1.0,1.0,1.0])[:testDim]
     condVals = np.array([100,1e-6])
 
     sigma = defineBlockConductivity(p0,p1,ccMesh,condVals)
@@ -167,6 +167,7 @@ if __name__ == '__main__':
     print sigma.shape
     M.plotImage(sigma)
     print 'Done with block! :)'
+    plt.show()
 
     # -----------------------------------------
     print('Testing the two layered model')
@@ -178,11 +179,17 @@ if __name__ == '__main__':
     M.plotImage(sigma)
     print sigma
     print 'layer model!'
+    plt.show()
 
     # -----------------------------------------
     print('Testing scalar conductivity')
 
-    pFunction = lambda x,y,z: np.exp(x+y+z)
+    if testDim == 1:
+        pFunction = lambda x: np.exp(x)
+    elif testDim == 2:
+        pFunction = lambda x,y: np.exp(x+y)
+    elif testDim == 3:
+        pFunction = lambda x,y,z: np.exp(x+y+z)
 
     sigma = scalarConductivity(ccMesh,pFunction)
 
@@ -190,5 +197,6 @@ if __name__ == '__main__':
     M.plotImage(sigma)
     print sigma
     print 'Scalar conductivity defined!'
+    plt.show()
 
     # -----------------------------------------