diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 376cb789..226b8720 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -115,9 +115,8 @@ class CylMesh(BaseTensorMesh): @property def vectorCCx(self): """Cell-centered grid vector (1D) in the x direction.""" - if self.nCy == 1: - return np.r_[-self.hx[0]*0.5, self.hx[:-1].cumsum()] + self.hx*0.5 # - self.hx[0]/2 - return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + firstEl = -self.hx[0]*0.5 if self.nCy == 1 else 0 + return np.r_[firstEl, self.hx[:-1].cumsum()] + self.hx*0.5 @property def vectorCCy(self): @@ -128,55 +127,18 @@ class CylMesh(BaseTensorMesh): def vectorNx(self): """Nodal grid vector (1D) in the x direction.""" if self.nCy == 1: - return self.hx.cumsum()# - self.hx[0]/2 + return self.hx.cumsum() return np.r_[0, self.hx].cumsum() @property def vectorNy(self): """Nodal grid vector (1D) in the y direction.""" if self.nCy == 1: + # There aren't really any nodes, but all the grids need + # somewhere to live, why not zero?! return np.r_[0] 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""" @@ -270,6 +232,8 @@ class CylMesh(BaseTensorMesh): @property def nodalGrad(self): """Construct gradient operator (nodes to edges).""" + if self.nCy == 1: + raise Exception('Nodal grad does not make sense for cylindrically symmetric mesh.') raise NotImplementedError('nodalGrad not yet implemented') @property diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 9a68d12a..44623883 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -95,6 +95,7 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridFx(self): + if self.nFx == 0: return """Face staggered grid in the x direction.""" if getattr(self, '_gridFx', None) is None: self._gridFx = Utils.ndgrid(self.getTensor('Fx')) @@ -102,6 +103,7 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridFy(self): + if self.nFy == 0: return """Face staggered grid in the y direction.""" if getattr(self, '_gridFy', None) is None and self.dim > 1: self._gridFy = Utils.ndgrid(self.getTensor('Fy')) @@ -109,6 +111,7 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridFz(self): + if self.nFz == 0: return """Face staggered grid in the z direction.""" if getattr(self, '_gridFz', None) is None and self.dim > 2: self._gridFz = Utils.ndgrid(self.getTensor('Fz')) @@ -116,6 +119,7 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridEx(self): + if self.nEx == 0: return """Edge staggered grid in the x direction.""" if getattr(self, '_gridEx', None) is None: self._gridEx = Utils.ndgrid(self.getTensor('Ex')) @@ -123,6 +127,7 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridEy(self): + if self.nEy == 0: return """Edge staggered grid in the y direction.""" if getattr(self, '_gridEy', None) is None and self.dim > 1: self._gridEy = Utils.ndgrid(self.getTensor('Ey')) @@ -130,6 +135,7 @@ class BaseTensorMesh(BaseRectangularMesh): @property def gridEz(self): + if self.nEz == 0: return """Edge staggered grid in the z direction.""" if getattr(self, '_gridEz', None) is None and self.dim > 2: self._gridEz = Utils.ndgrid(self.getTensor('Ez')) diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index ebde658c..9eaf0998 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -78,15 +78,44 @@ class TestCyl2DMesh(unittest.TestCase): 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.gridN.shape == (9, 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.gridFy is None) self.assertTrue(self.mesh.gridFz.shape == (self.mesh.nFz, 3)) - # self.assertTrue(self.mesh.gridEx.shape == (self.mesh.nEx, 3)) + self.assertTrue(self.mesh.gridEx is None) self.assertTrue(self.mesh.gridEy.shape == (self.mesh.nEy, 3)) - # self.assertTrue(self.mesh.gridEz.shape == (self.mesh.nEz, 3)) + self.assertTrue(self.mesh.gridEz is None) + + def test_gridCC(self): + x = np.r_[0,1.5,2.25,0,1.5,2.25] + y = np.zeros(6) + z = np.r_[1,1,1,2.5,2.5,2.5] + G = np.c_[x,y,z] + self.assertTrue(np.linalg.norm((G-self.mesh.gridCC).ravel()) == 0) + + def test_gridFx(self): + x = np.r_[1,2,2.5,1,2,2.5] + y = np.zeros(6) + z = np.r_[1,1,1,2.5,2.5,2.5] + G = np.c_[x,y,z] + self.assertTrue(np.linalg.norm((G-self.mesh.gridFx).ravel()) == 0) + + def test_gridFz(self): + x = np.r_[0,1.5,2.25,0,1.5,2.25,0,1.5,2.25] + y = np.zeros(9) + z = np.r_[0,0,0,2,2,2,3,3,3.] + G = np.c_[x,y,z] + self.assertTrue(np.linalg.norm((G-self.mesh.gridFz).ravel()) == 0) + + def test_gridEy(self): + x = np.r_[1,2,2.5,1,2,2.5,1,2,2.5] + y = np.zeros(9) + z = np.r_[0,0,0,2,2,2,3,3,3.] + G = np.c_[x,y,z] + self.assertTrue(np.linalg.norm((G-self.mesh.gridEy).ravel()) == 0) + MESHTYPES = ['uniformCylMesh']