wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-05 22:37:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Tested averaging operators.

Browse Source
This commit is contained in:
Rowan Cockett
2013-11-06 12:43:17 -08:00
parent 3f097a106d
commit 24b1307d8a
5 changed files with 1125 additions and 373 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/mesh/DiffOperators.py
+37 -16
View File
@@ -264,10 +264,10 @@ class DiffOperators(object):
if(self.dim == 1):
self._aveF2CC = av(n[0])
elif(self.dim == 2):
self._aveF2CC = sp.hstack((sp.kron(speye(n[1]), av(n[0])),
self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])),
sp.kron(av(n[1]), speye(n[0]))), format="csr")
elif(self.dim == 3):
self._aveF2CC = sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
kron3(speye(n[2]), av(n[1]), speye(n[0])),
kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
return self._aveF2CC
@@ -285,10 +285,10 @@ class DiffOperators(object):
if(self.dim == 1):
raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
elif(self.dim == 2):
self._aveE2CC = sp.hstack((sp.kron(av(n[1]), speye(n[0])),
self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])),
sp.kron(speye(n[1]), av(n[0]))), format="csr")
elif(self.dim == 3):
self._aveE2CC = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
kron3(av(n[2]), speye(n[1]), av(n[0])),
kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
return self._aveE2CC
@@ -314,26 +314,47 @@ class DiffOperators(object):
_aveN2CC = None
aveN2CC = property(**aveN2CC())
def aveN2CCv():
doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension separate."
def aveN2E():
doc = "Construct the averaging operator on cell nodes to cell edges, keeping each dimension separate."
def fget(self):
if(self._aveN2CCv is None):
if(self._aveN2E is None):
# The number of cell centers in each direction
n = self.n
if(self.dim == 1):
self._aveN2CCv = av(n[0])
self._aveN2E = av(n[0])
elif(self.dim == 2):
self._aveN2CCv = sp.block_diag((sp.kron(speye(n[1]), av(n[0])),
sp.kron(av(n[1]), speye(n[0]))), format="csr")
self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])),
sp.kron(av(n[1]), speye(n[0]+1))), format="csr")
elif(self.dim == 3):
self._aveN2CCv = sp.block_diag((kron3(speye(n[2]), speye(n[1]), av(n[0])),
kron3(speye(n[2]), av(n[1]), speye(n[0])),
kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
return self._aveN2CCv
self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), av(n[0])),
kron3(speye(n[2]+1), av(n[1]), speye(n[0]+1)),
kron3(av(n[2]), speye(n[1]+1), speye(n[0]+1))), format="csr")
return self._aveN2E
return locals()
_aveN2CCv = None
aveN2CCv = property(**aveN2CCv())
_aveN2E = None
aveN2E = property(**aveN2E())
def aveN2F():
doc = "Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate."
def fget(self):
if(self._aveN2F is None):
# The number of cell centers in each direction
n = self.n
if(self.dim == 1):
self._aveN2F = av(n[0])
elif(self.dim == 2):
self._aveN2F = sp.vstack((sp.kron(av(n[1]), speye(n[0]+1)),
sp.kron(speye(n[1]+1), av(n[0]))), format="csr")
elif(self.dim == 3):
self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), speye(n[0]+1)),
kron3(av(n[2]), speye(n[1]+1), av(n[0])),
kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr")
return self._aveN2F
return locals()
_aveN2F = None
aveN2F = property(**aveN2F())
def getMass(self, materialProp=None, loc='e'):
""" Produces mass matricies.
SimPEG/mesh/LogicallyOrthogonalMesh.py
+1 -2
View File
@@ -45,8 +45,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
def fget(self):
if self._gridCC is None:
ccV = (self.aveN2CCv*mkvc(self.gridN))
self._gridCC = ccV.reshape((-1, self.dim), order='F')
self._gridCC = np.concatenate([self.aveN2CC*self.gridN[:,i] for i in range(self.dim)]).reshape((-1,self.dim), order='F')
return self._gridCC
return locals()
_gridCC = None # Store grid by default
SimPEG/tests/test_LogicallyOrthogonalMesh.py
+4
View File
@@ -83,11 +83,15 @@ class BasicLOMTests(unittest.TestCase):
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LOM3.nFv[2])))
def test_grid(self):
self.assertTrue(np.all(self.LOM2.gridCC == self.TM2.gridCC))
self.assertTrue(np.all(self.LOM2.gridN == self.TM2.gridN))
self.assertTrue(np.all(self.LOM2.gridFx == self.TM2.gridFx))
self.assertTrue(np.all(self.LOM2.gridFy == self.TM2.gridFy))
self.assertTrue(np.all(self.LOM2.gridEx == self.TM2.gridEx))
self.assertTrue(np.all(self.LOM2.gridEy == self.TM2.gridEy))
self.assertTrue(np.all(self.LOM3.gridCC == self.TM3.gridCC))
self.assertTrue(np.all(self.LOM3.gridN == self.TM3.gridN))
self.assertTrue(np.all(self.LOM3.gridFx == self.TM3.gridFx))
self.assertTrue(np.all(self.LOM3.gridFy == self.TM3.gridFy))
self.assertTrue(np.all(self.LOM3.gridFz == self.TM3.gridFz))
SimPEG/tests/test_operators.py
+103
View File
@@ -155,6 +155,109 @@ class TestNodalGrad2D(OrderTest):
def test_order(self):
self.orderTest()
class TestAveraging2D(OrderTest):
name = "Averaging 2D"
meshTypes = MESHTYPES
meshDimension = 2
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 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, M.gridFx), call2(fun, M.gridFy)]
self.getThere = lambda M: call2(fun, M.gridCC)
self.getAve = lambda M: M.aveF2CC
self.orderTest()
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()
class TestAveraging3D(OrderTest):
name = "Averaging 3D"
meshTypes = MESHTYPES
meshDimension = 3
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_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()
if __name__ == '__main__':
unittest.main()
docs/api_TestResults.rst
+980 -355
View File
File diff suppressed because it is too large Load Diff