don't test 3D mesh2mesh (too big at the moment)

# Conflicts:
#	SimPEG/EM/Utils/testingUtils.py

SimPEG/DCIP/BaseDC.py

@@ -200,11 +200,11 @@ class ProblemDC_CC(Problem.BaseProblem):
 
         return F
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         """
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: Jv
 
@@ -225,11 +225,10 @@ class ProblemDC_CC(Problem.BaseProblem):
         # Set current model; clear dependent property $\mathbf{A(m)}$
         self.curModel = m
         sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
-        if u is None:
+        if f is None:
             # Run forward simulation if $u$ not provided
-            u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
-        else:
-            u = u[self.survey.srcList, 'phi_sol']
+            f = self.fields(self.curModel)
+        u = f[self.survey.srcList, 'phi_sol']
 
         D = self.mesh.faceDiv
         G = self.mesh.cellGrad
@@ -251,19 +250,18 @@ class ProblemDC_CC(Problem.BaseProblem):
         if self.Ainv is None:
             self.Ainv = self.Solver(dA_du, **self.solverOpts)
 
-        P        = self.survey.getP(self.mesh)
+        P     = self.survey.getP(self.mesh)
         Jv    = - P * mkvc( self.Ainv * dCdm_x_v )
         return Jv
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
 
         self.curModel = m
         sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
-        if u is None:
-            # Run forward simulation if $u$ not provided
-            u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
-        else:
-            u = u[self.survey.srcList, 'phi_sol']
+        if f is None:
+            # Run forward simulation if $f$ not provided
+            f = self.fields(self.curModel)
+        u = f[self.survey.srcList, 'phi_sol']
 
         shp = u.shape
         P = self.survey.getP(self.mesh)

SimPEG/DCIP/BaseIP.py

@@ -14,12 +14,12 @@ class SurveyIP(SurveyDC):
         Survey.BaseSurvey.__init__(self, **kwargs)
         self._Ps = {}
 
-    def dpred(self, m, u=None):
+    def dpred(self, m, f=None):
         """
             Predicted data.
 
             .. math::
-                d_\\text{pred} = Pu(m)
+                d_\\text{pred} = Pf(m)
         """
 
         return self.prob.forward(m)
@@ -143,10 +143,10 @@ class ProblemIP(Problem.BaseProblem):
         J_x_v    = - P * mkvc( self.Ainv * dCdm_x_v )
         return -J_x_v
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         return self.forward(v)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
 
         self.curModel = m
         # sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$

SimPEG/DataMisfit.py

@@ -22,11 +22,11 @@ class BaseDataMisfit(object):
         Utils.setKwargs(self,**kwargs)
 
     @Utils.timeIt
-    def eval(self, m, u=None):
-        """eval(m, u=None)
+    def eval(self, m, f=None):
+        """eval(m, f=None)
 
             :param numpy.array m: geophysical model
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: float
             :return: data misfit
 
@@ -34,11 +34,11 @@ class BaseDataMisfit(object):
         raise NotImplementedError('This method should be overwritten.')
 
     @Utils.timeIt
-    def evalDeriv(self, m, u=None):
-        """evalDeriv(m, u=None)
+    def evalDeriv(self, m, f=None):
+        """evalDeriv(m, f=None)
 
             :param numpy.array m: geophysical model
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: data misfit derivative
 
@@ -47,12 +47,12 @@ class BaseDataMisfit(object):
 
 
     @Utils.timeIt
-    def eval2Deriv(self, m, v, u=None):
-        """eval2Deriv(m, v, u=None)
+    def eval2Deriv(self, m, v, f=None):
+        """eval2Deriv(m, v, f=None)
 
             :param numpy.array m: geophysical model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: data misfit derivative
 
@@ -89,7 +89,7 @@ class l2_DataMisfit(BaseDataMisfit):
         """
 
         if getattr(self, '_Wd', None) is None:
-            
+
             survey = self.survey
 
             if getattr(survey,'std', None) is None:
@@ -108,24 +108,20 @@ class l2_DataMisfit(BaseDataMisfit):
         self._Wd = value
 
     @Utils.timeIt
-    def eval(self, m, u=None):
-        "eval(m, u=None)"
-        prob   = self.prob
-        survey = self.survey
-        R = self.Wd * survey.residual(m, u=u)
+    def eval(self, m, f=None):
+        "eval(m, f=None)"
+        if f is None: f = self.prob.fields(m)
+        R = self.Wd * self.survey.residual(m, f)
         return 0.5*np.vdot(R, R)
 
     @Utils.timeIt
-    def evalDeriv(self, m, u=None):
-        "evalDeriv(m, u=None)"
-        prob   = self.prob
-        survey = self.survey
-        if u is None: u = prob.fields(m)
-        return prob.Jtvec(m, self.Wd * (self.Wd * survey.residual(m, u=u)), u=u)
+    def evalDeriv(self, m, f=None):
+        "evalDeriv(m, f=None)"
+        if f is None: f = self.prob.fields(m)
+        return self.prob.Jtvec(m, self.Wd * (self.Wd * self.survey.residual(m, f=f)), f=f)
 
     @Utils.timeIt
-    def eval2Deriv(self, m, v, u=None):
-        "eval2Deriv(m, v, u=None)"
-        prob   = self.prob
-        if u is None: u = prob.fields(m)
-        return prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, u=u)), u=u)
+    def eval2Deriv(self, m, v, f=None):
+        "eval2Deriv(m, v, f=None)"
+        if f is None: f = self.prob.fields(m)
+        return self.prob.Jtvec_approx(m, self.Wd * (self.Wd * self.prob.Jvec_approx(m, v, f=f)), f=f)

SimPEG/Directives.py

@@ -123,10 +123,10 @@ class BetaEstimate_ByEig(InversionDirective):
         if self.debug: print 'Calculating the beta0 parameter.'
 
         m = self.invProb.curModel
-        u = self.invProb.getFields(m, store=True, deleteWarmstart=False)
+        f = self.invProb.getFields(m, store=True, deleteWarmstart=False)
 
         x0 = np.random.rand(*m.shape)
-        t = x0.dot(self.dmisfit.eval2Deriv(m,x0,u=u))
+        t = x0.dot(self.dmisfit.eval2Deriv(m,x0,f=f))
         b = x0.dot(self.reg.eval2Deriv(m, v=x0))
         self.beta0 = self.beta0_ratio*(t/b)
 

SimPEG/EM/Base.py

@@ -2,14 +2,14 @@ from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solve
 from scipy.constants import mu_0
 
 class EMPropMap(Maps.PropMap):
-    """ 
+    """
         Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m)
     """
 
     sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap))
-    mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap))
+    mu = Maps.Property("Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap))
 
-    rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap)) 
+    rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap))
     mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap))
 
 
@@ -21,7 +21,7 @@ class BaseEMProblem(Problem.BaseProblem):
 
     surveyPair = Survey.BaseSurvey
     dataPair = Survey.Data
-    
+
     PropMap = EMPropMap
 
     Solver = SimpegSolver
@@ -51,7 +51,7 @@ class BaseEMProblem(Problem.BaseProblem):
         if self.mapping.muMap is not None or self.mapping.muiMap is not None:
             toDelete += ['_MeMu', '_MeMuI','_MfMui','_MfMuiI']
         return toDelete
-    
+
     @property
     def Me(self):
         """
@@ -71,7 +71,7 @@ class BaseEMProblem(Problem.BaseProblem):
         return self._Mf
 
 
-    # ----- Magnetic Permeability ----- # 
+    # ----- Magnetic Permeability ----- #
     @property
     def MfMui(self):
         """
@@ -109,7 +109,7 @@ class BaseEMProblem(Problem.BaseProblem):
         return self._MeMuI
 
 
-    # ----- Electrical Conductivity ----- # 
+    # ----- Electrical Conductivity ----- #
     #TODO: hardcoded to sigma as the model
     @property
     def MeSigma(self):
@@ -120,18 +120,18 @@ class BaseEMProblem(Problem.BaseProblem):
             self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma)
         return self._MeSigma
 
-    # TODO: This should take a vector 
+    # TODO: This should take a vector
     def MeSigmaDeriv(self, u):
         """
             Derivative of MeSigma with respect to the model
-        """ 
+        """
         return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv
-    
+
 
     @property
     def MeSigmaI(self):
         """
-            Inverse of the edge inner product matrix for \\(\\sigma\\). 
+            Inverse of the edge inner product matrix for \\(\\sigma\\).
         """
         if getattr(self, '_MeSigmaI', None) is None:
             self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True)
@@ -140,8 +140,8 @@ class BaseEMProblem(Problem.BaseProblem):
     # TODO: This should take a vector
     def MeSigmaIDeriv(self, u):
         """
-            Derivative of :code:`MeSigma` with respect to the model 
-        """ 
+            Derivative of :code:`MeSigma` with respect to the model
+        """
         # TODO: only works for diagonal tensors. getEdgeInnerProductDeriv, invMat=True should be implemented in SimPEG
 
         dMeSigmaI_dI = -self.MeSigmaI**2
@@ -163,7 +163,7 @@ class BaseEMProblem(Problem.BaseProblem):
     # TODO: This should take a vector
     def MfRhoDeriv(self,u):
         """
-            Derivative of :code:`MfRho` with respect to the model. 
+            Derivative of :code:`MfRho` with respect to the model.
         """
         return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv)
         # self.curModel.rhoDeriv
@@ -181,6 +181,29 @@ class BaseEMProblem(Problem.BaseProblem):
     # TODO: This should take a vector
     def MfRhoIDeriv(self,u):
         """
-            Derivative of :code:`MfRhoI` with respect to the model. 
+            Derivative of :code:`MfRhoI` with respect to the model.
         """
         return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
+
+class BaseEMSurvey(Survey.BaseSurvey):
+
+    def __init__(self, srcList, **kwargs):
+        # Sort these by frequency
+        self.srcList = srcList
+        Survey.BaseSurvey.__init__(self, **kwargs)
+
+    def eval(self, u):
+        """
+        Project fields to receiver locations
+        :param Fields u: fields object
+        :rtype: numpy.ndarray
+        :return: data
+        """
+        data = Survey.Data(self)
+        for src in self.srcList:
+            for rx in src.rxList:
+                data[src, rx] = rx.eval(src, self.mesh, u)
+        return data
+
+    def evalDeriv(self, u):
+        raise Exception('Use Receivers to project fields deriv.')

SimPEG/EM/FDEM/FDEM.py

@@ -8,104 +8,104 @@ 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
     fieldsPair = Fields
 
-    def fields(self, m=None):
+    def fields(self, m):
         """
         Solve the forward problem for the fields.
-        
+
         :param numpy.array m: inversion model (nP,)
         :rtype numpy.array:
-        :return F: forward solution
+        :return f: forward solution
         """
 
         self.curModel = m
-        F = self.fieldsPair(self.mesh, self.survey)
+        f = self.fieldsPair(self.mesh, self.survey)
 
         for freq in self.survey.freqs:
             A = self.getA(freq)
             rhs = self.getRHS(freq)
             Ainv = self.Solver(A, **self.solverOpts)
-            sol = Ainv * rhs
+            u = Ainv * rhs
             Srcs = self.survey.getSrcByFreq(freq)
-            F[Srcs, self._solutionType] = sol
+            f[Srcs, self._solutionType] = u
             Ainv.clean()
-        return F
+        return f
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         """
         Sensitivity times a vector.
 
         :param numpy.array m: inversion model (nP,)
         :param numpy.array v: vector which we take sensitivity product with (nP,)
-        :param SimPEG.EM.FDEM.Fields u: fields object 
+        :param SimPEG.EM.FDEM.Fields u: fields object
         :rtype numpy.array:
