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

isScalar in Utils (tested) and support for scalar inner products and derivatives (tested).

Browse Source
This commit is contained in:
rowanc1
2014-03-04 12:19:57 -08:00
parent e59a039052
commit 90967ce7f9
6 changed files with 77 additions and 25 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/InnerProducts.py
+18 -10
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
from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor, invPropertyTensor, spzeros, isScalar
import numpy as np
@@ -183,25 +183,33 @@ class InnerProducts(object):
:rtype: scipy.csr_matrix
:return: dMdm, the derivative of the inner product matrix (n, nC*nA)
"""
if materialProperty is None:
return None
if v is None:
raise Exception('v must be supplied for this implementation.')
d = self.dim
Z = spzeros(self.nC, self.nC)
if d == 1:
dMdm = spzeros(n, self.nC)
if isScalar(materialProperty):
dMdm = spzeros(n, 1)
for i, p in enumerate(P):
dMdm = dMdm + p.T * sdiag( p * v )
dMdm = dMdm + sp.csr_matrix((p.T * (p * v), (range(n), np.zeros(n))), shape=(n,1))
if d == 1:
if materialProperty.size == self.nC:
dMdm = spzeros(n, self.nC)
for i, p in enumerate(P):
dMdm = dMdm + p.T * sdiag( p * v )
elif d == 2:
if materialProperty is None or materialProperty.size == self.nC:
if materialProperty.size == self.nC:
dMdm = spzeros(n, self.nC)
for i, p in enumerate(P):
Y = p * v
y1 = Y[:self.nC]
y2 = Y[self.nC:]
dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 )))
if materialProperty.size == self.nC*2:
elif materialProperty.size == self.nC*2:
dMdms = [spzeros(n, self.nC) for _ in range(2)]
for i, p in enumerate(P):
Y = p * v
@@ -210,7 +218,7 @@ class InnerProducts(object):
dMdms[0] = dMdms[0] + p.T * sp.vstack(( sdiag( y1 ), Z))
dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 )))
dMdm = sp.hstack(dMdms)
if materialProperty.size == self.nC*3:
elif materialProperty.size == self.nC*3:
dMdms = [spzeros(n, self.nC) for _ in range(3)]
for i, p in enumerate(P):
Y = p * v
@@ -221,7 +229,7 @@ class InnerProducts(object):
dMdms[2] = dMdms[2] + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 )))
dMdm = sp.hstack(dMdms)
elif d == 3:
if materialProperty is None or materialProperty.size == self.nC:
if materialProperty.size == self.nC:
dMdm = spzeros(n, self.nC)
for i, p in enumerate(P):
Y = p * v
@@ -229,7 +237,7 @@ class InnerProducts(object):
y2 = Y[self.nC:self.nC*2]
y3 = Y[self.nC*2:]
dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 ), sdiag( y3 )))
if materialProperty.size == self.nC*3:
elif materialProperty.size == self.nC*3:
dMdms = [spzeros(n, self.nC) for _ in range(3)]
for i, p in enumerate(P):
Y = p * v
@@ -240,7 +248,7 @@ class InnerProducts(object):
dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 ), Z))
dMdms[2] = dMdms[2] + p.T * sp.vstack(( Z, Z, sdiag( y3 )))
dMdm = sp.hstack(dMdms)
if materialProperty.size == self.nC*6:
elif materialProperty.size == self.nC*6:
dMdms = [spzeros(n, self.nC) for _ in range(6)]
for i, p in enumerate(P):
Y = p * v
SimPEG/Mesh/TensorMesh.py
+19 -10
View File
@@ -534,18 +534,18 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
"""
if materialProperty is None:
materialProperty = np.ones(self.nC)
elif type(materialProperty) in [float, int, long]:
materialProperty = materialProperty * np.ones(M.nC)
if invertProperty:
materialProperty = 1./materialProperty
if Utils.isScalar(materialProperty):
materialProperty = materialProperty*np.ones(self.nC)
if materialProperty.size == self.nC:
if invertProperty:
materialProperty = 1./materialProperty
Av = getattr(self, 'ave'+AvType+'2CC')
V = Utils.sdiag(self.vol)
return self.dim * Utils.sdiag(Av.T * V * materialProperty)
Vprop = self.vol * Utils.mkvc(materialProperty)
return self.dim * Utils.sdiag(Av.T * Vprop)
if materialProperty.size == self.nC*self.dim:
if invertProperty:
materialProperty = 1./materialProperty
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(materialProperty))
@@ -576,11 +576,20 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts):
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (nF, nF)
"""
if materialProperty is None or materialProperty.size == self.nC:
if materialProperty is None:
return None
if Utils.isScalar(materialProperty):
Av = getattr(self, 'ave'+AvType+'2CC')
V = Utils.sdiag(self.vol)
ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1))
if v is None:
return self.dim * Av.T * V * ones
return Utils.sdiag(v) * self.dim * Av.T * V * ones
if materialProperty.size == self.nC:
Av = getattr(self, 'ave'+AvType+'2CC')
V = Utils.sdiag(self.vol)
if v is None:
return self.dim * Av.T * Utils.sdiag(self.vol)
return self.dim * Av.T * V
return Utils.sdiag(v) * self.dim * Av.T * V
if materialProperty.size == self.nC*self.dim: # anisotropic
Av = getattr(self, 'ave'+AvType+'2CCV')
SimPEG/Tests/test_massMatrices.py
+1 -1
View File
@@ -122,7 +122,7 @@ class TestInnerProducts2D(OrderTest):
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)]
analytic = 189959./120 # Found using sympy. z=5
elif self.sigmaTest == 3:
sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
analytic = 781427./360 # Found using sympy. z=5
if self.location == 'edges':
SimPEG/Tests/test_massMatricesDerivs.py
+22 -2
View File
@@ -16,7 +16,7 @@ class TestInnerProductsDerivs(unittest.TestCase):
return M*v, Md
Md = mesh.getFaceInnerProductDeriv(sig, doFast=fast)
return M*v, Utils.sdiag(v)*Md
sig = np.random.rand(mesh.nC*rep)
sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
return checkDerivative(fun, sig, num=5, plotIt=False)
def doTestEdge(self, h, rep, vec, fast):
@@ -29,9 +29,15 @@ class TestInnerProductsDerivs(unittest.TestCase):
return M*v, Md
Md = mesh.getEdgeInnerProductDeriv(sig, doFast=fast)
return M*v, Utils.sdiag(v)*Md
sig = np.random.rand(mesh.nC*rep)
sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
return checkDerivative(fun, sig, num=5, plotIt=False)
def test_FaceIP_1D_float(self):
self.assertTrue(self.doTestFace([10],0,True, False))
def test_FaceIP_2D_float(self):
self.assertTrue(self.doTestFace([10, 4],0,True, False))
def test_FaceIP_3D_float(self):
self.assertTrue(self.doTestFace([10, 4, 5],0,True, False))
def test_FaceIP_1D_isotropic(self):
self.assertTrue(self.doTestFace([10],1,True, False))
def test_FaceIP_2D_isotropic(self):
@@ -47,6 +53,12 @@ class TestInnerProductsDerivs(unittest.TestCase):
def test_FaceIP_3D_tensor(self):
self.assertTrue(self.doTestFace([10, 4, 5],6,True, False))
def test_FaceIP_1D_float_fast(self):
self.assertTrue(self.doTestFace([10],0, False, True))
def test_FaceIP_2D_float_fast(self):
self.assertTrue(self.doTestFace([10, 4],0, False, True))
def test_FaceIP_3D_float_fast(self):
self.assertTrue(self.doTestFace([10, 4, 5],0, False, True))
def test_FaceIP_1D_isotropic_fast(self):
self.assertTrue(self.doTestFace([10],1, False, True))
def test_FaceIP_2D_isotropic_fast(self):
@@ -59,6 +71,10 @@ class TestInnerProductsDerivs(unittest.TestCase):
self.assertTrue(self.doTestFace([10, 4, 5],3, False, True))
def test_EdgeIP_2D_float(self):
self.assertTrue(self.doTestEdge([10, 4],0,True, False))
def test_EdgeIP_3D_float(self):
self.assertTrue(self.doTestEdge([10, 4, 5],0,True, False))
def test_EdgeIP_2D_isotropic(self):
self.assertTrue(self.doTestEdge([10, 4],1,True, False))
def test_EdgeIP_3D_isotropic(self):
@@ -72,6 +88,10 @@ class TestInnerProductsDerivs(unittest.TestCase):
def test_EdgeIP_3D_tensor(self):
self.assertTrue(self.doTestEdge([10, 4, 5],6,True, False))
def test_EdgeIP_2D_float_fast(self):
self.assertTrue(self.doTestEdge([10, 4],0, False, True))
def test_EdgeIP_3D_float_fast(self):
self.assertTrue(self.doTestEdge([10, 4, 5],0, False, True))
def test_EdgeIP_2D_isotropic_fast(self):
self.assertTrue(self.doTestEdge([10, 4],1, False, True))
def test_EdgeIP_3D_isotropic_fast(self):
SimPEG/Tests/test_utils.py
+6
View File
@@ -164,6 +164,12 @@ class TestSequenceFunctions(unittest.TestCase):
Z = B2*A - sp.identity(M.nC*3)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
def test_isFloat(self):
self.assertTrue(isScalar(1.))
self.assertTrue(isScalar(1))
self.assertTrue(isScalar(long(1)))
self.assertTrue(isScalar(np.r_[1.]))
self.assertTrue(isScalar(np.r_[1]))
if __name__ == '__main__':
unittest.main()
SimPEG/Utils/matutils.py
+11 -2
View File
@@ -1,6 +1,15 @@
import numpy as np
import scipy.sparse as sp
def isScalar(f):
scalarTypes = [float, int, long, np.float_, np.int_]
if type(f) in scalarTypes:
return True
elif type(f) == np.ndarray and f.size == 1 and type(f[0]) in scalarTypes:
return True
return False
def mkvc(x, numDims=1):
"""Creates a vector with the number of dimension specified
@@ -248,7 +257,7 @@ def makePropertyTensor(M, sigma):
if sigma is None: # default is ones
sigma = np.ones(M.nC)
if type(sigma) in [float, int, long]:
if isScalar(sigma):
sigma = sigma * np.ones(M.nC)
if M.dim == 1:
@@ -294,7 +303,7 @@ def invPropertyTensor(M, tensor, returnMatrix=False):
T = None
if type(tensor) in [float, int, long]:
if isScalar(tensor):
T = 1./tensor
elif tensor.size == M.nC: # Isotropic!