diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py new file mode 100644 index 00000000..f240a44f --- /dev/null +++ b/SimPEG/DiffOperators.py @@ -0,0 +1,156 @@ +import numpy as np +from scipy import sparse +from sputils import ddx, sdiag, speye, kron3, spzeros, av + + +class DiffOperators(object): + """ + Class creates the differential operators that you need! + """ + def __init__(self): + raise Exception('You should use a Mesh class.') + + def DIV(): + doc = "Construct the 3D divergence operator on Faces." + + def fget(self): + if(self._DIV is None): + # The number of cell centers in each direction + n = [x.size for x in self.h] + # Compute divergence operator on faces + dd = [ddx(x) for x in n] + if(self.dim == 1): + D = dd[0] + elif(self.dim == 2): + D1 = sparse.kron(speye(n[1]), dd[0]) + D2 = sparse.kron(dd[1], speye(n[0])) + D = sparse.hstack((D1, D2), format="csr") + elif(self.dim == 3): + D1 = kron3(speye(n[2]), speye(n[1]), dd[0]) + D2 = kron3(speye(n[2]), dd[1], speye(n[0])) + D3 = kron3(dd[2], speye(n[1]), speye(n[0])) + D = sparse.hstack((D1, D2, D3), format="csr") + # Compute areas of cell faces + S = self.area + # Compute cell volumes + V = self.vol + self._DIV = sdiag(1/V)*D*sdiag(S) + + return self._DIV + return locals() + _DIV = None + DIV = property(**DIV()) + + def GRAD(): + doc = "Construct the 3D nodal gradient operator." + + def fget(self): + if(self._GRAD is None): + # The number of cell centers in each direction + n1 = np.size(self.hx) + n2 = np.size(self.hy) + n3 = np.size(self.hz) + + # Compute lengths of cell edges + L = self.edge + + # Compute divergence operator on faces + d1 = ddx(n1) + d2 = ddx(n2) + d3 = ddx(n3) + D1 = kron3(speye(n3+1), speye(n2+1), d1) + D2 = kron3(speye(n3+1), d2, speye(n1+1)) + D3 = kron3(d3, speye(n2+1), speye(n1+1)) + + G = sparse.vstack((D1, D2, D3), format="csr") + self._GRAD = sdiag(1/L)*G + return self._GRAD + return locals() + _GRAD = None + GRAD = property(**GRAD()) + + def CURL(): + doc = "Construct the 3D curl operator." + + def fget(self): + if(self._CURL is None): + # The number of cell centers in each direction + n1 = np.size(self.hx) + n2 = np.size(self.hy) + n3 = np.size(self.hz) + + # Compute lengths of cell edges + L = self.edge + + # Compute areas of cell faces + S = self.area + + # Compute divergence operator on faces + d1 = ddx(n1) + d2 = ddx(n2) + d3 = ddx(n3) + + D32 = kron3(d3, speye(n2), speye(n1+1)) + D23 = kron3(speye(n3), d2, speye(n1+1)) + D31 = kron3(d3, speye(n2+1), speye(n1)) + D13 = kron3(speye(n3), speye(n2+1), d1) + D21 = kron3(speye(n3+1), d2, speye(n1)) + D12 = kron3(speye(n3+1), speye(n2), d1) + + O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1]) + O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) + O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) + + C = sparse.vstack((sparse.hstack((O1, -D32, D23)), + sparse.hstack((D31, O2, -D13)), + sparse.hstack((-D21, D12, O3))), format="csr") + + self._CURL = sdiag(1/S)*(C*sdiag(L)) + return self._CURL + return locals() + _CURL = None + CURL = property(**CURL()) + + def AVE_F(): + doc = "Construct the 3D averaging operator on cell faces." + + def fget(self): + if(self._AVE_F is None): + # The number of cell centers in each direction + n1 = np.size(self.hx) + n2 = np.size(self.hy) + n3 = np.size(self.hz) + + av1 = av(n1) + av2 = av(n2) + av3 = av(n3) + + self._AVE_F = sparse.hstack(kron3(speye(n3), speye(n2), av1), + kron3(speye(n3), av2, speye(n3)), + kron3(av3, speye(n2), speye(n3)), format="csr") + return self._AVE_F + return locals() + _AVE_F = None + AVE_F = property(**AVE_F()) + + def AVE_E(): + doc = "Construct the 3D averaging operator on cell edges." + + def fget(self): + if(self._AVE_E is None): + # The number of cell centers in each direction + n1 = np.size(self.hx) + n2 = np.size(self.hy) + n3 = np.size(self.hz) + + av1 = av(n1) + av2 = av(n2) + av3 = av(n3) + + self._AVE_E = sparse.hstack(kron3(av3, av2, speye(n1)), + kron3(av3, speye(n2), av1), + kron3(speye(n3), av2, av1), format="csr") + return self._AVE_E + return locals() + _AVE_E = None + AVE_E = property(**AVE_E()) diff --git a/SimPEG/TensorMesh.py b/SimPEG/TensorMesh.py index 5f5e44b3..023362ae 100644 --- a/SimPEG/TensorMesh.py +++ b/SimPEG/TensorMesh.py @@ -1,10 +1,11 @@ import numpy as np from BaseMesh import BaseMesh from TensorView import TensorView -from utils import ndgrid +from DiffOperators import DiffOperators +from utils import ndgrid, mkvc -class TensorMesh(BaseMesh, TensorView): +class TensorMesh(BaseMesh, TensorView, DiffOperators): """ TensorMesh is a mesh class that deals with tensor product meshes. @@ -186,6 +187,80 @@ class TensorMesh(BaseMesh, TensorView): def getCellNumbering(self): pass + # --------------- Geometries --------------------- + def vol(): + doc = "Construct cell volumes of the 3D model as 1d array." + + def fget(self): + if(self._vol is None): + vh = [mkvc(x, 2) for x in self.h] + # Compute cell volumes + if(self.dim == 1): + self._vol = mkvc(vh[0]) + elif(self.dim == 2): + # Cell sizes in each direction + self._vol = mkvc(vh[0]*vh[1].T) + elif(self.dim == 3): + # Cell sizes in each direction + v12 = vh[0]*vh[1].T + self._vol = mkvc(mkvc(v12, 2)*vh[2].T) + return self._vol + return locals() + _vol = None + vol = property(**vol()) + + def area(): + doc = "Construct face areas of the 3D model as 1d array." + + def fget(self): + if(self._area is None): + # Ensure that we are working with column vectors + vh = [mkvc(x, 2) for x in self.h] + # The number of cell centers in each direction + n = [x.size for x in self.h] + # Compute areas of cell faces + if(self.dim == 1): + self._area = np.ones((n[0]+1, 1)) + elif(self.dim == 2): + area1 = np.ones((n[0]+1, 1))*vh[1].T + area2 = vh[0]*np.ones((n[1]+1, 1)).T + self._area = np.r_[mkvc(area1), mkvc(area2)] + elif(self.dim == 3): + area1 = np.ones((n[0]+1, 1))*mkvc(vh[1]*vh[2].T) + area2 = vh[0]*mkvc(np.ones((n[1]+1, 1))*vh[2].T) + area3 = vh[0]*mkvc(vh[1]*np.ones((n[2]+1, 1)).T) + self._area = np.r_[mkvc(area1), mkvc(area2), mkvc(area3)] + return self._area + return locals() + _area = None + area = property(**area()) + + def edge(): + doc = "Construct edge legnths of the 3D model as 1d array." + + def fget(self): + if(self._area is None): + # Ensure that we are working with column vectors + vh = [mkvc(x, 2) for x in self.h] + # The number of cell centers in each direction + n = [x.size for x in self.h] + # Compute edge lengths + if(self.dim == 1): + self._edge = mkvc(vh[0]) + elif(self.dim == 2): + l1 = vh[0]*np.ones((n[1]+1, 1)).T + l2 = np.ones((n[0]+1, 1))*vh[1].T + self._edge = np.r_[mkvc(l1), mkvc(l2)] + elif(self.dim == 3): + l1 = vh[0]*mkvc(np.ones((n[1]+1, 1))*np.ones((n[2]+1, 1)).T) + l2 = np.ones((n[0]+1, 1))*mkvc(vh[1]*np.ones((n[2]+1, 1)).T) + l3 = np.ones((n[0]+1, 1))*mkvc(np.ones((n[1]+1, 1))*vh[2].T) + self._edge = np.r_[mkvc(l1), mkvc(l2), mkvc(l3)] + return self._edge + return locals() + _edge = None + edge = property(**edge()) + if __name__ == '__main__': print('Welcome to tensor mesh!') diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 2d1b1e99..bcbba681 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -1 +1,2 @@ from TensorMesh import TensorMesh +import utils diff --git a/SimPEG/getDiffop.py b/SimPEG/getDiffop.py deleted file mode 100644 index 6debc348..00000000 --- a/SimPEG/getDiffop.py +++ /dev/null @@ -1,198 +0,0 @@ -import numpy as np -from scipy import sparse -from utils import mkvc -from sputils import ddx, sdiag, speye, kron3, spzeros, av - -def getvol(h): - """Construct cell volumes of the 3D model as 1d array.""" - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # Compute cell volumes - v12 = h1.T*h2 - V = mkvc(v12.reshape(-1,1)*h3) - - return V - -def getarea(h): - """Construct face areas of the 3D model as 1d array.""" - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - # Compute areas of cell faces - area1 = np.ones((n1+1,1))*mkvc(h2.T*h3) - area2 = h1.T*mkvc(np.ones((n2+1,1))*h3) - area3 = h1.T*mkvc(h2.T*np.ones(n3+1)) - area = np.concatenate((mkvc(area1), mkvc(area2), mkvc(area3)), axis=0) - - return area - -def getlength_e(h): - """Construct edge legnths of the 3D model as 1d array.""" - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - # Compute areas of cell faces - l1 = h1.T*mkvc(np.ones((n2+1,1))*np.ones(n3+1)) - l2 = np.ones((n1+1,1))*mkvc(h2.T*np.ones(n3+1)) - l3 = np.ones((n1+1,1))*mkvc(np.ones((n2+1,1))*h3) - #l = np.hstack((np.hstack((mkvc(area1), mkvc(area2))), mkvc(area3))) - l = np.concatenate((mkvc(l1), mkvc(l2), mkvc(l3)), axis=0) - - return l - -def getDivMatrix(h): - """Construct the 3D divergence operator on Faces.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - # Compute areas of cell faces - S = getarea(h) - - # Compute cell volumes - V = getvol(h) - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - D1 = kron3(speye(n3), speye(n2), d1) - D2 = kron3(speye(n3), d2, speye(n1)) - D3 = kron3(d3, speye(n2), speye(n1)) - - D = sparse.hstack((D1, D2, D3), format="csr") - return sdiag(1/V)*D*sdiag(S) - -def getGradMatrix(h): - """Construct the 3D nodal gradient operator.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - # Compute lengths of cell edges - L = getlength_e(h) - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - D1 = kron3(speye(n3+1), speye(n2+1), d1) - D2 = kron3(speye(n3+1), d2, speye(n1+1)) - D3 = kron3(d3, speye(n2+1), speye(n1+1)) - - G = sparse.vstack((D1, D2, D3), format="csr") - return sdiag(1/L)*G - -def getCurlMatrix(h): - """Construct the 3D curl operator.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - # Compute lengths of cell edges - L = getlength_e(h) - - # Compute areas of cell faces - S = getarea(h) - - # Compute divergence operator on faces - d1 = ddx(n1) - d2 = ddx(n2) - d3 = ddx(n3) - - D32 = kron3(d3, speye(n2), speye(n1+1)) - D23 = kron3(speye(n3), d2, speye(n1+1)) - D31 = kron3(d3, speye(n2+1), speye(n1)) - D13 = kron3(speye(n3), speye(n2+1), d1) - D21 = kron3(speye(n3+1), d2, speye(n1)) - D12 = kron3(speye(n3+1), speye(n2), d1) - - O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1]) - O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1]) - O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1]) - - C = sparse.vstack((sparse.hstack((O1,-D32, D23)), - sparse.hstack((D31,O2, -D13)), - sparse.hstack((-D21,D12, O3))), format="csr") - - return sdiag(1/S)*(C*sdiag(L)) - -def getAverageMatrixF(h): - """Construct the 3D averaging operator on cell faces.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - av1 = av(n1) - av2 = av(n2) - av3 = av(n3) - - AvF = sparse.hstack(kron3(speye(n3), speye(n2), av1), - kron3(speye(n3), av2, speye(n3)), - kron3(av3, speye(n2), speye(n3)), format="csr") - return AvF - -def getAverageMatrixE(h): - """Construct the 3D averaging operator on cell edges.""" - - # Cell sizes in each direction - h1 = h[0] - h2 = h[1] - h3 = h[2] - - # The number of cell centers in each direction - n1 = np.size(h1) - n2 = np.size(h2) - n3 = np.size(h3) - - av1 = av(n1) - av2 = av(n2) - av3 = av(n3) - - AvE = sparse.hstack(kron3(av3, av2, speye(n1)), - kron3(av3, speye(n2), av1), - kron3(speye(n3), av2, av1), format="csr") - return AvE \ No newline at end of file diff --git a/SimPEG/tests/test_curl.py b/SimPEG/tests/test_curl.py index 8175ab2f..e4d9812c 100644 --- a/SimPEG/tests/test_curl.py +++ b/SimPEG/tests/test_curl.py @@ -3,7 +3,6 @@ import numpy as np import sys sys.path.append('../') from TensorMesh import TensorMesh -from getDiffop import getCurlMatrix err=0. @@ -11,7 +10,7 @@ print '>> Test Curl operator' for i in range(4): icount=i+1 nc = 2**icount - # Define the mesh + # Define the mesh h1 = np.ones((1,nc))/nc h2 = np.ones((1,nc))/nc h3 = np.ones((1,nc))/nc @@ -21,28 +20,28 @@ for i in range(4): #n = M.plotGrid() # Generate DIV matrix - CURL = getCurlMatrix(h) + CURL = M.CURL #Test function fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) - + Ex = fun(M.gridEx[:,1]) Ey = fun(M.gridEy[:,2]) Ez = fun(M.gridEz[:,0]) - E = np.concatenate((Ex,Ey,Ez)) + E = np.concatenate((Ex,Ey,Ez)) Fx = sol(M.gridFx[:,2]) Fy = sol(M.gridFy[:,0]) Fz = sol(M.gridFz[:,1]) - curlE_anal = np.concatenate((Fx,Fy,Fz)) + curlE_anal = np.concatenate((Fx,Fy,Fz)) - curlE = CURL*E + curlE = CURL*E err = np.linalg.norm((curlE-curlE_anal), np.inf) if icount == 1: - print 'h | inf norm | error ratio' - print '---------------------------------------' + print 'h | inf norm | error ratio' + print '---------------------------------------' print '%6.4f | %8.2e |'% (h1[0,0], err) else: print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) - err_old = err \ No newline at end of file + err_old = err \ No newline at end of file diff --git a/SimPEG/tests/test_div.py b/SimPEG/tests/test_div.py index 3ab8f6e2..35f4e408 100644 --- a/SimPEG/tests/test_div.py +++ b/SimPEG/tests/test_div.py @@ -3,7 +3,6 @@ import numpy as np import sys sys.path.append('../') from TensorMesh import TensorMesh -from getDIV import getDivMatrix, getarea, getvol err=0. @@ -11,7 +10,7 @@ print '>> Test face Divergence operator' for i in range(4): icount=i+1 nc = 2**icount - # Define the mesh + # Define the mesh h1 = np.ones((1,nc))/nc h2 = np.ones((1,nc))/nc h3 = np.ones((1,nc))/nc @@ -21,28 +20,26 @@ for i in range(4): #n = M.plotGrid() # Generate DIV matrix - DIV = getDivMatrix(h) - + DIV = M.DIV + #Test function fun = lambda x: np.sin(x) Fx = fun(M.gridFx[:,0]) Fy = fun(M.gridFy[:,1]) Fz = fun(M.gridFz[:,2]) - + F = np.concatenate((Fx,Fy,Fz)) divF = DIV*F sol = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) divF_anal = sol(M.gridCC[:,0], M.gridCC[:,1], M.gridCC[:,2]) - - area = getarea(h) - vol = getvol(h) + #err = np.linalg.norm((divF-divF_anal)*np.sqrt(vol), 2) err = np.linalg.norm((divF-divF_anal), np.inf) if icount == 1: - print 'h | inf norm | error ratio' - print '---------------------------------------' + print 'h | inf norm | error ratio' + print '---------------------------------------' print '%6.4f | %8.2e |'% (h1[0,0], err) else: print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) - err_old = err + err_old = err diff --git a/SimPEG/tests/test_grad.py b/SimPEG/tests/test_grad.py index a433c28d..245d9111 100644 --- a/SimPEG/tests/test_grad.py +++ b/SimPEG/tests/test_grad.py @@ -3,15 +3,14 @@ import numpy as np import sys sys.path.append('../') from TensorMesh import TensorMesh -from getDiffop import getGradMatrix err=0. -print '>> Test nodal Gradient operator' +print '>> Test nodal Gradient operator' for i in range(4): icount=i+1 nc = 2**icount - # Define the mesh + # Define the mesh h1 = np.ones((1,nc))/nc h2 = np.ones((1,nc))/nc h3 = np.ones((1,nc))/nc @@ -21,11 +20,11 @@ for i in range(4): #n = M.plotGrid() # Generate DIV matrix - GRAD = getGradMatrix(h) + GRAD = M.GRAD #Test function - fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) + fun = lambda x, y, z: (np.cos(x)+np.cos(y)+np.cos(z)) sol = lambda x: -np.sin(x) # i (sin(x)) + j (sin(y)) + k (sin(z)) - + phi = fun(M.gridN[:,0], M.gridN[:,1], M.gridN[:,2]) gradE = GRAD*phi @@ -33,14 +32,14 @@ for i in range(4): Ey = sol(M.gridEy[:,1]) Ez = sol(M.gridEz[:,2]) - gradE_anal = np.concatenate((Ex,Ey,Ez)) + gradE_anal = np.concatenate((Ex,Ey,Ez)) err = np.linalg.norm((gradE-gradE_anal), np.inf) if icount == 1: - print 'h | inf norm | error ratio' - print '---------------------------------------' + print 'h | inf norm | error ratio' + print '---------------------------------------' print '%6.4f | %8.2e |'% (h1[0,0], err) else: print '%6.4f | %8.2e | %6.4f' % (h1[0,0], err, err_old/err) - err_old = err - + err_old = err + diff --git a/SimPEG/tests/test_tensorMesh.py b/SimPEG/tests/test_tensorMesh.py index e1b74374..2a6cb0a1 100644 --- a/SimPEG/tests/test_tensorMesh.py +++ b/SimPEG/tests/test_tensorMesh.py @@ -10,8 +10,9 @@ class TestSequenceFunctions(unittest.TestCase): def setUp(self): a = np.array([1, 1, 1]) b = np.array([1, 2]) - x0 = np.array([3, 5]) - self.mesh2 = TensorMesh([a, b], x0) + c = np.array([1, 4]) + self.mesh2 = TensorMesh([a, b], np.array([3, 5])) + self.mesh3 = TensorMesh([a, b, c]) def test_vectorN_2D(self): testNx = np.array([3, 4, 5, 6]) @@ -29,6 +30,31 @@ class TestSequenceFunctions(unittest.TestCase): ytest = np.all(self.mesh2.vectorCCy == testNy) self.assertTrue(xtest and ytest) + def test_area_3D(self): + test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2]) + t1 = np.all(self.mesh3.area == test_area) + self.assertTrue(t1) + + def test_vol_3D(self): + test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8]) + t1 = np.all(self.mesh3.vol == test_vol) + self.assertTrue(t1) + + def test_vol_2D(self): + test_vol = np.array([1, 1, 1, 2, 2, 2]) + t1 = np.all(self.mesh2.vol == test_vol) + self.assertTrue(t1) + + def test_edge_3D(self): + test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]) + t1 = np.all(self.mesh3.edge == test_edge) + self.assertTrue(t1) + + def test_edge_2D(self): + test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2]) + t1 = np.all(self.mesh2.edge == test_edge) + self.assertTrue(t1) + if __name__ == '__main__': unittest.main()