-        :return: Jv (ndata,) 
+        :return: Jv (ndata,)
         """
 
-        if u is None:
-           u = self.fields(m)
+        if f is None:
+           f = self.fields(m)
 
         self.curModel = m
 
         Jv = self.dataPair(self.survey)
 
         for freq in self.survey.freqs:
-            A = self.getA(freq) 
-            Ainv = self.Solver(A, **self.solverOpts)
+            A = self.getA(freq)
+            Ainv = self.Solver(A, **self.solverOpts) # create the concept of Ainv (actually a solve)
 
             for src in self.survey.getSrcByFreq(freq):
-                u_src = u[src, self._solutionType]
+                u_src = f[src, self._solutionType]
                 dA_dm_v = self.getADeriv(freq, u_src, v)
-                dRHS_dm_v = self.getRHSDeriv(freq, src, v) 
+                dRHS_dm_v = self.getRHSDeriv(freq, src, v)
                 du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
-                
+
                 for rx in src.rxList:
-                    df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None)
+                    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, u, df_dm_v)
+                    Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
             Ainv.clean()
         return Utils.mkvc(Jv)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
         """
         Sensitivity transpose times a vector
 
         :param numpy.array m: inversion model (nP,)
         :param numpy.array v: vector which we take adjoint product with (nP,)
-        :param SimPEG.EM.FDEM.Fields u: fields object 
+        :param SimPEG.EM.FDEM.Fields u: fields object
         :rtype numpy.array:
-        :return: Jv (ndata,) 
+        :return: Jv (ndata,)
         """
 
-        if u is None:
-            u = self.fields(m)
+        if f is None:
+            f = self.fields(m)
 
         self.curModel = m
 
@@ -120,12 +120,12 @@ class BaseFDEMProblem(BaseEMProblem):
             ATinv = self.Solver(AT, **self.solverOpts)
 
             for src in self.survey.getSrcByFreq(freq):
-                u_src = u[src, self._solutionType]
+                u_src = f[src, self._solutionType]
 
                 for rx in src.rxList:
-                    PTv = rx.evalDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
+                    PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
 
-                    df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None)
+                    df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
                     df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
 
                     ATinvdf_duT = ATinv * df_duT
@@ -144,7 +144,7 @@ class BaseFDEMProblem(BaseEMProblem):
                         Jtv += - np.array(df_dmT, dtype=complex).real
                     else:
                         raise Exception('Must be real or imag')
-            
+
             ATinv.clean()
 
         return Utils.mkvc(Jtv)
@@ -154,23 +154,23 @@ class BaseFDEMProblem(BaseEMProblem):
         Evaluates the sources for a given frequency and puts them in matrix form
 
         :param float freq: Frequency
-        :rtype: (numpy.ndarray, numpy.ndarray)  
-        :return: S_m, S_e (nE or nF, nSrc)
+        :rtype: (numpy.ndarray, numpy.ndarray)
+        :return: s_m, s_e (nE or nF, nSrc)
         """
         Srcs = self.survey.getSrcByFreq(freq)
         if self._formulation is 'EB':
-            S_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
-            S_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
+            s_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
+            s_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
         elif self._formulation is 'HJ':
-            S_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
-            S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
+            s_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
+            s_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
 
         for i, src in enumerate(Srcs):
             smi, sei = src.eval(self)
-            S_m[:,i] = S_m[:,i] + smi
-            S_e[:,i] = S_e[:,i] + sei
+            s_m[:,i] = s_m[:,i] + smi
+            s_e[:,i] = s_e[:,i] + sei
 
-        return S_m, S_e
+        return s_m, s_e
 
 
 ##########################################################################################
@@ -204,10 +204,21 @@ 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
-        
+
         .. math ::
             \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}}
 
@@ -230,12 +241,12 @@ class Problem_e(BaseFDEMProblem):
         .. math ::
             \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = i \omega \\frac{d \mathbf{M^e_{\sigma}}\mathbf{v} }{d\mathbf{m}}
 
-        :param float freq: frequency 
-        :param numpy.ndarray u: solution vector (nE,) 
+        :param float freq: frequency
+        :param numpy.ndarray u: solution vector (nE,)
         :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
-        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) 
+        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
         """
 
         dsig_dm = self.curModel.sigmaDeriv
@@ -248,25 +259,25 @@ class Problem_e(BaseFDEMProblem):
 
     def getRHS(self, freq):
         """
-        Right hand side for the system 
+        Right hand side for the system
 
         .. math ::
             \mathbf{RHS} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e}
 
         :param float freq: Frequency
-        :rtype: numpy.ndarray 
+        :rtype: numpy.ndarray
         :return: RHS (nE, nSrc)
         """
 
-        S_m, S_e = self.getSourceTerm(freq)
+        s_m, s_e = self.getSourceTerm(freq)
         C = self.mesh.edgeCurl
         MfMui = self.MfMui
 
-        return C.T * (MfMui * S_m) -1j * omega(freq) * S_e
+        return C.T * (MfMui * s_m) -1j * omega(freq) * s_e
 
     def getRHSDeriv(self, freq, src, v, adjoint=False):
         """
-        Derivative of the right hand side with respect to the model 
+        Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
         :param SimPEG.EM.FDEM.Src src: FDEM source
@@ -278,14 +289,13 @@ class Problem_e(BaseFDEMProblem):
 
         C = self.mesh.edgeCurl
         MfMui = self.MfMui
-        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
+        s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
 
         if adjoint:
             dRHS = MfMui * (C * v)
-            return S_mDeriv(dRHS) - 1j * omega(freq) * S_eDeriv(v)
-
+            return s_mDeriv(dRHS) - 1j * omega(freq) * s_eDeriv(v)
         else:
-            return C.T * (MfMui * S_mDeriv(v)) -1j * omega(freq) * S_eDeriv(v)
+            return C.T * (MfMui * s_mDeriv(v)) -1j * omega(freq) * s_eDeriv(v)
 
 
 class Problem_b(BaseFDEMProblem):
@@ -315,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
@@ -346,12 +366,12 @@ class Problem_b(BaseFDEMProblem):
         .. math ::
             \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \\frac{\mathbf{M^e_{\sigma}} \mathbf{v}}{d\mathbf{m}}
 
-        :param float freq: frequency 
-        :param numpy.ndarray u: solution vector (nF,) 
+        :param float freq: frequency
+        :param numpy.ndarray u: solution vector (nF,)
         :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
-        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)  
+        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
         """
 
         MfMui = self.MfMui
@@ -373,21 +393,21 @@ class Problem_b(BaseFDEMProblem):
 
     def getRHS(self, freq):
         """
-        Right hand side for the system 
+        Right hand side for the system
 
         .. math ::
             \mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e}
 
         :param float freq: Frequency
-        :rtype: numpy.ndarray 
+        :rtype: numpy.ndarray
         :return: RHS (nE, nSrc)
         """
 
-        S_m, S_e = self.getSourceTerm(freq)
+        s_m, s_e = self.getSourceTerm(freq)
         C = self.mesh.edgeCurl
         MeSigmaI = self.MeSigmaI
 
-        RHS = S_m + C * ( MeSigmaI * S_e )
+        RHS = s_m + C * ( MeSigmaI * s_e )
 
         if self._makeASymmetric is True:
             MfMui = self.MfMui
@@ -408,21 +428,21 @@ class Problem_b(BaseFDEMProblem):
         """
 
         C = self.mesh.edgeCurl
-        S_m, S_e = src.eval(self)
+        s_m, s_e = src.eval(self)
         MfMui = self.MfMui
 
         if self._makeASymmetric and adjoint:
             v = self.MfMui * v
 
-        MeSigmaIDeriv = self.MeSigmaIDeriv(S_e)
-        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
+        MeSigmaIDeriv = self.MeSigmaIDeriv(s_e)
+        s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
 
         if not adjoint:
             RHSderiv = C * (MeSigmaIDeriv * v)
-            SrcDeriv = S_mDeriv(v) + C * (self.MeSigmaI * S_eDeriv(v))
+            SrcDeriv = s_mDeriv(v) + C * (self.MeSigmaI * s_eDeriv(v))
         elif adjoint:
             RHSderiv = MeSigmaIDeriv.T * (C.T * v)
-            SrcDeriv = S_mDeriv(v) + self.MeSigmaI.T * (C.T * S_eDeriv(v))
+            SrcDeriv = s_mDeriv(v) + self.MeSigmaI.T * (C.T * s_eDeriv(v))
 
         if self._makeASymmetric is True and not adjoint:
             return MfMui.T * (SrcDeriv + RHSderiv)
@@ -463,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
@@ -497,12 +527,12 @@ class Problem_j(BaseFDEMProblem):
 
             \\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \mathbf{m}} = \mathbf{C}  \mathbf{M^e_{mu^{-1}}} \mathbf{C^{\\top}} \\frac{d \mathbf{M^f_{\sigma^{-1}}}\mathbf{v} }{d \mathbf{m}}
 
-        :param float freq: frequency 
-        :param numpy.ndarray u: solution vector (nF,) 
+        :param float freq: frequency
+        :param numpy.ndarray u: solution vector (nF,)
         :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
-        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)  
+        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
         """
 
         MeMuI = self.MeMuI
@@ -522,7 +552,7 @@ class Problem_j(BaseFDEMProblem):
 
     def getRHS(self, freq):
         """
-        Right hand side for the system 
+        Right hand side for the system
 
         .. math ::
 
@@ -533,11 +563,11 @@ class Problem_j(BaseFDEMProblem):
         :return: RHS
         """
 
-        S_m, S_e = self.getSourceTerm(freq)
+        s_m, s_e = self.getSourceTerm(freq)
         C = self.mesh.edgeCurl
         MeMuI = self.MeMuI
 
-        RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e
+        RHS = C * (MeMuI * s_m) - 1j * omega(freq) * s_e
         if self._makeASymmetric is True:
             MfRho = self.MfRho
             return MfRho.T*RHS
@@ -546,7 +576,7 @@ class Problem_j(BaseFDEMProblem):
 
     def getRHSDeriv(self, freq, src, v, adjoint=False):
         """
-        Derivative of the right hand side with respect to the model 
+        Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
         :param SimPEG.EM.FDEM.Src src: FDEM source
@@ -558,16 +588,16 @@ class Problem_j(BaseFDEMProblem):
 
         C = self.mesh.edgeCurl
         MeMuI = self.MeMuI
-        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
+        s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
 
         if adjoint:
             if self._makeASymmetric:
                 MfRho = self.MfRho
                 v = MfRho*v
-            return S_mDeriv(MeMuI.T * (C.T * v)) - 1j * omega(freq) * S_eDeriv(v)
+            return s_mDeriv(MeMuI.T * (C.T * v)) - 1j * omega(freq) * s_eDeriv(v)
 
         else:
-            RHSDeriv = C * (MeMuI * S_mDeriv(v)) - 1j * omega(freq) * S_eDeriv(v)
+            RHSDeriv = C * (MeMuI * s_mDeriv(v)) - 1j * omega(freq) * s_eDeriv(v)
 
             if self._makeASymmetric:
                 MfRho = self.MfRho
@@ -601,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
@@ -626,12 +667,12 @@ class Problem_h(BaseFDEMProblem):
         .. math::
             \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C}^{\\top}\\frac{d \mathbf{M^f_{\\rho}}\mathbf{v} }{d\mathbf{m}}
 
-        :param float freq: frequency 
-        :param numpy.ndarray u: solution vector (nE,) 
+        :param float freq: frequency
+        :param numpy.ndarray u: solution vector (nE,)
         :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
-        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)  
+        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
         """
 
         MeMu = self.MeMu
@@ -644,26 +685,26 @@ class Problem_h(BaseFDEMProblem):
 
     def getRHS(self, freq):
         """
-        Right hand side for the system 
+        Right hand side for the system
 
         .. math ::
 
             \mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
 
         :param float freq: Frequency
-        :rtype: numpy.ndarray 
+        :rtype: numpy.ndarray
         :return: RHS (nE, nSrc)
         """
 
-        S_m, S_e = self.getSourceTerm(freq)
+        s_m, s_e = self.getSourceTerm(freq)
         C = self.mesh.edgeCurl
         MfRho  = self.MfRho
 
-        return S_m + C.T * ( MfRho * S_e )
+        return s_m + C.T * ( MfRho * s_e )
 
     def getRHSDeriv(self, freq, src, v, adjoint=False):
         """
-        Derivative of the right hand side with respect to the model 
+        Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
         :param SimPEG.EM.FDEM.Src src: FDEM source
