wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-29 20:17:38 +08:00
Code Issues Packages Projects Releases Wiki Activity

aveE2CCV for 3D

Browse Source
This commit is contained in:
Lindsey Heagy
2015-11-13 08:32:01 -08:00
parent d130b50f7f
commit 4fe9475ffc
2 changed files with 65 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/TreeMesh.py
+56 -28
View File
@@ -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
tests/mesh/test_TreeOperators.py
+9 -9
View File
@@ -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"