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

changed matutils to handle vector material properties not just columns

Browse Source
This commit is contained in:
rowanc1
2014-02-28 15:19:57 -08:00
parent e08467af3e
commit 1e5053085d
2 changed files with 21 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Tests/test_utils.py
+1
View File
@@ -71,6 +71,7 @@ class TestSequenceFunctions(unittest.TestCase):
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, [[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):
SimPEG/Utils/matutils.py
+20 -11
View File
@@ -146,7 +146,7 @@ def sub2ind(shape, subs):
"""From the given shape, returns the index of the given subscript"""
if type(subs) is not np.ndarray:
subs = np.array(subs)
if subs.size == len(shape):
if len(subs.shape) == 1:
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')
@@ -261,9 +261,11 @@ def makePropertyTensor(M, sigma):
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
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.shape[1] == 3: # Fully anisotropic
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))
@@ -273,15 +275,18 @@ def makePropertyTensor(M, sigma):
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
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.shape[1] == 6: # Fully anisotropic
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')
return Sigma
@@ -296,27 +301,31 @@ def invPropertyTensor(M, tensor, returnMatrix=False):
T = 1./mkvc(tensor) # ensure it is a vector.
elif M.dim == 2:
if tensor.shape[1] == 2: # Diagonal tensor
if tensor.size == M.nC*2: # Diagonal tensor
T = 1./tensor
elif tensor.shape[1] == 3: # Fully anisotropic
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.c_[b11, b22, b12]
T = np.r_[b11, b22, b12]
elif M.dim == 3:
if tensor.shape[1] == 3: # Diagonal tensor
if tensor.size == M.nC*3: # Diagonal tensor
T = 1./tensor
elif tensor.shape[1] == 6: # Fully anisotropic
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.c_[b11, b22, b33, b12, b13, b23]
T = np.r_[b11, b22, b33, b12, b13, b23]
if T is None:
raise Exception('Unexpected shape of tensor')
if returnMatrix:
return makePropertyTensor(M, T)
return T