@@ -673,17 +714,17 @@ class Problem_h(BaseFDEMProblem):
         :return: product of rhs deriv with a vector
         """
 
-        _, S_e = src.eval(self)
+        _, s_e = src.eval(self)
         C = self.mesh.edgeCurl
         MfRho  = self.MfRho
 
-        MfRhoDeriv = self.MfRhoDeriv(S_e)
+        MfRhoDeriv = self.MfRhoDeriv(s_e)
         if not adjoint:
             RHSDeriv = C.T * (MfRhoDeriv * v)
         elif adjoint:
             RHSDeriv = MfRhoDeriv.T * (C * v)
 
-        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
+        s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
 
-        return RHSDeriv + S_mDeriv(v) + C.T * (MfRho * S_eDeriv(v))
+        return RHSDeriv + s_mDeriv(v) + C.T * (MfRho * s_eDeriv(v))
 

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -8,7 +8,7 @@ from SimPEG.Utils import Zero, Identity, sdiag
 
 class Fields(SimPEG.Problem.Fields):
     """
-    
+
     Fancy Field Storage for a FDEM survey. Only one field type is stored for
     each problem, the rest are computed. The fields obejct acts like an array and is indexed by
 
@@ -34,56 +34,56 @@ class Fields(SimPEG.Problem.Fields):
 
     def _e(self, solution, srcList):
         """
-        Total electric field is sum of primary and secondary 
-        
+        Total electric field is sum of primary and secondary
+
         :param numpy.ndarray solution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
         :return: total electric field
         """
-        if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None: 
+        if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
             raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0])
 
         return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
 
     def _b(self, solution, srcList):
         """
-        Total magnetic flux density is sum of primary and secondary 
-        
+        Total magnetic flux density is sum of primary and secondary
+
         :param numpy.ndarray solution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
-        :return: total magnetic flux density 
+        :return: total magnetic flux density
         """
-        if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None: 
+        if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
             raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
 
         return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
 
     def _h(self, solution, srcList):
         """
-        Total magnetic field is sum of primary and secondary 
-        
+        Total magnetic field is sum of primary and secondary
+
         :param numpy.ndarray solution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
         :return: total magnetic field
         """
-        if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None: 
+        if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
             raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0])
 
         return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
 
     def _j(self, solution, srcList):
         """
-        Total current density is sum of primary and secondary 
-        
+        Total current density is sum of primary and secondary
+
         :param numpy.ndarray solution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
-        :return: total current density 
+        :return: total current density
         """
-        if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None: 
+        if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
             raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0])
 
         return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
@@ -99,7 +99,7 @@ class Fields(SimPEG.Problem.Fields):
         :rtype: numpy.ndarray
         :return: derivative times a vector (or tuple for adjoint)
         """
-        if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None: 
+        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:
@@ -117,12 +117,12 @@ class Fields(SimPEG.Problem.Fields):
         :rtype: numpy.ndarray
         :return: derivative times a vector (or tuple for adjoint)
         """
-        if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None: 
+        if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
             raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
         if adjoint:
             return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
-        return np.array(self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint), dtype = complex) 
+        return np.array(self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint), dtype = complex)
 
     def _hDeriv(self, src, du_dm_v, v, adjoint = False):
         """
@@ -135,10 +135,10 @@ class Fields(SimPEG.Problem.Fields):
         :rtype: numpy.ndarray
         :return: derivative times a vector (or tuple for adjoint)
         """
-        if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None: 
+        if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
             raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
-        if adjoint: 
+        if adjoint:
             return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
         return np.array(self._hDeriv_u(src, du_dm_v, adjoint) + self._hDeriv_m(src, v, adjoint), dtype = complex)
 
@@ -153,7 +153,7 @@ class Fields(SimPEG.Problem.Fields):
         :rtype: numpy.ndarray
         :return: derivative times a vector (or tuple for adjoint)
         """
-        if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None: 
+        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:
@@ -162,10 +162,10 @@ class Fields(SimPEG.Problem.Fields):
 
 class Fields_e(Fields):
     """
-    Fields object for Problem_e. 
+    Fields object for Problem_e.
 
     :param Mesh mesh: mesh
-    :param Survey survey: survey 
+    :param Survey survey: survey
     """
 
     knownFields = {'eSolution':'E'}
@@ -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):
         """
@@ -233,9 +223,9 @@ class Fields_e(Fields):
 
     def _eDeriv_u(self, src, v, adjoint = False):
         """
-        Partial derivative of the total electric field with respect to the thing we 
+        Partial derivative of the total electric field with respect to the thing we
         solved for.
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -247,8 +237,8 @@ class Fields_e(Fields):
 
     def _eDeriv_m(self, src, v, adjoint = False):
         """
-        Partial derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. 
-        
+        Partial derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -289,14 +279,14 @@ class Fields_e(Fields):
         b = (C * eSolution)
         for i, src in enumerate(srcList):
             b[:,i] *= - 1./(1j*omega(src.freq))
-            S_m, _ = src.eval(self.prob)
-            b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m
+            s_m, _ = src.eval(self.prob)
+            b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * s_m
         return b
 
     def _bDeriv_u(self, src, du_dm_v, adjoint = False):
         """
         Derivative of the magnetic flux density with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -312,8 +302,8 @@ class Fields_e(Fields):
 
     def _bDeriv_m(self, src, v, adjoint = False):
         """
-        Derivative of the magnetic flux density with respect to the inversion model. 
-        
+        Derivative of the magnetic flux density with respect to the inversion model.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -321,8 +311,8 @@ class Fields_e(Fields):
         :return: product of the magnetic flux density derivative with respect to the inversion model with a vector
         """
 
-        S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint)
-        return 1./(1j * omega(src.freq)) * S_mDeriv
+        s_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint)
+        return 1./(1j * omega(src.freq)) * s_mDeriv
 
     def _j(self,  eSolution, srcList):
         """
@@ -341,7 +331,7 @@ class Fields_e(Fields):
     def _jDeriv_u(self, src, du_dm_v, adjoint = False):
         """
         Derivative of the current density with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -351,15 +341,15 @@ class Fields_e(Fields):
         n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
-        if adjoint: 
+        if adjoint:
             return  self._eDeriv_u(src, self._MeSigma.T * (self._aveE2CCV.T * (VI.T * du_dm_v) ), adjoint = adjoint)
         return VI * (self._aveE2CCV * (self._MeSigma * (self._eDeriv_u(src, du_dm_v, adjoint=adjoint) ) ) )
-        
+
 
     def _jDeriv_m(self, src, v, adjoint = False):
         """
-        Derivative of the current density with respect to the inversion model. 
-        
+        Derivative of the current density with respect to the inversion model.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -373,7 +363,7 @@ class Fields_e(Fields):
         if adjoint:
             return self._MeSigmaDeriv(e).T * (self._aveE2CCV.T * (VI.T * v)) + self._eDeriv_m(src, self._aveE2CCV.T * (VI.T * v), adjoint=adjoint)
         return VI * (self._aveE2CCV * ( self._eDeriv_m(src, v, adjoint=adjoint) + self._MeSigmaDeriv(e) * v))
-        
+
 
 
     def _h(self, eSolution, srcList):
@@ -393,7 +383,7 @@ class Fields_e(Fields):
     def _hDeriv_u(self, src, du_dm_v, adjoint = False):
         """
         Derivative of the magnetic field with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -409,8 +399,8 @@ class Fields_e(Fields):
 
     def _hDeriv_m(self, src, v, adjoint = False):
         """
-        Derivative of the magnetic field with respect to the inversion model. 
-        
+        Derivative of the magnetic field with respect to the inversion model.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -428,10 +418,10 @@ class Fields_e(Fields):
 
 class Fields_b(Fields):
     """
-    Fields object for Problem_b. 
+    Fields object for Problem_b.
 
     :param Mesh mesh: mesh
-    :param Survey survey: survey 
+    :param Survey survey: survey
     """
 
     knownFields = {'bSolution':'F'}
@@ -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
@@ -506,9 +485,9 @@ class Fields_b(Fields):
 
     def _bDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the total magnetic flux density with respect to the thing we 
+        Partial derivative of the total magnetic flux density with respect to the thing we
         solved for.
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -520,8 +499,8 @@ class Fields_b(Fields):
 
     def _bDeriv_m(self, src, v, adjoint=False):
         """
-        Partial derivative of the total magnetic flux density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. 
-        
+        Partial derivative of the total magnetic flux density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -560,15 +539,15 @@ class Fields_b(Fields):
 
         e = ( self._edgeCurl.T * ( self._MfMui * bSolution))
         for i,src in enumerate(srcList):
-            _,S_e = src.eval(self.prob)
-            e[:,i] = e[:,i] + - S_e
+            _,s_e = src.eval(self.prob)
+            e[:,i] = e[:,i] + - s_e
 
         return self._MeSigmaI * e
 
     def _eDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the electric field with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -583,8 +562,8 @@ class Fields_b(Fields):
 
     def _eDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the electric field with respect to the inversion model 
-        
+        Derivative of the electric field with respect to the inversion model
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -593,15 +572,15 @@ class Fields_b(Fields):
         """
 
         bSolution = Utils.mkvc(self[src, 'bSolution'])
-        _,S_e = src.eval(self.prob)
+        _,s_e = src.eval(self.prob)
 
-        w = -S_e + self._edgeCurl.T * (self._MfMui * bSolution)
-        _, S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
+        w = -s_e + self._edgeCurl.T * (self._MfMui * bSolution)
+        _, s_eDeriv = src.evalDeriv(self.prob, v, adjoint)
 
 
         if adjoint:
-            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * S_eDeriv
-        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * S_eDeriv
+            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv
+        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv
 
     def _j(self, bSolution, srcList):
         """
@@ -617,13 +596,13 @@ class Fields_b(Fields):
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
         return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) )
-        
+
 
     def _jDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the current density with respect to the thing we 
+        Partial derivative of the current density with respect to the thing we
         solved for.
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -639,8 +618,8 @@ class Fields_b(Fields):
 
     def _jDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the current density with respect to the inversion model 
-        
+        Derivative of the current density with respect to the inversion model
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -664,9 +643,9 @@ class Fields_b(Fields):
 
     def _hDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the magnetic field with respect to the thing we 
+        Partial derivative of the magnetic field with respect to the thing we
         solved for.
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -682,8 +661,8 @@ class Fields_b(Fields):
 
     def _hDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the magnetic field with respect to the inversion model 
-        
+        Derivative of the magnetic field with respect to the inversion model
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -695,10 +674,10 @@ class Fields_b(Fields):
 
 class Fields_j(Fields):
     """
-    Fields object for Problem_j. 
+    Fields object for Problem_j.
 
     :param Mesh mesh: mesh
-    :param Survey survey: survey 
+    :param Survey survey: survey
     """
 
     knownFields = {'jSolution':'F'}
@@ -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
@@ -769,12 +738,12 @@ class Fields_j(Fields):
 
     def _j(self, jSolution, srcList):
         """
-        Total current density is sum of primary and secondary 
-        
+        Total current density is sum of primary and secondary
+
         :param numpy.ndarray jSolution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
-        :return: total current density 
+        :return: total current density
         """
 
         return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList)
@@ -782,9 +751,9 @@ class Fields_j(Fields):
 
     def _jDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the total current density with respect to the thing we 
+        Partial derivative of the total current density with respect to the thing we
         solved for.
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -797,8 +766,8 @@ class Fields_j(Fields):
 
     def _jDeriv_m(self, src, v, adjoint=False):
         """
-        Partial derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. 
-        
+        Partial derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -837,15 +806,15 @@ class Fields_j(Fields):
         h = (self._edgeCurl.T * (self._MfRho * jSolution) )
         for i, src in enumerate(srcList):
             h[:,i] *= -1./(1j*omega(src.freq))
-            S_m,_ = src.eval(self.prob)
-            h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (S_m)
+            s_m,_ = src.eval(self.prob)
+            h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (s_m)
         return self._MeMuI * h
 
 
     def _hDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the magnetic field with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -856,13 +825,13 @@ class Fields_j(Fields):
         if adjoint:
             return  -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm_v))
         return  -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm_v) )
-        
+
 
 
     def _hDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the magnetic field with respect to the inversion model 
