diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index 9ddb5332..34046a52 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -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): diff --git a/SimPEG/mesh/BaseMesh.py b/SimPEG/mesh/BaseMesh.py index 6a9a8032..2647cf44 100644 --- a/SimPEG/mesh/BaseMesh.py +++ b/SimPEG/mesh/BaseMesh.py @@ -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()) diff --git a/SimPEG/mesh/Cyl1DMesh.py b/SimPEG/mesh/Cyl1DMesh.py index 93b82b25..fa2bb729 100644 --- a/SimPEG/mesh/Cyl1DMesh.py +++ b/SimPEG/mesh/Cyl1DMesh.py @@ -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 diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 598b392a..fa3f42e0 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -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. diff --git a/SimPEG/mesh/InnerProducts.py b/SimPEG/mesh/InnerProducts.py index f0ac4ab0..5c2e7b5f 100644 --- a/SimPEG/mesh/InnerProducts.py +++ b/SimPEG/mesh/InnerProducts.py @@ -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() diff --git a/SimPEG/mesh/LogicallyOrthogonalMesh.py b/SimPEG/mesh/LogicallyOrthogonalMesh.py index b510a754..5c4a73db 100644 --- a/SimPEG/mesh/LogicallyOrthogonalMesh.py +++ b/SimPEG/mesh/LogicallyOrthogonalMesh.py @@ -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 diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index 6cc5d8a6..90f8fd10 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -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) diff --git a/SimPEG/tests/runTests.py b/SimPEG/tests/runTests.py index c44f9a1e..58d94303 100644 --- a/SimPEG/tests/runTests.py +++ b/SimPEG/tests/runTests.py @@ -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 == '