Merge pull request #67 from simpeg/ref/extractCoreMesh

Extract Core Mesh

discretize/utils/meshutils.py

@@ -173,17 +173,16 @@ def closestPoints(mesh, pts, gridLoc='CC'):
 
 
 def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
-    """Extracts Core Mesh from Global mesh
+    """
+    Extracts Core Mesh from Global mesh
 
     :param numpy.ndarray xyzlim: 2D array [ndim x 2]
     :param BaseMesh mesh: The mesh
 
     This function ouputs::
 
         - actind: corresponding boolean index from global to core
-        - meshcore: core mesh
-
-    Warning: 1D and 2D has not been tested
+        - meshcore: core sdiscretize mesh
     """
     import discretize
     if mesh.dim == 1:
@@ -196,18 +195,18 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
 
         hx = mesh.hx[xind]
 
-        x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
+        x0 = [xc[0]-hx[0]*0.5]
 
-        meshCore = discretize.TensorMesh([hx, hy], x0=x0)
+        meshCore = discretize.TensorMesh([hx], x0=x0)
 
-        actind = (mesh.gridCC[:, 0] > xmin) & (mesh.gridCC[:, 0] < xmax)
+        actind = (mesh.gridCC > xmin) & (mesh.gridCC < xmax)
 
     elif mesh.dim == 2:
         xmin, xmax = xyzlim[0, 0], xyzlim[0, 1]
         ymin, ymax = xyzlim[1, 0], xyzlim[1, 1]
 
+        xind = np.logical_and(mesh.vectorCCx > xmin, mesh.vectorCCx < xmax)
         yind = np.logical_and(mesh.vectorCCy > ymin, mesh.vectorCCy < ymax)
-        zind = np.logical_and(mesh.vectorCCz > zmin, mesh.vectorCCz < zmax)
 
         xc = mesh.vectorCCx[xind]
         yc = mesh.vectorCCy[yind]
@@ -219,8 +218,10 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
 
         meshCore = discretize.TensorMesh([hx, hy], x0=x0)
 
-        actind = (mesh.gridCC[:, 0] > xmin) & (mesh.gridCC[:, 0] < xmax) \
-               & (mesh.gridCC[:, 1] > ymin) & (mesh.gridCC[:, 1] < ymax) \
+        actind = (
+            (mesh.gridCC[:, 0] > xmin) & (mesh.gridCC[:, 0] < xmax) &
+            (mesh.gridCC[:, 1] > ymin) & (mesh.gridCC[:, 1] < ymax)
+        )
 
     elif mesh.dim == 3:
         xmin, xmax = xyzlim[0, 0], xyzlim[0, 1]

tests/base/test_utils.py

@@ -6,7 +6,8 @@
     sdiag, sub2ind, ndgrid, mkvc, isScalar,
     inv2X2BlockDiagonal, inv3X3BlockDiagonal,
     invPropertyTensor, makePropertyTensor, indexCube,
-    ind2sub, asArray_N_x_Dim, TensorType, Zero, Identity
+    ind2sub, asArray_N_x_Dim, TensorType, Zero, Identity,
+    ExtractCoreMesh
 )
 from discretize.Tests import checkDerivative
 import discretize
@@ -381,5 +382,47 @@ def test_both(self):
         assert o*z + o == 1
         assert o-z == 1
 
+
+class TestMeshUtils(unittest.TestCase):
+
+    def test_ExtractCoreMesh(self):
+
+        # 1D Test on TensorMesh
+        meshtest1d = discretize.TensorMesh([[(50., 10)]])
+        xzlim1d = np.r_[[[0., 250.]]]
+        actind1d, meshCore1d = ExtractCoreMesh(xzlim1d, meshtest1d)
+
+        assert len(actind1d) == meshtest1d.nC
+        assert meshCore1d.nC == np.count_nonzero(actind1d)
+        assert meshCore1d.vectorCCx.min() > xzlim1d[0, :].min()
+        assert meshCore1d.vectorCCx.max() < xzlim1d[0, :].max()
+
+        # 2D Test on TensorMesh
+        meshtest2d = discretize.TensorMesh([[(50., 10)], [(25., 10)]])
+        xzlim2d = xyzlim = np.r_[[[0., 200.], [0., 200.]]]
+        actind2d, meshCore2d = ExtractCoreMesh(xzlim2d, meshtest2d)
+
+        assert len(actind2d) == meshtest2d.nC
+        assert meshCore2d.nC == np.count_nonzero(actind2d)
+        assert meshCore2d.vectorCCx.min() > xzlim2d[0, :].min()
+        assert meshCore2d.vectorCCx.max() < xzlim2d[0, :].max()
+        assert meshCore2d.vectorCCy.min() > xzlim2d[1, :].min()
+        assert meshCore2d.vectorCCy.max() < xzlim2d[1, :].max()
+
+        # 3D Test on TensorMesh
+        meshtest3d = discretize.TensorMesh([[(50., 10)], [(25., 10)], [(5., 40)]])
+        xzlim3d = np.r_[[[0., 250.], [0., 200.], [0., 150]]]
+        actind3d, meshCore3d = ExtractCoreMesh(xzlim3d, meshtest3d)
+
+        assert len(actind3d) == meshtest3d.nC
+        assert meshCore3d.nC == np.count_nonzero(actind3d)
+        assert meshCore3d.vectorCCx.min() > xzlim3d[0, :].min()
+        assert meshCore3d.vectorCCx.max() < xzlim3d[0, :].max()
+        assert meshCore3d.vectorCCy.min() > xzlim3d[1, :].min()
+        assert meshCore3d.vectorCCy.max() < xzlim3d[1, :].max()
+        assert meshCore3d.vectorCCz.min() > xzlim3d[2, :].min()
+        assert meshCore3d.vectorCCz.max() < xzlim3d[2, :].max()
+
+
 if __name__ == '__main__':
     unittest.main()