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

Merged in Averaging (pull request #22)

Browse Source 
Averaging Testing and Fixes.
This commit is contained in:
rowanc1
2013-11-12 14:00:56 -08:00
parent cea1141dcd 7d3e41d591
commit 082583fe51
14 changed files with 1348 additions and 532 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/forward/DCProblem.py
-1
View File
@@ -5,7 +5,6 @@ from SimPEG.utils import ModelBuilder, sdiag, mkvc
from SimPEG import Solver
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as linalg
class DCProblem(ModelTransforms.LogModel, Problem):
SimPEG/mesh/BaseMesh.py
+38 -13
View File
@@ -102,7 +102,7 @@ class BaseMesh(object):
elif xType in ['F', 'E']:
# This will only deal with components of fields, not full 'F' or 'E'
xx = mkvc(xx) # unwrap it in case it is a matrix
nn = self.nF if xType == 'F' else self.nE
nn = self.nFv if xType == 'F' else self.nEv
nn = np.r_[0, nn]
nx = [0, 0, 0]
@@ -331,7 +331,7 @@ class BaseMesh(object):
return locals()
nEz = property(**nEz())
def nE():
def nEv():
doc = """
Total number of edges in each direction
@@ -346,6 +346,18 @@ class BaseMesh(object):
"""
fget = lambda self: np.array([np.prod(x) for x in [self.nEx, self.nEy, self.nEz] if not x is None])
return locals()
nEv = property(**nEv())
def nE():
doc = """
Total number of edges.
:rtype: int
:return: sum([prod(nEx), prod(nEy), prod(nEz)])
"""
fget = lambda self: np.sum(self.nEv)
return locals()
nE = property(**nE())
def nFx():
@@ -391,7 +403,7 @@ class BaseMesh(object):
return locals()
nFz = property(**nFz())
def nF():
def nFv():
doc = """
Total number of faces in each direction
@@ -406,6 +418,19 @@ class BaseMesh(object):
"""
fget = lambda self: np.array([np.prod(x) for x in [self.nFx, self.nFy, self.nFz] if not x is None])
return locals()
nFv = property(**nFv())
def nF():
doc = """
Total number of faces.
:rtype: int
:return: sum([prod(nFx), prod(nFy), prod(nFz)])
"""
fget = lambda self: np.sum(self.nFv)
return locals()
nF = property(**nF())
def normals():
@@ -418,13 +443,13 @@ class BaseMesh(object):
def fget(self):
if self.dim == 2:
nX = np.c_[np.ones(self.nF[0]), np.zeros(self.nF[0])]
nY = np.c_[np.zeros(self.nF[1]), np.ones(self.nF[1])]
nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0])]
nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1])]
return np.r_[nX, nY]
elif self.dim == 3:
nX = np.c_[np.ones(self.nF[0]), np.zeros(self.nF[0]), np.zeros(self.nF[0])]
nY = np.c_[np.zeros(self.nF[1]), np.ones(self.nF[1]), np.zeros(self.nF[1])]
nZ = np.c_[np.zeros(self.nF[2]), np.zeros(self.nF[2]), np.ones(self.nF[2])]
nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0]), np.zeros(self.nFv[0])]
nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1]), np.zeros(self.nFv[1])]
nZ = np.c_[np.zeros(self.nFv[2]), np.zeros(self.nFv[2]), np.ones(self.nFv[2])]
return np.r_[nX, nY, nZ]
return locals()
normals = property(**normals())
@@ -439,13 +464,13 @@ class BaseMesh(object):
def fget(self):
if self.dim == 2:
tX = np.c_[np.ones(self.nE[0]), np.zeros(self.nE[0])]
tY = np.c_[np.zeros(self.nE[1]), np.ones(self.nE[1])]
tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0])]
tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1])]
return np.r_[tX, tY]
elif self.dim == 3:
tX = np.c_[np.ones(self.nE[0]), np.zeros(self.nE[0]), np.zeros(self.nE[0])]
tY = np.c_[np.zeros(self.nE[1]), np.ones(self.nE[1]), np.zeros(self.nE[1])]
tZ = np.c_[np.zeros(self.nE[2]), np.zeros(self.nE[2]), np.ones(self.nE[2])]
tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0]), np.zeros(self.nEv[0])]
tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1]), np.zeros(self.nEv[1])]
tZ = np.c_[np.zeros(self.nEv[2]), np.zeros(self.nEv[2]), np.ones(self.nEv[2])]
return np.r_[tX, tY, tZ]
return locals()
tangents = property(**tangents())
SimPEG/mesh/Cyl1DMesh.py
+6
View File
@@ -110,6 +110,12 @@ class Cyl1DMesh(object):
return locals()
nFz = property(**nFz())
def nFv():
doc = "Total number of faces in each direction"
fget = lambda self: np.array([self.nFr, self.nFz])
return locals()
nFv = property(**nFv())
def nF():
doc = "Total number of faces"
fget = lambda self: self.nFr + self.nFz
SimPEG/mesh/DiffOperators.py
+76 -36
View File
@@ -1,11 +1,6 @@
import numpy as np
from scipy import sparse as sp
from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros
def ddx(n):
"""Define 1D derivatives, inner, this means we go from n+1 to n+1"""
return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr")
from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av
def checkBC(bc):
@@ -39,11 +34,6 @@ def ddxCellGrad(n, bc):
return D
def av(n):
"""Define 1D averaging operator from cell-centres to nodes."""
return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
class DiffOperators(object):
"""
Class creates the differential operators that you need!
@@ -107,6 +97,38 @@ class DiffOperators(object):
_nodalGrad = None
nodalGrad = property(**nodalGrad())
def nodalLaplacian():
doc = "Construct laplacian operator (nodes to edges)."
def fget(self):
if(self._nodalLaplacian is None):
print 'Warning: Laplacian has not been tested rigorously.'
# The number of cell centers in each direction
n = self.n
# Compute divergence operator on faces
if(self.dim == 1):
D1 = sdiag(1./self.hx) * ddx(mesh.nCx)
L = - D1.T*D1
elif(self.dim == 2):
D1 = sdiag(1./self.hx) * ddx(n[0])
D2 = sdiag(1./self.hy) * ddx(n[1])
L1 = sp.kron(speye(n[1]+1), - D1.T * D1)
L2 = sp.kron(- D2.T * D2, speye(n[0]+1))
L = L1 + L2
elif(self.dim == 3):
D1 = sdiag(1./self.hx) * ddx(n[0])
D2 = sdiag(1./self.hy) * ddx(n[1])
D3 = sdiag(1./self.hz) * ddx(n[2])
L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1)
L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1))
L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1))
L = L1 + L2 + L3
self._nodalLaplacian = L
return self._nodalLaplacian
return locals()
_nodalLaplacian = None
nodalLaplacian = property(**nodalLaplacian())
def setCellGradBC(self, BC):
"""
Function that sets the boundary conditions for cell-centred derivative operators.
@@ -182,7 +204,6 @@ class DiffOperators(object):
return locals()
cellGradx = property(**cellGradx())
def cellGrady():
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
def fget(self):
@@ -203,8 +224,6 @@ class DiffOperators(object):
return locals()
cellGrady = property(**cellGrady())
def cellGradz():
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
def fget(self):
@@ -222,7 +241,6 @@ class DiffOperators(object):
return locals()
cellGradz = property(**cellGradz())
def edgeCurl():
doc = "Construct the 3D curl operator."
@@ -265,6 +283,8 @@ class DiffOperators(object):
_edgeCurl = None
edgeCurl = property(**edgeCurl())
# --------------- Averaging ---------------------
def aveF2CC():
doc = "Construct the averaging operator on cell faces to cell centers."
@@ -274,10 +294,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
@@ -295,10 +315,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
@@ -316,37 +336,57 @@ class DiffOperators(object):
if(self.dim == 1):
self._aveN2CC = av(n[0])
elif(self.dim == 2):
self._aveN2CC = sp.hstack((sp.kron(av(n[1]), av(n[0])),
sp.kron(av(n[1]), av(n[0]))), format="csr")
self._aveN2CC = sp.kron(av(n[1]), av(n[0])).tocsr()
elif(self.dim == 3):
self._aveN2CC = sp.hstack((kron3(av(n[2]), av(n[1]), av(n[0])),
kron3(av(n[2]), av(n[1]), av(n[0])),
kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr")
self._aveN2CC = kron3(av(n[2]), av(n[1]), av(n[0])).tocsr()
return self._aveN2CC
return locals()
_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(av(n[1]), av(n[0])),
sp.kron(av(n[1]), av(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(av(n[2]), av(n[1]), av(n[0])),
kron3(av(n[2]), av(n[1]), av(n[0])),
kron3(av(n[2]), av(n[1]), av(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())
# --------------- Methods ---------------------
def getMass(self, materialProp=None, loc='e'):
""" Produces mass matricies.
SimPEG/mesh/InnerProducts.py
+6 -6
View File
@@ -174,8 +174,8 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False):
def Pxxx(pos):
ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]])
ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nF[0]
ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nF[0] + mesh.nF[1]
ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nFv[0]
ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nFv[0] + mesh.nFv[1]
IND = np.r_[ind1, ind2, ind3].flatten()
@@ -286,7 +286,7 @@ def getFaceInnerProduct2D(mesh, mu=None, returnP=False):
def Pxx(pos):
ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1]])
ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nF[0]
ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nFv[0]
IND = np.r_[ind1, ind2].flatten()
@@ -387,8 +387,8 @@ def getEdgeInnerProduct(mesh, sigma=None, returnP=False):
def Pxxx(pos):
ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]])
ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nE[0]
ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nE[0] + mesh.nE[1]
ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nEv[0]
ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nEv[0] + mesh.nEv[1]
IND = np.r_[ind1, ind2, ind3].flatten()
@@ -499,7 +499,7 @@ def getEdgeInnerProduct2D(mesh, sigma=None, returnP=False):
def Pxx(pos):
ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1]])
ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nE[0]
ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nEv[0]
IND = np.r_[ind1, ind2].flatten()
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/mesh/TensorMesh.py
+1 -1
View File
@@ -396,7 +396,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1
if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind:
nF_nE = self.nF if 'F' in locType else self.nE
nF_nE = self.nFv if 'F' in locType else self.nEv
components = [spzeros(loc.shape[0], n) for n in nF_nE]
components[ind] = interpmat(loc, *self.getTensor(locType))
Q = sp.hstack(components)
SimPEG/tests/runTests.py
+43 -36
View File
@@ -4,47 +4,54 @@ import unittest
import HTMLTestRunner
# This code will run all tests in directory named test_*.py
def main(html=False):
TITLE = 'Test Results'
test_file_strings = glob.glob('test_*.py')
module_strings = [str[0:len(str)-3] for str in test_file_strings]
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
in module_strings]
testSuite = unittest.TestSuite(suites)
TITLE = 'Test Results'
test_file_strings = glob.glob('test_*.py')
module_strings = [str[0:len(str)-3] for str in test_file_strings]
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
in module_strings]
testSuite = unittest.TestSuite(suites)
unittest.TextTestRunner(verbosity=2).run(testSuite)
if not html:
unittest.TextTestRunner(verbosity=2).run(testSuite)
return
outfile = open("report.html", "w")
runner = HTMLTestRunner.HTMLTestRunner(
stream=outfile,
title=TITLE,
description='SimPEG Test Report was automatically generated.'
)
outfile = open("report.html", "w")
runner = HTMLTestRunner.HTMLTestRunner(
stream=outfile,
title=TITLE,
description='SimPEG Test Report was automatically generated.',
verbosity=2
)
runner.run(testSuite)
outfile.close()
runner.run(testSuite)
outfile.close()
reader = open("report.html", "r")
writer = open("../../docs/api_TestResults.rst", "w")
reader = open("report.html", "r")
writer = open("../../docs/api_TestResults.rst", "w")
writer.write('.. _api_TestResults:\n\nTest Results\n============\n\n.. raw:: html\n\n')
writer.write('.. _api_TestResults:\n\nTest Results\n============\n\n.. raw:: html\n\n')
go = False
for line in reader:
skip = False
if line == '<style type="text/css" media="screen">\n':
go = True
elif line == "<div id='ending'>&nbsp;</div>\n":
go = False
elif line == '</head>\n':
skip = True
elif line == '<h1>'+TITLE+'</h1>\n':
skip = True
elif line == '<body>\n':
skip = True
if go and not skip:
writer.write(' '+line)
go = False
for line in reader:
skip = False
if line == '<style type="text/css" media="screen">\n':
go = True
elif line == "<div id='ending'>&nbsp;</div>\n":
go = False
elif line == '</head>\n':
skip = True
elif line == '<h1>'+TITLE+'</h1>\n':
skip = True
elif line == '<body>\n':
skip = True
if go and not skip:
writer.write(' '+line)
writer.close()
reader.close()
os.remove("report.html")
writer.close()
reader.close()
os.remove("report.html")
if __name__ == '__main__':
main(True)
SimPEG/tests/test_LogicallyOrthogonalMesh.py
+30 -26
View File
@@ -44,50 +44,54 @@ class BasicLOMTests(unittest.TestCase):
def test_tangents(self):
T = self.LOM2.tangents
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM2.nE[0])))
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM2.nE[0])))
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM2.nE[1])))
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM2.nE[1])))
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM2.nEv[0])))
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM2.nEv[0])))
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM2.nEv[1])))
self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM2.nEv[1])))
T = self.LOM3.tangents
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM3.nE[0])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM3.nE[0])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LOM3.nE[0])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM3.nEv[0])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM3.nEv[0])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LOM3.nEv[0])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM3.nE[1])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM3.nE[1])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LOM3.nE[1])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM3.nEv[1])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM3.nEv[1])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LOM3.nEv[1])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LOM3.nE[2])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LOM3.nE[2])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LOM3.nE[2])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LOM3.nEv[2])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LOM3.nEv[2])))
self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LOM3.nEv[2])))
def test_normals(self):
N = self.LOM2.normals
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM2.nF[0])))
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM2.nF[0])))
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM2.nF[1])))
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM2.nF[1])))
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM2.nFv[0])))
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM2.nFv[0])))
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM2.nFv[1])))
self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM2.nFv[1])))
N = self.LOM3.normals
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM3.nF[0])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM3.nF[0])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LOM3.nF[0])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM3.nFv[0])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM3.nFv[0])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LOM3.nFv[0])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM3.nF[1])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM3.nF[1])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LOM3.nF[1])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM3.nFv[1])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM3.nFv[1])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LOM3.nFv[1])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LOM3.nF[2])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LOM3.nF[2])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LOM3.nF[2])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LOM3.nFv[2])))
self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LOM3.nFv[2])))
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_basemesh.py
+20 -16
View File
@@ -38,15 +38,17 @@ class TestBaseMesh(unittest.TestCase):
def test_mesh_numbers(self):
c = self.mesh.nC == 36
f = np.all(self.mesh.nF == [42, 54, 48])
e = np.all(self.mesh.nE == [72, 56, 63])
fv = np.all(self.mesh.nFv == [42, 54, 48])
ev = np.all(self.mesh.nEv == [72, 56, 63])
f = np.all(self.mesh.nF == np.sum([42, 54, 48]))
e = np.all(self.mesh.nE == np.sum([72, 56, 63]))
self.assertTrue(np.all([c, f, e]))
self.assertTrue(np.all([c, fv, ev, f, e]))
def test_mesh_r_E_V(self):
ex = np.ones(self.mesh.nE[0])
ey = np.ones(self.mesh.nE[1])*2
ez = np.ones(self.mesh.nE[2])*3
ex = np.ones(self.mesh.nEv[0])
ey = np.ones(self.mesh.nEv[1])*2
ez = np.ones(self.mesh.nEv[2])*3
e = np.r_[ex, ey, ez]
tex = self.mesh.r(e, 'E', 'Ex', 'V')
tey = self.mesh.r(e, 'E', 'Ey', 'V')
@@ -60,9 +62,9 @@ class TestBaseMesh(unittest.TestCase):
self.assertTrue(np.all(tez == ez))
def test_mesh_r_F_V(self):
fx = np.ones(self.mesh.nF[0])
fy = np.ones(self.mesh.nF[1])*2
fz = np.ones(self.mesh.nF[2])*3
fx = np.ones(self.mesh.nFv[0])
fy = np.ones(self.mesh.nFv[1])*2
fz = np.ones(self.mesh.nFv[2])*3
f = np.r_[fx, fy, fz]
tfx = self.mesh.r(f, 'F', 'Fx', 'V')
tfy = self.mesh.r(f, 'F', 'Fy', 'V')
@@ -146,14 +148,16 @@ class TestMeshNumbers2D(unittest.TestCase):
def test_mesh_numbers(self):
c = self.mesh.nC == 12
f = np.all(self.mesh.nF == [14, 18])
e = np.all(self.mesh.nE == [18, 14])
fv = np.all(self.mesh.nFv == [14, 18])
ev = np.all(self.mesh.nEv == [18, 14])
f = np.all(self.mesh.nF == np.sum([14, 18]))
e = np.all(self.mesh.nE == np.sum([18, 14]))
self.assertTrue(np.all([c, f, e]))
self.assertTrue(np.all([c, fv, ev, f, e]))
def test_mesh_r_E_V(self):
ex = np.ones(self.mesh.nE[0])
ey = np.ones(self.mesh.nE[1])*2
ex = np.ones(self.mesh.nEv[0])
ey = np.ones(self.mesh.nEv[1])*2
e = np.r_[ex, ey]
tex = self.mesh.r(e, 'E', 'Ex', 'V')
tey = self.mesh.r(e, 'E', 'Ey', 'V')
@@ -165,8 +169,8 @@ class TestMeshNumbers2D(unittest.TestCase):
self.assertRaises(AssertionError, self.mesh.r, e, 'E', 'Ez', 'V')
def test_mesh_r_F_V(self):
fx = np.ones(self.mesh.nF[0])
fy = np.ones(self.mesh.nF[1])*2
fx = np.ones(self.mesh.nFv[0])
fy = np.ones(self.mesh.nFv[1])*2
f = np.r_[fx, fy]
tfx = self.mesh.r(f, 'F', 'Fx', 'V')
tfy = self.mesh.r(f, 'F', 'Fy', 'V')
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()
SimPEG/utils/__init__.py
+1 -1
View File
@@ -4,7 +4,7 @@ import lomutils
import interputils
import ModelBuilder
from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
from sputils import spzeros, kron3, speye, sdiag
from sputils import spzeros, kron3, speye, sdiag, ddx, av
from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird
from interputils import interpmat
from ipythonUtils import easyAnimate as animate
SimPEG/utils/sputils.py
+11
View File
@@ -1,5 +1,6 @@
from scipy import sparse as sp
from matutils import mkvc
import numpy as np
def sdiag(h):
@@ -20,3 +21,13 @@ def kron3(A, B, C):
def spzeros(n1, n2):
"""spzeros"""
return sp.coo_matrix((n1, n2)).tocsr()
def ddx(n):
"""Define 1D derivatives, inner, this means we go from n+1 to n"""
return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr")
def av(n):
"""Define 1D averaging operator from cell-centres to nodes."""
return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
docs/api_TestResults.rst
+1012 -394
View File
File diff suppressed because it is too large Load Diff