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

invMat option in getInnerProducts

Browse Source
This commit is contained in:
rowanc1
2014-04-25 17:21:31 -07:00
parent c3333225ec
commit d08da1a3b0
6 changed files with 200 additions and 83 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/InnerProducts.py
+19 -5
View File
@@ -1,5 +1,5 @@
from scipy import sparse as sp
from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor, invPropertyTensor, spzeros, isScalar
from SimPEG.Utils import *
import numpy as np
@@ -11,11 +11,12 @@ class InnerProducts(object):
raise Exception('InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class.')
def getFaceInnerProduct(self, prop=None, returnP=False,
invProp=False, doFast=True):
invProp=False, invMat=False, doFast=True):
"""
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
:param bool returnP: returns the projection matrices
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:param bool doFast: do a faster implementation if available.
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (nF, nF)
@@ -23,7 +24,7 @@ class InnerProducts(object):
fast = None
if returnP is False and hasattr(self, '_fastFaceInnerProduct') and doFast:
fast = self._fastFaceInnerProduct(prop=prop, invProp=invProp)
fast = self._fastFaceInnerProduct(prop=prop, invProp=invProp, invMat=invMat)
if fast is not None:
return fast
@@ -67,6 +68,12 @@ class InnerProducts(object):
if d > 2:
A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111
P += [P001, P101, P011, P111]
if invMat and tensorType(self, prop) < 3:
A = sdInv(A)
elif invMat and tensorType(self, prop) == 3:
raise Exception('Solver needed to invert A.')
if returnP:
return A, P
else:
@@ -95,11 +102,12 @@ class InnerProducts(object):
return self._getInnerProductDeriv(prop, v, P, self.nF)
def getEdgeInnerProduct(self, prop=None, returnP=False,
invProp=False, doFast=True):
invProp=False, invMat=False, doFast=True):
"""
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
:param bool returnP: returns the projection matrices
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:param bool doFast: do a faster implementation if available.
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (nE, nE)
@@ -107,7 +115,7 @@ class InnerProducts(object):
fast = None
if returnP is False and hasattr(self, '_fastEdgeInnerProduct') and doFast:
fast = self._fastEdgeInnerProduct(prop=prop, invProp=invProp)
fast = self._fastEdgeInnerProduct(prop=prop, invProp=invProp, invMat=invMat)
if fast is not None:
return fast
@@ -146,6 +154,12 @@ class InnerProducts(object):
if d == 3:
A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111
P += [P001, P101, P011, P111]
if invMat and tensorType(self, prop) < 3:
A = sdInv(A)
elif invMat and tensorType(self, prop) == 3:
raise Exception('Solver needed to invert A.')
if returnP:
return A, P
else:
SimPEG/Mesh/TensorMesh.py
+17 -8
View File
@@ -241,7 +241,7 @@ class BaseTensorMesh(BaseRectangularMesh):
return Q.tocsr()
def _fastFaceInnerProduct(self, prop=None, invProp=False):
def _fastFaceInnerProduct(self, prop=None, invProp=False, invMat=False):
"""
Fast version of getFaceInnerProduct.
This does not handle the case of a full tensor prop.
@@ -249,13 +249,14 @@ class BaseTensorMesh(BaseRectangularMesh):
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
:param bool returnP: returns the projection matrices
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (nF, nF)
"""
return self._fastInnerProduct('F', prop=prop, invProp=invProp)
return self._fastInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat)
def _fastEdgeInnerProduct(self, prop=None, invProp=False):
def _fastEdgeInnerProduct(self, prop=None, invProp=False, invMat=False):
"""
Fast version of getEdgeInnerProduct.
This does not handle the case of a full tensor prop.
@@ -263,13 +264,14 @@ class BaseTensorMesh(BaseRectangularMesh):
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
:param bool returnP: returns the projection matrices
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (nE, nE)
"""
return self._fastInnerProduct('E', prop=prop, invProp=invProp)
return self._fastInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat)
def _fastInnerProduct(self, AvType, prop=None, invProp=False):
def _fastInnerProduct(self, AvType, prop=None, invProp=False, invMat=False):
"""
Fast version of getFaceInnerProduct.
This does not handle the case of a full tensor prop.
@@ -278,6 +280,7 @@ class BaseTensorMesh(BaseRectangularMesh):
:param str AvType: 'E' or 'F'
:param bool returnP: returns the projection matrices
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (nF, nF)
"""
@@ -293,12 +296,18 @@ class BaseTensorMesh(BaseRectangularMesh):
if prop.size == self.nC:
Av = getattr(self, 'ave'+AvType+'2CC')
Vprop = self.vol * Utils.mkvc(prop)
return self.dim * Utils.sdiag(Av.T * Vprop)
if prop.size == self.nC*self.dim:
M = self.dim * Utils.sdiag(Av.T * Vprop)
elif prop.size == self.nC*self.dim:
Av = getattr(self, 'ave'+AvType+'2CCV')
V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
return Utils.sdiag(Av.T * V * Utils.mkvc(prop))
M = Utils.sdiag(Av.T * V * Utils.mkvc(prop))
else:
return None
if invMat:
return Utils.sdInv(M)
else:
return M
def _fastFaceInnerProductDeriv(self, prop=None, v=None):
"""
SimPEG/Tests/test_massMatrices.py
+42 -3
View File
@@ -1,6 +1,7 @@
import numpy as np
import unittest
from TestUtils import OrderTest
from SimPEG import Utils
class TestInnerProducts(OrderTest):
@@ -130,14 +131,21 @@ class TestInnerProducts2D(OrderTest):
Ec = np.vstack((cart(self.M.gridEx),
cart(self.M.gridEy)))
E = self.M.projectEdgeVector(Ec)
A = self.M.getEdgeInnerProduct(sigma)
if self.invProp:
A = self.M.getEdgeInnerProduct(Utils.invPropertyTensor(self.M, sigma), invProp=True)
else:
A = self.M.getEdgeInnerProduct(sigma)
numeric = E.T.dot(A.dot(E))
elif self.location == 'faces':
cart = lambda g: np.c_[call(ex, g), call(ey, g)]
Fc = np.vstack((cart(self.M.gridFx),
cart(self.M.gridFy)))
F = self.M.projectFaceVector(Fc)
A = self.M.getFaceInnerProduct(sigma)
if self.invProp:
A = self.M.getFaceInnerProduct(Utils.invPropertyTensor(self.M, sigma), invProp=True)
else:
A = self.M.getFaceInnerProduct(sigma)
numeric = F.T.dot(A.dot(F))
err = np.abs(numeric - analytic)
@@ -147,36 +155,60 @@ class TestInnerProducts2D(OrderTest):
self.name = "2D Edge Inner Product - Isotropic"
self.location = 'edges'
self.sigmaTest = 1
self.invProp = True
self.orderTest()
self.name += " - invProp"
self.invProp = False
self.orderTest()
def test_order3_edges(self):
self.name = "2D Edge Inner Product - Anisotropic"
self.location = 'edges'
self.sigmaTest = 2
self.invProp = True
self.orderTest()
self.name += " - invProp"
self.invProp = False
self.orderTest()
def test_order6_edges(self):
self.name = "2D Edge Inner Product - Full Tensor"
self.location = 'edges'
self.sigmaTest = 3
self.invProp = True
self.orderTest()
self.name += " - invProp"
self.invProp = False
self.orderTest()
def test_order1_faces(self):
self.name = "2D Face Inner Product - Isotropic"
self.location = 'faces'
self.sigmaTest = 1
self.invProp = True
self.orderTest()
self.name += " - invProp"
self.invProp = False
self.orderTest()
def test_order2_faces(self):
self.name = "2D Face Inner Product - Anisotropic"
self.location = 'faces'
self.sigmaTest = 2
self.invProp = True
self.orderTest()
self.name += " - invProp"
self.invProp = False
self.orderTest()
def test_order3_faces(self):
self.name = "2D Face Inner Product - Full Tensor"
self.location = 'faces'
self.sigmaTest = 3
self.invProp = True
self.orderTest()
self.name += " - invProp"
self.invProp = False
self.orderTest()
@@ -205,7 +237,10 @@ class TestInnerProducts1D(OrderTest):
if self.location == 'faces':
F = call(ex, self.M.gridFx)
A = self.M.getFaceInnerProduct(sigma)
if self.invProp:
A = self.M.getFaceInnerProduct(1/sigma, invProp=True)
else:
A = self.M.getFaceInnerProduct(sigma)
numeric = F.T.dot(A.dot(F))
err = np.abs(numeric - analytic)
@@ -215,6 +250,10 @@ class TestInnerProducts1D(OrderTest):
self.name = "1D Face Inner Product"
self.location = 'faces'
self.sigmaTest = 1
self.invProp = True
self.orderTest()
self.name += " - invProp"
self.invProp = False
self.orderTest()
SimPEG/Tests/test_massMatricesDerivs.py
+28
View File
@@ -70,6 +70,22 @@ class TestInnerProductsDerivs(unittest.TestCase):
def test_FaceIP_3D_anisotropic_fast(self):
self.assertTrue(self.doTestFace([10, 4, 5],3, False, True))
def test_FaceIP_1D_float_fast_vec(self):
self.assertTrue(self.doTestFace([10],0, True, True))
def test_FaceIP_2D_float_fast_vec(self):
self.assertTrue(self.doTestFace([10, 4],0, True, True))
def test_FaceIP_3D_float_fast_vec(self):
self.assertTrue(self.doTestFace([10, 4, 5],0, True, True))
def test_FaceIP_1D_isotropic_fast_vec(self):
self.assertTrue(self.doTestFace([10],1, True, True))
def test_FaceIP_2D_isotropic_fast_vec(self):
self.assertTrue(self.doTestFace([10, 4],1, True, True))
def test_FaceIP_3D_isotropic_fast_vec(self):
self.assertTrue(self.doTestFace([10, 4, 5],1, True, True))
def test_FaceIP_2D_anisotropic_fast_vec(self):
self.assertTrue(self.doTestFace([10, 4],2, True, True))
def test_FaceIP_3D_anisotropic_fast_vec(self):
self.assertTrue(self.doTestFace([10, 4, 5],3, True, True))
def test_EdgeIP_2D_float(self):
self.assertTrue(self.doTestEdge([10, 4],0,True, False))
@@ -101,6 +117,18 @@ class TestInnerProductsDerivs(unittest.TestCase):
def test_EdgeIP_3D_anisotropic_fast(self):
self.assertTrue(self.doTestEdge([10, 4, 5],3, False, True))
def test_EdgeIP_2D_float_fast_vec(self):
self.assertTrue(self.doTestEdge([10, 4],0, True, True))
def test_EdgeIP_3D_float_fast_vec(self):
self.assertTrue(self.doTestEdge([10, 4, 5],0, True, True))
def test_EdgeIP_2D_isotropic_fast_vec(self):
self.assertTrue(self.doTestEdge([10, 4],1, True, True))
def test_EdgeIP_3D_isotropic_fast_vec(self):
self.assertTrue(self.doTestEdge([10, 4, 5],1, True, True))
def test_EdgeIP_2D_anisotropic_fast_vec(self):
self.assertTrue(self.doTestEdge([10, 4],2, True, True))
def test_EdgeIP_3D_anisotropic_fast_vec(self):
self.assertTrue(self.doTestEdge([10, 4, 5],3, True, True))
SimPEG/Tests/test_utils.py
+33
View File
@@ -140,6 +140,39 @@ class TestSequenceFunctions(unittest.TestCase):
Z = B2*A - sp.identity(M.nC*2)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
def test_tensorType2D(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 ii, prop in enumerate([4, prop1, prop2, prop3]):
self.assertTrue(tensorType(M, prop) == ii)
self.assertRaises(Exception, tensorType, M, np.c_[a1, a2, a3, a3])
self.assertTrue(tensorType(M, None) == -1)
def test_tensorType3D(self):
M = Mesh.TensorMesh([6, 6, 7])
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 ii, prop in enumerate([4, prop1, prop2, prop3]):
self.assertTrue(tensorType(M, prop) == ii)
self.assertRaises(Exception, tensorType, M, np.c_[a1, a2, a3, a3])
self.assertTrue(tensorType(M, None) == -1)
def test_invPropertyTensor3D(self):
M = Mesh.TensorMesh([6, 6, 6])
SimPEG/Utils/matutils.py
+61 -67
View File
@@ -251,85 +251,79 @@ def inv2X2BlockDiagonal(a11, a12, a21, a22, returnMatrix=True):
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)
def tensorType(M, tensor):
if tensor is None: # default is ones
return -1
if isScalar(sigma):
sigma = sigma * np.ones(M.nC)
if isScalar(tensor):
return 0
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.size == M.nC*2: # Diagonal tensor
sigma = sigma.reshape((M.nC,2), order='F')
Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]])
elif sigma.size == M.nC*3: # Fully anisotropic
sigma = sigma.reshape((M.nC,3), order='F')
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.size == M.nC*3: # Diagonal tensor
sigma = sigma.reshape((M.nC,3), order='F')
Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]])
elif sigma.size == M.nC*6: # Fully anisotropic
sigma = sigma.reshape((M.nC,6), order='F')
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')
if tensor.size == M.nC:
return 1
if ((M.dim == 2 and tensor.size == M.nC*2) or
(M.dim == 3 and tensor.size == M.nC*3)):
return 2
if ((M.dim == 2 and tensor.size == M.nC*3) or
(M.dim == 3 and tensor.size == M.nC*6)):
return 3
raise Exception('Unexpected shape of tensor')
def makePropertyTensor(M, tensor):
if tensor is None: # default is ones
tensor = np.ones(M.nC)
if isScalar(tensor):
tensor = tensor * np.ones(M.nC)
propType = tensorType(M, tensor)
if propType == 1: # Isotropic!
Sigma = sp.kron(sp.identity(M.dim), sdiag(mkvc(tensor)))
elif propType == 2: # Diagonal tensor
Sigma = sdiag(mkvc(tensor))
elif M.dim == 2 and tensor.size == M.nC*3: # Fully anisotropic, 2D
tensor = tensor.reshape((M.nC,3), order='F')
row1 = sp.hstack((sdiag(tensor[:, 0]), sdiag(tensor[:, 2])))
row2 = sp.hstack((sdiag(tensor[:, 2]), sdiag(tensor[:, 1])))
Sigma = sp.vstack((row1, row2))
elif M.dim == 3 and tensor.size == M.nC*6: # Fully anisotropic, 3D
tensor = tensor.reshape((M.nC,6), order='F')
row1 = sp.hstack((sdiag(tensor[:, 0]), sdiag(tensor[:, 3]), sdiag(tensor[:, 4])))
row2 = sp.hstack((sdiag(tensor[:, 3]), sdiag(tensor[:, 1]), sdiag(tensor[:, 5])))
row3 = sp.hstack((sdiag(tensor[:, 4]), sdiag(tensor[:, 5]), sdiag(tensor[:, 2])))
Sigma = sp.vstack((row1, row2, row3))
else:
raise Exception('Unexpected shape of tensor')
return Sigma
def invPropertyTensor(M, tensor, returnMatrix=False):
T = None
propType = tensorType(M, tensor)
if isScalar(tensor):
T = 1./tensor
elif tensor.size == M.nC: # Isotropic!
elif propType < 3: # Isotropic or Diagonal
T = 1./mkvc(tensor) # ensure it is a vector.
elif M.dim == 2:
if tensor.size == M.nC*2: # Diagonal tensor
T = 1./tensor
elif tensor.size == M.nC*3: # Fully anisotropic
tensor = tensor.reshape((M.nC,3), order='F')
B = inv2X2BlockDiagonal(tensor[:,0], tensor[:,2],
tensor[:,2], tensor[:,1],
returnMatrix=False)
b11, b12, b21, b22 = B
T = np.r_[b11, b22, b12]
elif M.dim == 3:
if tensor.size == M.nC*3: # Diagonal tensor
T = 1./tensor
elif tensor.size == M.nC*6: # Fully anisotropic
tensor = tensor.reshape((M.nC,6), order='F')
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.r_[b11, b22, b33, b12, b13, b23]
if T is None:
elif M.dim == 2 and tensor.size == M.nC*3: # Fully anisotropic, 2D
tensor = tensor.reshape((M.nC,3), order='F')
B = inv2X2BlockDiagonal(tensor[:,0], tensor[:,2],
tensor[:,2], tensor[:,1],
returnMatrix=False)
b11, b12, b21, b22 = B
T = np.r_[b11, b22, b12]
elif M.dim == 3 and tensor.size == M.nC*6: # Fully anisotropic, 3D
tensor = tensor.reshape((M.nC,6), order='F')
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.r_[b11, b22, b33, b12, b13, b23]
else:
raise Exception('Unexpected shape of tensor')
if returnMatrix: