diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py
index c204d76a..e2c2ff80 100644
--- a/SimPEG/EM/FDEM/FDEM.py
+++ b/SimPEG/EM/FDEM/FDEM.py
@@ -8,29 +8,29 @@ from SimPEG.EM.Utils import omega
 
 class BaseFDEMProblem(BaseEMProblem):
     """
-        We start by looking at Maxwell's equations in the electric
-        field \\\(\\\mathbf{e}\\\) and the magnetic flux
-        density \\\(\\\mathbf{b}\\\)
+    We start by looking at Maxwell's equations in the electric
+    field \\\(\\\mathbf{e}\\\) and the magnetic flux
+    density \\\(\\\mathbf{b}\\\)
 
-        .. math ::
+    .. math ::
 
-            \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
-            {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
+        \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
+        {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
 
-        if using the E-B formulation (:code:`Problem_e`
-        or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
+    if using the E-B formulation (:code:`Problem_e`
+    or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
 
-        If we write Maxwell's equations in terms of
-        \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
+    If we write Maxwell's equations in terms of
+    \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
 
-        .. math ::
+    .. math ::
 
-            \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
-            \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
+        \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
+        \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
 
-        if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
+    if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
 
-        The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
+    The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
     """
 
     surveyPair = SurveyFDEM
@@ -204,6 +204,17 @@ class Problem_e(BaseFDEMProblem):
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
 
+    def _GLoc(self, fieldType):
+        if fieldType == 'e':
+            return 'E'
+        elif fieldType == 'b':
+            return 'F'
+        elif (fieldType == 'h') or (fieldType == 'j'):
+            return 'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
+
+
     def getA(self, freq):
         """
         System matrix
@@ -314,6 +325,16 @@ class Problem_b(BaseFDEMProblem):
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
 
+    def _GLoc(self, fieldType):
+        if fieldType == 'e':
+            return 'E'
+        elif fieldType == 'b':
+            return 'F'
+        elif (fieldType == 'h') or (fieldType == 'j'):
+            return'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
+
     def getA(self, freq):
         """
         System matrix
@@ -462,6 +483,16 @@ class Problem_j(BaseFDEMProblem):
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
 
+    def _GLoc(self, fieldType):
+        if fieldType == 'h':
+            return 'E'
+        elif fieldType == 'j':
+            return 'F'
+        elif (fieldType == 'e') or (fieldType == 'b'):
+            return 'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
+
     def getA(self, freq):
         """
         System matrix
@@ -600,6 +631,17 @@ class Problem_h(BaseFDEMProblem):
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
 
+    def _GLoc(self, fieldType):
+        if fieldType == 'h':
+            return 'E'
+        elif fieldType == 'j':
+            return 'F'
+        elif (fieldType == 'e') or (fieldType == 'b'):
+            return 'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
+
+
     def getA(self, freq):
         """
         System matrix
diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py
index e2193973..bfe7ae9e 100644
--- a/SimPEG/EM/FDEM/FieldsFDEM.py
+++ b/SimPEG/EM/FDEM/FieldsFDEM.py
@@ -193,16 +193,6 @@ class Fields_e(Fields):
         self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv
         self._MfMui = self.survey.prob.MfMui
 
-    def _GLoc(self, fieldType):
-        if fieldType == 'e':
-            return 'E'
-        elif fieldType == 'b':
-            return 'F'
-        elif (fieldType == 'h') or (fieldType == 'j'):
-            return 'CCV'
-        else:
-            raise Exception('Field type must be e, b, h, j')
-
 
     def _ePrimary(self, eSolution, srcList):
         """
@@ -465,17 +455,6 @@ class Fields_b(Fields):
         self._nC = self.survey.prob.mesh.nC
 
 
-
-    def _GLoc(self,fieldType):
-        if fieldType == 'e':
-            return 'E'
-        elif fieldType == 'b':
-            return 'F'
-        elif (fieldType == 'h') or (fieldType == 'j'):
-            return'CCV'
-        else:
-            raise Exception('Field type must be e, b, h, j')
-
     def _bPrimary(self, bSolution, srcList):
         """
         Primary magnetic flux density from source
@@ -729,16 +708,6 @@ class Fields_j(Fields):
         self._aveE2CCV = self.survey.prob.mesh.aveE2CCV
         self._nC = self.survey.prob.mesh.nC
 