-        
+        Derivative of the magnetic field with respect to the inversion model
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -875,19 +844,19 @@ class Fields_j(Fields):
         C = self._edgeCurl
         MfRho = self._MfRho
         MfRhoDeriv = self._MfRhoDeriv
-        S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
+        s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
 
         if not adjoint:
             hDeriv_m =  -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) )
-            S_mDeriv = S_mDeriv(v)
-            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * ( S_mDeriv)
+            s_mDeriv = s_mDeriv(v)
+            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * ( s_mDeriv)
 
         elif adjoint:
             hDeriv_m =  -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) )
 
-            S_mDeriv = S_mDeriv(MeMuI.T * v)
-            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv
-        
+            s_mDeriv = s_mDeriv(MeMuI.T * v)
+            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * s_mDeriv
+
         return hDeriv_m
 
     def _e(self, jSolution, srcList):
@@ -901,12 +870,12 @@ class Fields_j(Fields):
         """
         n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
-        return VI * (self._aveF2CCV * (self._MfRho *  self._j(jSolution, srcList))) 
+        return VI * (self._aveF2CCV * (self._MfRho *  self._j(jSolution, srcList)))
 
     def _eDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the electric field with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -921,8 +890,8 @@ class Fields_j(Fields):
 
     def _eDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the electric field with respect to the inversion model 
-        
+        Derivative of the electric field with respect to the inversion model
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -943,17 +912,17 @@ class Fields_j(Fields):
         :param numpy.ndarray hSolution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
-        :return: secondary magnetic flux density 
+        :return: secondary magnetic flux density
         """
         n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
-        return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) ) 
+        return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) )
 
     def _bDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the magnetic flux density with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -969,8 +938,8 @@ class Fields_j(Fields):
 
     def _bDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the magnetic flux density with respect to the inversion model 
-        
+        Derivative of the magnetic flux density with respect to the inversion model
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -980,20 +949,20 @@ class Fields_j(Fields):
         jSolution = self[src,'jSolution']
         n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
-        S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
+        s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
 
         if adjoint:
             v = self._aveE2CCV.T * ( VI.T * v)
-            return 1./(1j * omega(src.freq)) * ( S_mDeriv(v) - self._MfRhoDeriv(jSolution).T * (self._edgeCurl * v ))
-        return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( S_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
+            return 1./(1j * omega(src.freq)) * ( s_mDeriv(v) - self._MfRhoDeriv(jSolution).T * (self._edgeCurl * v ))
+        return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
 
 
 class Fields_h(Fields):
     """
-    Fields object for Problem_h. 
+    Fields object for Problem_h.
 
     :param Mesh mesh: mesh
-    :param Survey survey: survey 
+    :param Survey survey: survey
     """
 
     knownFields = {'hSolution':'E'}
@@ -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
@@ -1065,9 +1024,9 @@ class Fields_h(Fields):
 
     def _hDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the total magnetic field with respect to the thing we 
+        Partial derivative of the total magnetic field with respect to the thing we
         solved for.
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -1079,8 +1038,8 @@ class Fields_h(Fields):
 
     def _hDeriv_m(self, src, v, adjoint=False):
         """
-        Partial derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. 
-        
+        Partial derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -1119,14 +1078,14 @@ class Fields_h(Fields):
 
         j = self._edgeCurl*hSolution
         for i, src in enumerate(srcList):
-            _,S_e = src.eval(self.prob)
-            j[:,i] = j[:,i]+ -S_e
+            _,s_e = src.eval(self.prob)
+            j[:,i] = j[:,i]+ -s_e
         return j
 
     def _jDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the current density with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -1142,8 +1101,8 @@ class Fields_h(Fields):
 
     def _jDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the current density with respect to the inversion model. 
-        
+        Derivative of the current density with respect to the inversion model.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -1151,9 +1110,9 @@ class Fields_h(Fields):
         :return: product of the current density derivative with respect to the inversion model with a vector
         """
 
-        _,S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
-        return -S_eDeriv
-    
+        _,s_eDeriv = src.evalDeriv(self.prob, v, adjoint)
+        return -s_eDeriv
+
     def _e(self, hSolution, srcList):
         """
         Electric field from hSolution
@@ -1165,12 +1124,12 @@ class Fields_h(Fields):
         """
         n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
-        return VI * (self._aveF2CCV * (self._MfRho *  self._j(hSolution, srcList))) 
+        return VI * (self._aveF2CCV * (self._MfRho *  self._j(hSolution, srcList)))
 
     def _eDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the electric field with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -1181,12 +1140,12 @@ class Fields_h(Fields):
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         if adjoint:
             return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) )
-        return VI * (self._aveF2CCV * (self._MfRho *  self._edgeCurl * du_dm_v )) 
+        return VI * (self._aveF2CCV * (self._MfRho *  self._edgeCurl * du_dm_v ))
 
     def _eDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the electric field with respect to the inversion model. 
-        
+        Derivative of the electric field with respect to the inversion model.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -1196,7 +1155,7 @@ class Fields_h(Fields):
         hSolution = Utils.mkvc(self[src,'hSolution'])
         n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
-        if adjoint: 
+        if adjoint:
             return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) )
         return VI * (self._aveF2CCV * (self._MfRhoDeriv(self._edgeCurl * hSolution) * v ))
 
@@ -1207,10 +1166,10 @@ class Fields_h(Fields):
         :param numpy.ndarray hSolution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
-        :return: magnetic flux density 
+        :return: magnetic flux density
         """
         h = self._h(hSolution, srcList)
-        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components 
+        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
         return VI * (self._aveE2CCV * (self._MeMu * h))
@@ -1218,14 +1177,14 @@ class Fields_h(Fields):
     def _bDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the magnetic flux density with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
         :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
         """
-        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components 
+        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         if adjoint:
             return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v ))
@@ -1233,8 +1192,8 @@ class Fields_h(Fields):
 
     def _bDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the magnetic flux density with respect to the inversion model. 
-        
+        Derivative of the magnetic flux density with respect to the inversion model.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?

SimPEG/EM/FDEM/SrcFDEM.py

@@ -9,28 +9,30 @@ class BaseSrc(Survey.BaseSrc):
     """
 
     freq = None
-    # rxPair = RxFDEM
-    integrate = True
+    integrate = False
+
+    def __init__(self, rxList, **kwargs):
+        Survey.BaseSrc.__init__(self, rxList, **kwargs)
 
     def eval(self, prob):
         """
         Evaluate the source terms.
-        - :math:`S_m` : magnetic source term
-        - :math:`S_e` : electric source term
+        - :math:`s_m` : magnetic source term
+        - :math:`s_e` : electric source term
 
         :param Problem prob: FDEM Problem
         :rtype: (numpy.ndarray, numpy.ndarray)
         :return: tuple with magnetic source term and electric source term
         """
-        S_m = self.S_m(prob)
-        S_e = self.S_e(prob)
-        return S_m, S_e
+        s_m = self.s_m(prob)
+        s_e = self.s_e(prob)
+        return s_m, s_e
 
     def evalDeriv(self, prob, v=None, adjoint=False):
         """
         Derivatives of the source terms with respect to the inversion model
-        - :code:`S_mDeriv` : derivative of the magnetic source term
-        - :code:`S_eDeriv` : derivative of the electric source term
+        - :code:`s_mDeriv` : derivative of the magnetic source term
+        - :code:`s_eDeriv` : derivative of the electric source term
 
         :param Problem prob: FDEM Problem
         :param numpy.ndarray v: vector to take product with
@@ -39,9 +41,9 @@ class BaseSrc(Survey.BaseSrc):
         :return: tuple with magnetic source term and electric source term derivatives times a vector
         """
         if v is not None:
-            return self.S_mDeriv(prob, v, adjoint), self.S_eDeriv(prob, v, adjoint)
+            return self.s_mDeriv(prob, v, adjoint), self.s_eDeriv(prob, v, adjoint)
         else:
-            return lambda v: self.S_mDeriv(prob, v, adjoint), lambda v: self.S_eDeriv(prob, v, adjoint)
+            return lambda v: self.s_mDeriv(prob, v, adjoint), lambda v: self.s_eDeriv(prob, v, adjoint)
 
     def bPrimary(self, prob):
         """
@@ -83,7 +85,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         return Zero()
 
-    def S_m(self, prob):
+    def s_m(self, prob):
         """
         Magnetic source term
 
@@ -93,7 +95,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         return Zero()
 
-    def S_e(self, prob):
+    def s_e(self, prob):
         """
         Electric source term
 
@@ -103,7 +105,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 +118,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
 
@@ -131,22 +133,21 @@ class BaseSrc(Survey.BaseSrc):
 
 class RawVec_e(BaseSrc):
     """
-    RawVec electric source. It is defined by the user provided vector S_e
+    RawVec electric source. It is defined by the user provided vector s_e
 
     :param list rxList: receiver list
     :param float freq: frequency
-    :param numpy.array S_e: electric source term
-    :param bool integrate: Integrate the source term (multiply by Me) [True]
+    :param numpy.array s_e: electric source term
+    :param bool integrate: Integrate the source term (multiply by Me) [False]
     """
 
-    def __init__(self, rxList, freq, S_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
-        self._S_e = np.array(S_e, dtype=complex)
+    def __init__(self, rxList, freq, s_e):
+        self._s_e = np.array(s_e, dtype=complex)
         self.freq = float(freq)
-        self.integrate = integrate
 
         BaseSrc.__init__(self, rxList)
 
-    def S_e(self, prob):
+    def s_e(self, prob):
         """
         Electric source term
 
@@ -155,28 +156,27 @@ class RawVec_e(BaseSrc):
         :return: electric source term on mesh
         """
         if prob._formulation is 'EB' and self.integrate is True:
-            return prob.Me * self._S_e
-        return self._S_e
+            return prob.Me * self._s_e
+        return self._s_e
 
 
 class RawVec_m(BaseSrc):
     """
-    RawVec magnetic source. It is defined by the user provided vector S_m
+    RawVec magnetic source. It is defined by the user provided vector s_m
 
     :param float freq: frequency
     :param rxList: receiver list
-    :param numpy.array S_m: magnetic source term
-    :param bool integrate: Integrate the source term (multiply by Me) [True]
+    :param numpy.array s_m: magnetic source term
+    :param bool integrate: Integrate the source term (multiply by Me) [False]
     """
 
-    def __init__(self, rxList, freq, S_m, integrate=True):  #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
-        self._S_m = np.array(S_m, dtype=complex)
+    def __init__(self, rxList, freq, s_m, integrate=True):  #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
+        self._s_m = np.array(s_m, dtype=complex)
         self.freq = float(freq)
-        self.integrate = integrate
 
         BaseSrc.__init__(self, rxList)
 
-    def S_m(self, prob):
+    def s_m(self, prob):
         """
         Magnetic source term
 
@@ -185,28 +185,27 @@ class RawVec_m(BaseSrc):
         :return: magnetic source term on mesh
         """
         if prob._formulation is 'HJ' and self.integrate is True:
-            return prob.Me * self._S_m
-        return self._S_m
+            return prob.Me * self._s_m
+        return self._s_m
 
 
 class RawVec(BaseSrc):
     """
-    RawVec source. It is defined by the user provided vectors S_m, S_e
+    RawVec source. It is defined by the user provided vectors s_m, s_e
 
     :param rxList: receiver list
     :param float freq: frequency
-    :param numpy.array S_m: magnetic source term
-    :param numpy.array S_e: electric source term
-    :param bool integrate: Integrate the source term (multiply by Me) [True]
+    :param numpy.array s_m: magnetic source term
+    :param numpy.array s_e: electric source term
+    :param bool integrate: Integrate the source term (multiply by Me) [False]
     """
-    def __init__(self, rxList, freq, S_m, S_e, integrate=True):
-        self._S_m = np.array(S_m, dtype=complex)
-        self._S_e = np.array(S_e, dtype=complex)
+    def __init__(self, rxList, freq, s_m, s_e, **kwargs):
+        self._s_m = np.array(s_m, dtype=complex)
+        self._s_e = np.array(s_e, dtype=complex)
         self.freq = float(freq)
-        self.integrate = integrate
-        BaseSrc.__init__(self, rxList)
+        BaseSrc.__init__(self, rxList, **kwargs)
 
-    def S_m(self, prob):
+    def s_m(self, prob):
         """
         Magnetic source term
 
@@ -215,10 +214,10 @@ class RawVec(BaseSrc):
         :return: magnetic source term on mesh
         """
         if prob._formulation is 'HJ' and self.integrate is True:
-            return prob.Me * self._S_m
-        return self._S_m
+            return prob.Me * self._s_m
+        return self._s_m
 
-    def S_e(self, prob):
+    def s_e(self, prob):
         """
         Electric source term
 
@@ -227,8 +226,8 @@ class RawVec(BaseSrc):
         :return: electric source term on mesh
         """
         if prob._formulation is 'EB' and self.integrate is True:
-            return prob.Me * self._S_e
-        return self._S_e
+            return prob.Me * self._s_e
+        return self._s_e
 
 
 class MagDipole(BaseSrc):
@@ -278,14 +277,13 @@ class MagDipole(BaseSrc):
     :param float mu: background magnetic permeability
     """
 
-    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0):
+    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0, **kwargs):
         self.freq = float(freq)
         self.loc = loc
         self.orientation = orientation
         assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
         self.moment = moment
         self.mu = mu
-        self.integrate = False
         BaseSrc.__init__(self, rxList)
 
     def bPrimary(self, prob):
@@ -335,9 +333,9 @@ class MagDipole(BaseSrc):
         :return: primary magnetic field
         """
         b = self.bPrimary(prob)
-        return 1./self.mu * b 
+        return 1./self.mu * b
 
-    def S_m(self, prob):
+    def s_m(self, prob):
         """
         The magnetic source term
 
@@ -348,10 +346,10 @@ class MagDipole(BaseSrc):
 
         b_p = self.bPrimary(prob)
         if prob._formulation is 'HJ':
-            b_p = prob.Me * b_p 
+            b_p = prob.Me * b_p
         return -1j*omega(self.freq)*b_p
 
-    def S_e(self, prob):
+    def s_e(self, prob):
         """
         The electric source term
 
@@ -453,7 +451,7 @@ class MagDipole_Bfield(BaseSrc):
         b = self.bPrimary(prob)
         return 1/self.mu * b
 
-    def S_m(self, prob):
+    def s_m(self, prob):
         """
         The magnetic source term
 
@@ -466,7 +464,7 @@ class MagDipole_Bfield(BaseSrc):
             b = prob.Me * b
         return -1j*omega(self.freq)*b
 
-    def S_e(self, prob):
+    def s_e(self, prob):
         """
         The electric source term
 
@@ -543,10 +541,10 @@ class CircularLoop(BaseSrc):
             if not prob.mesh.isSymmetric:
                 # TODO ?
                 raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
+            a = MagneticLoopVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
 
         else:
-            srcfct = MagneticDipoleVectorPotential
+            srcfct = MagneticLoopVectorPotential
             ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
             ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
             az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)
