From 19bcdbfbf5361f6d9cb48199a8b086349b230e0e Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 13 Nov 2015 10:39:08 -0800 Subject: [PATCH] EM.FDEM.SrcFDEM_XXX --> EM.FDEM.Src.XXX --- SimPEG/EM/FDEM/SrcFDEM.py | 348 ++++++++++++++++++ SimPEG/EM/FDEM/SurveyFDEM.py | 621 ++++++++++++++++---------------- SimPEG/EM/FDEM/__init__.py | 2 +- SimPEG/EM/Utils/testingUtils.py | 10 +- 4 files changed, 665 insertions(+), 316 deletions(-) create mode 100644 SimPEG/EM/FDEM/SrcFDEM.py diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py new file mode 100644 index 00000000..f7b95925 --- /dev/null +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -0,0 +1,348 @@ +from SimPEG import Survey, Problem, Utils, np, sp +# import SimPEG.EM as EM +from SimPEG.EM.Utils import * +from scipy.constants import mu_0 +# from SurveyFDEM import RxFDEM + + +class BaseSrcFDEM(Survey.BaseSrc): + freq = None + # rxPair = EM.FDEM.RxFDEM + integrate = True + + def eval(self, prob): + S_m = self.S_m(prob) + S_e = self.S_e(prob) + return S_m, S_e + + def evalDeriv(self, prob, v, adjoint=False): + return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint) + + def bPrimary(self, prob): + return None + + def hPrimary(self, prob): + return None + + def ePrimary(self, prob): + return None + + def jPrimary(self, prob): + return None + + def S_m(self, prob): + return None + + def S_e(self, prob): + return None + + def S_mDeriv(self, prob, v, adjoint = False): + return None + + def S_eDeriv(self, prob, v, adjoint = False): + return None + + +class RawVec_e(BaseSrcFDEM): + """ + RawVec electric source. It is defined by the user provided vector S_e + + :param numpy.array S_e: electric source term + :param float freq: frequency + :param rxList: receiver list + """ + + def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): + self._S_e = np.array(S_e,dtype=complex) + self._ePrimary = ePrimary + self._bPrimary = bPrimary + self._hPrimary = hPrimary + self._jPrimary = jPrimary + self.freq = float(freq) + BaseSrcFDEM.__init__(self, rxList) + + def S_e(self, prob): + return self._S_e + + def ePrimary(self, prob): + return self._ePrimary + + def bPrimary(self, prob): + return self._bPrimary + + def hPrimary(self, prob): + return self._hPrimary + + def jPrimary(self, prob): + return self._jPrimary + + +class RawVec_m(BaseSrcFDEM): + """ + RawVec magnetic source. It is defined by the user provided vector S_m + + :param numpy.array S_m: magnetic source term + :param float freq: frequency + :param rxList: receiver list + """ + + def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): + self._S_m = np.array(S_m,dtype=complex) + self.freq = float(freq) + self.integrate = integrate + self._ePrimary = np.array(ePrimary,dtype=complex) + self._bPrimary = np.array(bPrimary,dtype=complex) + self._hPrimary = np.array(hPrimary,dtype=complex) + self._jPrimary = np.array(jPrimary,dtype=complex) + + BaseSrcFDEM.__init__(self, rxList) + + def S_m(self, prob): + return self._S_m + + def ePrimary(self, prob): + return self._ePrimary + + def bPrimary(self, prob): + return self._bPrimary + + def hPrimary(self, prob): + return self._hPrimary + + def jPrimary(self, prob): + return self._jPrimary + + +class RawVec(BaseSrcFDEM): + """ + RawVec source. It is defined by the user provided vectors S_m, S_e + + :param numpy.array S_m: magnetic source term + :param numpy.array S_e: electric source term + :param float freq: frequency + :param rxList: receiver list + """ + 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) + self.freq = float(freq) + self.integrate = integrate + BaseSrcFDEM.__init__(self, rxList) + + def S_m(self, prob): + if prob._eqLocs is 'EF' and self.integrate is True: + return prob.Me * self._S_m + return self._S_m + + def S_e(self, prob): + if prob._eqLocs is 'FE' and self.integrate is True: + return prob.Me * self._S_e + return self._S_e + + +class MagDipole(BaseSrcFDEM): + + #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): + self.freq = float(freq) + self.loc = loc + self.orientation = orientation + self.moment = moment + self.mu = mu + self.integrate = False + BaseSrcFDEM.__init__(self, rxList) + + def bPrimary(self, prob): + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl.T + + + if prob.mesh._meshType is 'CYL': + if not prob.mesh.isSymmetric: + # TODO ? + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + + else: + srcfct = MagneticDipoleVectorPotential + ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + a = np.concatenate((ax, ay, az)) + + return C*a + + def hPrimary(self, prob): + b = self.bPrimary(prob) + return h_from_b(prob,b) + + def S_m(self, prob): + b_p = self.bPrimary(prob) + return -1j*omega(self.freq)*b_p + + def S_e(self, prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) + + +class MagDipole_Bfield(BaseSrcFDEM): + + #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that + #TODO: neither does moment + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): + self.freq = float(freq) + self.loc = loc + self.orientation = orientation + self.moment = moment + self.mu = mu + BaseSrcFDEM.__init__(self, rxList) + + def bPrimary(self, prob): + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl.T + + srcfct = MagneticDipoleFields + if prob.mesh._meshType is 'CYL': + if not prob.mesh.isSymmetric: + # TODO ? + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + b = np.concatenate((bx,bz)) + else: + bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + b = np.concatenate((bx,by,bz)) + + return b + + def hPrimary(self, prob): + b = self.bPrimary(prob) + return h_from_b(prob, b) + + def S_m(self, prob): + b = self.bPrimary(prob) + return -1j*omega(self.freq)*b + + def S_e(self, prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) + + +class CircularLoop(BaseSrcFDEM): + + #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that + def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0): + self.freq = float(freq) + self.orientation = orientation + self.radius = radius + self.mu = mu + self.loc = loc + self.integrate = False + BaseSrcFDEM.__init__(self, rxList) + + def bPrimary(self, prob): + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl.T + + if prob.mesh._meshType is 'CYL': + 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) + + else: + srcfct = MagneticDipoleVectorPotential + 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) + a = np.concatenate((ax, ay, az)) + + return C*a + + def hPrimary(self, prob): + b = self.bPrimary(prob) + return 1./self.mu*b + + def S_m(self, prob): + b = self.bPrimary(prob) + return -1j*omega(self.freq)*b + + def S_e(self, prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) + + diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 94cac712..1cdc104c 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -1,6 +1,7 @@ from SimPEG import Survey, Problem, Utils, np, sp from SimPEG.EM.Utils import * from scipy.constants import mu_0 +import SrcFDEM as Src #################################################### # Receivers @@ -90,345 +91,345 @@ class RxFDEM(Survey.BaseRx): # Sources #################################################### -class SrcFDEM(Survey.BaseSrc): - freq = None - rxPair = RxFDEM - integrate = True +# class SrcFDEM(Survey.BaseSrc): +# freq = None +# rxPair = RxFDEM +# integrate = True - def eval(self, prob): - S_m = self.S_m(prob) - S_e = self.S_e(prob) - return S_m, S_e +# def eval(self, prob): +# S_m = self.S_m(prob) +# S_e = self.S_e(prob) +# return S_m, S_e - def evalDeriv(self, prob, v, adjoint=False): - return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint) +# def evalDeriv(self, prob, v, adjoint=False): +# return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint) - def bPrimary(self, prob): - return None +# def bPrimary(self, prob): +# return None - def hPrimary(self, prob): - return None +# def hPrimary(self, prob): +# return None - def ePrimary(self, prob): - return None +# def ePrimary(self, prob): +# return None - def jPrimary(self, prob): - return None +# def jPrimary(self, prob): +# return None - def S_m(self, prob): - return None +# def S_m(self, prob): +# return None - def S_e(self, prob): - return None +# def S_e(self, prob): +# return None - def S_mDeriv(self, prob, v, adjoint = False): - return None +# def S_mDeriv(self, prob, v, adjoint = False): +# return None - def S_eDeriv(self, prob, v, adjoint = False): - return None +# def S_eDeriv(self, prob, v, adjoint = False): +# return None -class SrcFDEM_RawVec_e(SrcFDEM): - """ - RawVec electric source. It is defined by the user provided vector S_e +# class SrcFDEM_RawVec_e(SrcFDEM): +# """ +# RawVec electric source. It is defined by the user provided vector S_e - :param numpy.array S_e: electric source term - :param float freq: frequency - :param rxList: receiver list - """ +# :param numpy.array S_e: electric source term +# :param float freq: frequency +# :param rxList: receiver list +# """ - def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): - self._S_e = np.array(S_e,dtype=complex) - self._ePrimary = ePrimary - self._bPrimary = bPrimary - self._hPrimary = hPrimary - self._jPrimary = jPrimary - self.freq = float(freq) - SrcFDEM.__init__(self, rxList) +# def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): +# self._S_e = np.array(S_e,dtype=complex) +# self._ePrimary = ePrimary +# self._bPrimary = bPrimary +# self._hPrimary = hPrimary +# self._jPrimary = jPrimary +# self.freq = float(freq) +# SrcFDEM.__init__(self, rxList) - def S_e(self, prob): - return self._S_e +# def S_e(self, prob): +# return self._S_e - def ePrimary(self, prob): - return self._ePrimary +# def ePrimary(self, prob): +# return self._ePrimary - def bPrimary(self, prob): - return self._bPrimary +# def bPrimary(self, prob): +# return self._bPrimary - def hPrimary(self, prob): - return self._hPrimary +# def hPrimary(self, prob): +# return self._hPrimary - def jPrimary(self, prob): - return self._jPrimary +# def jPrimary(self, prob): +# return self._jPrimary -class SrcFDEM_RawVec_m(SrcFDEM): - """ - RawVec magnetic source. It is defined by the user provided vector S_m +# class SrcFDEM_RawVec_m(SrcFDEM): +# """ +# RawVec magnetic source. It is defined by the user provided vector S_m - :param numpy.array S_m: magnetic source term - :param float freq: frequency - :param rxList: receiver list - """ +# :param numpy.array S_m: magnetic source term +# :param float freq: frequency +# :param rxList: receiver list +# """ - def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): - self._S_m = np.array(S_m,dtype=complex) - self.freq = float(freq) - self.integrate = integrate - self._ePrimary = np.array(ePrimary,dtype=complex) - self._bPrimary = np.array(bPrimary,dtype=complex) - self._hPrimary = np.array(hPrimary,dtype=complex) - self._jPrimary = np.array(jPrimary,dtype=complex) +# def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): +# self._S_m = np.array(S_m,dtype=complex) +# self.freq = float(freq) +# self.integrate = integrate +# self._ePrimary = np.array(ePrimary,dtype=complex) +# self._bPrimary = np.array(bPrimary,dtype=complex) +# self._hPrimary = np.array(hPrimary,dtype=complex) +# self._jPrimary = np.array(jPrimary,dtype=complex) - SrcFDEM.__init__(self, rxList) +# SrcFDEM.__init__(self, rxList) - def S_m(self, prob): - return self._S_m +# def S_m(self, prob): +# return self._S_m - def ePrimary(self, prob): - return self._ePrimary - - def bPrimary(self, prob): - return self._bPrimary - - def hPrimary(self, prob): - return self._hPrimary +# def ePrimary(self, prob): +# return self._ePrimary + +# def bPrimary(self, prob): +# return self._bPrimary + +# def hPrimary(self, prob): +# return self._hPrimary - def jPrimary(self, prob): - return self._jPrimary +# def jPrimary(self, prob): +# return self._jPrimary -class SrcFDEM_RawVec(SrcFDEM): - """ - RawVec source. It is defined by the user provided vectors S_m, S_e +# class SrcFDEM_RawVec(SrcFDEM): +# """ +# RawVec source. It is defined by the user provided vectors S_m, S_e - :param numpy.array S_m: magnetic source term - :param numpy.array S_e: electric source term - :param float freq: frequency - :param rxList: receiver list - """ - 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) - self.freq = float(freq) - self.integrate = integrate - SrcFDEM.__init__(self, rxList) - - def S_m(self, prob): - if prob._eqLocs is 'EF' and self.integrate is True: - return prob.Me * self._S_m - return self._S_m - - def S_e(self, prob): - if prob._eqLocs is 'FE' and self.integrate is True: - return prob.Me * self._S_e - return self._S_e - - -class SrcFDEM_MagDipole(SrcFDEM): - - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): - self.freq = float(freq) - self.loc = loc - self.orientation = orientation - self.moment = moment - self.mu = mu - self.integrate = False - SrcFDEM.__init__(self, rxList) - - def bPrimary(self, prob): - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - gridX = prob.mesh.gridEx - gridY = prob.mesh.gridEy - gridZ = prob.mesh.gridEz - C = prob.mesh.edgeCurl - - elif eqLocs is 'EF': - gridX = prob.mesh.gridFx - gridY = prob.mesh.gridFy - gridZ = prob.mesh.gridFz - C = prob.mesh.edgeCurl.T - - - if prob.mesh._meshType is 'CYL': - if not prob.mesh.isSymmetric: - # TODO ? - raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) - - else: - srcfct = MagneticDipoleVectorPotential - ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) - ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) - az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) - a = np.concatenate((ax, ay, az)) - - return C*a - - def hPrimary(self, prob): - b = self.bPrimary(prob) - return h_from_b(prob,b) - - def S_m(self, prob): - b_p = self.bPrimary(prob) - return -1j*omega(self.freq)*b_p - - def S_e(self, prob): - if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): - return None - else: - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - mui_s = prob.curModel.mui - 1./self.mu - MMui_s = prob.mesh.getFaceInnerProduct(mui_s) - C = prob.mesh.edgeCurl - elif eqLocs is 'EF': - mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) - C = prob.mesh.edgeCurl.T - - return -C.T * (MMui_s * self.bPrimary(prob)) - - -class SrcFDEM_MagDipole_Bfield(SrcFDEM): - - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - #TODO: neither does moment - def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): - self.freq = float(freq) - self.loc = loc - self.orientation = orientation - self.moment = moment - self.mu = mu - SrcFDEM.__init__(self, rxList) - - def bPrimary(self, prob): - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - gridX = prob.mesh.gridFx - gridY = prob.mesh.gridFy - gridZ = prob.mesh.gridFz - C = prob.mesh.edgeCurl - - elif eqLocs is 'EF': - gridX = prob.mesh.gridEx - gridY = prob.mesh.gridEy - gridZ = prob.mesh.gridEz - C = prob.mesh.edgeCurl.T - - srcfct = MagneticDipoleFields - if prob.mesh._meshType is 'CYL': - if not prob.mesh.isSymmetric: - # TODO ? - raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) - bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) - b = np.concatenate((bx,bz)) - else: - bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) - by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) - bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) - b = np.concatenate((bx,by,bz)) - - return b - - def hPrimary(self, prob): - b = self.bPrimary(prob) - return h_from_b(prob, b) - - def S_m(self, prob): - b = self.bPrimary(prob) - return -1j*omega(self.freq)*b - - def S_e(self, prob): - if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): - return None - else: - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - mui_s = prob.curModel.mui - 1./self.mu - MMui_s = prob.mesh.getFaceInnerProduct(mui_s) - C = prob.mesh.edgeCurl - elif eqLocs is 'EF': - mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) - C = prob.mesh.edgeCurl.T - - return -C.T * (MMui_s * self.bPrimary(prob)) - - -class SrcFDEM_CircularLoop(SrcFDEM): - - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0): - self.freq = float(freq) - self.orientation = orientation - self.radius = radius - self.mu = mu - self.loc = loc - self.integrate = False - SrcFDEM.__init__(self, rxList) - - def bPrimary(self, prob): - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - gridX = prob.mesh.gridEx - gridY = prob.mesh.gridEy - gridZ = prob.mesh.gridEz - C = prob.mesh.edgeCurl - - elif eqLocs is 'EF': - gridX = prob.mesh.gridFx - gridY = prob.mesh.gridFy - gridZ = prob.mesh.gridFz - C = prob.mesh.edgeCurl.T - - if prob.mesh._meshType is 'CYL': - 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) - - else: - srcfct = MagneticDipoleVectorPotential - 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) - a = np.concatenate((ax, ay, az)) - - return C*a - - def hPrimary(self, prob): - b = self.bPrimary(prob) - return 1./self.mu*b - - def S_m(self, prob): - b = self.bPrimary(prob) - return -1j*omega(self.freq)*b - - def S_e(self, prob): - if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): - return None - else: - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - mui_s = prob.curModel.mui - 1./self.mu - MMui_s = prob.mesh.getFaceInnerProduct(mui_s) - C = prob.mesh.edgeCurl - elif eqLocs is 'EF': - mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) - C = prob.mesh.edgeCurl.T - - return -C.T * (MMui_s * self.bPrimary(prob)) +# :param numpy.array S_m: magnetic source term +# :param numpy.array S_e: electric source term +# :param float freq: frequency +# :param rxList: receiver list +# """ +# 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) +# self.freq = float(freq) +# self.integrate = integrate +# SrcFDEM.__init__(self, rxList) + +# def S_m(self, prob): +# if prob._eqLocs is 'EF' and self.integrate is True: +# return prob.Me * self._S_m +# return self._S_m + +# def S_e(self, prob): +# if prob._eqLocs is 'FE' and self.integrate is True: +# return prob.Me * self._S_e +# return self._S_e + + +# class SrcFDEM_MagDipole(SrcFDEM): + +# #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that +# def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): +# self.freq = float(freq) +# self.loc = loc +# self.orientation = orientation +# self.moment = moment +# self.mu = mu +# self.integrate = False +# SrcFDEM.__init__(self, rxList) + +# def bPrimary(self, prob): +# eqLocs = prob._eqLocs + +# if eqLocs is 'FE': +# gridX = prob.mesh.gridEx +# gridY = prob.mesh.gridEy +# gridZ = prob.mesh.gridEz +# C = prob.mesh.edgeCurl + +# elif eqLocs is 'EF': +# gridX = prob.mesh.gridFx +# gridY = prob.mesh.gridFy +# gridZ = prob.mesh.gridFz +# C = prob.mesh.edgeCurl.T + + +# if prob.mesh._meshType is 'CYL': +# if not prob.mesh.isSymmetric: +# # TODO ? +# raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') +# a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + +# else: +# srcfct = MagneticDipoleVectorPotential +# ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) +# ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) +# az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) +# a = np.concatenate((ax, ay, az)) + +# return C*a + +# def hPrimary(self, prob): +# b = self.bPrimary(prob) +# return h_from_b(prob,b) + +# def S_m(self, prob): +# b_p = self.bPrimary(prob) +# return -1j*omega(self.freq)*b_p + +# def S_e(self, prob): +# if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): +# return None +# else: +# eqLocs = prob._eqLocs + +# if eqLocs is 'FE': +# mui_s = prob.curModel.mui - 1./self.mu +# MMui_s = prob.mesh.getFaceInnerProduct(mui_s) +# C = prob.mesh.edgeCurl +# elif eqLocs is 'EF': +# mu_s = prob.curModel.mu - self.mu +# MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) +# C = prob.mesh.edgeCurl.T + +# return -C.T * (MMui_s * self.bPrimary(prob)) + + +# class SrcFDEM_MagDipole_Bfield(SrcFDEM): + +# #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that +# #TODO: neither does moment +# def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): +# self.freq = float(freq) +# self.loc = loc +# self.orientation = orientation +# self.moment = moment +# self.mu = mu +# SrcFDEM.__init__(self, rxList) + +# def bPrimary(self, prob): +# eqLocs = prob._eqLocs + +# if eqLocs is 'FE': +# gridX = prob.mesh.gridFx +# gridY = prob.mesh.gridFy +# gridZ = prob.mesh.gridFz +# C = prob.mesh.edgeCurl + +# elif eqLocs is 'EF': +# gridX = prob.mesh.gridEx +# gridY = prob.mesh.gridEy +# gridZ = prob.mesh.gridEz +# C = prob.mesh.edgeCurl.T + +# srcfct = MagneticDipoleFields +# if prob.mesh._meshType is 'CYL': +# if not prob.mesh.isSymmetric: +# # TODO ? +# raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') +# bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) +# bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) +# b = np.concatenate((bx,bz)) +# else: +# bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) +# by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) +# bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) +# b = np.concatenate((bx,by,bz)) + +# return b + +# def hPrimary(self, prob): +# b = self.bPrimary(prob) +# return h_from_b(prob, b) + +# def S_m(self, prob): +# b = self.bPrimary(prob) +# return -1j*omega(self.freq)*b + +# def S_e(self, prob): +# if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): +# return None +# else: +# eqLocs = prob._eqLocs + +# if eqLocs is 'FE': +# mui_s = prob.curModel.mui - 1./self.mu +# MMui_s = prob.mesh.getFaceInnerProduct(mui_s) +# C = prob.mesh.edgeCurl +# elif eqLocs is 'EF': +# mu_s = prob.curModel.mu - self.mu +# MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) +# C = prob.mesh.edgeCurl.T + +# return -C.T * (MMui_s * self.bPrimary(prob)) + + +# class SrcFDEM_CircularLoop(SrcFDEM): + +# #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that +# def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0): +# self.freq = float(freq) +# self.orientation = orientation +# self.radius = radius +# self.mu = mu +# self.loc = loc +# self.integrate = False +# SrcFDEM.__init__(self, rxList) + +# def bPrimary(self, prob): +# eqLocs = prob._eqLocs + +# if eqLocs is 'FE': +# gridX = prob.mesh.gridEx +# gridY = prob.mesh.gridEy +# gridZ = prob.mesh.gridEz +# C = prob.mesh.edgeCurl + +# elif eqLocs is 'EF': +# gridX = prob.mesh.gridFx +# gridY = prob.mesh.gridFy +# gridZ = prob.mesh.gridFz +# C = prob.mesh.edgeCurl.T + +# if prob.mesh._meshType is 'CYL': +# 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) + +# else: +# srcfct = MagneticDipoleVectorPotential +# 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) +# a = np.concatenate((ax, ay, az)) + +# return C*a + +# def hPrimary(self, prob): +# b = self.bPrimary(prob) +# return 1./self.mu*b + +# def S_m(self, prob): +# b = self.bPrimary(prob) +# return -1j*omega(self.freq)*b + +# def S_e(self, prob): +# if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): +# return None +# else: +# eqLocs = prob._eqLocs + +# if eqLocs is 'FE': +# mui_s = prob.curModel.mui - 1./self.mu +# MMui_s = prob.mesh.getFaceInnerProduct(mui_s) +# C = prob.mesh.edgeCurl +# elif eqLocs is 'EF': +# mu_s = prob.curModel.mu - self.mu +# MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) +# C = prob.mesh.edgeCurl.T + +# return -C.T * (MMui_s * self.bPrimary(prob)) #################################################### @@ -440,7 +441,7 @@ class SurveyFDEM(Survey.BaseSurvey): docstring for SurveyFDEM """ - srcPair = SrcFDEM + srcPair = Src.BaseSrcFDEM def __init__(self, srcList, **kwargs): # Sort these by frequency diff --git a/SimPEG/EM/FDEM/__init__.py b/SimPEG/EM/FDEM/__init__.py index 110b4d1e..9a262ffb 100644 --- a/SimPEG/EM/FDEM/__init__.py +++ b/SimPEG/EM/FDEM/__init__.py @@ -1,3 +1,3 @@ -from SurveyFDEM import * +from SurveyFDEM import RxFDEM, Src, SurveyFDEM from FDEM import BaseFDEMProblem, ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h from FieldsFDEM import * \ No newline at end of file diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index b4570c0d..cb1e1d26 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -23,25 +23,25 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): for SrcType in SrcList: if SrcType is 'MagDipole': - Src.append(EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.])) + Src.append(EM.FDEM.Src.MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.])) elif SrcType is 'MagDipole_Bfield': - Src.append(EM.FDEM.SrcFDEM_MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.])) + Src.append(EM.FDEM.Src.MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.])) elif SrcType is 'CircularLoop': - Src.append(EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.])) + 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])] = 1. S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - Src.append(EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)) + Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) 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])] = 1. S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - Src.append(EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)) + Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) if verbose: print ' Fetching %s problem' % (fdemType)