mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-02 01:55:51 +08:00
Added code from KrisDavis for sources written for SimPEG.
This commit is contained in:
rowanc1
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -0,0 +1,219 @@
|
||||
import numpy as np
|
||||
|
||||
class Sources(object):
|
||||
"""
|
||||
Class creates base sources
|
||||
"""
|
||||
def __init__(self):
|
||||
raise Exception('Sources is a base class that requires a mesh. Inherit to your favorite Mesh class.')
|
||||
|
||||
|
||||
def magneticDipoleSource(self,dipoleLoc,dimDip=3):
|
||||
"""
|
||||
Calculates the response from magnetic dipole source(s).
|
||||
|
||||
Inputs:
|
||||
mesh - mesh object
|
||||
dipoleLoc - an array of [n x 3] dipole locations (xyz)
|
||||
dimDip - dipole dimension (e.g., x = 1,y = 2,z = 3)
|
||||
|
||||
Outputs:
|
||||
b - the magnetic field on the face grid
|
||||
"""
|
||||
|
||||
def getMagneticDipolePotential(loc,m,grid,dimAxis):
|
||||
# Get the potential of the dipole
|
||||
# Returns A(# faces x # transmitters)
|
||||
|
||||
nL = grid.shape[0]
|
||||
nT = loc.shape[0]
|
||||
fullLoc = np.ones((nL,3))
|
||||
A = np.zeros((nL,nT))
|
||||
for ii in range(nT):
|
||||
fullM = np.ones((nL,3))*m[ii,:]
|
||||
br = grid - fullLoc*loc[ii,:]
|
||||
cp = np.cross(fullM,br)
|
||||
r = np.sqrt(br[:,0]**2 + br[:,1]**2 + br[:,2]**2)
|
||||
A[:,ii] = ((1e-6)*cp[:,dimAxis-1])/r**3
|
||||
return A
|
||||
|
||||
dipoleLoc = np.atleast_2d(dipoleLoc)
|
||||
dimDip = np.atleast_2d(dimDip)
|
||||
m = np.zeros((dipoleLoc.shape[0],3))
|
||||
if (dimDip.shape[0] == 1):
|
||||
try:
|
||||
m[:,dimDip-1] = 1
|
||||
except IndexError:
|
||||
print "magneticSource:Error: Dipole dimension should be 1 (x), 2 (y), or 3 (z)."
|
||||
raise
|
||||
elif (dimDip.shape[0] == dipoleLoc.shape[0]):
|
||||
try:
|
||||
for jj in range(dipoleLoc.shape[0]):
|
||||
m[jj,dimDip[jj] - 1] = 1
|
||||
|
||||
except IndexError:
|
||||
print "magneticSource:Error: Dipole dimension should be 1 (x), 2 (y), or 3 (z)."
|
||||
raise
|
||||
else:
|
||||
print "magneticSource:Error: Dipole direction should also be vector of same length as dipole locations."
|
||||
raise
|
||||
|
||||
# Get magnetic potential at each set of orthogonal faces:
|
||||
Ax = getMagneticDipolePotential(dipoleLoc,m,self.gridEx,1)
|
||||
Ay = getMagneticDipolePotential(dipoleLoc,m,self.gridEy,2)
|
||||
Az = getMagneticDipolePotential(dipoleLoc,m,self.gridEz,3)
|
||||
|
||||
# Combine potential
|
||||
A = np.concatenate((Ax, Ay, Az),axis=0)
|
||||
|
||||
# B = curl A
|
||||
CURL = self.edgeCurl
|
||||
self._src = CURL*A
|
||||
|
||||
return self._src
|
||||
|
||||
def simpleLoopSource(self,loc,Lx,Ly=None):
|
||||
"""
|
||||
Returns unit values for a simple loop source(s)
|
||||
|
||||
Inputs:
|
||||
mesh - mesh object
|
||||
loc - an array of [n x 3] loop centre location (xyz)
|
||||
Lx - length of loop in x-direction (meters)
|
||||
Ly - length of loop in y-direction (default is Lx)
|
||||
|
||||
Outputs:
|
||||
e - Unit values of current on the edge grid
|
||||
"""
|
||||
|
||||
# Check for default value of Ly:
|
||||
Lx = np.atleast_2d(Lx)
|
||||
if Ly is None:
|
||||
Ly = Lx
|
||||
else:
|
||||
Ly = np.atleast_2d(Ly)
|
||||
|
||||
# Number of loops
|
||||
loc = np.atleast_2d(loc)
|
||||
nL = loc.shape[0]
|
||||
normLoc = np.zeros((nL,3))
|
||||
sub = np.zeros((nL,3),'i')
|
||||
|
||||
# Check number of values for loop sizes:
|
||||
if (Lx.shape[0] != Ly.shape[0]):
|
||||
print "simpleLoopSource:Error: Lx and Ly differ in lengths"
|
||||
raise
|
||||
|
||||
if (Lx.shape[0] == 1):
|
||||
Lx = np.ones((nL,1))*Lx
|
||||
Ly = np.ones((nL,1))*Ly
|
||||
|
||||
# widths of cells:
|
||||
hx = self.x0[0] + np.cumsum(self.hx)
|
||||
hy = self.x0[1] + np.cumsum(self.hy)
|
||||
hz = self.x0[2] + np.cumsum(self.hz)
|
||||
# Normalize the location to the grid:
|
||||
for ii in range(nL):
|
||||
# x
|
||||
diff = abs(loc[ii,0] - hx)
|
||||
iInd = np.argmin(diff)
|
||||
normLoc[ii,0] = iInd + diff[iInd]/hx[iInd]
|
||||
# y
|
||||
diff = abs(loc[ii,1] - hy)
|
||||
jInd = np.argmin(diff)
|
||||
normLoc[ii,1] = jInd + diff[jInd]/hy[jInd]
|
||||
# z
|
||||
diff = abs(loc[ii,2] - hz)
|
||||
kInd = np.argmin(diff)
|
||||
normLoc[ii,2] = kInd - diff[kInd]/hz[kInd]
|
||||
sub[ii,0] = iInd
|
||||
sub[ii,1] = jInd
|
||||
sub[ii,2] = kInd
|
||||
|
||||
# Get node and edge grids
|
||||
self._src = np.zeros((self.edge.shape[0],nL))
|
||||
# Edge grid matrix
|
||||
[Ex,Ey,Ez] = self.r(self.edge,'E','E','M')
|
||||
# No topo for now, but this could be changed
|
||||
EzSource = np.zeros((Ez.shape))
|
||||
ezF = self.r(EzSource,'Ez','Ez','V')
|
||||
# Loops for each loc
|
||||
for ii in range(nL):
|
||||
ExSource = np.zeros((Ex.shape))
|
||||
EySource = np.zeros((Ey.shape))
|
||||
# Find the node closest to loc (round)
|
||||
nodeLoc = np.around(normLoc[ii,:])
|
||||
# Put at bottom of cell
|
||||
nodeLoc[2] = np.floor(normLoc[ii,2])
|
||||
|
||||
# Fill edges in counter-clockwise format
|
||||
# First x:
|
||||
i = sub[ii,0]
|
||||
j = sub[ii,1]
|
||||
k = sub[ii,2]
|
||||
eDist1 = 0.0; eDist2 = 0.0
|
||||
end1 = False; end2 = False
|
||||
i1 = i-1; i2 = i+1
|
||||
#### GET I1,I2,J1,J2 FIRST AND THEN FILL EX AND EY. NEED THEM FOR BOTH
|
||||
while True:
|
||||
try:
|
||||
eDist1 = eDist1 + Ex[i1,j,k]
|
||||
eDist2 = eDist2 + Ex[i2,j,k]
|
||||
except IndexError:
|
||||
print "simpleLoopSource:Error: Loop goes off of x-edge of mesh."
|
||||
raise
|
||||
|
||||
if (np.around(abs(eDist1 - Lx[ii,0]/2)/Ex[i1,j,k]) < 1 and end1 == False):
|
||||
end1 = True
|
||||
else:
|
||||
i1 = i1 - 1
|
||||
|
||||
if (np.around(abs(eDist2 - Lx[ii,0]/2)/Ex[i2,j,k]) < 1 and end2 == False):
|
||||
end2 = True
|
||||
else:
|
||||
i2 = i2 + 1
|
||||
|
||||
if (end1 == True and end2 == True):
|
||||
break
|
||||
|
||||
eDist1 = 0.0; eDist2 = 0.0
|
||||
end1 = False; end2 = False
|
||||
j1 = j-1; j2 = j+1
|
||||
while True:
|
||||
try:
|
||||
eDist1 = eDist1 + Ey[i1,j1,k]
|
||||
eDist2 = eDist2 + Ey[i2,j2,k]
|
||||
except IndexError:
|
||||
print "simpleLoopSource:Error: Loop goes off of y-edge of mesh."
|
||||
raise
|
||||
|
||||
if (np.around(abs(eDist1 - Ly[ii,0]/2)/Ey[i1,j1,k]) < 1 and end1 == False):
|
||||
end1 = True
|
||||
else:
|
||||
j1 = j1 - 1
|
||||
|
||||
|
||||
if (np.around(abs(eDist2 - Ly[ii,0]/2)/Ey[i2,j2,k]) < 1 and end2 == False):
|
||||
end2 = True
|
||||
else:
|
||||
j2 = j2 + 1
|
||||
|
||||
if (end1 == True and end2 == True):
|
||||
break
|
||||
|
||||
# Fill values counter-clockwise:
|
||||
for jj in range(i1,i2):
|
||||
ExSource[jj,j1,k] = 1.0
|
||||
ExSource[jj,j2,k] = -1.0
|
||||
|
||||
for jj in range(j1,j2):
|
||||
EySource[i1,jj,k] = -1.0
|
||||
EySource[i2,jj,k] = 1.0
|
||||
|
||||
exF = self.r(ExSource,'Ex','Ex','V')
|
||||
eyF = self.r(ExSource,'Ey','Ey','V')
|
||||
# Set final e:
|
||||
self._src[:,ii] = np.concatenate((exF, eyF, ezF),axis=0)
|
||||
|
||||
#Finished: return e
|
||||
return self._src
|
||||
Reference in New Issue
Block a user