diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index ea1089ee..a7987cd6 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -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,46 +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)) - elif type(sigma) is float: - sigma = np.ones(self.nC)*sigma - - 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' return _getFacePx_Rectangular(M) diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index 60ca5558..6c47be27 100644 --- a/SimPEG/Solver.py +++ b/SimPEG/Solver.py @@ -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'} diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py index 2732c311..7e2660c2 100644 --- a/SimPEG/Tests/test_utils.py +++ b/SimPEG/Tests/test_utils.py @@ -1,6 +1,6 @@ -import numpy as np import unittest -from SimPEG.Utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal,sub2ind,ind2sub +from SimPEG.Utils import * +from SimPEG import Mesh, np, sp from SimPEG.Tests import checkDerivative @@ -96,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) @@ -120,6 +118,50 @@ class TestSequenceFunctions(unittest.TestCase): self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < 1e-10) + 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() diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index b2d0d56a..f45a1039 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -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 diff --git a/SimPEG/Utils/interputils.py b/SimPEG/Utils/interputils.py index 3f9b7e71..f1dc28b9 100644 --- a/SimPEG/Utils/interputils.py +++ b/SimPEG/Utils/interputils.py @@ -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): """ diff --git a/SimPEG/Utils/lomutils.py b/SimPEG/Utils/lomutils.py index 653cfdde..b641ef7f 100644 --- a/SimPEG/Utils/lomutils.py +++ b/SimPEG/Utils/lomutils.py @@ -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 diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index b2e1daf2..9e350957 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -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): """ @@ -97,6 +133,7 @@ def ndgrid(*args, **kwargs): else: return XYZ[2], XYZ[1], XYZ[0] + def ind2sub(shape, inds): """From the given shape, returns the subscripts of the given index""" if type(inds) is not np.ndarray: @@ -104,6 +141,7 @@ def ind2sub(shape, 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""" if type(subs) is not np.ndarray: @@ -114,6 +152,7 @@ def sub2ind(shape, subs): inds = np.ravel_multi_index(subs.T, shape, order='F') return mkvc(inds) + def getSubArray(A, ind): """subArray""" assert type(ind) == list, "ind must be a list of vectors" @@ -125,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 diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index d7fe0836..f8add82b 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -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" diff --git a/SimPEG/Utils/sputils.py b/SimPEG/Utils/sputils.py deleted file mode 100644 index f671be43..00000000 --- a/SimPEG/Utils/sputils.py +++ /dev/null @@ -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