From ef557beb84ce3d6f39f216e011d2a0fdf0d24cc8 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 18 Jul 2013 10:41:58 -0700 Subject: [PATCH] Merged Seogi's code into the Tensor Mesh code base. This is under the DiffOperators class, which can be inherited with BaseMesh to any Mesh Object. Cell area and other dimension calculations are included in the TensorMesh class. Wrote unit tests for cell vol, area, and edges. --- SimPEG/DiffOperators.py | 156 +++++++++++++++++++++++++ SimPEG/TensorMesh.py | 79 ++++++++++++- SimPEG/__init__.py | 1 + SimPEG/getDiffop.py | 198 -------------------------------- SimPEG/tests/test_curl.py | 19 ++- SimPEG/tests/test_div.py | 19 ++- SimPEG/tests/test_grad.py | 21 ++-- SimPEG/tests/test_tensorMesh.py | 30 ++++- 8 files changed, 289 insertions(+), 234 deletions(-) create mode 100644 SimPEG/DiffOperators.py delete mode 100644 SimPEG/getDiffop.py 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()