wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-29 11:27:23 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' of https://github.com/simpeg/simpeg into TreeMesh

Browse Source
This commit is contained in:
Rowan Cockett
2015-04-04 14:45:42 -07:00
parent 6ffeb5cc1a 905e877e6b
commit 6e4aac140b
9 changed files with 296 additions and 77 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/CylMesh.py
+61 -2
View File
@@ -1,7 +1,7 @@
import numpy as np
import scipy.sparse as sp
from scipy.constants import pi
from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, ddx, av, avExtrap
from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, spzeros, ddx, av, avExtrap
from TensorMesh import BaseTensorMesh
from InnerProducts import InnerProducts
from View import CylView
@@ -27,11 +27,14 @@ class CylMesh(BaseTensorMesh, InnerProducts, CylView):
_unitDimensions = [1, 2*np.pi, 1]
def __init__(self, h, x0=None):
def __init__(self, h, x0=None, cartesianOrigin=None):
BaseTensorMesh.__init__(self, h, x0)
assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi"
if self.dim == 2:
print 'Warning, a disk mesh has not been tested thoroughly.'
cartesianOrigin = np.zeros(self.dim) if cartesianOrigin is None else cartesianOrigin
assert len(cartesianOrigin) == self.dim, "cartesianOrigin must be the same length as the dimension of the mesh."
self.cartesianOrigin = np.array(cartesianOrigin, dtype=float)
@property
@@ -327,6 +330,62 @@ class CylMesh(BaseTensorMesh, InnerProducts, CylView):
raise NotImplementedError('wrapping in the averaging is not yet implemented')
return self._aveF2CCV
def getInterpolationMatCartMesh(self, Mrect, locType='CC'):
"""
Takes a cartesian mesh and returns a projection to translate onto the cartesian grid.
"""
assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes"
if locType == 'F':
# do this three times for each component
X = self.getInterpolationMatCartMesh(Mrect, locType='Fx')
Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy')
Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz')
return sp.vstack((X,Y,Z))
if locType == 'E':
X = self.getInterpolationMatCartMesh(Mrect, locType='Ex')
Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey')
Z = spzeros(Mrect.nEz, self.nE)
return sp.vstack((X,Y,Z))
grid = getattr(Mrect, 'grid' + locType)
# This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice
theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2
theta[theta < 0] += np.pi*2.0
r = ((grid[:,0] - self.cartesianOrigin[0])**2 + (grid[:,1] - self.cartesianOrigin[1])**2)**0.5
if locType in ['CC', 'N', 'Fz', 'Ez']:
G, proj = np.c_[r, theta, grid[:,2]], np.ones(r.size)
else:
dotMe = {
'Fx': Mrect.normals[:Mrect.nFx,:],
'Fy': Mrect.normals[Mrect.nFx:(Mrect.nFx+Mrect.nFy),:],
'Fz': Mrect.normals[-Mrect.nFz:,:],
'Ex': Mrect.tangents[:Mrect.nEx,:],
'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:],
'Ez': Mrect.tangents[-Mrect.nEz:,:],
}[locType]
if 'F' in locType:
normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)]
proj = ( normals * dotMe ).sum(axis=1)
if 'E' in locType:
tangents = np.c_[-np.sin(theta), np.cos(theta), np.zeros(theta.size)]
proj = ( tangents * dotMe ).sum(axis=1)
G = np.c_[r, theta, grid[:,2]]
interpType = locType
if interpType == 'Fy':
interpType = 'Fx'
elif interpType == 'Ex':
interpType = 'Ey'
Pc2r = self.getInterpolationMat(G, interpType)
Proj = sdiag(proj)
return Proj * Pc2r
if __name__ == '__main__':
SimPEG/Mesh/DiffOperators.py
+89 -42
View File
@@ -571,35 +571,58 @@ class DiffOperators(object):
@property
def aveF2CC(self):
"Construct the averaging operator on cell faces to cell centers."
if getattr(self, '_aveF2CC', None) is None:
n = self.vnC
if(self.dim == 1):
self._aveF2CC = av(n[0])
elif(self.dim == 2):
self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])),
sp.kron(av(n[1]), speye(n[0]))), format="csr")
elif(self.dim == 3):
self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
kron3(speye(n[2]), av(n[1]), speye(n[0])),
kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
return self._aveF2CC
if(self.dim == 1):
return self.aveFx2CC
elif(self.dim == 2):
return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC), format="csr")
elif(self.dim == 3):
return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr")
@property
def aveF2CCV(self):
"Construct the averaging operator on cell faces to cell centers."
if getattr(self, '_aveF2CCV', None) is None:
if(self.dim == 1):
return self.aveFx2CC
elif(self.dim == 2):
return sp.block_diag((self.aveFx2CC, self.aveFy2CC), format="csr")
elif(self.dim == 3):
return sp.block_diag((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr")
@property
def aveFx2CC(self):
"Construct the averaging operator on cell faces in the x direction to cell centers."
if getattr(self, '_aveFx2CC', None) is None:
n = self.vnC
if(self.dim == 1):
self._aveF2CCV = av(n[0])
self._aveFx2CC = av(n[0])
elif(self.dim == 2):
self._aveF2CCV = sp.block_diag((sp.kron(speye(n[1]), av(n[0])),
sp.kron(av(n[1]), speye(n[0]))), format="csr")
self._aveFx2CC = sp.kron(speye(n[1]), av(n[0]))
elif(self.dim == 3):
self._aveF2CCV = sp.block_diag((kron3(speye(n[2]), speye(n[1]), av(n[0])),
kron3(speye(n[2]), av(n[1]), speye(n[0])),
kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
return self._aveF2CCV
self._aveFx2CC = kron3(speye(n[2]), speye(n[1]), av(n[0]))
return self._aveFx2CC
@property
def aveFy2CC(self):
"Construct the averaging operator on cell faces in the y direction to cell centers."
if self.dim < 2: return None
if getattr(self, '_aveFy2CC', None) is None:
n = self.vnC
if(self.dim == 2):
self._aveFy2CC = sp.kron(av(n[1]), speye(n[0]))
elif(self.dim == 3):
self._aveFy2CC = kron3(speye(n[2]), av(n[1]), speye(n[0]))
return self._aveFy2CC
@property
def aveFz2CC(self):
"Construct the averaging operator on cell faces in the z direction to cell centers."
if self.dim < 3: return None
if getattr(self, '_aveFz2CC', None) is None:
n = self.vnC
if(self.dim == 3):
self._aveFz2CC = kron3(av(n[2]), speye(n[1]), speye(n[0]))
return self._aveFz2CC
@property
def aveCC2F(self):
@@ -620,36 +643,60 @@ class DiffOperators(object):
@property
def aveE2CC(self):
"Construct the averaging operator on cell edges to cell centers."
if getattr(self, '_aveE2CC', None) is None:
# The number of cell centers in each direction
n = self.vnC
if(self.dim == 1):
self._aveE2CC = speye(n[0])
elif(self.dim == 2):
self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])),
sp.kron(speye(n[1]), av(n[0]))), format="csr")
elif(self.dim == 3):
self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
kron3(av(n[2]), speye(n[1]), av(n[0])),
kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
return self._aveE2CC
if(self.dim == 1):
return self.aveEx2CC
elif(self.dim == 2):
return 0.5*sp.hstack((self.aveEx2CC, self.aveEy2CC), format="csr")
elif(self.dim == 3):
return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr")
@property
def aveE2CCV(self):
"Construct the averaging operator on cell edges to cell centers."
if getattr(self, '_aveE2CCV', None) is None:
if(self.dim == 1):
return self.aveEx2CC
elif(self.dim == 2):
return sp.block_diag((self.aveEx2CC, self.aveEy2CC), format="csr")
elif(self.dim == 3):
return sp.block_diag((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr")
@property
def aveEx2CC(self):
"Construct the averaging operator on cell edges in the x direction to cell centers."
if getattr(self, '_aveEx2CC', None) is None:
# The number of cell centers in each direction
n = self.vnC
if(self.dim == 1):
raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
self._aveEx2CC = speye(n[0])
elif(self.dim == 2):
self._aveE2CCV = sp.block_diag((sp.kron(av(n[1]), speye(n[0])),
sp.kron(speye(n[1]), av(n[0]))), format="csr")
self._aveEx2CC = sp.kron(av(n[1]), speye(n[0]))
elif(self.dim == 3):
self._aveE2CCV = sp.block_diag((kron3(av(n[2]), av(n[1]), speye(n[0])),
kron3(av(n[2]), speye(n[1]), av(n[0])),
kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
return self._aveE2CCV
self._aveEx2CC = kron3(av(n[2]), av(n[1]), speye(n[0]))
return self._aveEx2CC
@property
def aveEy2CC(self):
"Construct the averaging operator on cell edges in the y direction to cell centers."
if self.dim < 2: return None
if getattr(self, '_aveEy2CC', None) is None:
# The number of cell centers in each direction
n = self.vnC
if(self.dim == 2):
self._aveEy2CC = sp.kron(speye(n[1]), av(n[0]))
elif(self.dim == 3):
self._aveEy2CC = kron3(av(n[2]), speye(n[1]), av(n[0]))
return self._aveEy2CC
@property
def aveEz2CC(self):
"Construct the averaging operator on cell edges in the z direction to cell centers."
if self.dim < 3: return None
if getattr(self, '_aveEz2CC', None) is None:
# The number of cell centers in each direction
n = self.vnC
if(self.dim == 3):
self._aveEz2CC = kron3(speye(n[2]), av(n[1]), av(n[0]))
return self._aveEz2CC
@property
def aveN2CC(self):
SimPEG/Mesh/TensorMesh.py
+11 -9
View File
@@ -168,21 +168,21 @@ class BaseTensorMesh(BaseRectangularMesh):
"""
if key is 'Fx':
if key == 'Fx':
ten = [self.vectorNx , self.vectorCCy, self.vectorCCz]
elif key is 'Fy':
elif key == 'Fy':
ten = [self.vectorCCx, self.vectorNy , self.vectorCCz]
elif key is 'Fz':
elif key == 'Fz':
ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ]
elif key is 'Ex':
elif key == 'Ex':
ten = [self.vectorCCx, self.vectorNy , self.vectorNz ]
elif key is 'Ey':
elif key == 'Ey':
ten = [self.vectorNx , self.vectorCCy, self.vectorNz ]
elif key is 'Ez':
elif key == 'Ez':
ten = [self.vectorNx , self.vectorNy , self.vectorCCz]
elif key is 'CC':
elif key == 'CC':
ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz]
elif key is 'N':
elif key == 'N':
ten = [self.vectorNx , self.vectorNy , self.vectorNz ]
return [t for t in ten if t is not None]
@@ -204,10 +204,12 @@ class BaseTensorMesh(BaseRectangularMesh):
if locType == 'N' and self._meshType == 'CYL':
#NOTE: for a CYL mesh we add a node to check if we are inside in the radial direction!
tensors[0] = np.r_[0.,tensors[0]]
tensors[1] = np.r_[tensors[1], 2.0*np.pi]
inside = np.ones(pts.shape[0],dtype=bool)
for i, tensor in enumerate(tensors):
inside = inside & (pts[:,i] >= tensor.min()) & (pts[:,i] <= tensor.max())
TOL = np.diff(tensor).min() * 1.0e-10
inside = inside & (pts[:,i] >= tensor.min()-TOL) & (pts[:,i] <= tensor.max()+TOL)
return inside
def getInterpolationMat(self, loc, locType, zerosOutside=False):
SimPEG/Mesh/View.py
+20 -3
View File
@@ -109,6 +109,10 @@ class TensorView(object):
vc = (self.aveN2CC*v).reshape(self.vnC, order='F')
elif vType in ['Fx', 'Fy', 'Fz', 'Ex', 'Ey', 'Ez']:
aveOp = 'ave' + vType[0] + '2CCV'
# n = getattr(self,'vn'+vType[0])
# if 'x' in vType: v = np.r_[v,np.zeros(n[1]),np.zeros(n[2])]
# if 'y' in vType: v = np.r_[np.zeros(n[0]),v,np.zeros(n[2])]
# if 'z' in vType: v = np.r_[np.zeros(n[0]),np.zeros(n[1]),v]
v = getattr(self,aveOp)*v # average to cell centers
ind_xyz = {'x':0,'y':1,'z':2}[vType[1]]
vc = self.r(v.reshape((self.nC,-1),order='F'), 'CC','CC','M')[ind_xyz]
@@ -190,9 +194,16 @@ class TensorView(object):
M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
"""
if type(vType) in [list, tuple]:
assert ax is None, "cannot specify an axis to plot on with this function."
fig, axs = plt.subplots(1,len(vType))
out = []
for vTypeI, ax in zip(vType, axs):
out += [self.plotSlice(v,vType=vTypeI, normal=normal, ind=ind, grid=grid, view=view, ax=ax, clim=clim, showIt=False, pcolorOpts=pcolorOpts, streamOpts=streamOpts, gridOpts=gridOpts)]
return out
viewOpts = ['real','imag','abs','vec']
normalOpts = ['X', 'Y', 'Z']
vTypeOpts = ['CC', 'CCv','F','E']
vTypeOpts = ['CC', 'CCv','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez']
# Some user error checking
assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
@@ -219,9 +230,15 @@ class TensorView(object):
elif vType == 'CCv':
assert view == 'vec', 'Other types for CCv not supported'
else:
# Now just deal with 'F' and 'E'
# Now just deal with 'F' and 'E' (x,y,z, maybe...)
aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC')
v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
Av = getattr(self,aveOp)
if v.size == Av.shape[1]:
v = Av * v
else:
v = self.r(v,vType[0],vType) # get specific component
v = Av * v
# we should now be averaged to cell centers (might be a vector)
v = self.r(v.reshape((self.nC,-1),order='F'),'CC','CC','M')
if view == 'vec':
outSlice = []
SimPEG/Problem.py
+1 -13
View File
@@ -437,19 +437,7 @@ class BaseTimeProblem(BaseProblem):
del self.timeMesh
return
if type(value) is not list:
raise Exception('timeSteps must be a np.ndarray or a list of scalars and tuples.')
proposed = []
for v in value:
if Utils.isScalar(v):
proposed += [float(v)]
elif type(v) is tuple and len(v) == 2:
proposed += [float(v[0])]*int(v[1])
else:
raise Exception('timeSteps list must contain only scalars and len(2) tuples.')
self._timeSteps = np.array(proposed)
self._timeSteps = Utils.meshTensor(value)
del self.timeMesh
@property
SimPEG/Tests/TestUtils.py
+2 -2
View File
@@ -1,7 +1,7 @@
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm
from SimPEG.Utils import mkvc, sdiag
from SimPEG.Utils import mkvc, sdiag, diagEst
from SimPEG import Utils
from SimPEG.Mesh import TensorMesh, LogicallyRectMesh, CylMesh
import numpy as np
@@ -267,7 +267,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
# Ensure we are about precision
order0 = order0[E0[1:] > eps]
order1 = order1[E1[1:] > eps]
belowTol = order1.size == 0 and order0.size > 0
belowTol = (order1.size == 0 and order0.size >= 0)
# Make sure we get the correct order
correctOrder = order1.size > 0 and np.mean(order1) > tolerance * expectedOrder
SimPEG/Tests/test_cylMesh.py
+74
View File
@@ -133,6 +133,80 @@ class TestCyl2DMesh(unittest.TestCase):
def test_lightOperators(self):
self.assertTrue(self.mesh.nodalGrad is None)
def test_getInterpMatCartMesh_Cells(self):
Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
mc = np.arange(Mc.nC)
xr = np.linspace(0,0.4,50)
xc = np.linspace(0,0.4,50) + 0.2
Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'CC')
Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC')
Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC')
assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
def test_getInterpMatCartMesh_Faces(self):
Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
Pf = Mc.getInterpolationMatCartMesh(Mr, 'F')
mf = np.ones(Mc.nF)
frect = Pf * mf
fxcc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx')
fycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy')
fzcc = Mr.r(frect, 'F', 'Fz')
indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
TOL = 1e-2
assert np.abs(float(fxcc[indX]) - 1) < TOL
assert np.abs(float(fxcc[indY]) - 0) < TOL
assert np.abs(float(fycc[indX]) - 0) < TOL
assert np.abs(float(fycc[indY]) - 1) < TOL
assert np.abs((fzcc - 1).sum()) < TOL
mag = (fxcc**2 + fycc**2)**0.5
dist = ((Mr.gridCC[:,0] + 0.2)**2 + (Mr.gridCC[:,1] + 0.2)**2)**0.5
assert np.abs(mag[dist > 0.1].max() - 1) < TOL
assert np.abs(mag[dist > 0.1].min() - 1) < TOL
def test_getInterpMatCartMesh_Edges(self):
Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
Pe = Mc.getInterpolationMatCartMesh(Mr, 'E')
me = np.ones(Mc.nE)
erect = Pe * me
excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex')
eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey')
ezcc = Mr.r(erect, 'E', 'Ez')
indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
TOL = 1e-2
assert np.abs(float(excc[indX]) - 0) < TOL
assert np.abs(float(excc[indY]) + 1) < TOL
assert np.abs(float(eycc[indX]) - 1) < TOL
assert np.abs(float(eycc[indY]) - 0) < TOL
assert np.abs(ezcc.sum()) < TOL
mag = (excc**2 + eycc**2)**0.5
dist = ((Mr.gridCC[:,0] + 0.2)**2 + (Mr.gridCC[:,1] + 0.2)**2)**0.5
assert np.abs(mag[dist > 0.1].max() - 1) < TOL
assert np.abs(mag[dist > 0.1].min() - 1) < TOL
MESHTYPES = ['uniformCylMesh']
SimPEG/Tests/test_utils.py
+13
View File
@@ -246,5 +246,18 @@ class TestSequenceFunctions(unittest.TestCase):
self.assertTrue(np.all(true == listArray))
self.assertTrue(true.shape == listArray.shape)
class TestDiagEst(unittest.TestCase):
def setUp(self):
self.n = 10
self.A = np.random.rand(self.n,self.n)
self.Adiag = np.diagonal(self.A)
def testOnes(self):
Adiagtest = diagEst(self.A,self.n,self.n)
r = np.abs(Adiagtest-self.Adiag)
self.assertTrue(r.dot(r) < TOL)
if __name__ == '__main__':
unittest.main()
SimPEG/Utils/matutils.py
+25 -6
View File
@@ -341,30 +341,49 @@ def invPropertyTensor(M, tensor, returnMatrix=False):
return T
def diagEst(matFun, n, k=None, type='Probing'):
""" Based on Saad http://www-users.cs.umn.edu/~saad/PDF/umsi-2005-082.pdf, and http://www.cita.utoronto.ca/~niels/diagonal.pdf"""
def diagEst(matFun, n, k=None, approach='Probing'):
"""
Estimate the diagonal of a matrix, A. Note that the matrix may be a function which returns A times a vector.
Three different approaches have been implemented,
1. Probing : uses cyclic permutations of vectors with ones and zeros (default)
2. Ones : random +/- 1 entries
3. Random : random vectors
:param lambda (numpy.array) matFun: matrix to estimate the diagonal of
:param int64 n: size of the vector that should be used to compute matFun(v)
:param int64 k: number of vectors to be used to estimate the diagonal
:param str approach: approach to be used for getting vectors
:rtype: numpy.array
:return: est_diag(A)
Based on Saad http://www-users.cs.umn.edu/~saad/PDF/umsi-2005-082.pdf, and http://www.cita.utoronto.ca/~niels/diagonal.pdf
"""
if type(matFun).__name__=='ndarray':
A = matFun
matFun = lambda v: A.dot(v)
if k is None:
k = np.floor(n/10.)
if type =='Ones':
if approach =='Ones':
def getv(n,i=None):
v = np.random.randn(n)
v[v<0] = -1.
v[v>=0] = 1.
return v
elif type == 'Random':
elif approach == 'Random':
def getv(n,i=None):
return np.random.randn(n)
else: #if type == 'Probing':
else: #if approach == 'Probing':
def getv(n,i):
v = np.zeros(n)
v[i:n:k] = 1.
return v
#d = np.zeros(n)
Mv = np.zeros(n)
vv = np.zeros(n)