From e83b76a70c2c4d3546ec6f4d0c171e06c52f9055 Mon Sep 17 00:00:00 2001 From: seogi Date: Fri, 21 Nov 2014 09:37:37 -0800 Subject: [PATCH] Modification for Cylinderical mesh --- simpegEM/FDEM/FDEM.py | 70 ++++++++++++++++++++++++------ simpegEM/FDEM/SurveyFDEM.py | 5 ++- simpegEM/Sources/__init__.py | 1 + simpegEM/Sources/magneticDipole.py | 46 ++++++++++++++++++++ simpegEM/Utils/__init__.py | 3 ++ 5 files changed, 111 insertions(+), 14 deletions(-) create mode 100644 simpegEM/Utils/__init__.py diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 9e31002a..c4137486 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -162,9 +162,18 @@ class ProblemFDEM_e(BaseFDEMProblem): for i, tx in enumerate(Txs): if tx.txType == 'VMD': src = Sources.MagneticDipoleVectorPotential + SRCx = src(tx.loc, self.mesh.gridEx, 'x') + SRCy = src(tx.loc, self.mesh.gridEy, 'y') + SRCz = src(tx.loc, self.mesh.gridEz, 'z') + + elif tx.txType == 'CircularLoop': + src = Sources.MagneticLoopVectorPotential + SRCx = src(tx.loc, self.mesh.gridEx, 'x', tx.radius) + SRCy = src(tx.loc, self.mesh.gridEy, 'y', tx.radius) + SRCz = src(tx.loc, self.mesh.gridEz, 'z', tx.radius) else: raise NotImplemented('%s txType is not implemented' % tx.txType) - rhs[i] = src(tx.loc, self.mesh, ['Ex','Ey','Ez']) + rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') mui = self.MfMui @@ -241,20 +250,55 @@ class ProblemFDEM_b(BaseFDEMProblem): Txs = self.survey.getTransmitters(freq) rhs = range(len(Txs)) for i, tx in enumerate(Txs): - if tx.txType == 'VMD': - src = Sources.MagneticDipoleVectorPotential + + if self.mesh._meshType is 'CYL': + if self.mesh.isSymmetric: + if tx.txType == 'VMD': + SRC = Sources.MagneticDipoleVectorPotential(tx.loc, self.mesh.gridEy, 'y') + elif tx.txType =='CircularLoop': + SRC = Sources.MagneticLoopVectorPotential(tx.loc, self.mesh.gridEy, 'y', tx.radius) + else: + raise NotImplementedError('Only VMD and CircularLoop') + else: + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + + elif self.mesh._meshType is 'TENSOR': + + if tx.txType == 'VMD': + src = Sources.MagneticDipoleVectorPotential + SRCx = src(tx.loc, self.mesh.gridEx, 'x') + SRCy = src(tx.loc, self.mesh.gridEy, 'y') + SRCz = src(tx.loc, self.mesh.gridEz, 'z') + + elif tx.txType == 'VMD_B': + src = Sources.MagneticDipoleFields + SRCx = src(tx.loc, self.mesh.gridFx, 'x') + SRCy = src(tx.loc, self.mesh.gridFy, 'y') + SRCz = src(tx.loc, self.mesh.gridFz, 'z') + + elif tx.txType == 'CircularLoop': + src = Sources.MagneticLoopVectorPotential + SRCx = src(tx.loc, self.mesh.gridEx, 'x', tx.radius) + SRCy = src(tx.loc, self.mesh.gridEy, 'y', tx.radius) + SRCz = src(tx.loc, self.mesh.gridEz, 'z', tx.radius) + else: + + raise NotImplemented('%s txType is not implemented' % tx.txType) + SRC = np.concatenate((SRCx, SRCy, SRCz)) + 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)) + raise Exception('Unknown mesh for VMD') + + rhs[i] = SRC + + mui = self.MfMui + if tx.txType == 'VMD_B': + b_0 = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') + else: + a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') + C = self.mesh.edgeCurl + b_0 = C*a - a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') - mui = self.MfMui - C = self.mesh.edgeCurl - - b_0 = C*a return -1j*omega(freq)*mui*b_0 def calcFields(self, sol, freq, fieldType, adjoint=False): diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index dce06ba5..9e87465f 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -17,6 +17,7 @@ class RxFDEM(Survey.BaseRx): 'byi':['b', 'Fy', 'imag'], 'bzi':['b', 'Fz', 'imag'], } + radius = None def __init__(self, locs, rxType): Survey.BaseRx.__init__(self, locs, rxType) @@ -71,7 +72,9 @@ class TxFDEM(Survey.BaseTx): rxPair = RxFDEM - knownTxTypes = ['VMD'] + knownTxTypes = ['VMD', 'VMD_B', 'CircularLoop'] + + radius = None def __init__(self, loc, txType, freq, rxList): self.freq = float(freq) diff --git a/simpegEM/Sources/__init__.py b/simpegEM/Sources/__init__.py index 6e2fedc6..1f04379e 100644 --- a/simpegEM/Sources/__init__.py +++ b/simpegEM/Sources/__init__.py @@ -1,2 +1,3 @@ from magneticDipole import MagneticDipoleVectorPotential from CircularLoop import MagneticLoopVectorPotential +from magneticDipole import MagneticDipoleFields diff --git a/simpegEM/Sources/magneticDipole.py b/simpegEM/Sources/magneticDipole.py index 9f6166c8..3a102e74 100644 --- a/simpegEM/Sources/magneticDipole.py +++ b/simpegEM/Sources/magneticDipole.py @@ -52,3 +52,49 @@ def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0. if nTx == 1: return A.flatten() return A + +def MagneticDipoleFields(txLoc, obsLoc, component, dipoleMoment=1.): + """ + Calculate the vector potential of a set of magnetic dipoles + at given locations 'ref. ' + + :param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z) + :param numpy.ndarray obsLoc: Where the potentials will be calculated (x, y, z) + :param str component: The component to calculate - 'x', 'y', or 'z' + :param numpy.ndarray dipoleMoment: The vector dipole moment (vertical) + :rtype: numpy.ndarray + :return: The vector potential each dipole at each observation location + """ + + if component=='x': + dimInd = 0 + elif component=='y': + dimInd = 1 + elif component=='z': + dimInd = 2 + else: + raise ValueError('Invalid component') + + txLoc = np.atleast_2d(txLoc) + obsLoc = np.atleast_2d(obsLoc) + dipoleMoment = np.atleast_2d(dipoleMoment) + + nFaces = obsLoc.shape[0] + nTx = txLoc.shape[0] + + m = np.array(dipoleMoment).repeat(nFaces, axis=0) + B = np.empty((nFaces, nTx)) + for i in range(nTx): + dR = obsLoc - txLoc[i, np.newaxis].repeat(nFaces, axis=0) + r = np.sqrt((dR**2).sum(axis=1)) + if dimInd == 0: + B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2) + elif dimInd == 1: + B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2) + elif dimInd == 2: + B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1) + else: + raise Exception("Not Implemented") + if nTx == 1: + return B.flatten() + return B diff --git a/simpegEM/Utils/__init__.py b/simpegEM/Utils/__init__.py new file mode 100644 index 00000000..21404ecc --- /dev/null +++ b/simpegEM/Utils/__init__.py @@ -0,0 +1,3 @@ +import Sources +import Ana +import Solver \ No newline at end of file