From 26b334585fb30fbfd6b7141968ad0ea5340828f7 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 3 Aug 2013 16:39:34 -0700 Subject: [PATCH] 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 | 6 +++--- SimPEG/InnerProducts.py | 19 ++++++++++++++--- SimPEG/sputils.py | 34 +++++++++++++++++++++++++++++++ SimPEG/tests/OrderTest.py | 4 ++-- SimPEG/tests/test_massMatrices.py | 2 ++ 5 files changed, 57 insertions(+), 8 deletions(-) diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index df321e5f..d02b26bf 100644 --- a/SimPEG/DiffOperators.py +++ b/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 diff --git a/SimPEG/InnerProducts.py b/SimPEG/InnerProducts.py index b8819e6a..d1d0d623 100644 --- a/SimPEG/InnerProducts.py +++ b/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 diff --git a/SimPEG/sputils.py b/SimPEG/sputils.py index d9aab184..ae41a210 100644 --- a/SimPEG/sputils.py +++ b/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 diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/OrderTest.py index 1b7ba7de..5bddf118 100644 --- a/SimPEG/tests/OrderTest.py +++ 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]) diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index c304dda5..da9f0ca6 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/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])