-    def _GLoc(self,fieldType):
-        if fieldType == 'h':
-            return 'E'
-        elif fieldType == 'j':
-            return 'F'
-        elif (fieldType == 'e') or (fieldType == 'b'):
-            return 'CCV'
-        else:
-            raise Exception('Field type must be e, b, h, j')
-
     def _jPrimary(self, jSolution, srcList):
         """
         Primary current density from source
@@ -1024,16 +993,6 @@ class Fields_h(Fields):
         self._aveE2CCV = self.survey.prob.mesh.aveE2CCV
         self._nC = self.survey.prob.mesh.nC
 
-    def _GLoc(self,fieldType):
-        if fieldType == 'h':
-            return 'E'
-        elif fieldType == 'j':
-            return 'F'
-        elif (fieldType == 'e') or (fieldType == 'b'):
-            return 'CCV'
-        else:
-            raise Exception('Field type must be e, b, h, j')
-
     def _hPrimary(self, hSolution, srcList):
         """
         Primary magnetic field from source
diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py
index e8622eb2..874afba2 100644
--- a/SimPEG/EM/FDEM/SrcFDEM.py
+++ b/SimPEG/EM/FDEM/SrcFDEM.py
@@ -103,7 +103,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         return Zero()
 
-    def s_mDeriv(self, prob, v, adjoint = False):
+    def s_mDeriv(self, prob, v, adjoint=False):
         """
         Derivative of magnetic source term with respect to the inversion model
 
@@ -116,7 +116,7 @@ class BaseSrc(Survey.BaseSrc):
 
         return Zero()
 
-    def s_eDeriv(self, prob, v, adjoint = False):
+    def s_eDeriv(self, prob, v, adjoint=False):
         """
         Derivative of electric source term with respect to the inversion model
 
@@ -605,92 +605,231 @@ class CircularLoop(BaseSrc):
             return -C.T * (MMui_s * self.bPrimary(prob))
 
 
-class PrimSecSigma(BaseSrc):
+class PrimSec(BaseSrc):
+    """
+    Primary-Secondary source in the physical properties. A primary problem is
+    first solved, and the fields from this problem are used to construct a
+    source term for the secondary problem. Either a mesh and
+    fields need to be provided or a prob and a survey.
 
-    # TODO: This will only work for E-B formulation
-    def __init__(self, rxList, freq, mPrimary, primaryProblem=None, primarySurvey=None, primaryFields=None):
+    For the EB formulation, we start the derivation from Maxwell's equations:
+
+    .. math::
+        \\nabla \\times \\vec{E} + i \omega \\vec{B} = \\vec{s_m} \\\\
+        \\nabla \\times \\mu^{-1} \\vec{B} - \sigma \\vec{E} = \\vec{s_e}
+
+    we consider the physical properties, fields, and fluxes to be composed of
+    two parts, a primary and a secondary:
+
+        - :math:`\sigma   = \sigma_p + \sigma_s`
+        - :math:`\mu^{-1} = \mu^{-1}_p + \mu^{-1}_s`
+        - :math:`\\vec{E} = \\vec{E_p} + \\vec{E_s}`
+        - :math:`\\vec{B} = \\vec{B_p} + \\vec{B_s}`
+
+    and choose our primary such that
+
+    .. math::
+        \\nabla \\times \\vec{E}_p + i \omega \\vec{B}_p = \\vec{s_m} \\\\
+        \\nabla \\times \\mu^{-1}_p \\vec{B}_p - \sigma_p \\vec{E}_p = \\vec{s_e}_p
+
+    so the secondary problem is then
+
+    .. math::
+        \\nabla \\times \\vec{E}_s + i \omega \\vec{B}_s = 0 \\\\
+        \\nabla \\times \\mu^{-1} \\vec{B}_s - \sigma \\vec{E}_s = - \\nabla \\times \\mu^{-1}_s \\vec{B}_p + \sigma_s \\vec{E}_p
+
+
+    If instead, HJ formulation is considered, then we start off with
+
+    .. math::
+        \\nabla \\times \\rho \\vec{J} + i \omega \\mu \\vec{H} = \\vec{s_m} \\\\
+        \\nabla \\times \\vec{H} - \\vec{J} = \\vec{s_e}
+
+    and we define the primary secondary problem in terms of
+
+        - :math:`\\rho  = \\rho_p + \\rho_s`
+        - :math:`\mu = \mu_p + \mu_s`
+        - :math:`\\vec{J} = \\vec{J_p} + \\vec{J_s}`
+        - :math:`\\vec{H} = \\vec{H_p} + \\vec{H_s}`
+
+    with the primary being defined by
+
+    .. math::
+        \\nabla \\times \\rho_p \\vec{J}_p + i \omega \\mu_p \\vec{H}_p = \\vec{s_m} \\\\
+        \\nabla \\times \\vec{H}_p - \\vec{J}_p = \\vec{s_e}
+
+    so the secondary problem is given by
+
+    .. math::
+        \\nabla \\times \\rho \\vec{J}_s + i \omega \\mu \\vec{H} = - \\nabla \\times \\rho_s \\vec{J}_p - i \omega \\mu_s \\vec{H}_p \\
+        \\nabla \\times \\vec{H}_p - \\vec{J}_p = 0
+
+    Note: if different meshes are employed for the primary and secondary
+    problems, then we need to interpolate the fields from the primary mesh to
+    the secondary mesh. We do this by always interpolating the field and
+    computing a flux if need be in order to ensure that fluxes remain
+    numerically divergence free.
+
+    :param list rxList: Receiver list
+    :param float freq: frequency
+    :param numpy.array m: primary model
+    :param Problem prob: primary problem
+    :param Survey survey: primary survey
+    """
+
+
+    def __init__(self, rxList, freq, m, prob, survey):
         self.freq = float(freq)
