From 3d936e07a1328cffa8cde7dee7da108c9e60fcfc Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 5 Aug 2013 12:16:55 -0700 Subject: [PATCH] Edge Inner products in 2D NOTE: still has bugs. --- SimPEG/InnerProducts.py | 71 ++++++++++++++++++++++++++++++- SimPEG/exampleGrid.py | 2 +- SimPEG/tests/OrderTest.py | 6 ++- SimPEG/tests/test_massMatrices.py | 34 +++++++-------- 4 files changed, 92 insertions(+), 21 deletions(-) diff --git a/SimPEG/InnerProducts.py b/SimPEG/InnerProducts.py index 5cc391c4..5105295b 100644 --- a/SimPEG/InnerProducts.py +++ b/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])] diff --git a/SimPEG/exampleGrid.py b/SimPEG/exampleGrid.py index a49326df..2179a17c 100644 --- a/SimPEG/exampleGrid.py +++ b/SimPEG/exampleGrid.py @@ -1,5 +1,5 @@ import numpy as np -from utils import mkvc, ndgrid +from utils import ndgrid def exampleLomGird(nC, exType): diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/OrderTest.py index 107a8d9f..238c3216 100644 --- a/SimPEG/tests/OrderTest.py +++ b/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__': diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index 83d8bf3e..c1276403 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/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