From 1e5053085d269e2c9bfe8a10625627943b76f02d Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Fri, 28 Feb 2014 15:19:57 -0800
Subject: [PATCH] changed matutils to handle vector material properties not
 just columns

---
 SimPEG/Tests/test_utils.py |  1 +
 SimPEG/Utils/matutils.py   | 31 ++++++++++++++++++++-----------
 2 files changed, 21 insertions(+), 11 deletions(-)

diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py
index 429c797c..e69c8473 100644
--- a/SimPEG/Tests/test_utils.py
+++ b/SimPEG/Tests/test_utils.py
@@ -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):
diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py
index 9e350957..c15f6234 100644
--- a/SimPEG/Utils/matutils.py
+++ b/SimPEG/Utils/matutils.py
@@ -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