From c259a6651a7760ff2cf93aca20642516f48673bd Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 30 Jul 2013 22:12:53 -0700 Subject: [PATCH] Integrated getEdge/FaceInnerProduct into the tensor mesh class. --- SimPEG/DiffOperators.py | 2 +- .../{getInnerProducts.py => InnerProducts.py} | 42 ++++++--- SimPEG/TensorMesh.py | 5 +- SimPEG/getEdgeInnerProducts.py | 87 ------------------- SimPEG/getFaceInnerProducts.py | 87 ------------------- SimPEG/tests/OrderTest.py | 2 +- SimPEG/tests/test_massMatrices.py | 5 +- 7 files changed, 39 insertions(+), 191 deletions(-) rename SimPEG/{getInnerProducts.py => InnerProducts.py} (82%) delete mode 100644 SimPEG/getEdgeInnerProducts.py delete mode 100644 SimPEG/getFaceInnerProducts.py diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index c5daa96d..6688f254 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -50,7 +50,7 @@ class DiffOperators(object): Class creates the differential operators that you need! """ def __init__(self): - raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.') + raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.') def faceDiv(): doc = "Construct divergence operator (face-stg to cell-centres)." diff --git a/SimPEG/getInnerProducts.py b/SimPEG/InnerProducts.py similarity index 82% rename from SimPEG/getInnerProducts.py rename to SimPEG/InnerProducts.py index 94c02826..bc482de4 100644 --- a/SimPEG/getInnerProducts.py +++ b/SimPEG/InnerProducts.py @@ -4,6 +4,28 @@ from utils import sub2ind, ndgrid, mkvc import numpy as np +class InnerProducts(object): + """ + Class creates the inner product matrices that you need! + """ + def __init__(self): + raise Exception('InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class.') + + def getFaceInnerProduct(self, mu): + if self._meshType == 'TENSOR': + pass + elif self._meshType == 'LOM': + pass # todo: we should be doing something slightly different here! + return getFaceInnerProduct(self, mu) + + def getEdgeInnerProduct(self, sigma): + if self._meshType == 'TENSOR': + pass + elif self._meshType == 'LOM': + pass # todo: we should be doing something slightly different here! + return getEdgeInnerProduct(self, sigma) + + def getFaceInnerProduct(mesh, mu): m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) @@ -39,15 +61,15 @@ def getFaceInnerProduct(mesh, mu): # | |/ # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) - # no | node | e1 | e2 | e3 - # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k - # 100 | i+1,j ,k | i+1 ,j ,k | i,j ,k | i,j ,k - # 010 | i ,j+1,k | i ,j,k | i ,j+1 ,k | i ,j,k - # 110 | i+1,j+1,k | i+1 ,j,k | i,j+1 ,k | i,j,k - # 001 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k+1 - # 101 | i+1,j ,k | i+1 ,j ,k | i,j ,k | i,j ,k+1 - # 011 | i ,j+1,k | i ,j,k | i ,j+1 ,k | i ,j,k+1 - # 111 | i+1,j+1,k | i+1 ,j,k | i,j+1 ,k | i,j,k+1 + # no | node | f1 | f2 | f3 + # 000 | i ,j ,k | i , j, k | i, j , k | i, j, k + # 100 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k + # 010 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k + # 110 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k + # 001 | i ,j ,k | i , j, k | i, j , k | i, j, k+1 + # 101 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k+1 + # 011 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k+1 + # 111 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k+1 P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) P100 = Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 0]]) P010 = Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 0]]) @@ -162,4 +184,4 @@ if __name__ == '__main__': h = [np.array([1, 2, 3, 4]), np.array([1, 2, 1, 4, 2]), np.array([1, 1, 4, 1])] mesh = TensorMesh(h) mu = np.ones((mesh.nC, 6)) - A = getFaceInnerProduct(mesh, mu) + A = mesh.getFaceInnerProduct(mu) diff --git a/SimPEG/TensorMesh.py b/SimPEG/TensorMesh.py index b247138b..acfdfbac 100644 --- a/SimPEG/TensorMesh.py +++ b/SimPEG/TensorMesh.py @@ -2,10 +2,11 @@ import numpy as np from BaseMesh import BaseMesh from TensorView import TensorView from DiffOperators import DiffOperators +from InnerProducts import InnerProducts from utils import ndgrid, mkvc -class TensorMesh(BaseMesh, TensorView, DiffOperators): +class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): """ TensorMesh is a mesh class that deals with tensor product meshes. @@ -21,6 +22,8 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators): mesh = TensorMesh([hx, hy, hz]) """ + _meshType = 'TENSOR' + def __init__(self, h, x0=None): super(TensorMesh, self).__init__(np.array([x.size for x in h]), x0) diff --git a/SimPEG/getEdgeInnerProducts.py b/SimPEG/getEdgeInnerProducts.py deleted file mode 100644 index 68423a3f..00000000 --- a/SimPEG/getEdgeInnerProducts.py +++ /dev/null @@ -1,87 +0,0 @@ -from scipy import sparse as sp -from sputils import sdiag -from utils import sub2ind, ndgrid, mkvc -import numpy as np - - -def getEdgeInnerProduct(mesh, sigma): - - m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) - nc = mesh.nC - - i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2])) - - iijjkk = ndgrid(i, j, k) - ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - - 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] - ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nE[0] + mesh.nE[1] - - 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() - - # node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1) - # / / - # / / | - # edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k) - # / / | - # / / | - # node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k) - # | | | - # | | node(i+1,j+1,k+1) - # | | / - # edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k) - # | | / - # | | / - # | |/ - # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) - - # no | node | e1 | e2 | e3 - # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k - # 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k - # 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k - # 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k - # 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k - # 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k - # 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k - # 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k - P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) - P100 = Pxxx([[0, 0, 0], [1, 0, 0], [1, 0, 0]]) - P010 = Pxxx([[0, 1, 0], [0, 0, 0], [0, 1, 0]]) - P110 = Pxxx([[0, 1, 0], [1, 0, 0], [1, 1, 0]]) - P001 = Pxxx([[0, 0, 1], [0, 0, 1], [0, 0, 0]]) - P101 = Pxxx([[0, 0, 1], [1, 0, 1], [1, 0, 0]]) - P011 = Pxxx([[0, 1, 1], [0, 0, 1], [0, 1, 0]]) - P111 = Pxxx([[0, 1, 1], [1, 0, 1], [1, 1, 0]]) - - if sigma.size == mesh.nC: # Isotropic! - sigma = mkvc(sigma) # ensure it is a vector. - Sigma = sdiag(np.r_[sigma, sigma, sigma]) - elif sigma.shape[1] == 3: # Diagonal tensor - Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]]) - elif sigma.shape[1] == 6: # Fully anisotropic - 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)) - - # Cell volume - v = np.sqrt(mesh.vol) - v3 = np.r_[v, v, v] - V = sdiag(v3)*Sigma*sdiag(v3) # to keep symmetry - - A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111 - - A = 0.125*A - - 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])] - mesh = TensorMesh(h) - sigma = np.ones((mesh.nC, 6)) - A = getEdgeInnerProduct(mesh, sigma) diff --git a/SimPEG/getFaceInnerProducts.py b/SimPEG/getFaceInnerProducts.py deleted file mode 100644 index 961013c8..00000000 --- a/SimPEG/getFaceInnerProducts.py +++ /dev/null @@ -1,87 +0,0 @@ -from scipy import sparse as sp -from sputils import sdiag -from utils import sub2ind, ndgrid, mkvc -import numpy as np - - -def getFaceInnerProduct(mesh, mu): - - m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) - nc = mesh.nC - - i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2])) - - iijjkk = ndgrid(i, j, k) - ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - - def Pxxx(pos): - ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) - ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nF[0] - ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nF[0] + mesh.nF[1] - - 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.nF))).tocsr() - - # node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1) - # / / - # / / | - # edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k) - # / / | - # / / | - # node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k) - # | | | - # | | node(i+1,j+1,k+1) - # | | / - # edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k) - # | | / - # | | / - # | |/ - # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) - - # no | node | e1 | e2 | e3 - # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k - # 100 | i+1,j ,k | i+1 ,j ,k | i,j ,k | i,j ,k - # 010 | i ,j+1,k | i ,j,k | i ,j+1 ,k | i ,j,k - # 110 | i+1,j+1,k | i+1 ,j,k | i,j+1 ,k | i,j,k - # 001 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k+1 - # 101 | i+1,j ,k | i+1 ,j ,k | i,j ,k | i,j ,k+1 - # 011 | i ,j+1,k | i ,j,k | i ,j+1 ,k | i ,j,k+1 - # 111 | i+1,j+1,k | i+1 ,j,k | i,j+1 ,k | i,j,k+1 - P000 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) - P100 = Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 0]]) - P010 = Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 0]]) - P110 = Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 0]]) - P001 = Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 1]]) - P101 = Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 1]]) - P011 = Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 1]]) - P111 = Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) - - if mu.size == mesh.nC: # Isotropic! - mu = mkvc(mu) # ensure it is a vector. - mu = sdiag(np.r_[mu, mu, mu]) - elif mu.shape[1] == 3: # Diagonal tensor - mu = sdiag(np.r_[mu[:, 0], mu[:, 1], mu[:, 2]]) - elif mu.shape[1] == 6: # Fully anisotropic - row1 = sp.hstack((sdiag(mu[:, 0]), sdiag(mu[:, 3]), sdiag(mu[:, 4]))) - row2 = sp.hstack((sdiag(mu[:, 3]), sdiag(mu[:, 1]), sdiag(mu[:, 5]))) - row3 = sp.hstack((sdiag(mu[:, 4]), sdiag(mu[:, 5]), sdiag(mu[:, 2]))) - mu = sp.vstack((row1, row2, row3)) - - # Cell volume - v = np.sqrt(mesh.vol) - v3 = np.r_[v, v, v] - V = sdiag(v3)*mu*sdiag(v3) # to keep symmetry - - A = P000.T*V*P000 + P001.T*V*P001 + P010.T*V*P010 + P011.T*V*P011 + P100.T*V*P100 + P101.T*V*P101 + P110.T*V*P110 + P111.T*V*P111 - - A = 0.125*A - - 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])] - mesh = TensorMesh(h) - mu = np.ones((mesh.nC, 6)) - A = getFaceInnerProduct(mesh, mu) diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/OrderTest.py index 59dce6d6..55ac7d34 100644 --- a/SimPEG/tests/OrderTest.py +++ b/SimPEG/tests/OrderTest.py @@ -64,7 +64,7 @@ class OrderTest(unittest.TestCase): name = "Order Test" expectedOrder = 2 tolerance = 0.85 - meshSizes = [4, 8, 16, 32, 64] + meshSizes = [4, 8, 16, 32] meshType = 'uniformTensorMesh' meshDimension = 3 diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index 0aae0615..4840abed 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/SimPEG/tests/test_massMatrices.py @@ -1,9 +1,6 @@ import numpy as np import unittest from OrderTest import OrderTest -import sys -sys.path.append('../') -from getEdgeInnerProducts import * class TestEdgeInnerProduct(OrderTest): @@ -35,7 +32,7 @@ class TestEdgeInnerProduct(OrderTest): sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc), call(sigma4, Gc), call(sigma5, Gc), call(sigma6, Gc)] - A = getEdgeInnerProduct(self.M, sigma) + A = self.M.getEdgeInnerProduct(sigma) numeric = E.T*A*E analytic = 69881./21600 # Found using matlab symbolic toolbox. err = np.abs(numeric - analytic)