@@ -565,7 +563,7 @@ class CircularLoop(BaseSrc):
         b = self.bPrimary(prob)
         return 1./self.mu*b
 
-    def S_m(self, prob):
+    def s_m(self, prob):
         """
         The magnetic source term
 
@@ -578,7 +576,7 @@ class CircularLoop(BaseSrc):
             b =  prob.Me *  b
         return -1j*omega(self.freq)*b
 
-    def S_e(self, prob):
+    def s_e(self, prob):
         """
         The electric source term
 
@@ -604,6 +602,232 @@ class CircularLoop(BaseSrc):
 
             return -C.T * (MMui_s * self.bPrimary(prob))
 
-            
+
+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.
+
+    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.m = m
+        self.prob = prob
+        self.survey = survey
+        self.fields = None
+
+        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.mesh = self.prob.mesh
+        self.prob.curModel = self.m
+
+        BaseSrc.__init__(self, rxList)
+
+    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._MeSigma = prob.mesh.getEdgeInnerProduct(sigmaprimary)
+        return self._MeSigma
+
+    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._MfRho = prob.mesh.getFaceInnerProduct(rhoprimary)
+        return self._MfRho
+
+    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):
+        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:
+                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
+
+    # 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):
+        if prob._formulation == 'EB':
+            # - \\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):
+        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()
+
 
 

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -1,5 +1,6 @@
 import SimPEG
 from SimPEG.EM.Utils import *
+from SimPEG.EM.Base import BaseEMSurvey
 from scipy.constants import mu_0
 from SimPEG.Utils import Zero, Identity
 import SrcFDEM as Src
@@ -64,9 +65,9 @@ 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):
+    def eval(self, src, mesh, f):
         """
         Project fields to recievers to get data.
 
@@ -76,30 +77,28 @@ 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]
+        P = self.getP(mesh, self.projGLoc(f))
+        f_part_complex = f[src, self.projField]
         # get the real or imag component
         real_or_imag = self.projComp
-        u_part = getattr(u_part_complex, real_or_imag)
-        
-        return P*u_part
+        f_part = getattr(f_part_complex, real_or_imag)
 
-    def evalDeriv(self, src, mesh, u, v, adjoint=False):
+        return P*f_part
+
+    def evalDeriv(self, src, mesh, f, v, adjoint=False):
         """
         Derivative of projected fields with respect to the inversion model times a vector.
 
         :param Source src: FDEM source
         :param Mesh mesh: mesh used
-        :param Fields u: fields object
+        :param Fields f: fields object
         :param numpy.ndarray v: vector to multiply
         :rtype: numpy.ndarray
         :return: fields projected to recievers
         """
 
-        P = self.getP(mesh, self.projGLoc(u))
+        P = self.getP(mesh, self.projGLoc(f))
 
         if not adjoint:
             Pv_complex = P * v
@@ -123,7 +122,7 @@ class Rx(SimPEG.Survey.BaseRx):
 # Survey
 ####################################################
 
-class Survey(SimPEG.Survey.BaseSurvey):
+class Survey(BaseEMSurvey):
     """
     Frequency domain electromagnetic survey
 
@@ -131,12 +130,12 @@ class Survey(SimPEG.Survey.BaseSurvey):
     """
 
     srcPair = Src.BaseSrc
-    rxPaair = Rx 
+    rxPair = Rx
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency
         self.srcList = srcList
-        SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
+        BaseEMSurvey.__init__(self, srcList, **kwargs)
 
         _freqDict = {}
         for src in srcList:
@@ -171,24 +170,8 @@ class Survey(SimPEG.Survey.BaseSurvey):
         Returns the sources associated with a specific frequency.
         :param float freq: frequency for which we look up sources
         :rtype: dictionary
-        :return: sources at the sepcified frequency 
+        :return: sources at the sepcified frequency
         """
         assert freq in self._freqDict, "The requested frequency is not in this survey."
         return self._freqDict[freq]
 
-    def eval(self, u):
-        """
-        Project fields to receiver locations
-        :param Fields u: fields object
-        :rtype: numpy.ndarray
-        :return: data
-        """
-        data = SimPEG.Survey.Data(self)
-        for src in self.srcList:
-            for rx in src.rxList:
-                data[src, rx] = rx.eval(src, self.mesh, u)
-        return data
-
-    def evalDeriv(self, u):
-        raise Exception('Use Receivers to project fields deriv.')
-

SimPEG/EM/TDEM/BaseTDEM.py

@@ -108,11 +108,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         Ainv.clean()
         return F
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         """
             :param numpy.array m: Conductivity model
             :param numpy.ndarray v: vector (model object)
-            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
+            :param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m
             :rtype: numpy.ndarray
             :return: w (data object)
 
@@ -125,15 +125,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         """
         if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
         self.curModel = m
-        if u is None:
-            u = self.fields(m)
-        p = self.Gvec(m, v, u)
+        if f is None:
+            f = self.fields(m)
+        p = self.Gvec(m, v, f)
         y = self.solveAh(m, p)
-        Jv = self.survey.evalDeriv(u, v=y)
+        Jv = self.survey.evalDeriv(f, v=y)
         if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
         return - mkvc(Jv)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
         """
             :param numpy.array m: Conductivity model
             :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object)
@@ -150,15 +150,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         """
         if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50)
         self.curModel = m
-        if u is None:
-            u = self.fields(m)
+        if f is None:
+            f = self.fields(m)
 
         if not isinstance(v, self.dataPair):
             v = self.dataPair(self.survey, v)
 
-        p = self.survey.evalDeriv(u, v=v, adjoint=True)
+        p = self.survey.evalDeriv(f, v=v, adjoint=True)
         y = self.solveAht(m, p)
-        w = self.Gtvec(m, y, u)
+        w = self.Gtvec(m, y, f)
         if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
         return - mkvc(w)
 

SimPEG/EM/Utils/testingUtils.py

@@ -20,7 +20,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
     mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
 
     if useMu is True:
-        mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] 
+        mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))]
     else:
         mapping = Maps.ExpMap(mesh)
 
@@ -37,20 +37,39 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
             Src.append(EM.FDEM.Src.MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
         elif SrcType is 'CircularLoop':
             Src.append(EM.FDEM.Src.CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
+
         elif SrcType is 'RawVec':
             if fdemType is 'e' or fdemType is 'b':
                 S_m = np.zeros(mesh.nF)
                 S_e = np.zeros(mesh.nE)
                 S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
                 S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
-                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
+                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e, integrate=True))
 
             elif fdemType is 'h' or fdemType is 'j':
                 S_m = np.zeros(mesh.nE)
                 S_e = np.zeros(mesh.nF)
                 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))
+                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e, integrate=True))
+
+        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.PrimSec([Rx0], freq, mPrimary, prob=primaryProblem, survey=primarySurvey))
+
+        elif SrcType is 'PrimSecCyl':
+                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')
+
+                primSrc = EM.FDEM.Src.MagDipole([], freq, np.r_[0.,0.,0.])
+                primarySurvey = EM.FDEM.Survey([primSrc])
+                primaryProblem = EM.FDEM.Problem_e(primmesh)
+                mPrimary = np.ones(primmesh.nC)*CONDUCTIVITY
+                Src.append(EM.FDEM.Src.PrimSec([Rx0], freq, mPrimary, prob=primaryProblem, survey=primarySurvey))
 
     if verbose:
         print '  Fetching %s problem' % (fdemType)
@@ -90,7 +109,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
     prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
     mesh = prb1.mesh
     print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
-    
+
     logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
     mu = np.ones(mesh.nC)*MU
 

SimPEG/FLOW/Richards/RichardsProblem.py

@@ -45,19 +45,19 @@ class RichardsSurvey(Survey.BaseSurvey):
 
     @Utils.count
     @Utils.requires('prob')
-    def dpred(self, m, u=None):
+    def dpred(self, m, f=None):
         """
             Create the projected data from a model.
-            The field, u, (if provided) will be used for the predicted data
+            The field, f, (if provided) will be used for the predicted data
             instead of recalculating the fields (which may be expensive!).
 
             .. math::
-                d_\\text{pred} = P(u(m), m)
+                d_\\text{pred} = P(f(m), m)
 
             Where P is a projection of the fields onto the data space.
         """
-        if u is None: u = self.prob.fields(m)
-        return Utils.mkvc(self.eval(u, m))
+        if f is None: f = self.prob.fields(m)
+        return Utils.mkvc(self.eval(f, m))
 
     @Utils.requires('prob')
     def eval(self, U, m):
@@ -233,16 +233,16 @@ class RichardsProblem(Problem.BaseTimeProblem):
         return r, J
 
     @Utils.timeIt
-    def Jfull(self, m, u=None):
-        if u is None:
-            u = self.fields(m)
+    def Jfull(self, m, f=None):
+        if f is None:
+            f = self.fields(m)
 
-        nn = len(u)-1
+        nn = len(f)-1
         Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
         for ii in range(nn):
             dt = self.timeSteps[ii]
-            bc = self.getBoundaryConditions(ii, u[ii])
-            Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, u[ii], u[ii+1], dt, bc)
+            bc = self.getBoundaryConditions(ii, f[ii])
+            Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, f[ii], f[ii+1], dt, bc)
         Ad = sp.block_diag(Adiags)
         zRight = Utils.spzeros((len(Asubs)-1)*Asubs[0].shape[0],Adiags[0].shape[1])
         zTop = Utils.spzeros(Adiags[0].shape[0], len(Adiags)*Adiags[0].shape[1])