-        self.mPrimary = mPrimary
+        self.m = m
+        self.prob = prob
+        self.survey = survey
+        self.fields = None
 
-        if primaryFields is None:
-            assert primaryProblem is not None and primarySurvey is not None, 'If no primary fields are provided, a primaryProblem and primarySurvey must be provided'
+        if self.survey.ispaired:
+            if self.survey.prob is not self.prob:
+                raise Exception('The survey object is already paired to a problem. Use survey.unpair()')
         else:
-            self._fields = primaryFields
-
-        if primaryProblem is not None and primarySurvey is not None:
-            self.prob = primaryProblem
-            self.survey = primarySurvey
-            if self.survey.ispaired:
-                if self.survey.prob is not self.prob:
-                    raise Exception('The survey object is already paired to a problem. Use survey.unpair()')
-            else:
-                self.prob.pair(self.survey)
-
+            self.prob.pair(self.survey)
 
         self.mesh = self.prob.mesh
-        self.integrate = False
+        self.prob.curModel = self.m
+
         BaseSrc.__init__(self, rxList)
 
-    # @property
-    def MeSigmaPrimary(self, prob):
-        if getattr(self, '_MeSigmaPrimary', None) is None:
-            sigmaprimary = self.prob.mapping.sigmaMap * self.mPrimary
+    def MeSigma(self, prob):
+        if getattr(self, '_MeSigma', None) is None:
+            sigmaprimary = self.prob.curModel.sigma
             if self.mesh != prob.mesh:
                 P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='CC')
                 sigmaprimary = P * sigmaprimary
-            self._MeSigmaPrimary = prob.mesh.getEdgeInnerProduct(sigmaprimary)
-        return self._MeSigmaPrimary
+            self._MeSigma = prob.mesh.getEdgeInnerProduct(sigmaprimary)
+        return self._MeSigma
 
-    # @property
-    def MfRhoIPrimary(self, prob):
-        if getattr(self, '_MfRhoIPrimary', None) is None:
-            rhoprimary = self.prob.mapping.rhoMap * self.mPrimary
+    def MfMui(self, prob):
+        if getattr(self, '_MfMui', None) is None:
+            muiprimary = self.prob.curModel.mui
+            if self.mesh != prob.mesh and not isinstance(muiprimary,float): # if different meshes and mu is a vector --> need to interpolate
+                P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='CC')
+                muiprimary = P * muiprimary
+            self._MfMui = prob.mesh.getFaceInnerProduct(muiprimary)
+        return self._MfMui
+
+    def MfRho(self, prob):
+        if getattr(self, '_MfRho', None) is None:
+            rhoprimary = self.prob.curModel.rho
             if self.mesh != prob.mesh:
                 P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='CC')
                 rhoprimary = P * rhoprimary
