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

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

SimPEG/EM/Base.py

@@ -7,7 +7,7 @@ class EMPropMap(Maps.PropMap):
     """
 
     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))
     mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap))

SimPEG/EM/FDEM/FDEM.py

@@ -8,29 +8,29 @@ from SimPEG.EM.Utils import omega
 
 class BaseFDEMProblem(BaseEMProblem):
     """
-        We start by looking at Maxwell's equations in the electric
-        field \\\(\\\mathbf{e}\\\) and the magnetic flux
-        density \\\(\\\mathbf{b}\\\)
+    We start by looking at Maxwell's equations in the electric
+    field \\\(\\\mathbf{e}\\\) and the magnetic flux
+    density \\\(\\\mathbf{b}\\\)
 
-        .. math ::
+    .. math ::
 
-            \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
-            {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
+        \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
+        {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
 
-        if using the E-B formulation (:code:`Problem_e`
-        or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
+    if using the E-B formulation (:code:`Problem_e`
+    or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
 
-        If we write Maxwell's equations in terms of
-        \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
+    If we write Maxwell's equations in terms of
+    \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
 
-        .. math ::
+    .. math ::
 
-            \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
-            \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
+        \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
+        \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
 
-        if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
+    if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
 
-        The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
+    The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
     """
 
     surveyPair = SurveyFDEM
@@ -204,6 +204,17 @@ class Problem_e(BaseFDEMProblem):
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
 
+    def _GLoc(self, fieldType):
+        if fieldType == 'e':
+            return 'E'
+        elif fieldType == 'b':
+            return 'F'
+        elif (fieldType == 'h') or (fieldType == 'j'):
+            return 'CCV'
+        else:
+            raise Exception('Field type must be e, b, h, j')
+
+
     def getA(self, freq):
         """
         System matrix
@@ -283,7 +294,6 @@ class Problem_e(BaseFDEMProblem):
         if adjoint:
             dRHS = MfMui * (C * 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)
 
@@ -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
@@ -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
@@ -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

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -193,16 +193,6 @@ class Fields_e(Fields):
         self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv
         self._MfMui = self.survey.prob.MfMui
 
-    def _GLoc(self, fieldType):
-        if fieldType == 'e':
-            return 'E'
-        elif fieldType == 'b':
-            return 'F'
-        elif (fieldType == 'h') or (fieldType == 'j'):
-            return 'CCV'
-        else:
-            raise Exception('Field type must be e, b, h, j')
-
 
     def _ePrimary(self, eSolution, srcList):
         """
@@ -465,17 +455,6 @@ class Fields_b(Fields):
         self._nC = self.survey.prob.mesh.nC
 
 
-
-    def _GLoc(self,fieldType):
-        if fieldType == 'e':
-            return 'E'
-        elif fieldType == 'b':
-            return 'F'
-        elif (fieldType == 'h') or (fieldType == 'j'):
-            return'CCV'
-        else:
-            raise Exception('Field type must be e, b, h, j')
-
     def _bPrimary(self, bSolution, srcList):
         """
         Primary magnetic flux density from source
@@ -729,16 +708,6 @@ class Fields_j(Fields):
         self._aveE2CCV = self.survey.prob.mesh.aveE2CCV
         self._nC = self.survey.prob.mesh.nC
 
-    def _GLoc(self,fieldType):
-        if fieldType == 'h':
-            return 'E'
-        elif fieldType == 'j':
-            return 'F'
-        elif (fieldType == 'e') or (fieldType == 'b'):
-            return 'CCV'
-        else:
-            raise Exception('Field type must be e, b, h, j')
-
     def _jPrimary(self, jSolution, srcList):
         """
         Primary current density from source
@@ -1024,16 +993,6 @@ class Fields_h(Fields):
         self._aveE2CCV = self.survey.prob.mesh.aveE2CCV
         self._nC = self.survey.prob.mesh.nC
 
-    def _GLoc(self,fieldType):
-        if fieldType == 'h':
-            return 'E'
-        elif fieldType == 'j':
-            return 'F'
-        elif (fieldType == 'e') or (fieldType == 'b'):
-            return 'CCV'
-        else:
-            raise Exception('Field type must be e, b, h, j')
-
     def _hPrimary(self, hSolution, srcList):
         """
         Primary magnetic field from source

SimPEG/EM/FDEM/SrcFDEM.py

@@ -9,8 +9,10 @@ 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):
         """
@@ -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
 
@@ -136,13 +138,12 @@ class RawVec_e(BaseSrc):
     :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 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):
+    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)
 
@@ -166,13 +167,12 @@ class RawVec_m(BaseSrc):
     :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 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)
         self.freq = float(freq)
-        self.integrate = integrate
 
         BaseSrc.__init__(self, rxList)
 
@@ -197,14 +197,13 @@ class RawVec(BaseSrc):
     :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 bool integrate: Integrate the source term (multiply by Me) [False]
     """
-    def __init__(self, rxList, freq, s_m, s_e, integrate=True):
+    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):
         """
@@ -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):
@@ -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)
@@ -605,5 +603,231 @@ 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

@@ -65,7 +65,7 @@ class Rx(SimPEG.Survey.BaseRx):
 
     def projGLoc(self, u):
         """Grid Location projection (e.g. Ex Fy ...)"""
-        return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
+        return u.prob._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
 
     def eval(self, src, mesh, f):
         """
@@ -77,8 +77,6 @@ class Rx(SimPEG.Survey.BaseRx):
         :rtype: numpy.ndarray
         :return: fields projected to recievers
         """
-        # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
-        # projGLoc += self.knownRxTypes[self.rxType][1]
 
         P = self.getP(mesh, self.projGLoc(f))
         f_part_complex = f[src, self.projField]

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/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/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/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()