diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 0cd0cca3..b35dac1a 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -887,6 +887,28 @@ class TreeMesh(InnerProducts, BaseMesh): self._faceDiv = Utils.sdiag(1/VOL)*D*Utils.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 = [], [], [] + for face in self.faces: + for edge in face.edges: + j = face.edges[edge].index + I += [face.num]*len(j) + J += j + isNeg = lambda e: {'0':True,'1':False,'2':True,'3':False}[e[1]] + V += [-1 if isNeg(edge) else 1]*len(j) + C = sp.csr_matrix((V,(I,J)), shape=(self.nF, self.nE)) + S = self.area + L = self.edge + self._edgeCurl = Utils.sdiag(1/S)*C*Utils.sdiag(L) + return self._edgeCurl + @property def nodalGrad(self): if getattr(self, '_nodalGrad', None) is None: diff --git a/SimPEG/Tests/test_TreeMesh.py b/SimPEG/Tests/test_TreeMesh.py index acdc4308..7a744d0a 100644 --- a/SimPEG/Tests/test_TreeMesh.py +++ b/SimPEG/Tests/test_TreeMesh.py @@ -508,6 +508,10 @@ class SimpleOctreeOperatorTests(unittest.TestCase): self.assertTrue((self.tM.nodalGrad - self.oM.nodalGrad).toarray().sum() == 0) self.assertTrue((self.tM2.nodalGrad - self.oM2.nodalGrad).toarray().sum() == 0) + def test_edgeCurl(self): + self.assertTrue((self.tM.edgeCurl - self.oM.edgeCurl).toarray().sum() == 0) + # self.assertTrue((self.tM2.edgeCurl - self.oM2.edgeCurl).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)