Merge pull request #77 from simpeg/gridHfunction

GridH function

Author: lheagy

.bumpversion.cfg

@@ -1,4 +1,4 @@
 [bumpversion]
-current_version = 0.1.11
+current_version = 0.1.12
 files = setup.py discretize/__init__.py docs/conf.py
 

discretize/TensorMesh.py

@@ -159,6 +159,24 @@ def gridN(self):
         """Nodal grid."""
         return self._getTensorGrid('N')
 
+    @property
+    def h_gridded(self):
+        """
+        Returns an (nC, dim) numpy array with the widths of all cells in order
+        """
+
+        if self.dim == 1:
+            return np.reshape(self.h, (self.nC, 1))
+        elif self.dim == 2:
+            hx = np.kron(np.ones(self.nCy), self.h[0])
+            hy = np.kron(self.h[1], np.ones(self.nCx))
+            return np.c_[hx, hy]
+        elif self.dim == 3:
+            hx = np.kron(np.ones(self.nCy*self.nCz), self.h[0])
+            hy = np.kron(np.ones(self.nCz), np.kron(self.h[1], np.ones(self.nCx)))
+            hz = np.kron(self.h[2], np.ones(self.nCx*self.nCy))
+            return np.c_[hx, hy, hz]
+
     @property
     def gridFx(self):
         """Face staggered grid in the x direction."""

discretize/TreeMesh.py

@@ -170,7 +170,7 @@ def __dirty__(self, val):
             '_faceDiv', '_edgeCurl', '_nodalGrad',
             '_aveFx2CC', '_aveFy2CC', '_aveFz2CC', '_aveF2CC', '_aveF2CCV',
             '_aveEx2CC', '_aveEy2CC', '_aveEz2CC', '_aveE2CC', '_aveE2CCV',
-            '_aveN2CC',
+            '_aveN2CC', '_h_gridded'
         ]
         for p in deleteThese:
             if hasattr(self, p): delattr(self, p)
@@ -715,6 +715,12 @@ def balance(self, recursive=True, cells=None, verbose=False, _inRecursion=False)
 
     @property
     def gridCC(self):
+        """
+        Returns an M by N numpy array with the center locations of all cells
+        in order. M is the number of cells and N=1,2,3 is the dimension of the
+        mesh.
+        """
+
         if getattr(self, '_gridCC', None) is None:
             self._gridCC = np.zeros((len(self._cells),self.dim))
             for ii, ind in enumerate(self._sortedCells):
@@ -724,10 +730,28 @@ def gridCC(self):
 
     @property
     def gridN(self):
+        """
+        Returns an M by N numpy array with the widths of all cells in order.
+        M is the number of nodes and N=1,2,3 is the dimension of the mesh.
+        """
+
         self.number()
         R = self._deflationMatrix('N', withHanging=False)
         return R.T * self._gridN + np.repeat([self.x0],self.nN,axis=0)
 
+    @property
+    def h_gridded(self):
+        """
+        Returns an (nC, dim) numpy array with the widths of all cells in order
+        """
+
+        if getattr(self, '_h_gridded', None) is None:
+            self._h_gridded = np.zeros((len(self._cells), self.dim))
+            for ii, ind in enumerate(self._sortedCells):
+                p = self._asPointer(ind)
+                self._h_gridded[ii, :] = self._cellH(p)
+        return self._h_gridded
+
     @property
     def gridFx(self):
         self.number()

discretize/__init__.py

