diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index fc0e7265..77d5c9de 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -203,7 +203,9 @@ class TreeMesh(BaseMesh, InnerProducts): '_gridCC', '_gridN', '_gridFx', '_gridFy', '_gridFz', '_gridEx', '_gridEy', '_gridEz', '_area', '_edge', '_vol', '_faceDiv', '_edgeCurl', '_nodalGrad', - '_aveF2CC', '_aveF2CCV', '_aveE2CC', '_aveE2CCV','_aveN2CC' + '_aveFx2CC', '_aveFy2CC', '_aveFz2CC', '_aveF2CC', '_aveF2CCV', + '_aveEx2CC', '_aveEy2CC', '_aveEz2CC', '_aveE2CC', '_aveE2CCV', + '_aveN2CC', ] for p in deleteThese: if hasattr(self, p): delattr(self, p) @@ -1407,207 +1409,241 @@ class TreeMesh(BaseMesh, InnerProducts): # return self._nodalGrad @property - def aveE2CC(self): - "Construct the averaging operator on cell edges to cell centers." - if getattr(self, '_aveE2CC', None) is None: - - # TODO: preallocate + def aveEx2CC(self): + if getattr(self, '_aveEx2CC', None) is None: I, J, V = [], [], [] 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] + raise Exception('aveEx2CC not implemented in 2D') 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 + PM = [1./4.]*4 for ii, ind in enumerate(self._sortedCells): p = self._pointer(ind) w = self._levelWidth(p[-1]) - edges = [ + 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]])], + ] + + for pm, edge in zip(PM,edgesx): + I += [ii] + J += [edge] + V += [pm] + + Av = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntEx)) + Re = self._deflationMatrix('Ex',asOnes=False,withHanging=True) + + self._aveEx2CC = Av*Re + return self._aveEx2CC + + @property + def aveEy2CC(self): + "Construct the averaging operator on cell edges to cell centers." + if getattr(self, '_aveEy2CC', None) is None: + I, J, V = [], [], [] + + if self.dim == 2: + raise NotImplementedError('aveEy2CC not implemented in 2D') + + if self.dim == 3: + PM = [1./4.]*4 # plus / plus + + for ii, ind in enumerate(self._sortedCells): + p = self._pointer(ind) + w = self._levelWidth(p[-1]) + + 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]])], + ] + + for pm, edge in zip(PM,edgesy): + I += [ii] + J += [edge] + V += [pm] + + Av = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntEy)) + Re = self._deflationMatrix('Ey',asOnes=False,withHanging=True) + + self._aveEy2CC = Av*Re + return self._aveEy2CC + + @property + def aveEz2CC(self): + "Construct the averaging operator on cell edges to cell centers." + # raise Exception('Not yet implemented!') + if getattr(self, '_aveEz2CC', None) is None: + I, J, V = [], [], [] + + if self.dim == 2: + raise Exception('There are no z edges in 2D') + + if self.dim == 3: + PM = [1./4.]*4 # plus / plus + + for ii, ind in enumerate(self._sortedCells): + p = self._pointer(ind) + w = self._levelWidth(p[-1]) + + 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]])] + self._ez2i[self._index([ p[0] + w, p[1] + w, p[2] , p[3]])], ] - for off, pm, edge in zip(offset,PM,edges): + for pm, edge in zip(PM,edgesz): I += [ii] - J += [edge + off] + J += [edge] V += [pm] - Av = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntE)) - Re = self._deflationMatrix('E',asOnes=False,withHanging=True) + Av = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntEz)) + Re = self._deflationMatrix('Ez',asOnes=False,withHanging=True) - self._aveE2CC = Av*Re + self._aveEz2CC = Av*Re + return self._aveEz2CC + + @property + def aveE2CC(self): + "Construct the averaging operator on cell edges to cell centers." + if getattr(self, '_aveE2CC', None) is None: + if self.dim == 2: + raise Exception('aveE2CC not implemented in 2D') + elif self.dim == 3: + self._aveE2CC = 1./self.dim*sp.hstack([self.aveEx2CC, self.aveEy2CC, self.aveEz2CC]) return self._aveE2CC @property def aveE2CCV(self): "Construct the averaging operator on cell edges to cell centers." - raise Exception('Not yet implemented!') + # raise Exception('Not yet implemented!') + if getattr(self, '_aveE2CCV', None) is None: + if self.dim == 2: + raise Exception('aveE2CC not implemented in 2D') + elif self.dim == 3: + self._aveE2CCV = sp.block_diag([self.aveEx2CC, self.aveEy2CC, self.aveEz2CC]) + return self._aveE2CCV + + @property + def aveFx2CC(self): + if getattr(self, '_aveFx2CC', None) is None: + I, J, V = [], [], [] + PM = [1./2.]*self.dim # 0.5, 0.5 + + for ii, ind in enumerate(self._sortedCells): + p = self._pointer(ind) + w = self._levelWidth(p[-1]) + + if self.dim == 2: + facesx = [ + self._fx2i[self._index([ p[0] , p[1] , p[2]])], + self._fx2i[self._index([ p[0] + w, p[1] , p[2]])], + ] + + elif self.dim == 3: + 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]])], + ] + + for pm, face in zip(PM,facesx): + I += [ii] + J += [face] + V += [pm] + + Av = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntFx)) + Rf = self._deflationMatrix('Fx',asOnes=True,withHanging=True) + + self._aveFx2CC = Av*Rf + return self._aveFx2CC + + @property + def aveFy2CC(self): + if getattr(self, '_aveFy2CC', None) is None: + I, J, V = [], [], [] + PM = [1./2.]*2 # 0.5, 0.5 + + for ii, ind in enumerate(self._sortedCells): + p = self._pointer(ind) + w = self._levelWidth(p[-1]) + + if self.dim == 2: + facesy = [ + self._fy2i[self._index([ p[0] , p[1] , p[2]])], + self._fy2i[self._index([ p[0] , p[1] + w, p[2]])], + ] + elif self.dim == 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]])], + ] + + for pm, face in zip(PM,facesy): + I += [ii] + J += [face] + V += [pm] + + Av = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntFy)) + Rf = self._deflationMatrix('Fy',asOnes=True,withHanging=True) + + self._aveFy2CC = Av*Rf + return self._aveFy2CC + + @property + def aveFz2CC(self): + if getattr(self, '_aveFz2CC', None) is None: + I, J, V = [], [], [] + PM = [1./2.]*2 # 0.5, 0.5 + + for ii, ind in enumerate(self._sortedCells): + p = self._pointer(ind) + w = self._levelWidth(p[-1]) + + if self.dim == 2: + raise Exception('There are no z-faces in 2D') + elif self.dim == 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 pm, face in zip(PM,facesz): + I += [ii] + J += [face] + V += [pm] + + Av = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntFz)) + Rf = self._deflationMatrix('Fz',asOnes=True,withHanging=True) + self._aveFz2CC = Av*Rf + return self._aveFz2CC @property def aveF2CC(self): "Construct the averaging operator on cell faces to cell centers." if getattr(self, '_aveF2CC', None) is None: - # TODO: Preallocate! - I, J, V = [], [], [] - PM = [1./(2.*self.dim)]*2*self.dim # plus / plus - - # TODO total number of faces? - offset = [0]*2 + [self.ntFx]*2 + [self.ntFx+self.ntFy]*2 - - for ii, ind in enumerate(self._sortedCells): - - p = self._pointer(ind) - w = self._levelWidth(p[-1]) - - if self.dim == 2: - faces = [ - self._fx2i[self._index([ p[0] , p[1] , p[2]])], - self._fx2i[self._index([ p[0] + w, p[1] , p[2]])], - self._fy2i[self._index([ p[0] , p[1] , p[2]])], - self._fy2i[self._index([ p[0] , p[1] + w, p[2]])] - ] - elif self.dim == 3: - faces = [ - 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]])], - 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]])], - 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(offset,PM,faces): - I += [ii] - J += [face + off] - V += [pm] - - - Av = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntF)) - Rf = self._deflationMatrix('F',asOnes=True,withHanging=True) - - self._aveF2CC = Av*Rf + if self.dim == 2: + self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC]) + elif self.dim == 3: + self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]) return self._aveF2CC - @property def aveF2CCV(self): "Construct the averaging operator on cell faces to cell centers." 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 + self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC]) + elif self.dim == 3: + self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]) return self._aveF2CCV