diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 3f01be59..935c8ce0 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -6,7 +6,7 @@ from TensorMesh import TensorMesh class CylMesh(TensorMesh): """ - CylMesh is a mesh class for cylindrically problems + CylMesh is a mesh class for cylindrical problems """ _meshType = 'CYL' @@ -15,7 +15,8 @@ class CylMesh(TensorMesh): assert len(h) == 3, "len(h) must equal 3, for a cylindrically symmetric mesh use [hx, 1, hz]" if x0 is not None: - assert x0.size == 3, "x0.size must equal 1" + assert type(x0) == np.ndarray, "x0 must be an ndarray" + assert x0.size == 3, "x0 must have 3 elements" else: x0 = np.r_[0, 0, 0] @@ -39,7 +40,9 @@ class CylMesh(TensorMesh): :rtype: int :return: nNx """ - return self.nCx + if self.nCy == 1: + return self.nCx + return self.nCx + 1 @property def nNy(self): @@ -49,7 +52,9 @@ class CylMesh(TensorMesh): :rtype: int :return: nNy """ - return self.nCy - 1 + if self.nCy == 1: + return self.nCy - 1 + return self.nCy @property def nN(self): @@ -121,16 +126,38 @@ class CylMesh(TensorMesh): """ return (self._n + np.r_[0,-1,0]).prod() + @property + def vectorCCx(self): + """Cell-centered grid vector (1D) in the x direction.""" + if self.nCy == 1: + return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 - self.hx[0]/2 + return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + + @property + def vectorCCy(self): + """Cell-centered grid vector (1D) in the y direction.""" + return np.r_[0, self.hy[:-1]] + @property def vectorNx(self): - """Nodal grid vector (1D) in the r direction""" - return self.hr.cumsum() + """Nodal grid vector (1D) in the x direction.""" + if self.nCy == 1: + return self.hx.cumsum() - self.hx[0]/2 + return np.r_[0, self.hx].cumsum() + + @property + def vectorNy(self): + """Nodal grid vector (1D) in the y direction.""" + return np.r_[0, self.hy[:-1].cumsum()] + self.hy[0]*0.5 @property def edge(self): """Edge lengths""" if getattr(self, '_edge', None) is None: - self._edge = 2*pi*self.gridN[:,0] + if self.nCy == 1: + self._edge = 2*pi*self.gridN[:,0] + else: + raise NotImplementedError('edges not implemented for 3D cyl mesh') return self._edge @property diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index aa1af471..c0db7982 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -3,20 +3,18 @@ import sys from SimPEG import * -class TestCyl1DMesh(unittest.TestCase): +class TestCyl2DMesh(unittest.TestCase): def setUp(self): - hx = np.ones(3) - hz = np.ones(2) + hx = np.r_[1,1,0.5] + hz = np.r_[2,1] self.mesh = Mesh.CylMesh([hx, 1,hz]) def test_cylMeshInheritance(self): self.assertTrue(isinstance(self.mesh, Mesh.BaseMesh)) - def test_cylMeshDimensions(self): - self.assertTrue(self.mesh.dim == 3) - def test_cylMesh_numbers(self): + self.assertTrue(self.mesh.dim == 3) self.assertTrue(self.mesh.nCx == 3) self.assertTrue(self.mesh.nCy == 1) self.assertTrue(self.mesh.nCz == 2) @@ -45,7 +43,62 @@ class TestCyl1DMesh(unittest.TestCase): self.assertTrue(np.all(self.mesh.vnEz == [3, 0, 2])) self.assertTrue(self.mesh.nE == 9) + def test_vectorsCC(self): + v = np.r_[0, 1, 1.75] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0) + v = np.r_[0] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCy)) == 0) + v = np.r_[1, 2.5] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCz)) == 0) + def test_vectorsN(self): + v = np.r_[0.5, 1.5, 2] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNx)) == 0) + v = np.r_[np.pi] #This is kinda a fake. But it is where it would be if there was a radial connection + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNy)) == 0) + v = np.r_[0, 2, 3.] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNz)) == 0) + + def test_dimensions(self): + v = np.r_[0.5, 1.5, 2, 0.5, 1.5, 2, 0.5, 1.5, 2] * 2 * np.pi + self.assertTrue(np.linalg.norm((v-self.mesh.edge)) == 0) + + +class TestCyl3DMesh(unittest.TestCase): + + def setUp(self): + hx = np.r_[1,1,0.5] + hy = np.r_[np.pi, np.pi] + hz = np.r_[2,1] + self.mesh = Mesh.CylMesh([hx, hy,hz]) + + def test_cylMesh_numbers(self): + self.assertTrue(self.mesh.nCx == 3) + self.assertTrue(self.mesh.nCy == 2) + self.assertTrue(self.mesh.nCz == 2) + self.assertTrue(np.all(self.mesh.vnC == [3, 2, 2])) + + self.assertTrue(self.mesh.nN == 24) + self.assertTrue(self.mesh.nNx == 4) + self.assertTrue(self.mesh.nNy == 2) + self.assertTrue(self.mesh.nNz == 3) + self.assertTrue(np.all(self.mesh.vnN == [4, 2, 3])) + + def test_vectorsCC(self): + v = np.r_[0.5, 1.5, 2.25] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCx)) == 0) + v = np.r_[0, np.pi] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCy)) == 0) + v = np.r_[1, 2.5] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorCCz)) == 0) + + def test_vectorsN(self): + v = np.r_[0, 1, 2, 2.5] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNx)) == 0) + v = np.r_[np.pi/2, 1.5*np.pi] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNy)) == 0) + v = np.r_[0, 2, 3] + self.assertTrue(np.linalg.norm((v-self.mesh.vectorNz)) == 0) if __name__ == '__main__':