diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 9ab89e64..da0927e5 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -136,6 +136,45 @@ class CylMesh(BaseTensorMesh): """Nodal grid vector (1D) in the y direction.""" return np.r_[0, self.hy[:-1].cumsum()] + self.hy[0]*0.5 + + def getTensor(self, locType): + """ Returns a tensor list. + + :param str locType: What tensor (see below) + :rtype: list + :return: list of the tensors that make up the mesh. + + locType can be:: + + 'Ex' -> x-component of field defined on edges + 'Ey' -> y-component of field defined on edges + 'Ez' -> z-component of field defined on edges + 'Fx' -> x-component of field defined on faces + 'Fy' -> y-component of field defined on faces + 'Fz' -> z-component of field defined on faces + 'N' -> scalar field defined on nodes + 'CC' -> scalar field defined on cell centers + """ + + if locType is 'Fx': + ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] + elif locType is 'Fy': + ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] + elif locType is 'Fz': + ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] + elif locType is 'Ex': + ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] + elif locType is 'Ey': + ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] + elif locType is 'Ez': + ten = [self.vectorNx , self.vectorNy , self.vectorCCz] + elif locType is 'CC': + ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] + elif locType is 'N': + ten = [self.vectorNx , self.vectorNy , self.vectorNz ] + + return [t for t in ten if t is not None] + @property def edge(self): """Edge lengths""" @@ -219,6 +258,13 @@ class CylMesh(BaseTensorMesh): self._faceDivz = sdiag(1/V)*D3*sdiag(S) return self._faceDivz + + @property + def cellGrad(self): + """The cell centered Gradient, takes you to cell faces.""" + raise NotImplementedError('Cell Grad is not yet implemented.') + + @property def nodalGrad(self): """Construct gradient operator (nodes to edges).""" diff --git a/SimPEG/Tests/TestUtils.py b/SimPEG/Tests/TestUtils.py index d84408d4..a980546e 100644 --- a/SimPEG/Tests/TestUtils.py +++ b/SimPEG/Tests/TestUtils.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt from numpy.linalg import norm from SimPEG.Utils import mkvc, sdiag from SimPEG import Utils -from SimPEG.Mesh import TensorMesh, LogicallyOrthogonalMesh +from SimPEG.Mesh import TensorMesh, LogicallyOrthogonalMesh, CylMesh import numpy as np import scipy.sparse as sp import unittest @@ -88,10 +88,7 @@ class OrderTest(unittest.TestCase): """ if 'TensorMesh' in self._meshType: if 'uniform' in self._meshType: - h1 = np.ones(nc)/nc - h2 = np.ones(nc)/nc - h3 = np.ones(nc)/nc - h = [h1, h2, h3] + h = [nc, nc, nc] elif 'random' in self._meshType: h1 = np.random.rand(nc) h2 = np.random.rand(nc) @@ -104,6 +101,20 @@ class OrderTest(unittest.TestCase): max_h = max([np.max(hi) for hi in self.M.h]) return max_h + elif 'CylMesh' in self._meshType: + if 'uniform' in self._meshType: + h = [nc, nc, nc] + else: + raise Exception('Unexpected meshType') + + if self.meshDimension == 2: + self.M = CylMesh([h[0], 1, h[2]]) + max_h = max([np.max(hi) for hi in [self.M.hx, self.M.hz]]) + elif self.meshDimension == 3: + self.M = CylMesh(h) + max_h = max([np.max(hi) for hi in self.M.h]) + return max_h + elif 'LOM' in self._meshType: if 'uniform' in self._meshType: kwrd = 'rect' diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index b6006bff..584bcc2d 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -1,6 +1,7 @@ import unittest import sys from SimPEG import * +from TestUtils import OrderTest class TestCyl2DMesh(unittest.TestCase): @@ -75,8 +76,60 @@ class TestCyl2DMesh(unittest.TestCase): vol = np.r_[2*a,a] self.assertTrue(np.linalg.norm((vol-self.mesh.vol)) == 0) - def test_faceDiv(self): - print self.mesh.faceDiv + def test_gridSizes(self): + self.assertTrue(self.mesh.gridCC.shape == (self.mesh.nC, 3)) + # self.assertTrue(self.mesh.gridN.shape == (self.mesh.nN, 3)) + + self.assertTrue(self.mesh.gridFx.shape == (self.mesh.nFx, 3)) + # self.assertTrue(self.mesh.gridFy.shape == (self.mesh.nFy, 3)) + self.assertTrue(self.mesh.gridFz.shape == (self.mesh.nFz, 3)) + + # self.assertTrue(self.mesh.gridEx.shape == (self.mesh.nEx, 3)) + self.assertTrue(self.mesh.gridEy.shape == (self.mesh.nEy, 3)) + # self.assertTrue(self.mesh.gridEz.shape == (self.mesh.nEz, 3)) + + +MESHTYPES = ['uniformCylMesh'] +MESHDIMENSION = 2 +call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) +call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) +cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] +cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)] +cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy))) +cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey))) +cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz))) +cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez))) + + +class TestFaceDiv(OrderTest): + name = "FaceDiv" + meshTypes = MESHTYPES + meshDimension = MESHDIMENSION + + def getError(self): + + funX = lambda x, y, z: np.cos(2*np.pi*y) + funY = lambda x, y, z: 1+x*0 + funZ = lambda x, y, z: np.cos(2*np.pi*x) + + solX = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*z) + solY = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x) + solZ = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*y) + + Ec = cartE3(self.M, funX, funY, funZ) + print Ec.shape, self.M.nE + print self.M + E = self.M.projectEdgeVector(Ec) + + Fc = cartF3(self.M, solX, solY, solZ) + curlE_anal = self.M.projectFaceVector(Fc) + + curlE = self.M.edgeCurl.dot(E) + err = np.linalg.norm((curlE - curlE_anal), np.inf) + return err + + # def test_order(self): + # self.orderTest() class TestCyl3DMesh(unittest.TestCase):