-            self._MfRhoIPrimary = prob.mesh.getFaceInnerProduct(rhoprimary, invMat=True)
-        return self._MfRhoIPrimary
+            self._MfRho = prob.mesh.getFaceInnerProduct(rhoprimary)
+        return self._MfRho
 
-    @property
-    def fields(self):
-        if getattr(self, '_fields', None) is None:
-            # check if I have fields if not, solve the primary to get fields object
-            self._fields = self.prob.fields(self.mPrimary)
-        return self._fields
+    def MeMu(self, prob):
+        if getattr(self, '_MeMu', None) is None:
+            muprimary = self.prob.curModel.mu
+            if self.mesh != prob.mesh and not isinstance(muiprimary,float): # if different meshes and mu is a vector --> need to interpolate
+                P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='CC')
+                muprimary = P * muprimary
+            self._MeMu = prob.mesh.getEdgeInnerProduct(muprimary)
+        return self._MeMu
 
+    # note if you switch from one formulation to another, but are using the same mesh, this will break
     def ePrimary(self,prob):
-        # check if a primary problem is defined
         if getattr(self, '_ePrimary', None) is None:
+            if self.fields is None:
+                self.fields = self.prob.fields(self.m)
+
             ePrimary = self.fields[:,'e']
 
             if self.mesh != prob.mesh:
-                P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='E')
+                if self.prob._formulation == 'HJ':
+                    P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType=prob._GLoc('e'), locTypeFrom='CCV')
+                else:
+                    P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType=prob._GLoc('e'))
                 ePrimary = Utils.mkvc(P * ePrimary)
             self._ePrimary = Utils.mkvc(ePrimary)
 
         return self._ePrimary
 
-    def jPrimary(self,prob):
-        # check if a primary problem is defined
-        # if getattr(self, '_jPrimary', None) is None:
-        # if self.prob is None or self.survey is None:
-        #     raise Exception('if Not specifying a jPrimary, a primarySurvey and primaryProblem must be provided.')
-        jPrimary = self.prob.fields(self.mPrimary)[:,'j']
-        if self.mesh != prob.mesh:
-            P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType='F')
-            jPrimary = P * jPrimary
-        return jPrimary
+    # note if you switch from one formulation to another, but are using the same mesh, this will break
+    def bPrimary(self, prob):
+        if getattr(self, '_bPrimary', None) is None:
+            if self.fields is None:
+                self.fields = self.prob.fields(self.m)
+
+            if self.mesh == prob.mesh:
+                bPrimary = self.fields[:,'b']
+            else:
+                bPrimary = prob.mesh.edgeCurl * self.ePrimary(prob)
+
+            self._bPrimary = Utils.mkvc(bPrimary)
+
+        return self._bPrimary
+
+    # note if you switch from one formulation to another, but are using the same mesh, this will break
+    def hPrimary(self, prob):
+        if getattr(self, '_hPrimary', None) is None:
+            if self.fields is None:
+                self.fields = self.prob.fields(self.m)
+
+            hPrimary = self.fields[:,'h']
+
+            if self.mesh != prob.mesh:
+                if self.prob._formulation == 'EB':
+                    P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType=prob._GLoc('h'), locTypeFrom='CCV')
+                else:
+                    P = self.mesh.getInterpolationMatMesh2Mesh(prob.mesh, locType=prob._GLoc('h'))
+                print P.shape, hPrimary.shape, prob._GLoc('h')
+                hPrimary = Utils.mkvc(P * hPrimary)
+            self._hPrimary = Utils.mkvc(hPrimary)
+
+        return self._hPrimary
+
+    # note if you switch from one formulation to another, but are using the same mesh, this will break
+    def jPrimary(self, prob):
+        if getattr(self, '_jPrimary', None) is None:
+            if self.fields is None:
+                self.fields = self.prob.fields(self.m)
+
+            if self.mesh == prob.mesh:
+                jPrimary = self.fields[:,'j']
+            else:
+                jPrimary = prob.mesh.edgeCurl * self.hPrimary(prob)
+
+            self._jPrimary = Utils.mkvc(jPrimary)
+
+        return self._jPrimary
 
     def s_e(self,prob):
-        # todo --> see if the primary and secondary are on the same mesh or not
         if prob._formulation == 'EB':
