mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-06 05:16:51 +08:00
EM.FDEM.SrcFDEM_XXX --> EM.FDEM.Src.XXX
This commit is contained in:
Lindsey Heagy
@@ -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))
|
||||
|
||||
|
||||
+311
-310
@@ -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
|
||||
|
||||
@@ -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 *
|
||||
@@ -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)
|
||||
|
||||
Reference in New Issue
Block a user