From 0b47ee501134b96d806a600461dfdb4af1eff48e Mon Sep 17 00:00:00 2001
From: SEOGI KANG <sgkang09@gmail.com>
Date: Thu, 18 Sep 2014 17:02:29 -0700
Subject: [PATCH] Add 2d curl and testing

---
 SimPEG/Mesh/DiffOperators.py   | 44 ++++++++++++++++++----------------
 SimPEG/Tests/test_operators.py | 21 ++++++++++++++++
 2 files changed, 45 insertions(+), 20 deletions(-)

diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py
index 966397b3..e46bed20 100644
--- a/SimPEG/Mesh/DiffOperators.py
+++ b/SimPEG/Mesh/DiffOperators.py
@@ -422,11 +422,9 @@ class DiffOperators(object):
 
         def fget(self):
             if(self._edgeCurl is None):
-                assert self.dim > 2, "Edge Curl only programed for 3D."
+                assert self.dim > 1, "Edge Curl only programed for 2 or 3D."
                 # The number of cell centers in each direction
-                n1 = self.nCx
-                n2 = self.nCy
-                n3 = self.nCz
+                n = self.vnC
 
                 # Compute lengths of cell edges
                 L = self.edge
@@ -435,26 +433,32 @@ class DiffOperators(object):
                 S = self.area
 
                 # Compute divergence operator on faces
-                d1 = ddx(n1)
-                d2 = ddx(n2)
-                d3 = ddx(n3)
+                if self.dim == 2:
 
-                D32 = kron3(d3, speye(n2), speye(n1+1))
-                D23 = kron3(speye(n3), d2, speye(n1+1))
-                D31 = kron3(d3, speye(n2+1), speye(n1))
-                D13 = kron3(speye(n3), speye(n2+1), d1)
-                D21 = kron3(speye(n3+1), d2, speye(n1))
-                D12 = kron3(speye(n3+1), speye(n2), d1)
+                    D21 = sp.kron(ddx(n[1]), speye(n[0]))
+                    D12 = sp.kron(speye(n[1]), ddx(n[0]))
+                    C = sp.hstack((-D21, D12), format="csr")
+                    self._edgeCurl = C*sdiag(1/S)
 
-                O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1])
-                O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
-                O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])
+                elif self.dim == 3:
 
-                C = sp.vstack((sp.hstack((O1, -D32, D23)),
-                               sp.hstack((D31, O2, -D13)),
-                               sp.hstack((-D21, D12, O3))), format="csr")
+                    D32 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]+1))
+                    D23 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]+1))
+                    D31 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]))
+                    D13 = kron3(speye(n[2]), speye(n[1]+1), ddx(n[0]))
+                    D21 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]))
+                    D12 = kron3(speye(n[2]+1), speye(n[1]), ddx(n[0]))
+
+                    O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1])
+                    O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
+                    O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])
+
+                    C = sp.vstack((sp.hstack((O1, -D32, D23)),
+                                   sp.hstack((D31, O2, -D13)),
+                                   sp.hstack((-D21, D12, O3))), format="csr")
+
+                    self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
 
-                self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
             return self._edgeCurl
         return locals()
     _edgeCurl = None
diff --git a/SimPEG/Tests/test_operators.py b/SimPEG/Tests/test_operators.py
index 1610f975..c15a44ac 100644
--- a/SimPEG/Tests/test_operators.py
+++ b/SimPEG/Tests/test_operators.py
@@ -49,6 +49,27 @@ class TestCurl(OrderTest):
     def test_order(self):
         self.orderTest()
 
+class TestCurl2D(OrderTest):
+    name = "Cell Grad 2D - Dirichlet"
+    meshTypes = ['uniformTensorMesh']
+    meshDimension = 2
+    meshSizes = [8, 16, 32, 64]
+
+    def getError(self):
+        #Test function
+        ex = lambda x, y: np.cos(y)
+        ey = lambda x, y: np.cos(x)
+        sol = lambda x, y: -np.sin(x)+np.sin(y)
+
+        sol_curl2d = call2(sol, self.M.gridCC)
+        Ec = cartE2(self.M, ex, ey)
+        sol_anal = self.M.edgeCurl*self.M.projectFaceVector(Ec)
+        err = np.linalg.norm((sol_curl2d-sol_anal), np.inf)
+
+        return err
+
+    def test_order(self):
+        self.orderTest()
 
 class TestCellGrad1D_InhomogeneousDirichlet(OrderTest):
     name = "Cell Grad 1D - Dirichlet"