@@ -251,7 +251,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
         B = np.array(sp.vstack(Bs).todense())
 
         Ainv = self.Solver(A, **self.solverOpts)
-        P = self.survey.evalDeriv(u, m)
+        P = self.survey.evalDeriv(f, m)
         AinvB = Ainv * B
         z = np.zeros((self.mesh.nC, B.shape[1]))
         zAinvB = np.vstack((z, AinvB))
@@ -259,41 +259,41 @@ class RichardsProblem(Problem.BaseTimeProblem):
         return J
 
     @Utils.timeIt
-    def Jvec(self, m, v, u=None):
-        if u is None:
-            u = self.fields(m)
+    def Jvec(self, m, v, f=None):
+        if f is None:
+            f = self.fields(m)
 
-        JvC = range(len(u)-1) # Cell to hold each row of the long vector.
+        JvC = range(len(f)-1) # Cell to hold each row of the long vector.
 
         # This is done via forward substitution.
-        bc = self.getBoundaryConditions(0, u[0])
-        temp, Adiag, B = self.diagsJacobian(m, u[0], u[1], self.timeSteps[0], bc)
+        bc = self.getBoundaryConditions(0, f[0])
+        temp, Adiag, B = self.diagsJacobian(m, f[0], f[1], self.timeSteps[0], bc)
         Adiaginv = self.Solver(Adiag, **self.solverOpts)
         JvC[0] = Adiaginv * (B*v)
 
-        for ii in range(1,len(u)-1):
-            bc = self.getBoundaryConditions(ii, u[ii])
-            Asub, Adiag, B = self.diagsJacobian(m, u[ii], u[ii+1], self.timeSteps[ii], bc)
+        for ii in range(1,len(f)-1):
+            bc = self.getBoundaryConditions(ii, f[ii])
+            Asub, Adiag, B = self.diagsJacobian(m, f[ii], f[ii+1], self.timeSteps[ii], bc)
             Adiaginv = self.Solver(Adiag, **self.solverOpts)
             JvC[ii] = Adiaginv * (B*v - Asub*JvC[ii-1])
 
-        P = self.survey.evalDeriv(u, m)
+        P = self.survey.evalDeriv(f, m)
         return P * np.concatenate([np.zeros(self.mesh.nC)] + JvC)
 
     @Utils.timeIt
-    def Jtvec(self, m, v, u=None):
-        if u is None:
-            u = self.field(m)
+    def Jtvec(self, m, v, f=None):
+        if f is None:
+            f = self.field(m)
 
-        P = self.survey.evalDeriv(u, m)
+        P = self.survey.evalDeriv(f, m)
         PTv = P.T*v
 
         # This is done via backward substitution.
         minus = 0
         BJtv = 0
-        for ii in range(len(u)-1,0,-1):
-            bc = self.getBoundaryConditions(ii-1, u[ii-1])
-            Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii], self.timeSteps[ii-1], bc)
+        for ii in range(len(f)-1,0,-1):
+            bc = self.getBoundaryConditions(ii-1, f[ii-1])
+            Asub, Adiag, B = self.diagsJacobian(m, f[ii-1], f[ii], self.timeSteps[ii-1], bc)
             #select the correct part of v
             vpart = range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0])
             AdiaginvT = self.Solver(Adiag.T, **self.solverOpts)

SimPEG/InvProblem.py

@@ -82,23 +82,23 @@ class BaseInvProblem(object):
         self._warmstart = value
 
     def getFields(self, m, store=False, deleteWarmstart=True):
-        u = None
+        f = None
 
         for mtest, u_ofmtest in self.warmstart:
             if m is mtest:
-                u = u_ofmtest
+                f = u_ofmtest
                 if self.debug: print 'InvProb is Warm Starting!'
                 break
 
-        if u is None:
-            u = self.prob.fields(m)
+        if f is None:
+            f = self.prob.fields(m)
 
         if deleteWarmstart:
             self.warmstart = []
         if store:
-            self.warmstart += [(m,u)]
+            self.warmstart += [(m,f)]
 
-        return u
+        return f
 
     @Utils.timeIt
     def evalFunction(self, m, return_g=True, return_H=True):
@@ -109,21 +109,21 @@ class BaseInvProblem(object):
         gc.collect()
 
         # Store fields if doing a line-search
-        u = self.getFields(m, store=(return_g==False and return_H==False))
+        f = self.getFields(m, store=(return_g==False and return_H==False))
 
-        phi_d = self.dmisfit.eval(m, u=u)
+        phi_d = self.dmisfit.eval(m, f=f)
         phi_m = self.reg.eval(m)
 
-        self.dpred = self.survey.dpred(m, u=u)  # This is a cheap matrix vector calculation.
+        self.dpred = self.survey.dpred(m, f=f)  # This is a cheap matrix vector calculation.
 
         self.phi_d, self.phi_d_last  = phi_d, self.phi_d
         self.phi_m, self.phi_m_last  = phi_m, self.phi_m
 
-        f = phi_d + self.beta * phi_m
+        phi = phi_d + self.beta * phi_m
 
-        out = (f,)
+        out = (phi,)
         if return_g:
-            phi_dDeriv = self.dmisfit.evalDeriv(m, u=u)
+            phi_dDeriv = self.dmisfit.evalDeriv(m, f=f)
             phi_mDeriv = self.reg.evalDeriv(m)
 
             g = phi_dDeriv + self.beta * phi_mDeriv
@@ -131,7 +131,7 @@ class BaseInvProblem(object):
 
         if return_H:
             def H_fun(v):
-                phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, u=u)
+                phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, f=f)
                 phi_m2Deriv = self.reg.eval2Deriv(m, v=v)
 
                 return phi_d2Deriv + self.beta * phi_m2Deriv

SimPEG/MT/BaseMT.py

@@ -27,7 +27,7 @@ class BaseMTProblem(BaseFDEMProblem):
     # Might need to add more stuff here.
 
     ## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         """
         Function to calculate the data sensitivities dD/dm times a vector.
 
@@ -39,8 +39,8 @@ class BaseMTProblem(BaseFDEMProblem):
         """
 
         # Calculate the fields
-        if u is None:
-           u = self.fields(m)
+        if f is None:
+           f= self.fields(m)
         # Set current model
         self.curModel = m
         # Initiate the Jv object
@@ -56,9 +56,9 @@ class BaseMTProblem(BaseFDEMProblem):
                 # We need fDeriv_m = df/du*du/dm + df/dm
                 # Construct du/dm, it requires a solve
                 # NOTE: need to account for the 2 polarizations in the derivatives.
-                u_src = u[src,:]
+                f_src = f[src,:]
                 # dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
-                dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns.
+                dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns.
                 dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
                 if dRHS_dm is None:
                     du_dm = dA_duI * ( -dA_dm )
@@ -68,13 +68,13 @@ class BaseMTProblem(BaseFDEMProblem):
                 for rx in src.rxList:
                     # Get the projection derivative
                     # v should be of size 2*nE (for 2 polarizations)
-                    PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m
+                    PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m
                     Jv[src, rx] = PDeriv_u(mkvc(du_dm))
             dA_duI.clean()
         # Return the vectorized sensitivities
         return mkvc(Jv)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
         """
         Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
 
@@ -85,8 +85,8 @@ class BaseMTProblem(BaseFDEMProblem):
             :return: Data sensitivities wrt m
         """
 
-        if u is None:
-            u = self.fields(m)
+        if f is None:
+            f = self.fields(m)
 
         self.curModel = m
 
@@ -103,15 +103,15 @@ class BaseMTProblem(BaseFDEMProblem):
 
             for src in self.survey.getSrcByFreq(freq):
                 ftype = self._fieldType + 'Solution'
-                u_src = u[src, :]
+                f_src = f[src, :]
 
                 for rx in src.rxList:
                     # Get the adjoint evalDeriv
                     # PTv needs to be nE,
-                    PTv = rx.evalDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
+                    PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
                     # Get the
                     dA_duIT = ATinv * PTv
-                    dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True)
+                    dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
                     dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
                     # Make du_dmT
                     if dRHS_dmT is None:
@@ -129,4 +129,4 @@ class BaseMTProblem(BaseFDEMProblem):
                         raise Exception('Must be real or imag')
             # Clean the factorization, clear memory.
             ATinv.clean()
-        return Jtv
+        return Jtv

SimPEG/MT/SurveyMT.py

@@ -427,15 +427,15 @@ class Survey(SimPEGsurvey.BaseSurvey):
         assert freq in self._freqDict, "The requested frequency is not in this survey."
         return self._freqDict[freq]
 
-    def eval(self, u):
+    def eval(self, f):
         data = Data(self)
         for src in self.srcList:
             sys.stdout.flush()
             for rx in src.rxList:
-                data[src, rx] = rx.eval(src, self.mesh, u)
+                data[src, rx] = rx.eval(src, self.mesh, f)
         return data
 
-    def evalDeriv(self, u):
+    def evalDeriv(self, f):
         raise Exception('Use Transmitters to project fields deriv.')
 
 #################

SimPEG/Mesh/BaseMesh.py

@@ -1,4 +1,5 @@
 import numpy as np
+import scipy.sparse as sp
 from SimPEG import Utils
 
 
@@ -594,3 +595,63 @@ class BaseRectangularMesh(BaseMesh):
             return out
         else:
             return switchKernal(x)
+
+
+    def getInterpolationMatMesh2Mesh(self, mesh2, locType='CC', locTypeFrom=None):
+        """
+        Interpolates variables from the current mesh to a new mesh (mesh2)
+
+        :param Mesh mesh2: SimPEG mesh which we interpolate values to
+        :param string locType: location of variables 'CC', 'E', 'F', 'N'
+        :rtype: scipy.sparse.csr_matrix
+        :return P: interpolation matrix
+        """
+
+        # import warnings
+        # warnings.warn(
+        #      "`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"
+        if mesh2._meshType == 'CYL':
+            assert self._meshType == 'CYL', "Interpolation from 3D mesh to Cyl mesh is not supported"
+
+        # if Cyl to cart call
+        if self._meshType == 'CYL' and mesh2._meshType != 'CYL':
+            return self.getInterpolationMatCartMesh(mesh2, locType)
+
+        # Scalars
+        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', 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', locTypeFrom=locTypeFrom+'y')
+
+            if self.dim == 1:
+                return self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType, locTypeFrom=locTypeFrom+'x')
+            elif self.dim == 2:
+                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, 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])
+
+
+
+

SimPEG/Problem.py

@@ -88,28 +88,28 @@ class BaseProblem(object):
         return self.survey is not None
 
     @Utils.timeIt
-    def Jvec(self, m, v, u=None):
-        """Jvec(m, v, u=None)
+    def Jvec(self, m, v, f=None):
+        """Jvec(m, v, f=None)
 
             Effect of J(m) on a vector v.
 
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: Jv
         """
         raise NotImplementedError('J is not yet implemented.')
 
     @Utils.timeIt
-    def Jtvec(self, m, v, u=None):
-        """Jtvec(m, v, u=None)
+    def Jtvec(self, m, v, f=None):
+        """Jtvec(m, v, f=None)
 
             Effect of transpose of J(m) on a vector v.
 
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: JTv
         """
@@ -117,32 +117,32 @@ class BaseProblem(object):
 
 
     @Utils.timeIt
-    def Jvec_approx(self, m, v, u=None):
-        """Jvec_approx(m, v, u=None)
+    def Jvec_approx(self, m, v, f=None):
+        """Jvec_approx(m, v, f=None)
 
             Approximate effect of J(m) on a vector v
 
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: approxJv
         """
-        return self.Jvec(m, v, u)
+        return self.Jvec(m, v, f)
 
     @Utils.timeIt
-    def Jtvec_approx(self, m, v, u=None):
-        """Jtvec_approx(m, v, u=None)
+    def Jtvec_approx(self, m, v, f=None):
+        """Jtvec_approx(m, v, f=None)
 
             Approximate effect of transpose of J(m) on a vector v.
 
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: JTv
         """
-        return self.Jtvec(m, v, u)
+        return self.Jtvec(m, v, f)
 
     def fields(self, m):
         """
