diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 02a1e19a..8030a20e 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -1326,33 +1326,7 @@ class TreeMesh(BaseMesh, InnerProducts): if self.dim == 2: raise NotImplementedError('aveE2CC not implemented yet') - # PM = [1./(2.*self.dim)]*self.dim # plus / plus - # offset = [0]*2 + [self.ntEx]*2 - # for ii, ind in enumerate(self._sortedCells): - # p = self._pointer(ind) - # w = self._levelWidth(p[-1]) - - # edges = [ - # self._ex2i[self._index([ p[0] , p[1] , p[2] , p[3]])], - # self._ex2i[self._index([ p[0] , p[1] + w, p[2] , p[3]])], - # self._ex2i[self._index([ p[0] , p[1] , p[2] + w, p[3]])], - # self._ex2i[self._index([ p[0] , p[1] + w, p[2] + w, p[3]])], - # self._ey2i[self._index([ p[0] , p[1] , p[2] , p[3]])], - # self._ey2i[self._index([ p[0] + w, p[1] , p[2] , p[3]])], - # self._ey2i[self._index([ p[0] , p[1] , p[2] + w, p[3]])], - # self._ey2i[self._index([ p[0] + w, p[1] , p[2] + w, p[3]])], - # self._ez2i[self._index([ p[0] , p[1] , p[2] , p[3]])], - # self._ez2i[self._index([ p[0] + w, p[1] , p[2] , p[3]])], - # self._ez2i[self._index([ p[0] , p[1] + w, p[2] , p[3]])], - # self._ez2i[self._index([ p[0] + w, p[1] + w, p[2] , p[3]])] - # ] - - # for off, pm, edge in zip(offset,PM,edges): - # I += [ii] - # J += [edge + off] - # V += [pm] - if self.dim == 3: PM = [1./(4.*self.dim)]*4*self.dim # plus / plus offset = [0]*4 + [self.ntEx]*4 + [self.ntEx+self.ntEy]*4 @@ -1392,7 +1366,61 @@ class TreeMesh(BaseMesh, InnerProducts): @property def aveE2CCV(self): "Construct the averaging operator on cell edges to cell centers." - raise Exception('Not yet implemented!') + # raise Exception('Not yet implemented!') + + I, J, V = [], [], [] + + if self.dim == 2: + raise NotImplementedError('aveE2CC not implemented yet') + + if self.dim == 3: + PM = [1./4.]*4 # plus / plus + offsetx,offsety,offsetz = [0]*4, [self.ntEx]*4 , [self.ntEx+self.ntEy]*4 + + for ii, ind in enumerate(self._sortedCells): + p = self._pointer(ind) + w = self._levelWidth(p[-1]) + + edgesx = [ + self._ex2i[self._index([ p[0] , p[1] , p[2] , p[3]])], + self._ex2i[self._index([ p[0] , p[1] + w, p[2] , p[3]])], + self._ex2i[self._index([ p[0] , p[1] , p[2] + w, p[3]])], + self._ex2i[self._index([ p[0] , p[1] + w, p[2] + w, p[3]])], + ] + edgesy = [ + self._ey2i[self._index([ p[0] , p[1] , p[2] , p[3]])], + self._ey2i[self._index([ p[0] + w, p[1] , p[2] , p[3]])], + self._ey2i[self._index([ p[0] , p[1] , p[2] + w, p[3]])], + self._ey2i[self._index([ p[0] + w, p[1] , p[2] + w, p[3]])], + ] + edgesz = [ + self._ez2i[self._index([ p[0] , p[1] , p[2] , p[3]])], + self._ez2i[self._index([ p[0] + w, p[1] , p[2] , p[3]])], + self._ez2i[self._index([ p[0] , p[1] + w, p[2] , p[3]])], + self._ez2i[self._index([ p[0] + w, p[1] + w, p[2] , p[3]])] + ] + + for off, pm, edge in zip(offsetx,PM,edgesx): + I += [ii] + J += [edge + off] + V += [pm] + + for off, pm, edge in zip(offsety,PM,edgesy): + I += [ii + self.nC] + J += [edge + off] + V += [pm] + + for off, pm, edge in zip(offsetz,PM,edgesz): + I += [ii + self.nC*2.] + J += [edge + off] + V += [pm] + + + Av = sp.csr_matrix((V,(I,J)), shape=(self.nC*self.dim, self.ntE)) + Re = self._deflationMatrix('E',asOnes=False,withHanging=True) + + self._aveE2CCV = Av*Re + return self._aveE2CCV @property def aveF2CC(self): @@ -1433,7 +1461,7 @@ class TreeMesh(BaseMesh, InnerProducts): V += [pm] - Av = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntF)) + Av = sp.csr_matrix((V,(I,J)), shape=(self.nC*self.dim, self.ntF)) Rf = self._deflationMatrix('F',asOnes=True,withHanging=True) self._aveF2CC = Av*Rf diff --git a/tests/mesh/test_TreeOperators.py b/tests/mesh/test_TreeOperators.py index 6be03d33..7e36b15c 100644 --- a/tests/mesh/test_TreeOperators.py +++ b/tests/mesh/test_TreeOperators.py @@ -547,15 +547,15 @@ class TestAveraging3D(Tests.OrderTest): self.getAve = lambda M: M.aveE2CC self.orderTest() -# def test_orderE2CCV(self): -# self.name = "Averaging 3D: E2CCV" -# 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.gridEx), call3(funY, M.gridEy), call3(funZ, M.gridEz)] -# self.getThere = lambda M: np.r_[call3(funX, M.gridCC), call3(funY, M.gridCC), call3(funZ, M.gridCC)] -# self.getAve = lambda M: M.aveE2CCV -# self.orderTest() + def test_orderE2CCV(self): + self.name = "Averaging 3D: E2CCV" + 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.gridEx), call3(funY, M.gridEy), call3(funZ, M.gridEz)] + self.getThere = lambda M: np.r_[call3(funX, M.gridCC), call3(funY, M.gridCC), call3(funZ, M.gridCC)] + self.getAve = lambda M: M.aveE2CCV + self.orderTest() # def test_orderCC2F(self): # self.name = "Averaging 3D: CC2F"