diff --git a/SimPEG/EM/Static/DC/FieldsDC_2D.py b/SimPEG/EM/Static/DC/FieldsDC_2D.py
new file mode 100644
index 00000000..77e3199e
--- /dev/null
+++ b/SimPEG/EM/Static/DC/FieldsDC_2D.py
@@ -0,0 +1,106 @@
+import SimPEG
+import Utils, numpy as np, scipy.sparse as sp
+
+class Fields_ky(SimPEG.Problem.TimeFields):
+
+    """
+
+    Fancy Field Storage for a 2.5D code.
+
+    u[:,'phi', kyInd] = phi
+    print u[src0,'phi']
+
+    Only one field type is stored for
+    each problem, the rest are computed. The fields obejct acts like an array and is indexed by
+    .. code-block:: python
+        f = problem.fields(m)
+        e = f[srcList,'e']
+        j = f[srcList,'j']
+
+    If accessing all sources for a given field, use the :code:`:`
+    .. code-block:: python
+        f = problem.fields(m)
+        phi = f[:,'phi']
+        e = f[:,'e']
+        b = f[:,'b']
+    The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies)
+    """
+
+    knownFields = {}
+    dtype = float
+
+    def _phiDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
+        if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
+            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
+
+        if adjoint:
+            return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint)
+
+        return np.array(self._phiDeriv_u(kyInd, src, du_dm_v, adjoint) + self._phiDeriv_m(kyInd, src, v, adjoint), dtype = float)
+
+    def _eDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
+        if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
+
+        if adjoint:
+            return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint)
+        return np.array(self._eDeriv_u(kyInd, src, du_dm_v, adjoint) + self._eDeriv_m(kyInd, src, v, adjoint), dtype = float)
+
+    def _jDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
+        if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
+
+        if adjoint:
+            return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)
+        return np.array(self._jDeriv_u(kyInd, src, du_dm_v, adjoint) + self._jDeriv_m(kyInd, src, v, adjoint), dtype = float)
+
+
+    # def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+    #     if adjoint is True:
+    #         return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint)
+    #     return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v)
+
+    # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+    #     if adjoint is True:
+    #         return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
+    #     return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
+
+
+class Fields_ky_CC(Fields_ky):
+    knownFields = {'phiSolution':'CC'}
+    aliasFields = {
+                    'phi': ['phiSolution','CC','_phi'],
+                    'j' : ['phiSolution','F','_j'],
+                    'e' : ['phiSolution','F','_e'],
+                  }
+                  # primary - secondary
+                  # CC variables
+
+    def __init__(self, mesh, survey, **kwargs):
+        Fields_ky.__init__(self, mesh, survey, **kwargs)
+
+    def startup(self):
+        self.prob = self.survey.prob
+
+    def _GLoc(self, fieldType):
+        if fieldType == 'phi':
+            return 'CC'
+        elif fieldType == 'e' or fieldType == 'j':
+            return 'F'
+        else:
+            raise Exception('Field type must be phi, e, j')
+
+    def _phi(self, phiSolution, src, kyInd):
+        return phiSolution
+
+    def _phiDeriv_u(self, kyInd, src, v, adjoint = False):
+        return Identity()*v
+
+    def _phiDeriv_m(self, kyInd, src, v, adjoint = False):
+        return Zero()
+
+    def _j(self, phiSolution, srcList):
+        raise NotImplementedError
+
+    def _e(self, phiSolution, srcList):
+        raise NotImplementedError
diff --git a/SimPEG/EM/Static/DC/ProblemDC_2D.py b/SimPEG/EM/Static/DC/ProblemDC_2D.py
new file mode 100644
index 00000000..ee04a560
--- /dev/null
+++ b/SimPEG/EM/Static/DC/ProblemDC_2D.py
@@ -0,0 +1,235 @@
+from SimPEG import Problem, Utils
+from SimPEG.EM.Base import BaseEMProblem
+from SurveyDC import Survey, Survey_ky
+from FieldsDC_2D import Fields_ky, Fields_ky_CC
+from SimPEG.Utils import sdiag
+import numpy as np
+from SimPEG.Utils import Zero
+from BoundaryUtils import getxBCyBC_CC
+
+class BaseDCProblem_2D(BaseEMProblem):
+
+    surveyPair = Survey_ky
+    fieldsPair = Fields_ky
+    nky = 15
+    ky = np.logspace(-4, 1, nky)
+    Ainv = [None for i in range(nky)]
+    nT = nky # Only for using TimeFields
+
+    def fields(self, m):
+        self.curModel = m
+
+        if not self.Ainv[0] == None:
+            for i in range(self.nky):
+                self.Ainv[i].clean()
+
+        f = self.fieldsPair(self.mesh, self.survey)
+        Srcs = self.survey.srcList
+        for iky in range(self.nky):
+            ky = self.ky[iky]
+            A = self.getA(ky)
+            self.Ainv[iky] = self.Solver(A, **self.solverOpts)
+            RHS = self.getRHS(ky)
+            u = self.Ainv[iky] * RHS
+            f[Srcs, self._solutionType, iky] = u
+        return f
+
+    # def Jvec(self, m, v, f=None):
+
+    #     if f is None:
+    #         f = self.fields(m)
+
+    #     self.curModel = m
+
+    #     Jv = self.dataPair(self.survey) #same size as the data
+
+    #     A = self.getA()
+
+    #     for src in self.survey.srcList:
+    #         u_src = f[src, self._solutionType] # solution vector
+    #         dA_dm_v = self.getADeriv(u_src, v)
+    #         dRHS_dm_v = self.getRHSDeriv(src, v)
+    #         du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
+
+    #         for rx in src.rxList:
+    #             df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+    #             df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
+    #             Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
+    #     return Utils.mkvc(Jv)
+
+    # def Jtvec(self, m, v, f=None):
+    #     if f is None:
+    #         f = self.fields(m)
+
+    #     self.curModel = m
+
+    #     # Ensure v is a data object.
+    #     if not isinstance(v, self.dataPair):
+    #         v = self.dataPair(self.survey, v)
+
+    #     Jtv = np.zeros(m.size)
+    #     AT = self.getA()
+
+
+    #     for src in self.survey.srcList:
+    #         u_src = f[src, self._solutionType]
+    #         for rx in src.rxList:
+    #             PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
+    #             df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+    #             df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
+
+    #             ATinvdf_duT = self.Ainv * df_duT
+
+    #             dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
+    #             dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True)
+    #             du_dmT = -dA_dmT + dRHS_dmT
+    #             Jtv += df_dmT + du_dmT
+
+    #     return Utils.mkvc(Jtv)
+
+    def getSourceTerm(self, ky):
+        """
+        takes concept of source and turns it into a matrix
+        """
+        """
+        Evaluates the sources, and puts them in matrix form
+
+        :rtype: (numpy.ndarray, numpy.ndarray)
+        :return: q (nC or nN, nSrc)
+        """
+
+        Srcs = self.survey.srcList
+
+        if self._formulation is 'EB':
+            n = self.mesh.nN
+            # return NotImplementedError
+
+        elif self._formulation is 'HJ':
+            n = self.mesh.nC
+
+        q = np.zeros((n, len(Srcs)))
+
+        for i, src in enumerate(Srcs):
+            q[:,i] = src.eval(self)
+        return q
+
+class Problem2D_CC(BaseDCProblem_2D):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'HJ' # CC potentials means J is on faces
+    fieldsPair    = Fields_ky_CC
+
+    def __init__(self, mesh, **kwargs):
+        BaseDCProblem_2D.__init__(self, mesh, **kwargs)
+        self.setBC()
+
+    def getA(self, ky):
+        """
+
+        Make the A matrix for the cell centered DC resistivity problem
+
+        A = D MfRhoI D^\\top V
+
+        """
+
+        D = self.Div
+        G = self.Grad
+        vol = self.mesh.vol
+        # TODO: this won't work for full anisotropy
+        MfRhoI = self.MfRhoI
+        # Get resistivity rho
+        rho = self.curModel.rho
+        A = D * MfRhoI * G + Utils.sdiag(ky**2*vol/rho)
+        return A
+
+    def getADeriv(self, ky, u, v, adjoint= False):
+
+        D = self.Div
+        G = self.Grad
+        MfRhoIDeriv = self.MfRhoIDeriv
+
+        if adjoint:
+            return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + Utils.sdiag(ky**2*mesh.vol)*v
+        return D * ((MfRhoIDeriv( G * u )) * v) + Utils.sdiag(ky**2*mesh.vol)*v
+
+    def getRHS(self, ky):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm(ky)
+        return RHS
+
+    def getRHSDeriv(self, ky, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+        """
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, ky, adjoint=adjoint)
+        # return qDeriv
+        return Zero()
+
+    def setBC(self):
+        if self.mesh.dim==3:
+            fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd
+            gBFxm = self.mesh.gridFx[fxm,:]
+            gBFxp = self.mesh.gridFx[fxp,:]
+            gBFym = self.mesh.gridFy[fym,:]
+            gBFyp = self.mesh.gridFy[fyp,:]
+            gBFzm = self.mesh.gridFz[fzm,:]
+            gBFzp = self.mesh.gridFz[fzp,:]
+
+            # Setup Mixed B.C (alpha, beta, gamma)
+            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+            temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
+
+            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
+            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
+            alpha_zm, alpha_zp = temp_zm*0., temp_zp*0.
+
+            beta_xm, beta_xp = temp_xm, temp_xp
+            beta_ym, beta_yp = temp_ym, temp_yp
+            beta_zm, beta_zp = temp_zm, temp_zp
+
+            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
+            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
+            gamma_zm, gamma_zp = temp_zm*0., temp_zp*0.
+
+            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
+            beta =  [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
+            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
+
+        elif self.mesh.dim==2:
+
+            fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd
+            gBFxm = self.mesh.gridFx[fxm,:]
+            gBFxp = self.mesh.gridFx[fxp,:]
+            gBFym = self.mesh.gridFy[fym,:]
+            gBFyp = self.mesh.gridFy[fyp,:]
+
+            # Setup Mixed B.C (alpha, beta, gamma)
+            temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+            temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+
+            alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
+            alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
+
+            beta_xm, beta_xp = temp_xm, temp_xp
+            beta_ym, beta_yp = temp_ym, temp_yp
+
+            gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
+            gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
+
+            alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
+            beta =  [beta_xm, beta_xp, beta_ym, beta_yp]
+            gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
+
+        x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma)
+        V = self.Vol
+        self.Div = V * self.mesh.faceDiv
+        P_BC, B = self.mesh.getBCProjWF_simple()
+        M = B*self.mesh.aveCC2F
+        self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M
diff --git a/SimPEG/EM/Static/DC/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py
index 7d50392e..f7c7d352 100644
--- a/SimPEG/EM/Static/DC/RxDC.py
+++ b/SimPEG/EM/Static/DC/RxDC.py
@@ -1,5 +1,5 @@
 import SimPEG
-# from SimPEG.EM.Base import BaseEMSurvey
+import numpy as np
 from SimPEG.Utils import Zero, closestPoints
 
 class BaseRx(SimPEG.Survey.BaseRx):
@@ -43,9 +43,6 @@ class BaseRx(SimPEG.Survey.BaseRx):
         elif adjoint:
             return P.T*v
 
-
-
-
 # DC.Rx.Dipole(locs)
 class Dipole(BaseRx):
 
@@ -77,6 +74,56 @@ class Dipole(BaseRx):
 
         return P
 
-# class Pole(BaseRx):
 
+class Dipole_ky(BaseRx):
+
+    def __init__(self, locsM, locsN, rxType = 'phi', **kwargs):
+        assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size'
+        locs = [locsM, locsN]
+        # We may not need this ...
+        BaseRx.__init__(self, locs, rxType)
+
+    @property
+    def nD(self):
+        """Number of data in the receiver."""
+        return self.locs[0].shape[0]
+
+        # Not sure why ...
+        # return int(self.locs[0].size / 2)
+
+    def getP(self, mesh, Gloc):
+        if mesh in self._Ps:
+            return self._Ps[mesh]
+
+        P0 = mesh.getInterpolationMat(self.locs[0], Gloc)
+        P1 = mesh.getInterpolationMat(self.locs[1], Gloc)
+        P = P0 - P1
+        if self.storeProjections:
+            self._Ps[mesh] = P
+        return P
+
+    def eval(self, ky, src, mesh, f):
+        P = self.getP(mesh, self.projGLoc(f))
+        Pf = P*f[src, self.projField,:]
+        return self.IntTrapezoidal(ky, Pf, y=0.)
+
+    def evalDeriv(self, ky, src, mesh, f, v, adjoint=False):
+        P = self.getP(mesh, self.projGLoc(f))
+        if not adjoint:
+            return P*v
+        elif adjoint:
+            return P.T*v
+
+    def IntTrapezoidal(self, ky, Pf, y=0.):
+        phi = np.zeros(Pf.shape[0])
+        nky = ky.size
+        dky = np.diff(ky)
+        dky = np.r_[dky[0], dky]
+        phi0 = Pf[:,0]
+        for iky in range(nky):
+            phi1 = 2./np.pi*Pf[:,iky]/2.
+            phi += phi1*dky[iky]/2.*np.cos(ky[iky]*y)
+            phi += phi0*dky[iky]/2.*np.cos(ky[iky]*y)
+            phi0 = phi1.copy()
+        return phi
 
diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py
index 7879ba0c..60a53dff 100644
--- a/SimPEG/EM/Static/DC/SrcDC.py
+++ b/SimPEG/EM/Static/DC/SrcDC.py
@@ -57,5 +57,3 @@ class Pole(BaseSrc):
             q[inds] = self.current * np.r_[1.]
         return q
 
-    # def bc_contribution
-
diff --git a/SimPEG/EM/Static/DC/SurveyDC.py b/SimPEG/EM/Static/DC/SurveyDC.py
index 3b631bef..62e2922a 100644
--- a/SimPEG/EM/Static/DC/SurveyDC.py
+++ b/SimPEG/EM/Static/DC/SurveyDC.py
@@ -1,6 +1,6 @@
 import SimPEG
 from SimPEG.EM.Base import BaseEMSurvey
-from SimPEG import sp
+from SimPEG import sp, Survey
 from SimPEG.Utils import Zero, Identity
 from RxDC import BaseRx
 from SrcDC import BaseSrc
@@ -13,6 +13,24 @@ class Survey(BaseEMSurvey):
         self.srcList = srcList
         BaseEMSurvey.__init__(self, srcList, **kwargs)
 
+class Survey_ky(BaseEMSurvey):
+    rxPair  = BaseRx
+    srcPair = BaseSrc
 
+    def __init__(self, srcList, **kwargs):
+        self.srcList = srcList
+        BaseEMSurvey.__init__(self, srcList, **kwargs)
 
-
+    def eval(self, f):
+        """
+        Project fields to receiver locations
+        :param Fields u: fields object
+        :rtype: numpy.ndarray
+        :return: data
+        """
+        data = SimPEG.Survey.Data(self)
+        ky = self.prob.ky
+        for src in self.srcList:
+            for rx in src.rxList:
+                data[src, rx] = rx.eval(ky, src, self.mesh, f)
+        return data
diff --git a/SimPEG/EM/Static/DC/__init__.py b/SimPEG/EM/Static/DC/__init__.py
index 82cf76d5..8c790084 100644
--- a/SimPEG/EM/Static/DC/__init__.py
+++ b/SimPEG/EM/Static/DC/__init__.py
@@ -1,5 +1,6 @@
 from ProblemDC import Problem3D_CC, Problem3D_N
-from SurveyDC import Survey
+from ProblemDC_2D import Problem2D_CC
+from SurveyDC import Survey, Survey_ky
 import SrcDC as Src #Pole
 import RxDC as Rx
 from FieldsDC import Fields_CC