broke up averaging so that we have averaging for components as well

This commit is contained in:
Lindsey Heagy
2015-11-15 00:08:26 -08:00
parent 4fe9475ffc
commit d7bcc0c074
2 changed files with 291 additions and 191 deletions
+226 -173
View File
@@ -1316,236 +1316,289 @@ class TreeMesh(BaseMesh, InnerProducts):
# self._nodalGrad = Utils.sdiag(1/L)*G
# return self._nodalGrad
@property
def aveE2CC(self):
"Construct the averaging operator on cell edges to cell centers."
if getattr(self, '_aveE2CC', None) is None:
# @property
# def aveE2CC(self):
# "Construct the averaging operator on cell edges to cell centers."
# if getattr(self, '_aveE2CC', None) is None:
# TODO: preallocate
# # TODO: preallocate
# I, J, V = [], [], []
# if self.dim == 2:
# raise NotImplementedError('aveE2CC not implemented yet')
# 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
# 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]
# Av = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.ntE))
# Re = self._deflationMatrix('E',asOnes=False,withHanging=True)
# self._aveE2CC = Av*Re
# return self._aveE2CC
@property
def aveEx2CC(self):
if getattr(self, '_aveEx2CC', None) is None:
I, J, V = [], [], []
if self.dim == 2:
raise NotImplementedError('aveE2CC not implemented yet')
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!')
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
I, J, V = [], [], []
if self.dim == 2:
raise NotImplementedError('aveE2CC not implemented yet')
@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 == 3:
PM = [1./4.]*4 # plus / plus
offsetx,offsety,offsetz = [0]*4, [self.ntEx]*4 , [self.ntEx+self.ntEy]*4
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])
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]])]
]
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 off, pm, edge in zip(offsetx,PM,edgesx):
for pm, face in zip(PM,facesz):
I += [ii]
J += [edge + off]
J += [face]
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
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.dim, 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
if self.dim == 2:
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
+65 -18
View File
@@ -445,6 +445,22 @@ class TestTreeAveraging2D(Tests.OrderTest):
self.getAve = lambda M: M.aveF2CC
self.orderTest()
def test_orderFx2CC(self):
self.name = "Averaging 2D: Fx2CC"
funX = lambda x, y: (np.cos(x)+np.sin(y))
self.getHere = lambda M: np.r_[call2(funX, M.gridFx)]
self.getThere = lambda M: np.r_[call2(funX, M.gridCC)]
self.getAve = lambda M: M.aveFx2CC
self.orderTest()
def test_orderFy2CC(self):
self.name = "Averaging 2D: Fy2CC"
funY = lambda x, y: (np.cos(y)*np.sin(x))
self.getHere = lambda M: np.r_[call2(funY, M.gridFy)]
self.getThere = lambda M: np.r_[call2(funY, M.gridCC)]
self.getAve = lambda M: M.aveFy2CC
self.orderTest()
def test_orderF2CCV(self):
self.name = "Averaging 2D: F2CCV"
funX = lambda x, y: (np.cos(x)+np.sin(y))
@@ -464,28 +480,11 @@ class TestTreeAveraging2D(Tests.OrderTest):
# self.orderTest()
# self.expectedOrders = 2
# 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, np.r_[M.gridEx, M.gridEy])]
# self.getThere = lambda M: call2(fun, M.gridCC)
# self.getAve = lambda M: M.aveE2CC
# self.orderTest()
# def test_orderE2CCV(self):
# self.name = "Averaging 2D: E2CCV"
# 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.gridEx), call2(funY, M.gridEy)]
# self.getThere = lambda M: np.r_[call2(funX, M.gridCC), call2(funY, M.gridCC)]
# self.getAve = lambda M: M.aveE2CCV
# self.orderTest()
class TestAveraging3D(Tests.OrderTest):
name = "Averaging 3D"
meshTypes = ['notatreeTree', 'uniformTree']#, 'randomTree']
meshDimension = 3
meshSizes = [8,16,32]
meshSizes = [8,16]
expectedOrders = [2,1]
def getError(self):
@@ -529,6 +528,30 @@ class TestAveraging3D(Tests.OrderTest):
self.getAve = lambda M: M.aveF2CC
self.orderTest()
def test_orderFx2CC(self):
self.name = "Averaging 3D: Fx2CC"
funX = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
self.getHere = lambda M: np.r_[call3(funX, M.gridFx)]
self.getThere = lambda M: np.r_[call3(funX, M.gridCC)]
self.getAve = lambda M: M.aveFx2CC
self.orderTest()
def test_orderFy2CC(self):
self.name = "Averaging 3D: Fy2CC"
funY = lambda x, y, z: (np.cos(x)+np.sin(y)*np.exp(z))
self.getHere = lambda M: np.r_[call3(funY, M.gridFy)]
self.getThere = lambda M: np.r_[call3(funY, M.gridCC)]
self.getAve = lambda M: M.aveFy2CC
self.orderTest()
def test_orderFz2CC(self):
self.name = "Averaging 3D: Fz2CC"
funZ = lambda x, y, z: (np.cos(x)+np.sin(y)*np.exp(z))
self.getHere = lambda M: np.r_[call3(funZ, M.gridFz)]
self.getThere = lambda M: np.r_[call3(funZ, M.gridCC)]
self.getAve = lambda M: M.aveFz2CC
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))
@@ -539,6 +562,30 @@ class TestAveraging3D(Tests.OrderTest):
self.getAve = lambda M: M.aveF2CCV
self.orderTest()
def test_orderEx2CC(self):
self.name = "Averaging 3D: Ex2CC"
funX = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
self.getHere = lambda M: np.r_[call3(funX, M.gridEx)]
self.getThere = lambda M: np.r_[call3(funX, M.gridCC)]
self.getAve = lambda M: M.aveEx2CC
self.orderTest()
def test_orderEy2CC(self):
self.name = "Averaging 3D: Ey2CC"
funY = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
self.getHere = lambda M: np.r_[call3(funY, M.gridEy)]
self.getThere = lambda M: np.r_[call3(funY, M.gridCC)]
self.getAve = lambda M: M.aveEy2CC
self.orderTest()
def test_orderEz2CC(self):
self.name = "Averaging 3D: Ez2CC"
funZ = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
self.getHere = lambda M: np.r_[call3(funZ, M.gridEz)]
self.getThere = lambda M: np.r_[call3(funZ, M.gridCC)]
self.getAve = lambda M: M.aveEz2CC
self.orderTest()
def test_orderE2CC(self):
self.name = "Averaging 3D: E2CC"
fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))