diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 2de403ce..4fc9451f 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -1608,6 +1608,48 @@ class TreeMesh(BaseMesh, InnerProducts): self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]) return self._aveF2CCV + @property + def aveN2CC(self): + if getattr(self, '_aveN2CC', None) is None: + I, J, V = [], [], [] + PM = [1./2.**self.dim] * 2**self.dim + + for ii, ind in enumerate(self._sortedCells): + p = self._pointer(ind) + w = self._levelWidth(p[-1]) + + if self.dim == 2: + nodes = [ + self._n2i[self._index([ p[0] , p[1] , p[2] ])], + self._n2i[self._index([ p[0] + w, p[1] , p[2] ])], + self._n2i[self._index([ p[0] , p[1] + w, p[2] ])], + self._n2i[self._index([ p[0] + w, p[1] + w, p[2] ])], + ] + + + if self.dim == 3: + nodes = [ + self._n2i[self._index([ p[0] , p[1] , p[2] , p[3] ])], + self._n2i[self._index([ p[0] + w, p[1] , p[2] , p[3] ])], + self._n2i[self._index([ p[0] , p[1] + w, p[2] , p[3] ])], + self._n2i[self._index([ p[0] + w, p[1] + w, p[2] , p[3] ])], + self._n2i[self._index([ p[0] , p[1] , p[2] + w, p[3] ])], + self._n2i[self._index([ p[0] + w, p[1] , p[2] + w, p[3] ])], + self._n2i[self._index([ p[0] , p[1] + w, p[2] + w, p[3] ])], + self._n2i[self._index([ p[0] + w, p[1] + w, p[2] + w, p[3] ])], + ] + + for pm, node in zip(PM,nodes): + I += [ii] + J += [node] + V += [pm] + + Av = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntN)) + Re = self._deflationMatrix('N',asOnes=False,withHanging=True) + + self._aveN2CC = Av*Re + return self._aveN2CC + def _getFaceP(self, xFace, yFace, zFace): ind1, ind2, ind3 = [], [], [] diff --git a/tests/mesh/test_TreeOperators.py b/tests/mesh/test_TreeOperators.py index a5eb1358..1a62673e 100644 --- a/tests/mesh/test_TreeOperators.py +++ b/tests/mesh/test_TreeOperators.py @@ -413,13 +413,13 @@ class TestTreeAveraging2D(Tests.OrderTest): return err - # def test_orderN2CC(self): - # self.name = "Averaging 2D: N2CC" - # fun = lambda x, y: (np.cos(x)+np.sin(y)) - # self.getHere = lambda M: call2(fun, M.gridN) - # self.getThere = lambda M: call2(fun, M.gridCC) - # self.getAve = lambda M: M.aveN2CC - # self.orderTest() + def test_orderN2CC(self): + self.name = "Averaging 2D: N2CC" + fun = lambda x, y: (np.cos(x)+np.sin(y)) + self.getHere = lambda M: call2(fun, M.gridN) + self.getThere = lambda M: call2(fun, M.gridCC) + self.getAve = lambda M: M.aveN2CC + self.orderTest() # def test_orderN2F(self): # self.name = "Averaging 2D: N2F" @@ -496,13 +496,13 @@ class TestAveraging3D(Tests.OrderTest): err = np.linalg.norm((self.getThere(self.M)-num), np.inf) return err -# def test_orderN2CC(self): -# self.name = "Averaging 3D: N2CC" -# fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z)) -# self.getHere = lambda M: call3(fun, M.gridN) -# self.getThere = lambda M: call3(fun, M.gridCC) -# self.getAve = lambda M: M.aveN2CC -# self.orderTest() + def test_orderN2CC(self): + self.name = "Averaging 3D: N2CC" + fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z)) + self.getHere = lambda M: call3(fun, M.gridN) + self.getThere = lambda M: call3(fun, M.gridCC) + self.getAve = lambda M: M.aveN2CC + self.orderTest() # def test_orderN2F(self): # self.name = "Averaging 3D: N2F"