fake the sources for FDEM

simpegEM/FDEM/FDEM.py

@@ -1,6 +1,7 @@
 from SimPEG import Problem, Solver, Utils, np, sp
 from scipy.constants import mu_0
 from SurveyFDEM import SurveyFDEM, DataFDEM, FieldsFDEM
+from simpegEM.Utils import Sources
 
 def omega(freq):
     """Change frequency to angular frequency, omega"""
@@ -60,16 +61,33 @@ class ProblemFDEM_e(Problem.BaseProblem):
 
     def getA(self, freq):
         """
-            :param int fInd: Frequency index
+            :param float freq: Frequency
             :rtype: scipy.sparse.csr_matrix
             :return: A
         """
         return self.mesh.edgeCurl.T*self.MfMui*self.mesh.edgeCurl + 1j*omega(freq)*self.MeSigma
 
     def getRHS(self, freq):
-        #TODO: this needs to also depend on your transmitter!
+        """
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nE, nTx)
+            :return: RHS
+        """
+        Txs = self.survey.getTransmitters(freq)
+        rhs = range(len(Txs))
+        for i, tx in enumerate(Txs):
+            if tx.txType == 'VMD':
+                src = Sources.MagneticDipoleVectorPotential
+            else:
+                raise NotImplemented('%s txType is not implemented' % tx.txType)
+            SRCx = src(tx.loc, self.mesh.gridEx, 'x')
+            SRCy = src(tx.loc, self.mesh.gridEy, 'y')
+            SRCz = src(tx.loc, self.mesh.gridEz, 'z')
+            rhs[i] = np.concatenate((SRCx, SRCy, SRCz))
 
-        return -1j*omega(freq)*self.Me*self.j_s
+        #TODO: this is completely wrong. b0 = self.mesh.edgeCurl*rhs but we are doing an e formulation...
+        j_s = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F')
+        return -1j*omega(freq)*self.Me*j_s
 
 
     def fields(self, m, useThisRhs=None):

simpegEM/Tests/test_FDEM.py

@@ -9,45 +9,31 @@ class FDEM_bDerivTests(unittest.TestCase):
     def setUp(self):
 
         cs = 5.
-        ncx = 2
-        ncy = 2
-        ncz = 2
+        ncx, ncy, ncz = 2, 2, 2
         npad = 3
         hx = Utils.meshTensors(((npad,cs), (ncx,cs), (npad,cs)))
         hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
         hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs)))
         mesh = Mesh.TensorMesh([hx,hy,hz])
 
-        XY = Utils.ndgrid(np.linspace(20,50,3), np.linspace(20,50,3))
-        rxLoc = np.c_[XY, np.ones(XY.shape[0])*40]
-
         model = Model.LogModel(mesh)
 
         x = np.linspace(5,10,3)
         XYZ = Utils.ndgrid(x,x,np.r_[0])
         rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex')
-        Tx0 = EM.FDEM.TxFDEM(None, 'VMD', 3, rxList)
+        Tx0 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-2, rxList)
 
         x = np.linspace(5,10,3)
         XYZ = Utils.ndgrid(x,x,np.r_[0])
         rxList = EM.FDEM.RxListFDEM(XYZ, 'Ey')
-        Tx1 = EM.FDEM.TxFDEM(None, 'VMD', 3, rxList)
+        Tx1 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-4, rxList)
 
         survey = EM.FDEM.SurveyFDEM([Tx0, Tx1])
 
         prb = EM.FDEM.ProblemFDEM_e(model)
         prb.pair(survey)
 
-        sigma = np.log(np.ones(mesh.nC)*1e-3)
-
-        j_sx = np.zeros(mesh.vnEx)
-        j_sx[4,4,4] = 1
-        j_s = np.r_[Utils.mkvc(j_sx),np.zeros(mesh.nEy+mesh.nEz)]
-
-        prb.j_s = np.c_[j_s, j_s]
-        f = prb.fields(sigma)
-
-        self.sigma = sigma
+        self.sigma = np.log(np.ones(mesh.nC)*1e-3)
         self.prb = prb
         self.survey = survey