mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-21 12:50:58 +08:00
Averaging from faces to cell centres (2nd order on non-refined mesh, 1st order when we refine... which I think is ok)
This commit is contained in:
Lindsey Heagy
@@ -1328,6 +1328,62 @@ class TreeMesh(BaseMesh, InnerProducts):
|
||||
|
||||
@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))
|
||||
R = 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
|
||||
return self._aveF2CC
|
||||
|
||||
|
||||
|
||||
@property
|
||||
def aveF2CCV(self):
|
||||
"Construct the averaging operator on cell faces to cell centers."
|
||||
raise Exception('Not yet implemented!')
|
||||
|
||||
|
||||
@@ -14,6 +14,8 @@ cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_
|
||||
cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez)))
|
||||
|
||||
|
||||
plotit = False
|
||||
|
||||
class TestFaceDiv2D(Tests.OrderTest):
|
||||
name = "Face Divergence 2D"
|
||||
meshTypes = MESHTYPES
|
||||
@@ -388,6 +390,176 @@ class TestTreeInnerProducts2D(Tests.OrderTest):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
class TestTreeAveraging2D(Tests.OrderTest):
|
||||
"""Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
meshTypes = ['uniformTree', 'randomTree']
|
||||
meshDimension = 2
|
||||
meshSizes = [4,8,16,32]
|
||||
|
||||
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)
|
||||
|
||||
if plotit:
|
||||
self.M.plotImage(self.getThere(self.M)-num)
|
||||
plt.show()
|
||||
plt.tight_layout
|
||||
|
||||
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_orderN2F(self):
|
||||
# self.name = "Averaging 2D: N2F"
|
||||
# fun = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
# self.getHere = lambda M: call2(fun, M.gridN)
|
||||
# self.getThere = lambda M: np.r_[call2(fun, M.gridFx), call2(fun, M.gridFy)]
|
||||
# self.getAve = lambda M: M.aveN2F
|
||||
# self.orderTest()
|
||||
|
||||
# def test_orderN2E(self):
|
||||
# self.name = "Averaging 2D: N2E"
|
||||
# fun = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
# self.getHere = lambda M: call2(fun, M.gridN)
|
||||
# self.getThere = lambda M: np.r_[call2(fun, M.gridEx), call2(fun, M.gridEy)]
|
||||
# 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
|
||||
# self.orderTest()
|
||||
|
||||
# def test_orderCC2F(self):
|
||||
# self.name = "Averaging 2D: CC2F"
|
||||
# fun = lambda x, y: (np.cos(x)+np.sin(y))
|
||||
# self.getHere = lambda M: call2(fun, M.gridCC)
|
||||
# self.getThere = lambda M: np.r_[call2(fun, M.gridFx), call2(fun, M.gridFy)]
|
||||
# self.getAve = lambda M: M.aveCC2F
|
||||
# self.expectedOrders = 1
|
||||
# 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, M.gridEx), call2(fun, 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 = ['uniformTree', 'randomTree']
|
||||
meshDimension = 3
|
||||
meshSizes = [8,16]
|
||||
|
||||
def getError(self):
|
||||
num = self.getAve(self.M) * self.getHere(self.M)
|
||||
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_orderN2F(self):
|
||||
# self.name = "Averaging 3D: N2F"
|
||||
# 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: np.r_[call3(fun, M.gridFx), call3(fun, M.gridFy), call3(fun, M.gridFz)]
|
||||
# self.getAve = lambda M: M.aveN2F
|
||||
# self.orderTest()
|
||||
|
||||
# def test_orderN2E(self):
|
||||
# self.name = "Averaging 3D: N2E"
|
||||
# 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: np.r_[call3(fun, M.gridEx), call3(fun, M.gridEy), call3(fun, M.gridEz)]
|
||||
# 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_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"
|
||||
# fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
# self.getHere = lambda M: np.r_[call3(fun, M.gridEx), call3(fun, M.gridEy), call3(fun, M.gridEz)]
|
||||
# self.getThere = lambda M: call3(fun, M.gridCC)
|
||||
# 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_orderCC2F(self):
|
||||
# self.name = "Averaging 3D: CC2F"
|
||||
# fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
|
||||
# self.getHere = lambda M: call3(fun, M.gridCC)
|
||||
# self.getThere = lambda M: np.r_[call3(fun, M.gridFx), call3(fun, M.gridFy), call3(fun, M.gridFz)]
|
||||
# self.getAve = lambda M: M.aveCC2F
|
||||
# self.expectedOrders = 1
|
||||
# self.orderTest()
|
||||
# self.expectedOrders = 2
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
Reference in New Issue
Block a user