@@ -224,9 +224,9 @@ class LinearProblem(BaseProblem):
     def fields(self, m):
         return self.G.dot(m)
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         return self.G.dot(v)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
         return self.G.T.dot(v)
 

SimPEG/Survey.py

@@ -295,38 +295,38 @@ class BaseSurvey(object):
 
     @Utils.count
     @Utils.requires('prob')
-    def dpred(self, m, u=None):
-        """dpred(m, u=None)
+    def dpred(self, m, f=None):
+        """dpred(m, f=None)
 
             Create the projected data from a model.
-            The field, u, (if provided) will be used for the predicted data
+            The fields, f, (if provided) will be used for the predicted data
             instead of recalculating the fields (which may be expensive!).
 
             .. math::
 
-                d_\\text{pred} = P(u(m))
+                d_\\text{pred} = P(f(m))
 
             Where P is a projection of the fields onto the data space.
         """
-        if u is None: u = self.prob.fields(m)
-        return Utils.mkvc(self.eval(u))
+        if f is None: f = self.prob.fields(m)
+        return Utils.mkvc(self.eval(f))
 
 
     @Utils.count
-    def eval(self, u):
-        """eval(u)
+    def eval(self, f):
+        """eval(f)
 
             This function projects the fields onto the data space.
 
             .. math::
 
-                d_\\text{pred} = \mathbf{P} u(m)
+                d_\\text{pred} = \mathbf{P} f(m)
         """
         raise NotImplemented('eval is not yet implemented.')
 
     @Utils.count
-    def evalDeriv(self, u):
-        """evalDeriv(u)
+    def evalDeriv(self, f):
+        """evalDeriv(f)
 
             This function s the derivative of projects the fields onto the data space.
 
@@ -337,11 +337,11 @@ class BaseSurvey(object):
         raise NotImplemented('eval is not yet implemented.')
 
     @Utils.count
-    def residual(self, m, u=None):
-        """residual(m, u=None)
+    def residual(self, m, f=None):
+        """residual(m, f=None)
 
             :param numpy.array m: geophysical model
-            :param numpy.array u: fields
+            :param numpy.array f: fields
             :rtype: numpy.array
             :return: data residual
 
@@ -352,14 +352,14 @@ class BaseSurvey(object):
                 \mu_\\text{data} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}
 
         """
-        return Utils.mkvc(self.dpred(m, u=u) - self.dobs)
+        return Utils.mkvc(self.dpred(m, f=f) - self.dobs)
 
     @property
     def isSynthetic(self):
         "Check if the data is synthetic."
         return self.mtrue is not None
 
-    def makeSyntheticData(self, m, std=0.05, u=None, force=False):
+    def makeSyntheticData(self, m, std=0.05, f=None, force=False):
         """
             Make synthetic data given a model, and a standard deviation.
 
@@ -372,16 +372,16 @@ class BaseSurvey(object):
         if getattr(self, 'dobs', None) is not None and not force:
             raise Exception('Survey already has dobs. You can use force=True to override this exception.')
         self.mtrue = m
-        self.dtrue = self.dpred(m, u=u)
+        self.dtrue = self.dpred(m, f=f)
         noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape)
         self.dobs = self.dtrue+noise
         self.std = self.dobs*0 + std
         return self.dobs
 
 class LinearSurvey(BaseSurvey):
-    def eval(self, u):
-        return u
-    
+    def eval(self, f):
+        return f
+
     @property
     def nD(self):
         return self.prob.G.shape[0]

SimPEG/Tests.py

@@ -82,14 +82,14 @@ class OrderTest(unittest.TestCase):
     _meshType = meshTypes[0]
     meshDimension = 3
 
-    def setupMesh(self, nc):
+    def makeMesh(self, nc, meshType=_meshType, meshDimension=meshDimension):
         """
         For a given number of cells nc, generate a TensorMesh with uniform cells with edge length h=1/nc.
         """
-        if 'TensorMesh' in self._meshType:
-            if 'uniform' in self._meshType:
+        if 'TensorMesh' in meshType:
+            if 'uniform' in meshType:
                 h = [nc, nc, nc]
-            elif 'random' in self._meshType:
+            elif 'random' in meshType:
                 h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
                 h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
                 h3 = np.random.rand(nc)*nc*0.5 + nc*0.5
@@ -97,46 +97,46 @@ class OrderTest(unittest.TestCase):
             else:
                 raise Exception('Unexpected meshType')
 
-            self.M = TensorMesh(h[:self.meshDimension])
-            max_h = max([np.max(hi) for hi in self.M.h])
-            return max_h
+            M = TensorMesh(h[:meshDimension])
+            max_h = max([np.max(hi) for hi in M.h])
+            return M, max_h
 
-        elif 'CylMesh' in self._meshType:
-            if 'uniform' in self._meshType:
+        elif 'CylMesh' in meshType:
+            if 'uniform' in meshType:
                 h = [nc, nc, nc]
             else:
                 raise Exception('Unexpected meshType')
 
-            if self.meshDimension == 2:
-                self.M = CylMesh([h[0], 1, h[2]])
-                max_h = max([np.max(hi) for hi in [self.M.hx, self.M.hz]])
-            elif self.meshDimension == 3:
-                self.M = CylMesh(h)
-                max_h = max([np.max(hi) for hi in self.M.h])
-            return max_h
+            if meshDimension == 2:
+                M = CylMesh([h[0], 1, h[2]])
+                max_h = max([np.max(hi) for hi in [M.hx, M.hz]])
+            elif meshDimension == 3:
+                M = CylMesh(h)
+                max_h = max([np.max(hi) for hi in M.h])
+            return M, max_h
 
-        elif 'Curv' in self._meshType:
-            if 'uniform' in self._meshType:
+        elif 'Curv' in meshType:
+            if 'uniform' in meshType:
                 kwrd = 'rect'
-            elif 'rotate' in self._meshType:
+            elif 'rotate' in meshType:
                 kwrd = 'rotate'
             else:
                 raise Exception('Unexpected meshType')
-            if self.meshDimension == 1:
+            if meshDimension == 1:
                 raise Exception('Lom not supported for 1D')
-            elif self.meshDimension == 2:
+            elif meshDimension == 2:
                 X, Y = Utils.exampleLrmGrid([nc, nc], kwrd)
-                self.M = CurvilinearMesh([X, Y])
-            elif self.meshDimension == 3:
+                M = CurvilinearMesh([X, Y])
+            elif meshDimension == 3:
                 X, Y, Z = Utils.exampleLrmGrid([nc, nc, nc], kwrd)
-                self.M = CurvilinearMesh([X, Y, Z])
-            return 1./nc
+                M = CurvilinearMesh([X, Y, Z])
+            return M, 1./nc
 
-        elif 'Tree' in self._meshType:
+        elif 'Tree' in meshType:
             nc *= 2
-            if 'uniform' in self._meshType or 'notatree' in self._meshType:
+            if 'uniform' in meshType or 'notatree' in meshType:
                 h = [nc, nc, nc]
-            elif 'random' in self._meshType:
+            elif 'random' in meshType:
                 h1 = np.random.rand(nc)*nc*0.5 + nc*0.5
                 h2 = np.random.rand(nc)*nc*0.5 + nc*0.5
                 h3 = np.random.rand(nc)*nc*0.5 + nc*0.5
@@ -145,20 +145,29 @@ class OrderTest(unittest.TestCase):
                 raise Exception('Unexpected meshType')
 
             levels = int(np.log(nc)/np.log(2))
-            self.M = Tree(h[:self.meshDimension], levels=levels)
+            M = Tree(h[:meshDimension], levels=levels)
             def function(cell):
-                if 'notatree' in self._meshType:
+                if 'notatree' in meshType:
                     return levels - 1
                 r = cell.center - np.array([0.5]*len(cell.center))
                 dist = np.sqrt(r.dot(r))
                 if dist < 0.2:
                     return levels
                 return levels - 1
-            self.M.refine(function,balance=False)
-            self.M.number(balance=False)
-            # self.M.plotGrid(showIt=True)
-            max_h = max([np.max(hi) for hi in self.M.h])
-            return max_h
+            M.refine(function,balance=False)
+            M.number(balance=False)
+            # M.plotGrid(showIt=True)
+            max_h = max([np.max(hi) for hi in M.h])
+            return M, max_h
+
+
+    def setupMesh(self, nc):
+        """
+        For a given number of cells nc, generate a TensorMesh with uniform cells with edge length h=1/nc.
+        """
+        M, h = self.makeMesh(nc, meshType=self._meshType, meshDimension=self.meshDimension)
+        self.M = M
+        return h
 
     def getError(self):
         """For given h, generate A[h], f and A(f) and return norm of error."""

SimPEG/Utils/matutils.py

@@ -27,7 +27,7 @@ def mkvc(x, numDims=1):
 
     if isinstance(x, Zero):
         return x
-    
+
     assert isinstance(x, np.ndarray), "Vector must be a numpy array"
 
     if numDims == 1:
@@ -422,9 +422,9 @@ class Zero(object):
     def __ge__(self, v):return 0 >= v
     def __gt__(self, v):return 0 > v
 
-    @property 
+    @property
     def transpose(self): return Zero()
-    
+
     @property
     def T(self): return Zero()
 

tests/em/fdem/inverse/derivs/test_FDEM_derivs.py

@@ -8,8 +8,8 @@ from SimPEG.EM.Utils.testingUtils import getFDEMProblem
 
 testE = True
 testB = True
-testH = True
-testJ = True
+testH = False
+testJ = False
 
 verbose = False
 
@@ -20,8 +20,8 @@ MU = mu_0
 freq = 1e-1
 addrandoms = True
 
-SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
-
+# SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
+SrcType = ['PrimSecCyl']
 
 def derivTest(fdemType, comp):
 

tests/flow/test_Richards.py

@@ -116,8 +116,8 @@ class RichardsTests1D(unittest.TestCase):
         v = np.random.rand(self.survey.nD)
         z = np.random.rand(self.M.nC)
         Hs = self.prob.fields(self.Ks)
-        vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
-        zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
+        vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
+        zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
         tol = TOL*(10**int(np.log10(np.abs(zJv))))
         passed = np.abs(vJz - zJv) < tol
         print 'Richards Adjoint Test - PressureHead'
@@ -188,8 +188,8 @@ class RichardsTests2D(unittest.TestCase):
         v = np.random.rand(self.survey.nD)
         z = np.random.rand(self.M.nC)
         Hs = self.prob.fields(self.Ks)
-        vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
-        zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
+        vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
+        zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
         tol = TOL*(10**int(np.log10(np.abs(zJv))))
         passed = np.abs(vJz - zJv) < tol
         print '2D: Richards Adjoint Test - PressureHead'
@@ -260,8 +260,8 @@ class RichardsTests3D(unittest.TestCase):
         v = np.random.rand(self.survey.nD)
         z = np.random.rand(self.M.nC)
         Hs = self.prob.fields(self.Ks)
-        vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
-        zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
+        vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
+        zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
         tol = TOL*(10**int(np.log10(np.abs(zJv))))
         passed = np.abs(vJz - zJv) < tol
         print '3D: Richards Adjoint Test - PressureHead'

tests/mesh/test_InterpolationMesh2Mesh.py

