From 81a36684ca236beb6ddae6465ece4c9c5b0e03b9 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 18 Nov 2015 15:34:55 -0800 Subject: [PATCH] NodalGrad --- SimPEG/Mesh/TreeMesh.py | 84 +++++++++++++++++++++++--------- tests/mesh/test_TreeOperators.py | 62 ++++++++++++++++++++++- 2 files changed, 122 insertions(+), 24 deletions(-) diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 88b2d26e..75f22773 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -1231,7 +1231,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts): self._hanging(force=force) def _deflationMatrix(self, location, withHanging=True, asOnes=False): - assert location in ['N','F','Fx','Fy'] + (['Fz','E','Ex','Ey','Ez'] if self.dim == 3 else []) + assert location in ['N','F','Fx','Fy','E','Ex','Ey'] + (['Fz','Ez'] if self.dim == 3 else []) args = dict() args['N'] = (self._nodes, self._hangingN, self._n2i ) @@ -1242,6 +1242,9 @@ class TreeMesh(BaseTensorMesh, InnerProducts): args['Ex'] = (self._edgesX, self._hangingEx, self._ex2i) args['Ey'] = (self._edgesY, self._hangingEy, self._ey2i) args['Ez'] = (self._edgesZ, self._hangingEz, self._ez2i) + elif self.dim == 2: + args['Ex'] = (self._facesY, self._hangingFy, self._fy2i) + args['Ey'] = (self._facesX, self._hangingFx, self._fx2i) if location in ['F', 'E']: Rlist = [self._deflationMatrix(location + subLoc, withHanging=withHanging, asOnes=asOnes) for subLoc in ['x','y','z'][:self.dim]] return sp.block_diag(Rlist) @@ -1401,27 +1404,64 @@ class TreeMesh(BaseTensorMesh, InnerProducts): @property def nodalGrad(self): - raise Exception('Not yet implemented!') - # if getattr(self, '_nodalGrad', None) is None: - # self.number() - # # TODO: Preallocate! - # I, J, V = [], [], [] - # # kinda a hack for the 2D gradient - # # because edges are not stored - # edges = self.faces if self.dim == 2 else self.edges - # for edge in edges: - # if self.dim == 3: - # I += [edge.num, edge.num] - # elif self.dim == 2 and edge.faceType == 'x': - # I += [edge.num + self.nFy, edge.num + self.nFy] - # elif self.dim == 2 and edge.faceType == 'y': - # I += [edge.num - self.nFx, edge.num - self.nFx] - # 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 + if getattr(self, '_nodalGrad', None) is None: + self.number() + # TODO: Preallocate! + I, J, V = [], [], [] + # kinda a hack for the 2D gradient + # because edges are not stored + edgesX = self._facesY if self.dim == 2 else self._edgesX + offset = 0 + for ex in edgesX: + p = self._pointer(ex) + w = self._levelWidth(p[-1]) + if self.dim == 2: + I += [self._fy2i[ex] + offset]*2 + nodePlus = self._index([ p[0] + w, p[1], p[2]]) + elif self.dim == 3: + I += [self._ex2i[ex] + offset]*2 + nodePlus = self._index([ p[0] + w, p[1], p[2], p[3]]) + J += [self._n2i[ex], self._n2i[nodePlus]] + V += [-1, 1] + + edgesY = self._facesX if self.dim == 2 else self._edgesY + offset = self.ntFy if self.dim == 2 else self.ntEx + for ey in edgesY: + p = self._pointer(ey) + w = self._levelWidth(p[-1]) + if self.dim == 2: + I += [self._fx2i[ey] + offset]*2 + nodePlus = self._index([ p[0], p[1] + w, p[2]]) + elif self.dim == 3: + I += [self._ey2i[ey] + offset]*2 + nodePlus = self._index([ p[0], p[1] + w, p[2], p[3]]) + J += [self._n2i[ey], self._n2i[nodePlus]] + V += [-1, 1] + if self.dim == 3: + + edgesZ = self._edgesZ + offset = self.ntEx + self.ntEy + for ez in edgesZ: + p = self._pointer(ez) + w = self._levelWidth(p[-1]) + I += [self._ez2i[ez] + offset]*2 + nodePlus = self._index([ p[0], p[1], p[2] + w, p[3]]) + J += [self._n2i[ez], self._n2i[nodePlus]] + V += [-1, 1] + + G = sp.csr_matrix((V,(I,J)), shape=(self.ntE, self.ntN)) + if self.dim == 2: + L = np.r_[self._areaFyFull, self._areaFxFull] + elif self.dim == 3: + L = np.r_[self._edgeExFull, self._edgeEyFull, self._edgeEzFull] + + Rn = self._deflationMatrix('N') + Re = self._deflationMatrix('E', withHanging=True, asOnes=False) + + Re_ave = Utils.sdiag(1./Re.sum(axis=0)) * Re.T + + self._nodalGrad = Re_ave*Utils.sdiag(1/L)*G*Rn + return self._nodalGrad @property def aveEx2CC(self): diff --git a/tests/mesh/test_TreeOperators.py b/tests/mesh/test_TreeOperators.py index 73385842..c16f6b23 100644 --- a/tests/mesh/test_TreeOperators.py +++ b/tests/mesh/test_TreeOperators.py @@ -69,9 +69,9 @@ class TestFaceDiv3D(Tests.OrderTest): class TestCurl(Tests.OrderTest): name = "Curl" - meshTypes = MESHTYPES + meshTypes = ['notatreeTree', 'uniformTree'] #, 'randomTree']#, 'uniformTree'] meshSizes = [8, 16]#, 32] - expectedOrders = 1 # This is due to linear interpolation in the Re projection + expectedOrders = [2,1] # This is due to linear interpolation in the Re projection def getError(self): # fun: i (cos(y)) + j (cos(z)) + k (cos(x)) @@ -102,6 +102,62 @@ class TestCurl(Tests.OrderTest): self.orderTest() +class TestNodalGrad(Tests.OrderTest): + name = "Nodal Gradient" + meshTypes = ['notatreeTree', 'uniformTree'] #['randomTree', 'uniformTree'] + meshSizes = [8, 16]#, 32] + expectedOrders = [2,1] + + def getError(self): + #Test function + fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) + # i (sin(x)) + j (sin(y)) + k (sin(z)) + solX = lambda x, y, z: -np.sin(x) + solY = lambda x, y, z: -np.sin(y) + solZ = lambda x, y, z: -np.sin(z) + + phi = call3(fun, self.M.gridN) + gradE = self.M.nodalGrad.dot(phi) + + Ec = cartE3(self.M, solX, solY, solZ) + gradE_ana = self.M.projectEdgeVector(Ec) + + err = np.linalg.norm((gradE-gradE_ana), np.inf) + + return err + + def test_order(self): + self.orderTest() + + +class TestNodalGrad2D(Tests.OrderTest): + name = "Nodal Gradient 2D" + meshTypes = ['notatreeTree', 'uniformTree'] #['randomTree', 'uniformTree'] + meshSizes = [8, 16]#, 32] + expectedOrders = [2,1] + meshDimension = 2 + + def getError(self): + #Test function + fun = lambda x, y: (np.cos(x)+np.cos(y)) + # i (sin(x)) + j (sin(y)) + k (sin(z)) + solX = lambda x, y: -np.sin(x) + solY = lambda x, y: -np.sin(y) + + phi = call2(fun, self.M.gridN) + gradE = self.M.nodalGrad.dot(phi) + + Ec = cartE2(self.M, solX, solY) + gradE_ana = self.M.projectEdgeVector(Ec) + + err = np.linalg.norm((gradE-gradE_ana), np.inf) + + return err + + def test_order(self): + self.orderTest() + + class TestTreeInnerProducts(Tests.OrderTest): """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.""" @@ -613,5 +669,7 @@ class TestAveraging3D(Tests.OrderTest): # self.expectedOrders = 1 # self.orderTest() # self.expectedOrders = 2 + + if __name__ == '__main__': unittest.main()