mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-03 10:36:43 +08:00
Seogi's FDEM initial commit from 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
File diff suppressed because one or more lines are too long
@@ -0,0 +1,73 @@
|
||||
import numpy as np
|
||||
|
||||
def ismember(a, b):
|
||||
tf = np.array([i in b for i in a])
|
||||
return tf
|
||||
|
||||
def path2edgeModel(mesh, pts):
|
||||
edm_x = np.zeros(np.prod(mesh.nEx))
|
||||
edm_y = np.zeros(np.prod(mesh.nEy))
|
||||
edm_z = np.zeros(np.prod(mesh.nEz))
|
||||
|
||||
for ii in range (pts.shape[0]-1):
|
||||
pt1 = pts[ii,:]
|
||||
pt2 = pts[ii+1,:]
|
||||
delta = pt2 - pt1
|
||||
deltaDim = np.argwhere(delta)
|
||||
#assert(np.size(deltaDim)==1), "Path must be orthoginal to mesh"
|
||||
if deltaDim == 0:
|
||||
xLoc = mesh.vectorCCx[(min(pt1[0],pt2[0]) < mesh.vectorCCx ) & (mesh.vectorCCx < max(pt1[0],pt2[0]))]
|
||||
yLoc = pts[ii,1]
|
||||
zLoc = pts[ii,2]
|
||||
delDir = np.sign(pt2[0]-pt1[0])
|
||||
xyz = np.c_[xLoc, np.ones(np.size(xLoc))*yLoc, np.ones(np.size(xLoc))*zLoc]
|
||||
edgeInd=ismember(map(tuple,mesh.gridEx),map(tuple,xyz))
|
||||
edm_x[edgeInd] = delDir
|
||||
# print '>> x-direction', ii
|
||||
# print mesh.gridEx[edgeInd]
|
||||
if deltaDim == 1:
|
||||
xLoc = pts[ii,0]
|
||||
yLoc = mesh.vectorCCy[(min(pt1[1],pt2[1]) < mesh.vectorCCy ) & (mesh.vectorCCy < max(pt1[1],pt2[1]))]
|
||||
zLoc = pts[ii,2]
|
||||
delDir = np.sign(pt2[1]-pt1[1])
|
||||
xyz = np.c_[np.ones(np.size(yLoc))*xLoc, yLoc, np.ones(np.size(yLoc))*zLoc]
|
||||
edgeInd=ismember(map(tuple,mesh.gridEy),map(tuple,xyz))
|
||||
edm_y[edgeInd] = delDir
|
||||
# print '>> y-direction', ii
|
||||
# print mesh.gridEy[edgeInd]
|
||||
if deltaDim == 2:
|
||||
xLoc = pts[ii,0]
|
||||
yLoc = pts[ii,1]
|
||||
zLoc = mesh.vectorCCz[(min(pt1[2],pt2[2]) < mesh.vectorCCz ) & (mesh.vectorCCz < max(pt1[2],pt2[2]))]
|
||||
delDir = np.sign(pt2[2]-pt1[2])
|
||||
xyz = np.c_[np.ones(np.size(zLoc))*xLoc, np.ones(np.size(zLoc))*yLoc, zLoc]
|
||||
edgeInd=ismember(map(tuple,mesh.gridEz),map(tuple,xyz))
|
||||
edm_z[edgeInd] = delDir
|
||||
# print '>> z-direction', ii
|
||||
# print mesh.gridEz[edgeInd]
|
||||
|
||||
|
||||
edgeModel = np.r_[edm_x, edm_y, edm_z]
|
||||
return edgeModel
|
||||
|
||||
def rho(x1, y1, x, y):
|
||||
r = np.sqrt((x-x1)**2+(y-y1)**2)
|
||||
return r
|
||||
|
||||
def MMRhalf(loc1, loc2, x, y):
|
||||
""" Anaytic function for MMR response (B^{1D})
|
||||
- loc1=(x1,y1): x, y location for (+) charge
|
||||
- loc2=(x2,y2): x, y1
|
||||
- x : observation points in x-direction
|
||||
- y : observation points in y-direction
|
||||
"""
|
||||
x1=loc1[0]
|
||||
x2=loc2[0]
|
||||
y1=loc1[1]
|
||||
y2=loc2[1]
|
||||
mu0 = 4*np.pi*1e-7
|
||||
I = 1
|
||||
By =mu0*I/(4*np.pi)*np.array((x-x1)/rho(x1,y1,x,y)**2-(x-x2)/rho(x2,y2,x,y)**2)
|
||||
Bx =mu0*I/(4*np.pi)*np.array(-(y-y1)/rho(x1,y1,x,y)**2+(y-y2)/rho(x2,y2,x,y)**2)
|
||||
|
||||
return Bx, By
|
||||
@@ -0,0 +1,156 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
from SimPEG import utils, TensorMesh
|
||||
from SimPEG.utils import spzeros, mkvc
|
||||
|
||||
def interpmat(x,y,z,xr,yr,zr):
|
||||
|
||||
""" Local nterpolation computed for each receiver point in turn """
|
||||
|
||||
nx = max(x.shape)
|
||||
ny = max(y.shape)
|
||||
nz = max(z.shape)
|
||||
npts = max(xr.shape)
|
||||
|
||||
Q = sp.lil_matrix((npts, nx*ny*nz))
|
||||
|
||||
for i in range(npts):
|
||||
# in x-direction
|
||||
im = np.argmin(abs(x-xr[i]))
|
||||
if xr[i] - x[im] >= 0: # Point on the left
|
||||
ind_x1 = im
|
||||
ind_x2 = im+1
|
||||
elif xr[i] - x[im] < 0: # Point on the right
|
||||
ind_x1 = im-1
|
||||
ind_x2 = im
|
||||
dx1 = xr[i] - x[ind_x1]
|
||||
dx2 = x[ind_x2] - xr[i]
|
||||
# in y-direction
|
||||
im = np.argmin(abs(y-yr[i]))
|
||||
if yr[i] - y[im] >= 0: # Point on the left
|
||||
ind_y1 = im
|
||||
ind_y2 = im+1
|
||||
elif yr[i] - y[im] < 0: # Point on the right
|
||||
ind_y1 = im-1
|
||||
ind_y2 = im
|
||||
dy1 = yr[i] - y[ind_y1]
|
||||
dy2 = y[ind_y2] - yr[i]
|
||||
# in z-direction
|
||||
im = np.argmin(abs(z-zr[i]))
|
||||
if zr[i] - z[im] >= 0: # Point on the left
|
||||
ind_z1 = im
|
||||
ind_z2 = im+1
|
||||
elif zr[i] - z[im] < 0: # Point on the right
|
||||
ind_z1 = im-1
|
||||
ind_z2 = im
|
||||
dz1 = zr[i] - z[ind_z1]
|
||||
dz2 = z[ind_z2] - zr[i]
|
||||
dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1])
|
||||
|
||||
Dx = x[ind_x2] - x[ind_x1]
|
||||
Dy = y[ind_y2] - y[ind_y1]
|
||||
Dz = z[ind_z2] - z[ind_z1]
|
||||
|
||||
# Get the row in the matrix
|
||||
|
||||
inds = utils.sub2ind((nx,ny,nz),[
|
||||
( ind_x1, ind_y2, ind_z1),
|
||||
( ind_x1, ind_y1, ind_z1),
|
||||
( ind_x2, ind_y1, ind_z1),
|
||||
( ind_x2, ind_y2, ind_z1),
|
||||
( ind_x1, ind_y1, ind_z2),
|
||||
( ind_x1, ind_y2, ind_z2),
|
||||
( ind_x2, ind_y1, ind_z2),
|
||||
( ind_x2, ind_y2, ind_z2)])
|
||||
|
||||
vals = [(1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz),
|
||||
(1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz),
|
||||
(1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz),
|
||||
(1-dx2/Dx)*(1-dy2/Dy)*(1-dz1/Dz),
|
||||
(1-dx1/Dx)*(1-dy1/Dy)*(1-dz2/Dz),
|
||||
(1-dx1/Dx)*(1-dy2/Dy)*(1-dz2/Dz),
|
||||
(1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz),
|
||||
(1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz)]
|
||||
|
||||
Q[i, mkvc(inds)] = vals
|
||||
Q = Q.tocsr()
|
||||
return Q
|
||||
|
||||
def getInterpmat(mesh, rxLoc, dataType):
|
||||
""" """
|
||||
xr = rxLoc[:,0]
|
||||
yr = rxLoc[:,1]
|
||||
zr = rxLoc[:,2]
|
||||
nrx = rxLoc.shape[0]
|
||||
if dataType == 'fx':
|
||||
Qx = interpmat(np.unique(mesh.gridFx[:,0]),
|
||||
np.unique(mesh.gridFx[:,1]),
|
||||
np.unique(mesh.gridFx[:,2]),
|
||||
xr,yr,zr)
|
||||
Q = sp.hstack([Qx,spzeros(nrx,mesh.nF[1]),spzeros(nrx,mesh.nF[2])])
|
||||
elif dataType == 'fy':
|
||||
Qy = interpmat(np.unique(mesh.gridFy[:,0]),
|
||||
np.unique(mesh.gridFy[:,1]),
|
||||
np.unique(mesh.gridFy[:,2]),
|
||||
xr,yr,zr)
|
||||
Q = sp.hstack([spzeros(nrx,mesh.nF[0]),Qy,spzeros(nrx,mesh.nF[2])])
|
||||
elif dataType == 'fz':
|
||||
Qz = interpmat(np.unique(mesh.gridFz[:,0]),
|
||||
np.unique(mesh.gridFz[:,1]),
|
||||
np.unique(mesh.gridFz[:,2]),
|
||||
xr,yr,zr)
|
||||
Q = sp.hstack([spzeros(nrx,mesh.nF[0]),spzeros(nrx,mesh.nF[1]),Qz])
|
||||
elif dataType == 'ex':
|
||||
Qx = interpmat(np.unique(mesh.gridEx[:,0]),
|
||||
np.unique(mesh.gridEx[:,1]),
|
||||
np.unique(mesh.gridEx[:,2]),
|
||||
xr, yr, zr)
|
||||
Q = sp.hstack([Qx,spzeros(nrx,mesh.nE[1]),spzeros(nrx,mesh.nE[2])])
|
||||
elif dataType == 'ey':
|
||||
Qy = interpmat(np.unique(mesh.gridEy[:,0]),
|
||||
np.unique(mesh.gridEy[:,1]),
|
||||
np.unique(mesh.gridEy[:,2]),
|
||||
xr, yr, zr)
|
||||
Q = sp.hstack([spzeros(nrx,mesh.nE[0]),Qy,spzeros(nrx,mesh.nE[2])])
|
||||
elif dataType == 'ez':
|
||||
Qz = interpmat(np.unique(mesh.gridEz[:,0]),
|
||||
np.unique(mesh.gridEz[:,1]),
|
||||
np.unique(mesh.gridEz[:,2]),
|
||||
xr,yr,zr)
|
||||
Q = sp.hstack([spzeros(nrx,mesh.nE[0]),spzeros(nrx,mesh.nE[1]),Qz])
|
||||
else:
|
||||
assert(True), "Input either face (fx, fy, fz) or edge (ex, ey, ez) option"
|
||||
return Q
|
||||
|
||||
if __name__ == '__main__':
|
||||
pad = 1
|
||||
padfactor = 1.5
|
||||
cs = 100
|
||||
xpad = cs*(np.ones(pad)*padfactor)**np.arange(pad)
|
||||
ypad = cs*(np.ones(pad)*padfactor)**np.arange(pad)
|
||||
zpad = cs*(np.ones(pad)*padfactor)**np.arange(pad)
|
||||
|
||||
core = 10
|
||||
xcore = cs*np.ones(core)
|
||||
ycore = cs*np.ones(core)
|
||||
zcore = cs*np.ones(core)
|
||||
|
||||
hx = np.r_[xpad[::-1],xcore, cs, xcore,xpad]
|
||||
hy = np.r_[ypad[::-1],ycore, cs, ycore, ypad]
|
||||
hz = np.r_[zpad[::-1],zcore,zcore, zpad]
|
||||
x0 = np.array([-np.sum(hx)/2, -np.sum(hy)/2, -np.sum(hz)/2], )
|
||||
mesh = TensorMesh([hx, hy, hz],x0)
|
||||
|
||||
xr1 = np.linspace(-500,500,5)
|
||||
yr1 = np.linspace(-500,500,5)
|
||||
zr1 = 0
|
||||
xr, yr = np.meshgrid(xr1, yr1, indexing='ij')
|
||||
zr = np.ones((xr.shape[0],xr.shape[1]))*zr1
|
||||
xr = mkvc(xr)
|
||||
yr = mkvc(yr)
|
||||
zr = mkvc(zr)
|
||||
rxLoc = np.c_[xr, yr, zr]
|
||||
Q = getInterpmat(mesh, rxLoc, 'ex')
|
||||
|
||||
print Q
|
||||
|
||||
Reference in New Issue
Block a user