vectors in CylMesh

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

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__':