From aa539462a1aa5f959088eb00da61d220223a1bff Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Tue, 11 Feb 2014 13:58:31 -0800
Subject: [PATCH] Nodal Gradient!

---
 SimPEG/Mesh/TreeMesh.py       | 43 +++++++++++++++++++++++++++++++----
 SimPEG/Tests/test_TreeMesh.py |  4 ++++
 2 files changed, 42 insertions(+), 5 deletions(-)

diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py
index 0cc349b6..63280076 100644
--- a/SimPEG/Mesh/TreeMesh.py
+++ b/SimPEG/Mesh/TreeMesh.py
@@ -97,7 +97,7 @@ class TreeEdge(TreeObject):
 
         self.node0 = node0 if isinstance(node0,TreeNode) else TreeNode(mesh, x0=self.x0)
         self.node1 = node1 if isinstance(node1,TreeNode) else TreeNode(mesh, x0=self.x0 + self.tangent*self.sz[0])
-
+        self.nodes = {'n0':node0, 'n1':node1}
 
     def refine(self):
         if not self.isleaf: return
@@ -131,6 +131,10 @@ class TreeEdge(TreeObject):
     def center(self):
         return 0.5*(self.node0.x0 + self.node1.x0)
 
+    @property
+    def length(self):
+        return np.sqrt(((self.node1.x0 - self.node0.x0)**2).sum())
+
     @property
     def index(self):
         if self.isleaf: return [self.num]
@@ -233,6 +237,11 @@ class TreeFace(TreeObject):
         """area of the face"""
         return self.sz.prod()
 
+    @property
+    def length(self):
+        if self.dim == 3: raise Exception('face.length is not defined for 2D face')
+        return np.sqrt(((self.node1.x0 - self.node0.x0)**2).sum())
+
     def refine(self):
         if not self.isleaf: return
         self.mesh.isNumbered = False
@@ -846,12 +855,20 @@ class TreeMesh(InnerProducts, BaseMesh):
     @property
     def area(self):
         self.number()
-        if self.dim == 2:
-            faces = self.sortedFaceX + self.sortedFaceY
-        elif self.dim == 3:
-            faces = self.sortedFaceX + self.sortedFaceY + self.sortedFaceZ
+        faces = self.sortedFaceX + self.sortedFaceY
+        if self.dim == 3:
+            faces += self.sortedFaceZ
         return np.array([face.area for face in faces], dtype=float)
 
+    @property
+    def edge(self):
+        self.number()
+        if self.dim == 2:
+            edges = self.sortedFaceY + self.sortedFaceX
+        elif self.dim == 3:
+            edges = self.sortedEdgeX + self.sortedEdgeY + self.sortedEdgeZ
+        return np.array([e.length for e in edges], dtype=float)
+
     @property
     def faceDiv(self):
         if getattr(self, '_faceDiv', None) is None:
@@ -870,6 +887,22 @@ class TreeMesh(InnerProducts, BaseMesh):
             self._faceDiv = Utils.sdiag(1/VOL)*D*Utils.sdiag(S)
         return self._faceDiv
 
+    @property
+    def nodalGrad(self):
+        if getattr(self, '_nodalGrad', None) is None:
+            self.number()
+            # TODO: Preallocate!
+            I, J, V = [], [], []
+            edges = self.faces if self.dim == 2 else self.edges
+            for edge in edges:
+                I += [edge.num, edge.num]
+                J += [edge.node0.num, edge.node1.num]
+                V += [-1, 1]
+            G = sp.csr_matrix((V,(I,J)), shape=(self.nE, self.nN))
+            L = self.edge
+            self._nodalGrad = Utils.sdiag(1/L)*G
+        return self._nodalGrad
+
     def _getFaceP(self, face0, face1, face2):
         I, J, V = [], [], []
         for cell in self.sortedCells:
diff --git a/SimPEG/Tests/test_TreeMesh.py b/SimPEG/Tests/test_TreeMesh.py
index 254cf127..acdc4308 100644
--- a/SimPEG/Tests/test_TreeMesh.py
+++ b/SimPEG/Tests/test_TreeMesh.py
@@ -504,6 +504,10 @@ class SimpleOctreeOperatorTests(unittest.TestCase):
         self.assertTrue((self.tM.faceDiv - self.oM.faceDiv).toarray().sum() == 0)
         self.assertTrue((self.tM2.faceDiv - self.oM2.faceDiv).toarray().sum() == 0)
 
+    def test_nodalGrad(self):
+        self.assertTrue((self.tM.nodalGrad - self.oM.nodalGrad).toarray().sum() == 0)
+        self.assertTrue((self.tM2.nodalGrad - self.oM2.nodalGrad).toarray().sum() == 0)
+
     def test_InnerProducts(self):
         self.assertTrue((self.tM.getFaceInnerProduct() - self.oM.getFaceInnerProduct()).toarray().sum() == 0)
         self.assertTrue((self.tM2.getFaceInnerProduct() - self.oM2.getFaceInnerProduct()).toarray().sum() == 0)