wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-21 12:50:58 +08:00
Code Issues Packages Projects Releases Wiki Activity

aveE2CC for 3D

Browse Source
This commit is contained in:
Lindsey Heagy
2015-11-12 09:11:39 -08:00
parent 6a30ea1d83
commit d130b50f7f
2 changed files with 89 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/TreeMesh.py
+68 -3
View File
@@ -1320,8 +1320,73 @@ class TreeMesh(BaseMesh, InnerProducts):
def aveE2CC(self):
"Construct the averaging operator on cell edges to cell centers."
if getattr(self, '_aveE2CC', None) is None:
# TODO: preallocate
I, J, V = [], [], []
if self.dim == 2:
self._aveE2CC = self.aveF2CC
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
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
@@ -1335,7 +1400,7 @@ class TreeMesh(BaseMesh, InnerProducts):
if getattr(self, '_aveF2CC', None) is None:
# TODO: Preallocate!
I, J, V = [], [], []
PM = [1./(2*self.dim)]*2*self.dim # plus / plus
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
@@ -1451,7 +1516,7 @@ class TreeMesh(BaseMesh, InnerProducts):
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
return self._aveF2CCV
tests/mesh/test_TreeOperators.py
+21 -21
View File
@@ -437,13 +437,13 @@ class TestTreeAveraging2D(Tests.OrderTest):
# 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_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"
@@ -521,13 +521,13 @@ class TestAveraging3D(Tests.OrderTest):
# 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_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"
@@ -539,13 +539,13 @@ class TestAveraging3D(Tests.OrderTest):
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_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"