EdgeInnerProducts now working for LOM

Fixed DiffOps bug in CURL

Fixed bug in OrderTest --> must remember to do cumSum, and not just supply dx to ndgrid

test_massMatrices now tests uniformLOM

SimPEG/DiffOperators.py

@@ -168,9 +168,9 @@ class DiffOperators(object):
         def fget(self):
             if(self._edgeCurl is None):
                 # The number of cell centers in each direction
-                n1 = np.size(self.hx)
-                n2 = np.size(self.hy)
-                n3 = np.size(self.hz)
+                n1 = self.nCx
+                n2 = self.nCy
+                n3 = self.nCy
 
                 # Compute lengths of cell edges
                 L = self.edge

SimPEG/InnerProducts.py

@@ -1,6 +1,6 @@
 from scipy import sparse as sp
-from sputils import sdiag
-from utils import sub2ind, ndgrid, mkvc
+from sputils import sdiag, inv3X3BlockDiagonal
+from utils import sub2ind, ndgrid, mkvc, getSubArray
 import numpy as np
 
 
@@ -123,6 +123,11 @@ def getEdgeInnerProduct(mesh, sigma=None, returnP=False):
     iijjkk = ndgrid(i, j, k)
     ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
 
+    if mesh._meshType == 'LOM':
+        eT1 = mesh.r(mesh.tangents, 'E', 'Ex', 'M')
+        eT2 = mesh.r(mesh.tangents, 'E', 'Ey', 'M')
+        eT3 = mesh.r(mesh.tangents, 'E', 'Ez', 'M')
+
     def Pxxx(pos):
         ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]])
         ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nE[0]
@@ -130,7 +135,15 @@ def getEdgeInnerProduct(mesh, sigma=None, returnP=False):
 
         IND = np.r_[ind1, ind2, ind3].flatten()
 
-        return sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nE))).tocsr()
+        PXXX = sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nE))).tocsr()
+
+        if mesh._meshType == 'LOM':
+            I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(eT1[1], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(eT1[2], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]),
+                                       getSubArray(eT2[0], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(eT2[1], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(eT2[2], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]),
+                                       getSubArray(eT3[0], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(eT3[1], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(eT3[2], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]))
+            PXXX = I3x3 * PXXX
+
+        return PXXX
 
     # no  | node        | e1          | e2          | e3
     # 000 | i  ,j  ,k   | i  ,j  ,k   | i  ,j  ,k   | i  ,j  ,k

SimPEG/sputils.py

@@ -1,4 +1,5 @@
 from scipy import sparse as sp
+from utils import mkvc
 
 
 def sdiag(h):
@@ -19,3 +20,36 @@ def kron3(A, B, C):
 def spzeros(n1, n2):
     """spzeros"""
     return sp.coo_matrix((n1, n2)).tocsr()
+
+
+def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33):
+
+    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

SimPEG/tests/OrderTest.py

@@ -92,8 +92,8 @@ class OrderTest(unittest.TestCase):
 
         elif 'LOM' in self.meshType:
             if 'uniform' in self.meshType:
-                xx = np.ones(nc)/nc
-                X, Y, Z = utils.ndgrid(xx, xx, xx, vector=False)
+                xx = np.cumsum(np.r_[0, np.ones(nc)/nc])
+                X, Y, Z = utils.ndgrid([xx, xx, xx], vector=False)
             else:
                 raise Exception('Unexpected meshType')
             self.M = LogicallyOrthogonalMesh([X, Y, Z])

SimPEG/tests/test_massMatrices.py

@@ -32,6 +32,8 @@ from OrderTest import OrderTest
 class TestInnerProducts(OrderTest):
     """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
 
+    meshType = 'uniformLOM'
+
     def getError(self):
 
         call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])