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

Merge branch 'develop' of https://github.com/simpeg/simpeg into boundaryConditions

Browse Source
This commit is contained in:
rowanc1
2014-02-24 10:11:24 -08:00
parent 9c881b5b5d 33431c4e5f
commit ab399a2a28
15 changed files with 527 additions and 245 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/InnerProducts.py
+3 -40
View File
@@ -1,5 +1,5 @@
from scipy import sparse as sp
from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor
import numpy as np
@@ -174,7 +174,7 @@ class InnerProducts(object):
P011 = V3*Pxxx('fXm', 'fYp', 'fZp')
P111 = V3*Pxxx('fXp', 'fYp', 'fZp')
Mu = _makeTensor(M, mu)
Mu = makePropertyTensor(M, mu)
A = P000.T*Mu*P000 + P100.T*Mu*P100
P = [P000, P100]
@@ -283,7 +283,7 @@ class InnerProducts(object):
P011 = V*eP('eX3', 'eY2', 'eZ2')
P111 = V*eP('eX3', 'eY3', 'eZ3')
Sigma = _makeTensor(M, sigma)
Sigma = makePropertyTensor(M, sigma)
A = P000.T*Sigma*P000 + P100.T*Sigma*P100 + P010.T*Sigma*P010 + P110.T*Sigma*P110
P = [P000, P100, P010, P110]
if M.dim == 3:
@@ -313,43 +313,6 @@ class InnerProducts(object):
# | |/
# node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
def _makeTensor(M, sigma):
if sigma is None: # default is ones
sigma = np.ones((M.nC, 1))
if M.dim == 1:
if sigma.size == M.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
Sigma = sdiag(sigma)
else:
raise Exception('Unexpected shape of sigma')
elif M.dim == 2:
if sigma.size == M.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
Sigma = sdiag(np.r_[sigma, sigma])
elif sigma.shape[1] == 2: # Diagonal tensor
Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]])
elif sigma.shape[1] == 3: # Fully anisotropic
row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2])))
row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1])))
Sigma = sp.vstack((row1, row2))
else:
raise Exception('Unexpected shape of sigma')
elif M.dim == 3:
if sigma.size == M.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
Sigma = sdiag(np.r_[sigma, sigma, sigma])
elif sigma.shape[1] == 3: # Diagonal tensor
Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]])
elif sigma.shape[1] == 6: # Fully anisotropic
row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4])))
row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5])))
row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2])))
Sigma = sp.vstack((row1, row2, row3))
else:
raise Exception('Unexpected shape of sigma')
return Sigma
def _getFacePx(M):
assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh'
SimPEG/Mesh/TensorMesh.py
+33 -10
View File
@@ -367,7 +367,7 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
return [t for t in ten if t is not None]
def isInside(self, pts):
def isInside(self, pts, locType='N'):
"""
Determines if a set of points are inside a mesh.
@@ -376,15 +376,23 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
:return inside, numpy array of booleans
"""
pts = np.atleast_2d(pts)
inside = (pts[:,0] >= self.vectorNx.min()) & (pts[:,0] <= self.vectorNx.max())
tensors = self.getTensor(locType)
if type(pts) == list:
pts = np.array(pts)
assert type(pts) == np.ndarray, "must be a numpy array"
if self.dim > 1:
inside = inside & ((pts[:,1] >= self.vectorNy.min()) & (pts[:,1] <= self.vectorNy.max()))
if self.dim > 2:
inside = inside & ((pts[:,2] >= self.vectorNz.min()) & (pts[:,2] <= self.vectorNz.max()))
assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)"
elif len(pts.shape) == 1:
pts = pts[:,np.newaxis]
else:
assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)"
inside = np.ones(pts.shape[0],dtype=bool)
for i, tensor in enumerate(tensors):
inside = inside & (pts[:,i] >= tensor.min()) & (pts[:,i] <= tensor.max())
return inside
def getInterpolationMat(self, loc, locType):
def getInterpolationMat(self, loc, locType, zerosOutside=False):
""" Produces interpolation matrix
:param numpy.ndarray loc: Location of points to interpolate to
@@ -404,8 +412,21 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
'CC' -> scalar field defined on cell centers
"""
loc = np.atleast_2d(loc)
assert np.all(self.isInside(loc)), "Points outside of mesh"
if type(loc) == list:
loc = np.array(loc)
assert type(loc) == np.ndarray, "must be a numpy array"
if self.dim > 1:
assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)"
elif len(loc.shape) == 1:
loc = loc[:,np.newaxis]
else:
assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)"
if zerosOutside is False:
assert np.all(self.isInside(loc)), "Points outside of mesh"
else:
indZeros = np.logical_not(self.isInside(loc))
loc[indZeros, :] = np.array([v.mean() for v in self.getTensor('CC')])
ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1
if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind:
@@ -417,7 +438,9 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
Q = Utils.interpmat(loc, *self.getTensor(locType))
else:
raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))
return Q
if zerosOutside:
Q[indZeros, :] = 0
return Q.tocsr()
@property
SimPEG/Mesh/TensorView.py
+62 -9
View File
@@ -243,11 +243,30 @@ class TensorView(object):
def plotSlice(self, v, vType='CC',
normal='Z', ind=None, grid=False, view='real',
ax=None, clim=None, showIt=False,
pcolorOpts={},
streamOpts={'color':'k'},
gridOpts={'color':'k'}
):
"""
Plots a slice of a 3D mesh.
.. plot::
from SimPEG import *
mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2)))
M = Mesh.TensorMesh(mT)
q = np.zeros(M.vnC)
q[[4,4],[4,4],[2,6]]=[-1,1]
q = Utils.mkvc(q)
A = M.faceDiv*M.cellGrad
b = Solver(A).solve(q)
M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
"""
viewOpts = ['real','imag','abs','vec']
normalOpts = ['X', 'Y', 'Z']
vTypeOpts = ['CC','F','E']
vTypeOpts = ['CC', 'CCv','F','E']
# Some user error checking
assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
@@ -279,11 +298,15 @@ class TensorView(object):
def doSlice(v):
if vType == 'CC':
return getIndSlice(self.r(v,'CC','CC','M'))
# Now just deal with 'F' and 'E'
aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC')
v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
if view == 'vec':
elif vType == 'CCv':
v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M')
assert view == 'vec', 'Other types for CCv not yet supported'
else:
# Now just deal with 'F' and 'E'
aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC')
v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M')
if view == 'vec':
outSlice = []
if 'X' not in normal: outSlice.append(getIndSlice(v[0]))
if 'Y' not in normal: outSlice.append(getIndSlice(v[1]))
@@ -304,13 +327,31 @@ class TensorView(object):
v = doSlice(v)
if clim is None:
clim = [v.min(),v.max()]
out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, v.T, vmin=clim[0], vmax=clim[1]),)
out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, v.T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
elif view in ['vec']:
U, V = doSlice(v)
if clim is None:
clim = [v.min(),v.max()]
out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, 0.5*(U+V).T, vmin=clim[0], vmax=clim[1]),)
out += (plt.streamplot(tM.vectorCCx, tM.vectorCCy, U.T, V.T),)
uv = np.r_[mkvc(U), mkvc(V)]
uv = np.sqrt(uv**2)
clim = [uv.min(),uv.max()]
# Matplotlib seems to not support irregular
# spaced vectors at the moment. So we will
# Interpolate down to a regular mesh at the
# smallest mesh size in this 2D slice.
nxi = int(tM.hx.sum()/tM.hx.min())
nyi = int(tM.hy.sum()/tM.hy.min())
tMi = self.__class__([np.ones(nxi)*tM.hx.sum()/nxi,
np.ones(nyi)*tM.hy.sum()/nyi])
P = tM.getInterpolationMat(tMi.gridCC,'CC',zerosOutside=True)
Ui = P*mkvc(U)
Vi = P*mkvc(V)
Ui = tMi.r(Ui, 'CC', 'CC', 'M')
Vi = tMi.r(Vi, 'CC', 'CC', 'M')
# End Interpolation
out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, np.sqrt(U**2+V**2).T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
out += (ax.streamplot(tMi.vectorCCx, tMi.vectorCCy, Ui.T, Vi.T, **streamOpts),)
if grid:
xXGrid = np.c_[tM.vectorNx,tM.vectorNx,np.nan*np.ones(tM.nNx)].flatten()
@@ -322,6 +363,8 @@ class TensorView(object):
ax.set_xlabel('y' if normal == 'X' else 'x')
ax.set_ylabel('y' if normal == 'Z' else 'z')
ax.set_title('Slice %d' % ind)
ax.set_xlim(*tM.vectorNx[[0,-1]])
ax.set_ylim(*tM.vectorNy[[0,-1]])
if showIt: plt.show()
return out
@@ -514,3 +557,13 @@ class TensorView(object):
return animate(fig, animateFrame, frames=len(frames))
if __name__ == '__main__':
from SimPEG import *
mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2)))
M = Mesh.TensorMesh(mT)
q = np.zeros(M.vnC)
q[[4,4],[4,4],[2,6]]=[-1,1]
q = Utils.mkvc(q)
A = M.faceDiv*M.cellGrad
b = Solver(A).solve(q)
M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
SimPEG/Model.py
+98 -1
View File
@@ -174,12 +174,109 @@ class Vertical1DModel(BaseModel):
), shape=(repNum, 1))
return sp.kron(sp.identity(self.nP), repVec)
class Mesh2Mesh(BaseModel):
"""
Takes a model on one mesh are translates it to another mesh.
.. plot::
from SimPEG import *
M = Mesh.TensorMesh([100,100])
h1 = Utils.meshTensors(((7,6,1.5),(10,6),(7,6,1.5)))
h1 = h1/h1.sum()
M2 = Mesh.TensorMesh([h1,h1])
V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
v = Utils.mkvc(V)
modh = Model.Mesh2Mesh([M,M2])
modH = Model.Mesh2Mesh([M2,M])
H = modH.transform(v)
h = modh.transform(H)
ax = plt.subplot(131)
M.plotImage(v, ax=ax)
ax.set_title('Fine Mesh (Original)')
ax = plt.subplot(132)
M2.plotImage(H,clim=[0,1],ax=ax)
ax.set_title('Course Mesh')
ax = plt.subplot(133)
M.plotImage(h,clim=[0,1],ax=ax)
ax.set_title('Fine Mesh (Interpolated)')
"""
def __init__(self, meshes, **kwargs):
Utils.setKwargs(self, **kwargs)
assert type(meshes) is list, "meshes must be a list of two meshes"
assert len(meshes) == 2, "meshes must be a list of two meshes"
assert meshes[0].dim == meshes[1].dim, """The two meshes must be the same dimension"""
self.mesh = meshes[0]
self.mesh2 = meshes[1]
self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC,'CC',zerosOutside=True)
@property
def nP(self):
"""Number of parameters in the model."""
return self.mesh2.nC
def transform(self, m):
return self.P*m
def transformDeriv(self, m):
return self.P
class ActiveModel(BaseModel):
"""
Active model parameters.
"""
indActive = None #: Active Cells
valInactive = None #: Values of inactive Cells
nC = None #: Number of cells in the full model
def __init__(self, mesh, indActive, valInactive, nC=None):
self.mesh = mesh
self.nC = nC or mesh.nC
if indActive.dtype is not bool:
z = np.zeros(self.nC,dtype=bool)
z[indActive] = True
indActive = z
self.indActive = indActive
self.indInactive = np.logical_not(indActive)
if type(valInactive) in [float, int, long]:
valInactive = np.ones(self.nC)*float(valInactive)
valInactive[self.indActive] = 0
self.valInactive = valInactive
inds = np.nonzero(self.indActive)[0]
self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
@property
def nP(self):
"""Number of parameters in the model."""
return self.indActive.sum()
def transform(self, m):
return self.P*m + self.valInactive
def transformDeriv(self, m):
return self.P
class ComboModel(BaseModel):
"""Combination of various models."""
def __init__(self, mesh, models, **kwargs):
BaseModel.__init__(self, mesh, **kwargs)
self.models = [m(mesh, **kwargs) for m in models]
self.models = []
for m in models:
if not isinstance(m, BaseModel):
self.models += [m(mesh, **kwargs)]
else:
self.models += [m]
@property
def nP(self):
SimPEG/Regularization.py
+2 -2
View File
@@ -59,7 +59,7 @@ class BaseRegularization(object):
@Utils.timeIt
def modelObj(self, m):
r = self.W * (m - self.mref)
r = self.W * self.model.transform(m - self.mref)
return 0.5*r.dot(r)
@Utils.timeIt
@@ -79,7 +79,7 @@ class BaseRegularization(object):
R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
"""
return self.W.T * ( self.W * (m - self.mref) )
return self.W.T * ( self.W * self.model.transform(m - self.mref) )
@Utils.timeIt
def modelObj2Deriv(self):
SimPEG/Solver.py
+1 -2
View File
@@ -1,8 +1,7 @@
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as linalg
from Utils.matutils import mkvc
from Utils.sputils import sdiag
from Utils.matutils import mkvc, sdiag
import warnings
DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'triangular':'fortran', 'diagonal':'python'}
SimPEG/Tests/test_interpolation.py
+16
View File
@@ -2,6 +2,9 @@ import numpy as np
import unittest
from TestUtils import OrderTest
from SimPEG.Utils import mkvc
from SimPEG import Mesh
import unittest
MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh']
TOLERANCES = [0.9, 0.5]
@@ -50,6 +53,19 @@ class TestInterpolation1D(OrderTest):
self.name = 'Interpolation 1D: N'
self.orderTest()
class TestOutliersInterp1D(unittest.TestCase):
def setUp(self):
pass
def test_outliers(self):
M = Mesh.TensorMesh([4])
Q = M.getInterpolationMat(np.array([[0],[0.126],[0.127]]),'CC',zerosOutside=True)
x = np.arange(4)+1
self.assertTrue(np.all(Q*x == [1,1.004,1.008]))
Q = M.getInterpolationMat(np.array([[-1],[0.126],[0.127]]),'CC',zerosOutside=True)
self.assertTrue(np.all(Q*x == [0,1.004,1.008]))
class TestInterpolation2d(OrderTest):
name = "Interpolation 2D"
LOCS = np.random.rand(50,2)*0.6+0.2
SimPEG/Tests/test_model.py
+6 -1
View File
@@ -11,7 +11,8 @@ class ModelTests(unittest.TestCase):
a = np.array([1, 1, 1])
b = np.array([1, 2])
self.mesh2 = Mesh.TensorMesh([a, b], np.array([3, 5]))
self.mesh2 = Mesh.TensorMesh([a, b], x0=np.array([3, 5]))
self.mesh22 = Mesh.TensorMesh([b, a], x0=np.array([3, 5]))
def test_modelTransforms(self):
for M in dir(Model):
@@ -22,6 +23,10 @@ class ModelTests(unittest.TestCase):
continue
self.assertTrue(model.test())
def test_Mesh2MeshModel(self):
model = Model.Mesh2Mesh([self.mesh22, self.mesh2])
self.assertTrue(model.test())
def test_comboModels(self):
combos = [(Model.LogModel, Model.Vertical1DModel)]
for combo in combos:
SimPEG/Tests/test_utils.py
+59 -4
View File
@@ -1,6 +1,6 @@
import numpy as np
import unittest
from SimPEG.Utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
from SimPEG.Utils import *
from SimPEG import Mesh, np, sp
from SimPEG.Tests import checkDerivative
@@ -64,6 +64,19 @@ class TestSequenceFunctions(unittest.TestCase):
self.assertTrue(np.all(XYZ[:, 1] == X2_test))
self.assertTrue(np.all(XYZ[:, 2] == X3_test))
def test_sub2ind(self):
x = np.ones((5,2))
self.assertTrue(np.all(sub2ind(x.shape, [0,0]) == [0]))
self.assertTrue(np.all(sub2ind(x.shape, [4,0]) == [4]))
self.assertTrue(np.all(sub2ind(x.shape, [0,1]) == [5]))
self.assertTrue(np.all(sub2ind(x.shape, [4,1]) == [9]))
self.assertTrue(np.all(sub2ind(x.shape, [[0,0],[4,0],[0,1],[4,1]]) == [0,4,5,9]))
def test_ind2sub(self):
x = np.ones((5,2))
self.assertTrue(np.all(ind2sub(x.shape, [0,4,5,9])[0] == [0,4,0,4]))
self.assertTrue(np.all(ind2sub(x.shape, [0,4,5,9])[1] == [0,0,1,1]))
def test_indexCube_2D(self):
nN = np.array([3, 3])
self.assertTrue(np.all(indexCube('A', nN) == np.array([0, 1, 3, 4])))
@@ -83,8 +96,6 @@ class TestSequenceFunctions(unittest.TestCase):
self.assertTrue(np.all(indexCube('H', nN) == np.array([10, 11, 13, 14, 19, 20, 22, 23])))
def test_invXXXBlockDiagonal(self):
import scipy.sparse as sp
a = [np.random.rand(5, 1) for i in range(4)]
B = inv2X2BlockDiagonal(*a)
@@ -107,6 +118,50 @@ class TestSequenceFunctions(unittest.TestCase):
self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < 1e-12)
def test_invPropertyTensor2D(self):
M = Mesh.TensorMesh([6, 6])
a1 = np.random.rand(M.nC)
a2 = np.random.rand(M.nC)
a3 = np.random.rand(M.nC)
prop1 = a1
prop2 = np.c_[a1, a2]
prop3 = np.c_[a1, a2, a3]
for prop in [4, prop1, prop2, prop3]:
b = invPropertyTensor(M, prop)
A = makePropertyTensor(M, prop)
B1 = makePropertyTensor(M, b)
B2 = invPropertyTensor(M, prop, returnMatrix=True)
Z = B1*A - sp.identity(M.nC*2)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12)
Z = B2*A - sp.identity(M.nC*2)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12)
def test_invPropertyTensor3D(self):
M = Mesh.TensorMesh([6, 6, 6])
a1 = np.random.rand(M.nC)
a2 = np.random.rand(M.nC)
a3 = np.random.rand(M.nC)
a4 = np.random.rand(M.nC)
a5 = np.random.rand(M.nC)
a6 = np.random.rand(M.nC)
prop1 = a1
prop2 = np.c_[a1, a2, a3]
prop3 = np.c_[a1, a2, a3, a4, a5, a6]
for prop in [4, prop1, prop2, prop3]:
b = invPropertyTensor(M, prop)
A = makePropertyTensor(M, prop)
B1 = makePropertyTensor(M, b)
B2 = invPropertyTensor(M, prop, returnMatrix=True)
Z = B1*A - sp.identity(M.nC*3)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12)
Z = B2*A - sp.identity(M.nC*3)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12)
if __name__ == '__main__':
unittest.main()
SimPEG/Utils/__init__.py
+2 -3
View File
@@ -1,7 +1,6 @@
from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
from sputils import spzeros, kron3, speye, sdiag, sdInv, ddx, av, avExtrap
from matutils import *
from meshutils import exampleLomGird, meshTensors
from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube
from lomutils import volTetra, faceInfo, indexCube
from interputils import interpmat
from ipythonutils import easyAnimate as animate
import ModelBuilder
SimPEG/Utils/interputils.py
+37 -29
View File
@@ -1,7 +1,6 @@
import numpy as np
import scipy.sparse as sp
from sputils import spzeros
from matutils import mkvc, sub2ind
from matutils import mkvc, sub2ind, spzeros
def _interp_point_1D(x, xr_i):
"""
@@ -20,9 +19,17 @@ def _interp_point_1D(x, xr_i):
elif xr_i - x[im] < 0: # Point on the right
ind_x1 = im-1
ind_x2 = im
ind_x1 = max(min(ind_x1, x.size-1), 0)
ind_x2 = max(min(ind_x2, x.size-1), 0)
dx1 = xr_i - x[ind_x1]
dx2 = x[ind_x2] - xr_i
return ind_x1, ind_x2, dx1, dx2
Dx = x[ind_x2] - x[ind_x1]
if ind_x1 == ind_x2:
dx1 = 0.5
dx2 = 0.5
Dx = 1
return ind_x1, ind_x2, dx1, dx2, Dx
def interpmat(locs, x, y=None, z=None):
@@ -70,14 +77,17 @@ def _interpmat1D(locs, x):
Q = sp.lil_matrix((npts, nx))
for i in range(npts):
ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i])
dv = (x[ind_x2] - x[ind_x1])
Dx = x[ind_x2] - x[ind_x1]
ind_x1, ind_x2, dx1, dx2, Dx = _interp_point_1D(x, locs[i])
# dv = (x[ind_x2] - x[ind_x1])
# Get the row in the matrix
inds = [ind_x1, ind_x2]
vals = [(1-dx1/Dx),(1-dx2/Dx)]
Q[i, inds] = vals
return Q.tocsr()
for I, V in zip(inds, vals):
Q[i, I] += V
# Q[i, mkvc(inds)] = vals
return Q
@@ -91,13 +101,10 @@ def _interpmat2D(locs, x, y):
for i in range(npts):
ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, dy1, dy2 = _interp_point_1D(y, locs[i, 1])
ind_x1, ind_x2, dx1, dx2, Dx = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, dy1, dy2, Dy = _interp_point_1D(y, locs[i, 1])
dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1])
Dx = x[ind_x2] - x[ind_x1]
Dy = y[ind_y2] - y[ind_y1]
# dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1])
# Get the row in the matrix
@@ -112,9 +119,13 @@ def _interpmat2D(locs, x, y):
(1-dx2/Dx)*(1-dy1/Dy),
(1-dx2/Dx)*(1-dy2/Dy)]
Q[i, mkvc(inds)] = vals
return Q.tocsr()
for I, V in zip(mkvc(inds), vals):
Q[i, I] += V
# Q[i, mkvc(inds)] = vals
return Q
@@ -129,15 +140,11 @@ def _interpmat3D(locs, x, y, z):
for i in range(npts):
ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, dy1, dy2 = _interp_point_1D(y, locs[i, 1])
ind_z1, ind_z2, dz1, dz2 = _interp_point_1D(z, locs[i, 2])
ind_x1, ind_x2, dx1, dx2, Dx = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, dy1, dy2, Dy = _interp_point_1D(y, locs[i, 1])
ind_z1, ind_z2, dz1, dz2, Dz = _interp_point_1D(z, locs[i, 2])
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]
# dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1])
# Get the row in the matrix
@@ -160,20 +167,21 @@ def _interpmat3D(locs, x, y, z):
(1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz),
(1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz)]
Q[i, mkvc(inds)] = vals
for I, V in zip(mkvc(inds), vals):
Q[i, I] += V
# Q[i, mkvc(inds)] = vals
return Q.tocsr()
return Q
if __name__ == '__main__':
import SimPEG
import numpy as np
from SimPEG import *
import matplotlib.pyplot as plt
locs = np.random.rand(50)*0.8+0.1
x = np.linspace(0,1,7)
dense = np.linspace(0,1,200)
fun = lambda x: np.cos(2*np.pi*x)
Q = SimPEG.Utils.interpmat(locs, x)
Q = Utils.interpmat(locs, x)
plt.plot(x, fun(x), 'bs-')
plt.plot(dense, fun(dense), 'y:')
plt.plot(locs, Q*fun(x), 'mo')
SimPEG/Utils/lomutils.py
+1 -77
View File
@@ -1,7 +1,6 @@
import numpy as np
from scipy import sparse as sp
from matutils import mkvc, ndgrid, sub2ind
from sputils import sdiag
from matutils import mkvc, ndgrid, sub2ind, sdiag
def volTetra(xyz, A, B, C, D):
@@ -188,78 +187,3 @@ def faceInfo(xyz, A, B, C, D, average=True, normalizeNormals=True):
return N, area
def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33):
""" B = inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33)
inverts a stack of 3x3 matrices
Input:
A - a11, a12, a13, a21, a22, a23, a31, a32, a33
Output:
B - inverse
"""
a11 = mkvc(a11)
a12 = mkvc(a12)
a13 = mkvc(a13)
a21 = mkvc(a21)
a22 = mkvc(a22)
a23 = mkvc(a23)
a31 = mkvc(a31)
a32 = mkvc(a32)
a33 = mkvc(a33)
detA = a31*a12*a23 - a31*a13*a22 - a21*a12*a33 + a21*a13*a32 + a11*a22*a33 - a11*a23*a32
b11 = +(a22*a33 - a23*a32)/detA
b12 = -(a12*a33 - a13*a32)/detA
b13 = +(a12*a23 - a13*a22)/detA
b21 = +(a31*a23 - a21*a33)/detA
b22 = -(a31*a13 - a11*a33)/detA
b23 = +(a21*a13 - a11*a23)/detA
b31 = -(a31*a22 - a21*a32)/detA
b32 = +(a31*a12 - a11*a32)/detA
b33 = -(a21*a12 - a11*a22)/detA
B = sp.vstack((sp.hstack((sdiag(b11), sdiag(b12), sdiag(b13))),
sp.hstack((sdiag(b21), sdiag(b22), sdiag(b23))),
sp.hstack((sdiag(b31), sdiag(b32), sdiag(b33)))))
return B
def inv2X2BlockDiagonal(a11, a12, a21, a22):
""" B = inv2X2BlockDiagonal(a11, a12, a21, a22)
Inverts a stack of 2x2 matrices by using the inversion formula
inv(A) = (1/det(A)) * cof(A)^T
Input:
A - a11, a12, a13, a21, a22, a23, a31, a32, a33
Output:
B - inverse
"""
a11 = mkvc(a11)
a12 = mkvc(a12)
a21 = mkvc(a21)
a22 = mkvc(a22)
# compute inverse of the determinant.
detAinv = 1./(a11*a22 - a21*a12)
b11 = +detAinv*a22
b12 = -detAinv*a12
b21 = -detAinv*a21
b22 = +detAinv*a11
B = sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))),
sp.hstack((sdiag(b21), sdiag(b22)))))
return B
SimPEG/Utils/matutils.py
+206 -24
View File
@@ -1,5 +1,5 @@
import numpy as np
import scipy.sparse as sp
def mkvc(x, numDims=1):
"""Creates a vector with the number of dimension specified
@@ -30,6 +30,42 @@ def mkvc(x, numDims=1):
elif numDims == 3:
return x.flatten(order='F')[:, np.newaxis, np.newaxis]
def sdiag(h):
"""Sparse diagonal matrix"""
return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr")
def sdInv(M):
"Inverse of a sparse diagonal matrix"
return sdiag(1/M.diagonal())
def speye(n):
"""Sparse identity"""
return sp.identity(n, format="csr")
def kron3(A, B, C):
"""Three kron prods"""
return sp.kron(sp.kron(A, B), C, format="csr")
def spzeros(n1, n2):
"""spzeros"""
return sp.coo_matrix((n1, n2)).tocsr()
def ddx(n):
"""Define 1D derivatives, inner, this means we go from n+1 to n"""
return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr")
def av(n):
"""Define 1D averaging operator from nodes to cell-centers."""
return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
def avExtrap(n):
"""Define 1D averaging operator from cell-centers to nodes."""
Av = sp.spdiags((0.5*np.ones((n, 1))*[1, 1]).T, [-1, 0], n+1, n, format="csr") + sp.csr_matrix(([0.5,0.5],([0,n],[0,n-1])),shape=(n+1,n))
return Av
def ndgrid(*args, **kwargs):
"""
@@ -98,33 +134,23 @@ def ndgrid(*args, **kwargs):
return XYZ[2], XYZ[1], XYZ[0]
def ind2sub(shape, ind):
"""From the given shape, returns the subscrips of the given index"""
revshp = []
revshp.extend(shape)
mult = [1]
for i in range(0, len(revshp)-1):
mult.extend([mult[i]*revshp[i]])
mult = np.array(mult).reshape(len(mult))
sub = []
for i in range(0, len(shape)):
sub.extend([np.math.floor(ind / mult[i])])
ind = ind - (np.math.floor(ind/mult[i]) * mult[i])
return sub
def ind2sub(shape, inds):
"""From the given shape, returns the subscripts of the given index"""
if type(inds) is not np.ndarray:
inds = np.array(inds)
assert len(inds.shape) == 1, 'Indexing must be done as a 1D row vector, e.g. [3,6,6,...]'
return np.unravel_index(inds, shape, order='F')
def sub2ind(shape, subs):
"""From the given shape, returns the index of the given subscript"""
revshp = list(shape)
mult = [1]
for i in range(0, len(revshp)-1):
mult.extend([mult[i]*revshp[i]])
mult = np.array(mult).reshape(len(mult), 1)
idx = np.dot((subs), (mult))
return idx
if type(subs) is not np.ndarray:
subs = np.array(subs)
if subs.size == len(shape):
subs = subs[np.newaxis,:]
assert subs.shape[1] == len(shape), 'Indexing must be done as a column vectors. e.g. [[3,6],[6,2],...]'
inds = np.ravel_multi_index(subs.T, shape, order='F')
return mkvc(inds)
def getSubArray(A, ind):
@@ -138,3 +164,159 @@ def getSubArray(A, ind):
return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]]
else:
raise Exception("getSubArray does not support dimension asked.")
def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33, returnMatrix=True):
""" B = inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33)
inverts a stack of 3x3 matrices
Input:
A - a11, a12, a13, a21, a22, a23, a31, a32, a33
Output:
B - inverse
"""
a11 = mkvc(a11)
a12 = mkvc(a12)
a13 = mkvc(a13)
a21 = mkvc(a21)
a22 = mkvc(a22)
a23 = mkvc(a23)
a31 = mkvc(a31)
a32 = mkvc(a32)
a33 = mkvc(a33)
detA = a31*a12*a23 - a31*a13*a22 - a21*a12*a33 + a21*a13*a32 + a11*a22*a33 - a11*a23*a32
b11 = +(a22*a33 - a23*a32)/detA
b12 = -(a12*a33 - a13*a32)/detA
b13 = +(a12*a23 - a13*a22)/detA
b21 = +(a31*a23 - a21*a33)/detA
b22 = -(a31*a13 - a11*a33)/detA
b23 = +(a21*a13 - a11*a23)/detA
b31 = -(a31*a22 - a21*a32)/detA
b32 = +(a31*a12 - a11*a32)/detA
b33 = -(a21*a12 - a11*a22)/detA
if not returnMatrix:
return b11, b12, b13, b21, b22, b23, b31, b32, b33
return sp.vstack((sp.hstack((sdiag(b11), sdiag(b12), sdiag(b13))),
sp.hstack((sdiag(b21), sdiag(b22), sdiag(b23))),
sp.hstack((sdiag(b31), sdiag(b32), sdiag(b33)))))
def inv2X2BlockDiagonal(a11, a12, a21, a22, returnMatrix=True):
""" B = inv2X2BlockDiagonal(a11, a12, a21, a22)
Inverts a stack of 2x2 matrices by using the inversion formula
inv(A) = (1/det(A)) * cof(A)^T
Input:
A - a11, a12, a21, a22
Output:
B - inverse
"""
a11 = mkvc(a11)
a12 = mkvc(a12)
a21 = mkvc(a21)
a22 = mkvc(a22)
# compute inverse of the determinant.
detAinv = 1./(a11*a22 - a21*a12)
b11 = +detAinv*a22
b12 = -detAinv*a12
b21 = -detAinv*a21
b22 = +detAinv*a11
if not returnMatrix:
return b11, b12, b21, b22
return sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))),
sp.hstack((sdiag(b21), sdiag(b22)))))
def makePropertyTensor(M, sigma):
if sigma is None: # default is ones
sigma = np.ones(M.nC)
if type(sigma) in [float, int, long]:
sigma = sigma * np.ones(M.nC)
if M.dim == 1:
if sigma.size == M.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
Sigma = sdiag(sigma)
else:
raise Exception('Unexpected shape of sigma')
elif M.dim == 2:
if sigma.size == M.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
Sigma = sdiag(np.r_[sigma, sigma])
elif sigma.shape[1] == 2: # Diagonal tensor
Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]])
elif sigma.shape[1] == 3: # Fully anisotropic
row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2])))
row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1])))
Sigma = sp.vstack((row1, row2))
else:
raise Exception('Unexpected shape of sigma')
elif M.dim == 3:
if sigma.size == M.nC: # Isotropic!
sigma = mkvc(sigma) # ensure it is a vector.
Sigma = sdiag(np.r_[sigma, sigma, sigma])
elif sigma.shape[1] == 3: # Diagonal tensor
Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]])
elif sigma.shape[1] == 6: # Fully anisotropic
row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4])))
row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5])))
row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2])))
Sigma = sp.vstack((row1, row2, row3))
else:
raise Exception('Unexpected shape of sigma')
return Sigma
def invPropertyTensor(M, tensor, returnMatrix=False):
T = None
if type(tensor) in [float, int, long]:
T = 1./tensor
elif tensor.size == M.nC: # Isotropic!
T = 1./mkvc(tensor) # ensure it is a vector.
elif M.dim == 2:
if tensor.shape[1] == 2: # Diagonal tensor
T = 1./tensor
elif tensor.shape[1] == 3: # Fully anisotropic
B = inv2X2BlockDiagonal(tensor[:,0], tensor[:,2],
tensor[:,2], tensor[:,1],
returnMatrix=False)
b11, b12, b21, b22 = B
T = np.c_[b11, b22, b12]
elif M.dim == 3:
if tensor.shape[1] == 3: # Diagonal tensor
T = 1./tensor
elif tensor.shape[1] == 6: # Fully anisotropic
B = inv3X3BlockDiagonal(tensor[:,0], tensor[:,3], tensor[:,4],
tensor[:,3], tensor[:,1], tensor[:,5],
tensor[:,4], tensor[:,5], tensor[:,2],
returnMatrix=False)
b11, b12, b13, b21, b22, b23, b31, b32, b33 = B
T = np.c_[b11, b22, b33, b12, b13, b23]
if T is None:
raise Exception('Unexpected shape of tensor')
if returnMatrix:
return makePropertyTensor(M, T)
return T
SimPEG/Utils/meshutils.py
+1 -2
View File
@@ -1,7 +1,6 @@
import numpy as np
from scipy import sparse as sp
from matutils import mkvc, ndgrid, sub2ind
from sputils import sdiag
from matutils import mkvc, ndgrid, sub2ind, sdiag
def exampleLomGird(nC, exType):
assert type(nC) == list, "nC must be a list containing the number of nodes"
SimPEG/Utils/sputils.py
-41
View File
@@ -1,41 +0,0 @@
from scipy import sparse as sp
from matutils import mkvc
import numpy as np
def sdiag(h):
"""Sparse diagonal matrix"""
return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr")
def sdInv(M):
"Inverse of a sparse diagonal matrix"
return sdiag(1/M.diagonal())
def speye(n):
"""Sparse identity"""
return sp.identity(n, format="csr")
def kron3(A, B, C):
"""Three kron prods"""
return sp.kron(sp.kron(A, B), C, format="csr")
def spzeros(n1, n2):
"""spzeros"""
return sp.coo_matrix((n1, n2)).tocsr()
def ddx(n):
"""Define 1D derivatives, inner, this means we go from n+1 to n"""
return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr")
def av(n):
"""Define 1D averaging operator from nodes to cell-centers."""
return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
def avExtrap(n):
"""Define 1D averaging operator from cell-centers to nodes."""
Av = sp.spdiags((0.5*np.ones((n, 1))*[1, 1]).T, [-1, 0], n+1, n, format="csr") + sp.csr_matrix(([0.5,0.5],([0,n],[0,n-1])),shape=(n+1,n))
return Av