@@ -0,0 +1,353 @@
+import numpy as np
+import unittest
+from SimPEG.Utils import mkvc
+from SimPEG import Mesh, Tests
+import unittest
+
+test1D = True
+test2D = True
+test3D = False
+
+call1 = lambda fun, xyz: fun(xyz)
+call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, -1])
+call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
+cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
+cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)]
+cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy)))
+cartF2Cyl = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFz, fx, fy)))
+cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey)))
+cartE2Cyl = lambda M, ex, ey: cart_row2(M.gridEy, ex, ey)
+cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz)))
+cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez)))
+
+TOL = 1e-7
+
+if test1D:
+    class TestInterpolationMesh2Mesh_Tensor1D(Tests.OrderTest):
+
+        name = 'Mesh2Mesh Tensor1D'
+        meshSizes = [8, 16, 32]
+        meshTypes = ['uniformTensorMesh']
+        meshDimension = 1
+
+        def getError(self):
+            funX = lambda x: np.cos(2*np.pi*x)
+
+            mesh2, _ = self.makeMesh(self.M.nC-1, meshType=self._meshType, meshDimension=self.meshDimension )
+            ana = call1(funX, getattr(mesh2, 'grid%s'%self.type))
+
+            v = call1(funX, getattr(self.M, 'grid%s'%self.type))
+            P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type)
+            num = P*v
+
+            return np.linalg.norm((num - ana), np.inf)
+
+        def test_orderCC_1D(self):
+            self.type = 'CC'
+            self.name = 'Mesh2Mesh Tensor1D: CC'
+            self.orderTest()
+
+        def test_orderN_1D(self):
+            self.type = 'N'
+            self.name = 'Mesh2Mesh Tensor1D: N'
+            self.orderTest()
+
+        def test_orderEx_1D(self):
+            self.type = 'Ex'
+            self.name = 'Mesh2Mesh Tensor1D: Ex'
+            self.orderTest()
+
+        def test_orderFx_1D(self):
+            self.type = 'Fx'
+            self.name = 'Mesh2Mesh Tensor1D: Fx'
+            self.orderTest()
+
+if test2D:
+    class TestInterpolationMesh2Mesh_Tensor2D(Tests.OrderTest):
+
+        name = 'Mesh2Mesh Tensor2D'
+        meshSizes = [4, 8, 16]
+        meshTypes = ['uniformTensorMesh']
+        meshDimension = 2
+
+        def getError(self):
+            funX = lambda x, y: np.cos(2*np.pi*y)
+            funY = lambda x, y: np.cos(2*np.pi*x)
+
+            mesh2, _ = self.makeMesh(self.M.nC-1, meshType=self._meshType, meshDimension=self.meshDimension )
+
+            if 'x' in self.type:
+                ana = call2(funX, getattr(mesh2, 'grid%s'%self.type))
+            elif 'y' in self.type:
+                ana = call2(funY, getattr(mesh2, 'grid%s'%self.type))
+            elif 'F' in self.type:
+                ana = cartF2(mesh2, funX, funY)
+                ana = mesh2.projectFaceVector(ana)
+            elif 'E' in self.type:
+                ana = cartE2(mesh2, funX, funY)
+                ana = mesh2.projectEdgeVector(ana)
+            else:
+                ana = call2(funX, getattr(mesh2, 'grid%s'%self.type))
+
+
+            if 'F' in self.type:
+                v = cartF2(self.M, funX, funY)
+                if 'x' in self.type or 'y' in self.type:
+                    v = self.M.projectFaceVector(v)
+                else:
+                    v = mkvc(v)
+            elif 'E' in self.type:
+                v = cartE2(self.M, funX, funY)
+                if 'x' in self.type or 'y' in self.type:
+                    v = self.M.projectEdgeVector(v)
+                else:
+                    v = mkvc(v)
+            elif 'CC' == self.type:
+                v = call2(funX, self.M.gridCC)
+            elif 'N' == self.type:
+                v = call2(funX, self.M.gridN)
+
+            P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type)
+            # print P.shape, v.shape
+            num = P*v
+
+            return np.linalg.norm((num - ana), np.inf)
+
+        def test_orderCC_2D(self):
+            self.type = 'CC'
+            self.name = 'Mesh2Mesh Tensor2D: CC'
+            self.orderTest()
+
+        def test_orderN_2D(self):
+            self.type = 'N'
+            self.name = 'Mesh2Mesh Tensor2D: N'
+            self.orderTest()
+
+        def test_orderE_2D(self):
+            self.type = 'E'
+            self.name = 'Mesh2Mesh Tensor2D: E'
+            self.orderTest()
+
+        def test_orderEx_2D(self):
+            self.type = 'Ex'
+            self.name = 'Mesh2Mesh Tensor2D: Ex'
+            self.orderTest()
+
+        def test_orderEy_2D(self):
+            self.type = 'Ey'
+            self.name = 'Mesh2Mesh Tensor2D: Ey'
+            self.orderTest()
+
+        def test_orderF_2D(self):
+            self.type = 'F'
+            self.name = 'Mesh2Mesh Tensor2D: F'
+            self.orderTest()
+
+        def test_orderFx_2D(self):
+            self.type = 'Fx'
+            self.name = 'Mesh2Mesh Tensor2D: Fx'
+            self.orderTest()
+
+        def test_orderFy_2D(self):
+            self.type = 'Fy'
+            self.name = 'Mesh2Mesh Tensor2D: Fy'
+            self.orderTest()
+
+    class TestInterpolationMesh2Mesh_Cyl(Tests.OrderTest):
+
+        name = 'Mesh2Mesh Cyl'
+        meshSizes = [4, 8, 16]
+        meshTypes = ['uniformCylMesh']
+        meshDimension = 2
+
+        def getError(self):
+            funX = lambda x, y: np.cos(2*np.pi*y)
+            funY = lambda x, y: np.cos(2*np.pi*x)
+
+            mesh2, _ = self.makeMesh(self.M.nC-1, meshType=self._meshType, meshDimension=self.meshDimension )
+
+            if 'x' in self.type:
+                ana = call2(funX, getattr(mesh2, 'grid%s'%self.type))
+            elif 'y' in self.type:
+                ana = call2(funY, getattr(mesh2, 'grid%s'%self.type))
+            elif 'z' in self.type:
+                ana = call2(funY, getattr(mesh2, 'grid%s'%self.type))
+            elif 'F' in self.type:
+                ana = cartF2Cyl(mesh2, funX, funY)
+                ana = np.c_[ana[:,0], np.zeros_like(ana[:,0]), ana[:,1]]
+                ana = mesh2.projectFaceVector(ana)
+            elif 'E' in self.type:
+                ana = cartE2Cyl(mesh2, funX, funY)
+                ana = np.c_[np.zeros_like(ana[:,1]),ana[:,1],np.zeros_like(ana[:,1])]
+                ana = mesh2.projectEdgeVector(ana)
+            else:
+                ana = call2(funX, getattr(mesh2, 'grid%s'%self.type))
+
+
+            if 'F' in self.type:
+                v = cartF2Cyl(self.M, funX, funY)
+                v = np.c_[v[:,0], np.zeros_like(v[:,0]),v[:,1]]
+                if 'x' in self.type or 'z' in self.type:
+                    v = self.M.projectFaceVector(v)
+                else:
+                    v = np.c_[v[:,0], v[:,2]]
+                    v = mkvc(v)
+            elif 'E' in self.type:
+                v = cartE2Cyl(self.M, funX, funY)
+                v = np.c_[np.zeros_like(v[:,1]), v[:,1],np.zeros_like(v[:,1])]
+                v = self.M.projectEdgeVector(v)
+
+            elif 'CC' == self.type:
+                v = call2(funX, self.M.gridCC)
+            elif 'N' == self.type:
+                v = call2(funX, self.M.gridN)
+
+            P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type)
+            num = P*v
+
+            return np.linalg.norm((num - ana), np.inf)
+
+        def test_orderCC_Cyl(self):
+            self.type = 'CC'
+            self.name = 'Mesh2Mesh Tensor2D: CC'
+            self.orderTest()
+
+        def test_orderN_Cyl(self):
+            self.type = 'N'
+            self.name = 'Mesh2Mesh Tensor2D: N'
+            self.orderTest()
+
+        def test_orderE_Cyl(self):
+            self.type = 'E'
+            self.name = 'Mesh2Mesh Tensor2D: E'
+            self.orderTest()
+
+        def test_orderEy_Cyl(self):
+            self.type = 'Ey'
+            self.name = 'Mesh2Mesh Tensor2D: Ey'
+            self.orderTest()
+
+        def test_orderF_Cyl(self):
+            self.type = 'F'
+            self.name = 'Mesh2Mesh Tensor2D: F'
+            self.orderTest()
+
+        def test_orderFx_Cyl(self):
+            self.type = 'Fx'
+            self.name = 'Mesh2Mesh Tensor2D: Fx'
+            self.orderTest()
+
+        def test_orderFz_Cyl(self):
+            self.type = 'Fz'
+            self.name = 'Mesh2Mesh Tensor2D: Fz'
+            self.orderTest()
+
+if test3D:
+    class TestInterpolationMesh2Mesh_Tensor3D(Tests.OrderTest):
+
+        name = 'Mesh2Mesh Tensor3D'
+        meshSizes = [4, 8, 16]
+        meshTypes = ['uniformTensorMesh']
+        meshDimension = 3
+
+        def getError(self):
+            funX = lambda x, y, z: np.cos(2*np.pi*y)
+            funY = lambda x, y, z: np.cos(2*np.pi*z)
+            funZ = lambda x, y, z: np.cos(2*np.pi*x)
+
+            mesh2, _ = self.makeMesh(self.M.nC-1, meshType=self._meshType, meshDimension=self.meshDimension )
+
+            if 'x' in self.type:
+                ana = call3(funX, getattr(mesh2, 'grid%s'%self.type))
+            elif 'y' in self.type:
+                ana = call3(funY, getattr(mesh2, 'grid%s'%self.type))
+            elif 'z' in self.type:
+                ana = call3(funZ, getattr(mesh2, 'grid%s'%self.type))
+            elif 'F' in self.type:
+                ana = cartF3(mesh2, funX, funY, funZ)
+                ana = mesh2.projectFaceVector(ana)
+            elif 'E' in self.type:
+                ana = cartE3(mesh2, funX, funY, funZ)
+                ana = mesh2.projectFaceVector(ana)
+            else:
+                ana = call3(funX, getattr(mesh2, 'grid%s'%self.type))
+
+
+            if 'F' in self.type:
+                v = cartF3(self.M, funX, funY, funZ)
+                if 'x' in self.type or 'y' in self.type or 'z' in self.type:
+                    v = self.M.projectFaceVector(v)
+                else:
+                    v = mkvc(v)
+            elif 'E' in self.type:
+                v = cartE3(self.M, funX, funY, funZ)
+                if 'x' in self.type or 'y' in self.type or 'z' in self.type:
+                    v = self.M.projectFaceVector(v)
+                else:
+                    v = mkvc(v)
+            elif 'CC' == self.type:
+                v = call3(funX, self.M.gridCC)
+            elif 'N' == self.type:
+                v = call3(funX, self.M.gridN)
+
+            P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type)
+            # print P.shape, v.shape
+            num = P*v
+
+            return np.linalg.norm((num - ana), np.inf)
+
+        def test_orderCC_3D(self):
+            self.type = 'CC'
+            self.name = 'Mesh2Mesh Tensor3D: CC'
+            self.orderTest()
+
+        def test_orderN_3D(self):
+            self.type = 'N'
+            self.name = 'Mesh2Mesh Tensor3D: N'
+            self.orderTest()
+
+        def test_orderE_3D(self):
+            self.type = 'E'
+            self.name = 'Mesh2Mesh Tensor3D: E'
+            self.orderTest()
+
+        def test_orderEx_3D(self):
+            self.type = 'Ex'
+            self.name = 'Mesh2Mesh Tensor3D: Ex'
+            self.orderTest()
+
+        def test_orderEy_3D(self):
+            self.type = 'Ey'
+            self.name = 'Mesh2Mesh Tensor3D: Ey'
+            self.orderTest()
+
+        def test_orderEz_3D(self):
+            self.type = 'Ez'
+            self.name = 'Mesh2Mesh Tensor3D: Ez'
+            self.orderTest()
+
+        def test_orderF_3D(self):
+            self.type = 'F'
+            self.name = 'Mesh2Mesh Tensor3D: F'
+            self.orderTest()
+
+        def test_orderFx_3D(self):
+            self.type = 'Fx'
+            self.name = 'Mesh2Mesh Tensor3D: Fx'
+            self.orderTest()
+
+        def test_orderFy_3D(self):
+            self.type = 'Fy'
+            self.name = 'Mesh2Mesh Tensor3D: Fy'
+            self.orderTest()
+
+        def test_orderFz_3D(self):
+            self.type = 'Fz'
+            self.name = 'Mesh2Mesh Tensor3D: Fz'
+            self.orderTest()
+
+
+
+
+if __name__ == '__main__':
+    unittest.main()