mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-12 06:37:46 +08:00
Merge branch 'develop' of https://github.com/simpeg/simpeg into cylClean
Conflicts: SimPEG/Mesh/InnerProducts.py
This commit is contained in:
rowanc1
@@ -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)
|
||||
|
||||
+1
-2
@@ -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'}
|
||||
|
||||
@@ -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()
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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):
|
||||
"""
|
||||
|
||||
@@ -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
|
||||
|
||||
+196
-1
@@ -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
|
||||
|
||||
@@ -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"
|
||||
|
||||
@@ -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
|
||||
Reference in New Issue
Block a user