From 5cdab15aa1ffb04ef2ca7c841523f48667c63a7f Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 11 Feb 2015 21:57:05 -0800 Subject: [PATCH] nodal gradient --- SimPEG/Mesh/NewTreeMesh.py | 20 +++++++++++++++++++- SimPEG/Tests/test_NewTreeMesh.py | 6 +++--- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/SimPEG/Mesh/NewTreeMesh.py b/SimPEG/Mesh/NewTreeMesh.py index 64bac93f..ca793065 100644 --- a/SimPEG/Mesh/NewTreeMesh.py +++ b/SimPEG/Mesh/NewTreeMesh.py @@ -909,9 +909,27 @@ class TreeMesh(BaseMesh): C = sp.csr_matrix((V,(I,J)), shape=(self.nF, self.nE)) S = self.area L = self.edge - self._edgeCurl = sdiag(1/S)*C*sdiag(L) + self._edgeCurl = sdiag(1.0/S)*C*sdiag(L) return self._edgeCurl + @property + def nodalGrad(self): + if getattr(self, '_nodalGrad', None) is None: + self.number() + # TODO: Preallocate! + I, J, V = [], [], [] + E = self._edges + N = self._nodes + activeEdges = E[:,ACTIVE] == 1 + for edge in E[activeEdges]: + I += [edge[NUM], edge[NUM]] + J += [N[edge[ENODE0], NUM], N[edge[ENODE1], NUM]] + V += [-1, 1] + G = sp.csr_matrix((V,(I,J)), shape=(self.nE, self.nN)) + L = self.edge + self._nodalGrad = sdiag(1.0/L)*G + return self._nodalGrad + def plotGrid(self, ax=None, text=True, showIt=False): import matplotlib.pyplot as plt diff --git a/SimPEG/Tests/test_NewTreeMesh.py b/SimPEG/Tests/test_NewTreeMesh.py index 78100df5..ed1d2864 100644 --- a/SimPEG/Tests/test_NewTreeMesh.py +++ b/SimPEG/Tests/test_NewTreeMesh.py @@ -86,9 +86,9 @@ class SimpleOctreeOperatorTests(unittest.TestCase): 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_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)