-            return Utils.mkvc(prob.MeSigma * self.ePrimary(prob) - self.MeSigmaPrimary(prob) * self.ePrimary(prob))
+            # - \\nabla \\times \\mu^{-1}_s \\vec{B}_p + \sigma_s \\vec{E}_p
+            s_e =  -prob.mesh.edgeCurl.T * ((prob.MfMui - self.MfMui(prob)) * self.bPrimary(prob)) +  (prob.MeSigma - self.MeSigma(prob)) * self.ePrimary(prob)
+            return Utils.mkvc(s_e)
+        else:
+            return Zero()
 
-    def s_eDeriv(self, prob, v, adjoint = False):
-        MeSigmaDeriv = prob.MeSigmaDeriv
-        if adjoint is not True:
-            return MeSigmaDeriv(self.ePrimary(prob)) * v
-        return MeSigmaDeriv(self.ePrimary(prob)) * v
+    def s_eDeriv(self, prob, v, adjoint=False):
+        if prob._formulation == 'EB':
+            if adjoint is True:
+                return prob.MeSigmaDeriv(self.ePrimary(prob)).T * v
+            return prob.MeSigmaDeriv(self.ePrimary(prob)) * v
+        else:
+            return Zero()
+
+    def s_m(self,prob):
+        if prob._formulation == 'HJ':
+            # - \\nabla \\times \\rho_s \\vec{J}_p - i \omega \\mu_s \\vec{H}_p
+            s_m = - prob.mesh.edgeCurl.T * (prob.MfRho - self.MfRho(prob)) * self.jPrimary(prob) - 1j * omega(self.freq) * ((prob.MeMu - self.MeMu(prob)) * self.hPrimary(prob))
+            return s_m
+        else:
+            return Zero()
+
+    def s_mDeriv(self, prob, v, adjoint=False):
+        if prob._formulation == 'HJ':
+            if adjoint is True:
+                return - prob.MfRhoDeriv(self.jPrimary(prob)).T * (prob.mesh.edgeCurl * v)
+            return - prob.mesh.edgeCurl.T * (prob.MfRhoDeriv(self.jPrimary(prob)) * v)
+        else:
+            return Zero()
 
 
 
diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py
index 163f6914..f9446df2 100644
--- a/SimPEG/EM/FDEM/SurveyFDEM.py
+++ b/SimPEG/EM/FDEM/SurveyFDEM.py
@@ -65,7 +65,7 @@ class Rx(SimPEG.Survey.BaseRx):
 
     def projGLoc(self, u):
         """Grid Location projection (e.g. Ex Fy ...)"""
-        return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
+        return u.prob._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
 
     def eval(self, src, mesh, u):
         """
@@ -77,8 +77,6 @@ class Rx(SimPEG.Survey.BaseRx):
         :rtype: numpy.ndarray
         :return: fields projected to recievers
         """
-        # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
-        # projGLoc += self.knownRxTypes[self.rxType][1]
 
         P = self.getP(mesh, self.projGLoc(u))
         u_part_complex = u[src, self.projField]
diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py
index b913c558..da05a5b4 100644
--- a/SimPEG/EM/Utils/testingUtils.py
+++ b/SimPEG/EM/Utils/testingUtils.py
@@ -51,17 +51,14 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
                 S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
                 S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
                 Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
-        elif SrcType is 'PrimSecSigma':
-            if fdemType is 'e' or fdemType is 'b':
+        elif SrcType is 'PrimSec':
                 primSrc = EM.FDEM.Src.MagDipole([], freq, np.r_[0.,0.,0.])
                 primarySurvey = EM.FDEM.Survey([primSrc])
                 primaryProblem = EM.FDEM.Problem_e(mesh,mapping=mapping)
                 mPrimary = np.ones(mapping.nP)*np.log(CONDUCTIVITY)
-                Src.append(EM.FDEM.Src.PrimSecSigma([Rx0], freq, mPrimary, primaryProblem=primaryProblem, primarySurvey=primarySurvey))
+                Src.append(EM.FDEM.Src.PrimSec([Rx0], freq, mPrimary, prob=primaryProblem, survey=primarySurvey))
 
         elif SrcType is 'PrimSecCyl':
-            if fdemType is 'e' or fdemType is 'b':
-
                 hx = [(cs,ncx + 2), (cs,npad + 2,1.3)]
                 hz = [(cs,npad + 2 ,-1.3), (cs,ncz+2), (cs,npad+2,1.3)]
                 primmesh = Mesh.CylMesh([hx,1,hz], '00C')