@@ -18,7 +18,7 @@
         """
     )
 
-__version__   = '0.1.11'
+__version__   = '0.1.12'
 __author__    = 'SimPEG Team'
 __license__   = 'MIT'
 __copyright__ = '2013 - 2017, SimPEG Developers, http://simpeg.xyz'

docs/conf.py

@@ -60,9 +60,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = '0.1.11'
+version = '0.1.12'
 # The full version, including alpha/beta/rc tags.
-release = '0.1.11'
+release = '0.1.12'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

setup.py

@@ -56,7 +56,7 @@ def configuration(parent_package='', top_path=None):
 
 setup(
     name="discretize",
-    version="0.1.11",
+    version="0.1.12",
     install_requires=[
         'numpy>=1.7',
         'scipy>=0.13',

tests/base/test_tensor.py

@@ -16,6 +16,19 @@ def setUp(self):
         self.mesh2 = discretize.TensorMesh([a, b], [3, 5])
         self.mesh3 = discretize.TensorMesh([a, b, c])
 
+    def test_gridded_2D(self):
+        H = self.mesh2.h_gridded
+        test_hx = np.all(H[:, 0] == np.r_[1., 1., 1., 1., 1., 1.])
+        test_hy = np.all(H[:, 1] == np.r_[1., 1., 1., 2., 2., 2.])
+        self.assertTrue(test_hx and test_hy)
+
+    def test_gridded_3D(self):
+        H = self.mesh3.h_gridded
+        test_hx = np.all(H[:, 0] == np.r_[1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
+        test_hy = np.all(H[:, 1] == np.r_[1., 1., 1., 2., 2., 2., 1., 1., 1., 2., 2., 2.])
+        test_hz = np.all(H[:, 2] == np.r_[1., 1., 1., 1., 1., 1., 4., 4., 4., 4., 4., 4.])
+        self.assertTrue(test_hx and test_hy and test_hz)
+
     def test_vectorN_2D(self):
         testNx = np.array([3, 4, 5, 6])
         testNy = np.array([5, 6, 8])

tests/tree/test_tree.py

@@ -107,6 +107,40 @@ def test_corsen(self):
         assert M._index([0, 2, 2]) not in M
         assert M._index([2, 2, 2]) not in M
 
+    def test_h_gridded_2D(self):
+        hx, hy = np.ones(4), np.r_[1., 2., 3., 4.]
+
+        M = discretize.TreeMesh([hx, hy])
+
+        def refinefcn(cell):
+            xyz = cell.center
+            d = (xyz**2).sum()**0.5
+            if d < 3:
+                return 2
+            return 1
+
+        M.refine(refinefcn)
+        H = M.h_gridded
+
+        test_hx = np.all(H[:, 0] == np.r_[1., 1., 1., 1., 2., 2., 2.])
+        test_hy = np.all(H[:, 1] == np.r_[1., 1., 2., 2., 3., 7., 7.])
+
+        self.assertTrue(test_hx and test_hy)
+
+    def test_h_gridded_updates(self):
+        mesh = discretize.TreeMesh([8, 8])
+        mesh.refine(1)
+
+        H = mesh.h_gridded
+        self.assertTrue(np.all(H[:, 0] == 0.5*np.ones(4)))
+        self.assertTrue(np.all(H[:, 1] == 0.5*np.ones(4)))
+
+        # refine the mesh and make sure h_gridded is updated
+        mesh.refine(2)
+        H = mesh.h_gridded
+        self.assertTrue(np.all(H[:, 0] == 0.25*np.ones(16)))
+        self.assertTrue(np.all(H[:, 1] == 0.25*np.ones(16)))
+
     def test_faceDiv(self):
 
         hx, hy = np.r_[1., 2, 3, 4], np.r_[5., 6, 7, 8]
@@ -233,6 +267,26 @@ def test_VectorIdenties(self):
         assert np.max(np.abs((M.faceDiv * M.edgeCurl).todense().flatten())) < TOL
         assert np.max(np.abs((Mr.faceDiv * Mr.edgeCurl).todense().flatten())) < TOL
 
+    def test_h_gridded_3D(self):
+        hx, hy, hz = np.ones(4), np.r_[1., 2., 3., 4.], 2*np.ones(4)
+
+        M = discretize.TreeMesh([hx, hy, hz])
+
+        def refinefcn(cell):
+            xyz = cell.center
+            d = (xyz**2).sum()**0.5
+            if d < 3:
+                return 2
+            return 1
+
+        M.refine(refinefcn)
+        H = M.h_gridded
+
+        test_hx = np.all(H[:, 0] == np.r_[1., 1., 1., 1., 1., 1., 1., 1., 2., 2., 2., 2., 2., 2., 2.])
+        test_hy = np.all(H[:, 1] == np.r_[1., 1., 2., 2., 1., 1., 2., 2., 3., 7., 7., 3., 3., 7., 7.])
+        test_hz = np.all(H[:, 2] == np.r_[2., 2., 2., 2., 2., 2., 2., 2., 4., 4., 4., 4., 4., 4., 4.])
+
+        self.assertTrue(test_hx and test_hy and test_hz)
 
 class Test2DInterpolation(unittest.TestCase):