diff --git a/SimPEG/Mesh/NewTreeMesh.py b/SimPEG/Mesh/NewTreeMesh.py index 6e068b87..38aa375c 100644 --- a/SimPEG/Mesh/NewTreeMesh.py +++ b/SimPEG/Mesh/NewTreeMesh.py @@ -1,5 +1,5 @@ import numpy as np, scipy.sparse as sp -from SimPEG.Utils import ndgrid, mkvc +from SimPEG.Utils import ndgrid, mkvc, sdiag NUM, ACTIVE, NX, NY, NZ = range(5) @@ -265,7 +265,7 @@ class TreeMesh(object): q = np.sum((n3 - n1)**2,axis=1)**0.5 # Area of an arbitrary quadrilateral (in a plane) - V = 0.25 * (4.0*(p**2)*(q**2) - (a**2 + c**2 - b**2 - d**2))**0.5 + V = 0.25 * (4.0*(p**2)*(q**2) - (a**2 + c**2 - b**2 - d**2)**2)**0.5 P = np.argsort(C[activeCells,NUM]) self._vol = V[P] @@ -391,7 +391,7 @@ class TreeMesh(object): VOL = self.vol D = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF)) S = self.area - self._faceDiv = Utils.sdiag(1/VOL)*D*Utils.sdiag(S) + self._faceDiv = sdiag(1/VOL)*D*sdiag(S) return self._faceDiv def number(self): @@ -452,10 +452,10 @@ if __name__ == '__main__': from SimPEG import Mesh, Utils import matplotlib.pyplot as plt - tM = TreeMesh(np.ones(3),np.ones(2)) + tM = TreeMesh([np.ones(3),np.ones(2)]) - tM.refineFace(0) - tM.refineFace(9) + # tM.refineFace(0) + # tM.refineFace(9) # print tM._faces # print tM._edges[0,:] diff --git a/SimPEG/Tests/test_NewTreeMesh.py b/SimPEG/Tests/test_NewTreeMesh.py index 21a44220..3d9807e5 100644 --- a/SimPEG/Tests/test_NewTreeMesh.py +++ b/SimPEG/Tests/test_NewTreeMesh.py @@ -69,5 +69,39 @@ class TestQuadTreeMesh(unittest.TestCase): self.assertTrue(np.linalg.norm((np.r_[ax,ay]-self.M.area)) < TOL) + +class SimpleOctreeOperatorTests(unittest.TestCase): + + def setUp(self): + h1 = np.random.rand(5) + h2 = np.random.rand(7) + h3 = np.random.rand(3) + # self.tM = TensorMesh([h1,h2,h3]) + # self.oM = TreeMesh([h1,h2,h3]) + self.tM2 = TensorMesh([h1,h2]) + self.oM2 = TreeMesh([h1,h2]) + # self.oM2.plotGrid(showIt=True) + + def test_faceDiv(self): + # self.assertAlmostEqual((self.tM.faceDiv - self.oM.faceDiv).toarray().sum(), 0) + self.assertAlmostEqual((self.tM2.faceDiv - self.oM2.faceDiv).toarray().sum(), 0) + + # def test_nodalGrad(self): + # self.assertAlmostEqual((self.tM.nodalGrad - self.oM.nodalGrad).toarray().sum(), 0) + # self.assertAlmostEqual((self.tM2.nodalGrad - self.oM2.nodalGrad).toarray().sum(), 0) + + # def test_edgeCurl(self): + # self.assertAlmostEqual((self.tM.edgeCurl - self.oM.edgeCurl).toarray().sum(), 0) + # # self.assertAlmostEqual((self.tM2.edgeCurl - self.oM2.edgeCurl).toarray().sum(), 0) + + # def test_InnerProducts(self): + # self.assertAlmostEqual((self.tM.getFaceInnerProduct() - self.oM.getFaceInnerProduct()).toarray().sum(), 0) + # self.assertAlmostEqual((self.tM2.getFaceInnerProduct() - self.oM2.getFaceInnerProduct()).toarray().sum(), 0) + # self.assertAlmostEqual((self.tM2.getEdgeInnerProduct() - self.oM2.getEdgeInnerProduct()).toarray().sum(), 0) + # self.assertAlmostEqual((self.tM.getEdgeInnerProduct() - self.oM.getEdgeInnerProduct()).toarray().sum(), 0) + + + + if __name__ == '__main__': unittest.main()