From 3369096126ec15466d51d0c7486f0deaa02e4629 Mon Sep 17 00:00:00 2001
From: Rowan Cockett <rowan@3ptscience.com>
Date: Wed, 11 Feb 2015 17:30:41 -0800
Subject: [PATCH] updates to edge Curlssssszzzz

---
 SimPEG/Mesh/NewTreeMesh.py       | 28 ++++++++++++++++++++++++++++
 SimPEG/Tests/test_NewTreeMesh.py |  6 +++---
 2 files changed, 31 insertions(+), 3 deletions(-)

diff --git a/SimPEG/Mesh/NewTreeMesh.py b/SimPEG/Mesh/NewTreeMesh.py
index 19951b9f..64bac93f 100644
--- a/SimPEG/Mesh/NewTreeMesh.py
+++ b/SimPEG/Mesh/NewTreeMesh.py
@@ -888,6 +888,30 @@ class TreeMesh(BaseMesh):
             self._faceDiv = sdiag(1/VOL)*D*sdiag(S)
         return self._faceDiv
 
+    @property
+    def edgeCurl(self):
+        """Construct the 3D curl operator."""
+        assert self.dim > 2, "Edge Curl only programed for 3D."
+
+        if getattr(self, '_edgeCurl', None) is None:
+            self.number()
+            # TODO: Preallocate!
+            I, J, V = [], [], []
+            F = self._faces
+            sign_edge = zip([-1,1,-1,1],[FEDGE0, FEDGE1, FEDGE2, FEDGE3])
+            activeFaces = F[:,ACTIVE] == 1
+            for face in F[activeFaces]:
+                for sign, edge in sign_edge:
+                    ij, jrow = self._index('_edges', face[edge])
+                    I += [face[NUM]]*len(ij)
+                    J += list(jrow[:,0])
+                    V += [sign]*len(ij)
+            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)
+        return self._edgeCurl
+
     def plotGrid(self, ax=None, text=True, showIt=False):
         import matplotlib.pyplot as plt
 
@@ -957,6 +981,10 @@ if __name__ == '__main__':
     print Mr.vol
 
 
+    tM = TreeMesh([100,100,100])
+    # print tM.vol
+
+
     # print tM._faces
     # print tM._edges[0,:]
     # print tM.vol
diff --git a/SimPEG/Tests/test_NewTreeMesh.py b/SimPEG/Tests/test_NewTreeMesh.py
index 36d9021b..78100df5 100644
--- a/SimPEG/Tests/test_NewTreeMesh.py
+++ b/SimPEG/Tests/test_NewTreeMesh.py
@@ -90,9 +90,9 @@ class SimpleOctreeOperatorTests(unittest.TestCase):
     #     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_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)