From 821932ed821184a93acdc0c4dcada94458edbbbc Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 12 Nov 2015 08:19:15 -0800 Subject: [PATCH] aveF2CCV in 2 and 3D --- SimPEG/Mesh/TreeMesh.py | 96 +++++++++++++++++++++++++++----- tests/mesh/test_TreeOperators.py | 73 ++++++++++++------------ 2 files changed, 120 insertions(+), 49 deletions(-) diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index b9bc7c74..057c0f59 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -1319,7 +1319,10 @@ class TreeMesh(BaseMesh, InnerProducts): @property def aveE2CC(self): "Construct the averaging operator on cell edges to cell centers." - raise Exception('Not yet implemented!') + if getattr(self, '_aveE2CC', None) is None: + if self.dim == 2: + self._aveE2CC = self.aveF2CC + return self._aveE2CC @property def aveE2CCV(self): @@ -1366,26 +1369,91 @@ class TreeMesh(BaseMesh, InnerProducts): Av = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntF)) - R = self._deflationMatrix('F',asOnes=True,withHanging=True) + Rf = self._deflationMatrix('F',asOnes=True,withHanging=True) - # VOL = self.vol - # if self.dim == 2: - # S = np.r_[self._areaFxFull, self._areaFyFull] - # elif self.dim == 3: - # S = np.r_[self._areaFxFull, self._areaFyFull, self._areaFzFull] - # self._faceDiv = Utils.sdiag(1.0/VOL)*D*Utils.sdiag(S)*R - # return self._faceDiv - - # raise Exception('Not yet implemented!') - self._aveF2CC = Av*R + self._aveF2CC = Av*Rf return self._aveF2CC - @property def aveF2CCV(self): "Construct the averaging operator on cell faces to cell centers." - raise Exception('Not yet implemented!') + if getattr(self, '_aveF2CCV', None) is None: + # TODO: Preallocate! + I, J, V = [], [], [] + PM = [1./2.]*2 # 0.5, 0.5 + + offsetx = [0]*2 + offsety = [self.ntFx]*2 + + if self.dim == 2: + for ii, ind in enumerate(self._sortedCells): + p = self._pointer(ind) + w = self._levelWidth(p[-1]) + + facesx = [ + self._fx2i[self._index([ p[0] , p[1] , p[2]])], + self._fx2i[self._index([ p[0] + w, p[1] , p[2]])], + ] + + facesy = [ + self._fy2i[self._index([ p[0] , p[1] , p[2]])], + self._fy2i[self._index([ p[0] , p[1] + w, p[2]])], + ] + + for off, pm, face in zip(offsetx,PM,facesx): + I += [ii] + J += [face + off] + V += [pm] + + for off, pm, face in zip(offsety,PM,facesy): + I += [ii + self.nC] + J += [face + off] + V += [pm] + + + + if self.dim == 3: + offsetz = [self.ntFx + self.ntFy]*2 + + for ii, ind in enumerate(self._sortedCells): + p = self._pointer(ind) + w = self._levelWidth(p[-1]) + + facesx = [ + self._fx2i[self._index([ p[0] , p[1] , p[2] , p[3]])], + self._fx2i[self._index([ p[0] + w, p[1] , p[2] , p[3]])], + ] + + facesy = [ + self._fy2i[self._index([ p[0] , p[1] , p[2] , p[3]])], + self._fy2i[self._index([ p[0] , p[1] + w, p[2] , p[3]])], + ] + facesz = [ + self._fz2i[self._index([ p[0] , p[1] , p[2] , p[3]])], + self._fz2i[self._index([ p[0] , p[1] , p[2] + w, p[3]])] + ] + + for off, pm, face in zip(offsetx,PM,facesx): + I += [ii] + J += [face + off] + V += [pm] + + for off, pm, face in zip(offsety,PM,facesy): + I += [ii + self.nC] + J += [face + off] + V += [pm] + + for off, pm, face in zip(offsetz,PM,facesz): + I += [ii + self.nC*2] + J += [face + off] + V += [pm] + + Av = sp.csr_matrix((V,(I,J)), shape=(self.nC*self.dim, self.ntF)) + Rf = self._deflationMatrix('F',asOnes=True,withHanging=True) + + self._aveF2CCV = Av*Rf + return self._aveF2CCV def _getFaceP(self, xFace, yFace, zFace): diff --git a/tests/mesh/test_TreeOperators.py b/tests/mesh/test_TreeOperators.py index b50b4000..2f0e066a 100644 --- a/tests/mesh/test_TreeOperators.py +++ b/tests/mesh/test_TreeOperators.py @@ -395,7 +395,7 @@ class TestTreeAveraging2D(Tests.OrderTest): meshTypes = ['notatreeTree', 'uniformTree']#, 'randomTree'] meshDimension = 2 - meshSizes = [4,8,16,32] + meshSizes = [4,8,16] expectedOrders = [2,1] def getError(self): @@ -437,24 +437,23 @@ class TestTreeAveraging2D(Tests.OrderTest): # self.getAve = lambda M: M.aveN2E # self.orderTest() - - def test_orderF2CC(self): - self.name = "Averaging 2D: F2CC" - fun = lambda x, y: (np.cos(x)+np.sin(y)) - self.getHere = lambda M: np.r_[call2(fun, np.r_[M.gridFx, M.gridFy])] - self.getThere = lambda M: call2(fun, M.gridCC) - self.getAve = lambda M: M.aveF2CC - self.orderTest() - - # def test_orderF2CCV(self): - # self.name = "Averaging 2D: F2CCV" - # funX = lambda x, y: (np.cos(x)+np.sin(y)) - # funY = lambda x, y: (np.cos(y)*np.sin(x)) - # self.getHere = lambda M: np.r_[call2(funX, M.gridFx), call2(funY, M.gridFy)] - # self.getThere = lambda M: np.r_[call2(funX, M.gridCC), call2(funY, M.gridCC)] - # self.getAve = lambda M: M.aveF2CCV + # def test_orderF2CC(self): + # self.name = "Averaging 2D: F2CC" + # fun = lambda x, y: (np.cos(x)+np.sin(y)) + # self.getHere = lambda M: np.r_[call2(fun, np.r_[M.gridFx, M.gridFy])] + # self.getThere = lambda M: call2(fun, M.gridCC) + # self.getAve = lambda M: M.aveF2CC # self.orderTest() + def test_orderF2CCV(self): + self.name = "Averaging 2D: F2CCV" + funX = lambda x, y: (np.cos(x)+np.sin(y)) + funY = lambda x, y: (np.cos(y)*np.sin(x)) + self.getHere = lambda M: np.r_[call2(funX, M.gridFx), call2(funY, M.gridFy)] + self.getThere = lambda M: np.r_[call2(funX, M.gridCC), call2(funY, M.gridCC)] + self.getAve = lambda M: M.aveF2CCV + self.orderTest() + # def test_orderCC2F(self): # self.name = "Averaging 2D: CC2F" # fun = lambda x, y: (np.cos(x)+np.sin(y)) @@ -468,7 +467,7 @@ class TestTreeAveraging2D(Tests.OrderTest): # def test_orderE2CC(self): # self.name = "Averaging 2D: E2CC" # fun = lambda x, y: (np.cos(x)+np.sin(y)) - # self.getHere = lambda M: np.r_[call2(fun, M.gridEx), call2(fun, M.gridEy)] + # self.getHere = lambda M: np.r_[call2(fun, np.r_[M.gridEx, M.gridEy])] # self.getThere = lambda M: call2(fun, M.gridCC) # self.getAve = lambda M: M.aveE2CC # self.orderTest() @@ -486,10 +485,14 @@ class TestAveraging3D(Tests.OrderTest): name = "Averaging 3D" meshTypes = ['notatreeTree', 'uniformTree']#, 'randomTree'] meshDimension = 3 - meshSizes = [8,16] + meshSizes = [4,8,16] expectedOrders = [2,1] def getError(self): + if plotIt: + plt.spy(self.getAve(self.M)) + plt.show() + num = self.getAve(self.M) * self.getHere(self.M) err = np.linalg.norm((self.getThere(self.M)-num), np.inf) return err @@ -518,23 +521,23 @@ class TestAveraging3D(Tests.OrderTest): # self.getAve = lambda M: M.aveN2E # self.orderTest() - def test_orderF2CC(self): - self.name = "Averaging 3D: F2CC" - fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z)) - self.getHere = lambda M: np.r_[call3(fun, M.gridFx), call3(fun, M.gridFy), call3(fun, M.gridFz)] - self.getThere = lambda M: call3(fun, M.gridCC) - self.getAve = lambda M: M.aveF2CC - self.orderTest() + # def test_orderF2CC(self): + # self.name = "Averaging 3D: F2CC" + # fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z)) + # self.getHere = lambda M: np.r_[call3(fun, M.gridFx), call3(fun, M.gridFy), call3(fun, M.gridFz)] + # self.getThere = lambda M: call3(fun, M.gridCC) + # self.getAve = lambda M: M.aveF2CC + # self.orderTest() -# def test_orderF2CCV(self): -# self.name = "Averaging 3D: F2CCV" -# funX = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z)) -# funY = lambda x, y, z: (np.cos(x)+np.sin(y)*np.exp(z)) -# funZ = lambda x, y, z: (np.cos(x)*np.sin(y)+np.exp(z)) -# self.getHere = lambda M: np.r_[call3(funX, M.gridFx), call3(funY, M.gridFy), call3(funZ, M.gridFz)] -# self.getThere = lambda M: np.r_[call3(funX, M.gridCC), call3(funY, M.gridCC), call3(funZ, M.gridCC)] -# self.getAve = lambda M: M.aveF2CCV -# self.orderTest() + def test_orderF2CCV(self): + self.name = "Averaging 3D: F2CCV" + funX = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z)) + funY = lambda x, y, z: (np.cos(x)+np.sin(y)*np.exp(z)) + funZ = lambda x, y, z: (np.cos(x)*np.sin(y)+np.exp(z)) + self.getHere = lambda M: np.r_[call3(funX, M.gridFx), call3(funY, M.gridFy), call3(funZ, M.gridFz)] + self.getThere = lambda M: np.r_[call3(funX, M.gridCC), call3(funY, M.gridCC), call3(funZ, M.gridCC)] + self.getAve = lambda M: M.aveF2CCV + self.orderTest() # def test_orderE2CC(self): # self.name = "Averaging 3D: E2CC"