Edge Inner products in 2D NOTE: still has bugs.

SimPEG/InnerProducts.py

@@ -26,8 +26,10 @@ class InnerProducts(object):
             pass
         elif self._meshType == 'LOM':
             pass  # todo: we should be doing something slightly different here!
-        return getEdgeInnerProduct(self, sigma, returnP)
-
+        if self.dim == 2:
+            return getEdgeInnerProduct2D(self, sigma, returnP)
+        elif self.dim == 3:
+            return getEdgeInnerProduct(self, sigma, returnP)
 
 # ------------------------ Geometries ------------------------------
 #
@@ -275,6 +277,71 @@ def getEdgeInnerProduct(mesh, sigma=None, returnP=False):
         return A
 
 
+def getEdgeInnerProduct2D(mesh, sigma=None, returnP=False):
+
+    if sigma is None:  # default is ones
+        sigma = np.ones((mesh.nC, 1))
+
+    m = np.array([mesh.nCx, mesh.nCy])
+    nc = mesh.nC
+
+    i, j = np.int64(range(m[0])), np.int64(range(m[1]))
+
+    iijj = ndgrid(i, j)
+    ii, jj = iijj[:, 0], iijj[:, 1]
+
+    if mesh._meshType == 'LOM':
+        eT1 = mesh.r(mesh.tangents, 'E', 'Ex', 'M')
+        eT2 = mesh.r(mesh.tangents, 'E', 'Ey', 'M')
+
+    def Pxx(pos):
+        ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1]])
+        ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nE[0]
+
+        IND = np.r_[ind1, ind2].flatten()
+
+        PXX = sp.coo_matrix((np.ones(2*nc), (range(2*nc), IND)), shape=(2*nc, np.sum(mesh.nE))).tocsr()
+
+        if mesh._meshType == 'LOM':
+            I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i + pos[0][0], j + pos[0][1]]), getSubArray(eT1[1], [i + pos[0][0], j + pos[0][1]]),
+                                       getSubArray(eT2[0], [i + pos[1][0], j + pos[1][1]]), getSubArray(eT2[1], [i + pos[1][0], j + pos[1][1]]))
+            PXX = I2x2 * PXX
+
+        return PXX
+
+    # no | node      | e1      | e2
+    # 00 | i  ,j     | i  ,j   | i  ,j
+    # 10 | i+1,j     | i  ,j   | i+1,j
+    # 01 | i  ,j+1   | i  ,j+1 | i  ,j
+    # 11 | i+1,j+1   | i  ,j+1 | i+1,j
+
+    # Square root of cell volume multiplied by 1/4
+    v = np.sqrt(0.25*mesh.vol)
+    V2 = sdiag(np.r_[v, v])  # We will multiply on each side to keep symmetry
+
+    P00 = V2*Pxx([[0, 0], [0, 0]])
+    P10 = V2*Pxx([[0, 0], [1, 0]])
+    P01 = V2*Pxx([[0, 1], [0, 0]])
+    P11 = V2*Pxx([[0, 1], [1, 0]])
+
+    if sigma.size == mesh.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))
+
+    A = P00.T*Sigma*P00 + P10.T*Sigma*P10 + P01.T*Sigma*P01 + P11.T*Sigma*P11
+    P = [P00, P10, P01, P11]
+    if returnP:
+        return A, P
+    else:
+        return A
+
+
 if __name__ == '__main__':
     from TensorMesh import TensorMesh
     h = [np.array([1, 2, 3, 4]), np.array([1, 2, 1, 4, 2]), np.array([1, 1, 4, 1])]

SimPEG/exampleGrid.py

@@ -1,5 +1,5 @@
 import numpy as np
-from utils import mkvc, ndgrid
+from utils import ndgrid
 
 
 def exampleLomGird(nC, exType):

SimPEG/tests/OrderTest.py

@@ -140,7 +140,11 @@ class OrderTest(unittest.TestCase):
                 max_h_old = max_h
             print '---------------------------------------------'
             passTest = np.mean(np.array(order)) > self.tolerance*self.expectedOrder
-            # passTest = len(np.where(np.array(order) > self.tolerance*self.expectedOrder)[0]) > np.floor(0.75*len(order))
+            if passTest:
+                print ['The test be workin!', 'You get a gold star!', 'Yay passed!', 'Happy little convergence test!', 'That was easy!'][np.random.randint(5)]
+            else:
+                print 'Failed to pass test on ' + self._meshType + '.'
+            print ''
             self.assertTrue(passTest)
 
 if __name__ == '__main__':

SimPEG/tests/test_massMatrices.py

@@ -165,23 +165,23 @@ class TestInnerProducts2D(OrderTest):
         err = np.abs(numeric - analytic)
         return err
 
-    # def test_order1_edges(self):
-    #     self.name = "2D Edge Inner Product - Isotropic"
-    #     self.location = 'edges'
-    #     self.sigmaTest = 1
-    #     self.orderTest()
+    def test_order1_edges(self):
+        self.name = "2D Edge Inner Product - Isotropic"
+        self.location = 'edges'
+        self.sigmaTest = 1
+        self.orderTest()
 
-    # def test_order3_edges(self):
-    #     self.name = "2D Edge Inner Product - Anisotropic"
-    #     self.location = 'edges'
-    #     self.sigmaTest = 2
-    #     self.orderTest()
+    def test_order3_edges(self):
+        self.name = "2D Edge Inner Product - Anisotropic"
+        self.location = 'edges'
+        self.sigmaTest = 2
+        self.orderTest()
 
-    # def test_order6_edges(self):
-    #     self.name = "2D Edge Inner Product - Full Tensor"
-    #     self.location = 'edges'
-    #     self.sigmaTest = 3
-    #     self.orderTest()
+    def test_order6_edges(self):
+        self.name = "2D Edge Inner Product - Full Tensor"
+        self.location = 'edges'
+        self.sigmaTest = 3
+        self.orderTest()
 
     def test_order1_faces(self):
         self.name = "2D Face Inner Product - Isotropic"
@@ -189,13 +189,13 @@ class TestInnerProducts2D(OrderTest):
         self.sigmaTest = 1
         self.orderTest()
 
-    def test_order3_faces(self):
+    def test_order2_faces(self):
         self.name = "2D Face Inner Product - Anisotropic"
         self.location = 'faces'
         self.sigmaTest = 2
         self.orderTest()
 
-    def test_order6_faces(self):
+    def test_order3_faces(self):
         self.name = "2D Face Inner Product - Full Tensor"
         self.location = 'faces'
         self.sigmaTest = 3