@@ -70,7 +67,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
                 primarySurvey = EM.FDEM.Survey([primSrc])
                 primaryProblem = EM.FDEM.Problem_e(primmesh)
                 mPrimary = np.ones(primmesh.nC)*CONDUCTIVITY
-                Src.append(EM.FDEM.Src.PrimSecSigma([Rx0], freq, mPrimary, primaryProblem=primaryProblem, primarySurvey=primarySurvey))
+                Src.append(EM.FDEM.Src.PrimSec([Rx0], freq, mPrimary, prob=primaryProblem, survey=primarySurvey))
 
     if verbose:
         print '  Fetching %s problem' % (fdemType)
diff --git a/SimPEG/Mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py
index bccbd4cc..d92fbd04 100644
--- a/SimPEG/Mesh/BaseMesh.py
+++ b/SimPEG/Mesh/BaseMesh.py
@@ -597,7 +597,7 @@ class BaseRectangularMesh(BaseMesh):
             return switchKernal(x)
 
 
-    def getInterpolationMatMesh2Mesh(self, mesh2, locType='CC'):
+    def getInterpolationMatMesh2Mesh(self, mesh2, locType='CC', locTypeFrom=None):
         """
         Interpolates variables from the current mesh to a new mesh (mesh2)
 
@@ -612,6 +612,9 @@ class BaseRectangularMesh(BaseMesh):
         #      "`getInterpolationMatMesh2Mesh` will be slow. If you want to interpolate a vector from one mesh to another, use `InterpolateVecMesh2Mesh`",
         #      RuntimeWarning)
 
+        if locTypeFrom is None:
+            locTypeFrom = locType # assume that we are interpolating to and from the same place
+
         # Error Checking
         if self._meshType == 'CYL':
             assert self.isSymmetric, "Currently, we do not support non-symmetric cyl meshes"
@@ -623,30 +626,30 @@ class BaseRectangularMesh(BaseMesh):
             return self.getInterpolationMatCartMesh(mesh2, locType)
 
         # Scalars
-        if locType in ['CC', 'N', 'Fx', 'Fy', 'Fz', 'Ex', 'Ey', 'Ez']:
-            grid = getattr(mesh2, 'grid%s'%locType)
+        if locType in ['CC', 'CCVx', 'CCVy', 'CCVz', 'N', 'Fx', 'Fy', 'Fz', 'Ex', 'Ey', 'Ez']:
+            grid = getattr(mesh2, 'grid%s'%locTypeFrom)
             return self.getInterpolationMat(grid, locType)
 
         # Vectors
         else:
             if self._meshType == 'CYL':
                 if locType == 'F':
-                    X = self.getInterpolationMatMesh2Mesh(mesh2, locType='Fx')
-                    Z = self.getInterpolationMatMesh2Mesh(mesh2, locType='Fz')
+                    X = self.getInterpolationMatMesh2Mesh(mesh2, locType='Fx', locTypeFrom=locTypeFrom+'x')
+                    Z = self.getInterpolationMatMesh2Mesh(mesh2, locType='Fz', locTypeFrom=locTypeFrom+'z')
                     return sp.block_diag([X, Z])
                 elif locType == 'E':
-                    return self.getInterpolationMatMesh2Mesh(mesh2, locType='Ey')
+                    return self.getInterpolationMatMesh2Mesh(mesh2, locType='Ey', locTypeFrom=locTypeFrom+'y')
 
             if self.dim == 1:
-                return self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType)
+                return self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType, locTypeFrom=locTypeFrom+'x')
             elif self.dim == 2:
-                X = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType)
-                Y = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sy'%locType)
+                X = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType, locTypeFrom=locTypeFrom+'x')
+                Y = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sy'%locType, locTypeFrom=locTypeFrom+'y')
                 return sp.block_diag([X, Y])
             elif self.dim == 3:
-                X = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType)
-                Y = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sy'%locType)
-                Z = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sz'%locType)
+                X = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType, locTypeFrom=locTypeFrom+'x')
+                Y = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sy'%locType, locTypeFrom=locTypeFrom+'y')
+                Z = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sz'%locType, locTypeFrom=locTypeFrom+'z')
                 return sp.block_diag([X, Y, Z])