CylMesh interpolation tests.

SimPEG/Mesh/TensorMesh.py

@@ -180,17 +180,13 @@ class BaseTensorMesh(BaseRectangularMesh):
         :rtype numpy.ndarray
         :return inside, numpy array of booleans
         """
+        pts = Utils.asArray_N_x_Dim(pts, self.dim)
 
         tensors = self.getTensor(locType)
-        if type(pts) == list:
-            pts = np.array(pts)
-        assert type(pts) == np.ndarray, "must be a numpy array"
-        if self.dim > 1:
-            assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)"
-        elif len(pts.shape) == 1:
-            pts = pts[:,np.newaxis]
-        else:
-            assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)"
+
+        if locType == 'N' and self._meshType == 'CYL':
+            #NOTE: for a CYL mesh we add a node to check if we are inside in the radial direction!
+            tensors[0] = np.r_[0.,tensors[0]]
 
         inside = np.ones(pts.shape[0],dtype=bool)
         for i, tensor in enumerate(tensors):
@@ -217,19 +213,7 @@ class BaseTensorMesh(BaseRectangularMesh):
             'CC'    -> scalar field defined on cell centers
         """
 
-        if type(loc) == list:
-            loc = np.array(loc)
-        assert type(loc) == np.ndarray, "must be a numpy array"
-
-        if loc.size == self.dim:
-            loc = np.atleast_2d(loc)
-
-        if self.dim > 1:
-            assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim) not (%d,%d)" % loc.shape
-        elif len(loc.shape) == 1:
-            loc = loc[:,np.newaxis]
-        else:
-            assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)"
+        loc = Utils.asArray_N_x_Dim(loc, self.dim)
 
         if zerosOutside is False:
             assert np.all(self.isInside(loc)), "Points outside of mesh"
@@ -342,7 +326,7 @@ class BaseTensorMesh(BaseRectangularMesh):
             :return: M, the inner product matrix (nF, nF)
         """
         if materialProperty is None:
-            return 0.0
+            return None
 
         if Utils.isScalar(materialProperty):
             Av = getattr(self, 'ave'+AvType+'2CC')

SimPEG/Tests/test_interpolation.py

@@ -7,14 +7,16 @@ import unittest
 
 
 MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh']
-TOLERANCES = [0.9, 0.5]
+TOLERANCES = [0.9, 0.5, 0.5]
 call1 = lambda fun, xyz: fun(xyz)
-call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
+call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, -1])
 call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
 cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)]
 cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)]
 cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy)))
+cartF2Cyl = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFz, fx, fy)))
 cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey)))
+cartE2Cyl = lambda M, ex, ey: cart_row2(M.gridEy, ex, ey)
 cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz)))
 cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez)))
 
@@ -133,6 +135,74 @@ class TestInterpolation2d(OrderTest):
 
 
 
+class TestInterpolation2dCyl(OrderTest):
+    name = "Interpolation 2D"
+    LOCS = np.c_[np.random.rand(4)*0.6+0.2, np.zeros(4), np.random.rand(4)*0.6+0.2]
+    meshTypes = ['uniformCylMesh'] # MESHTYPES +
+    tolerance = 0.6
+    meshDimension = 2
+    meshSizes = [16, 32, 64]
+
+    def getError(self):
+        funX = lambda x, y: np.cos(2*np.pi*y)
+        funY = lambda x, y: np.cos(2*np.pi*x)
+
+        if 'x' in self.type:
+            anal = call2(funX, self.LOCS)
+        elif 'y' in self.type:
+            anal = call2(funY, self.LOCS)
+        elif 'z' in self.type:
+            anal = call2(funY, self.LOCS)
+        else:
+            anal = call2(funX, self.LOCS)
+
+        if 'Fx' == self.type:
+            Fc = cartF2Cyl(self.M, funX, funY)
+            Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]]
+            grid = self.M.projectFaceVector(Fc)
+        elif 'Fz' == self.type:
+            Fc = cartF2Cyl(self.M, funX, funY)
+            Fc = np.c_[Fc[:,0],np.zeros(self.M.nF),Fc[:,1]]
+
+            grid = self.M.projectFaceVector(Fc)
+        elif 'E' in self.type:
+            Ec = cartE2Cyl(self.M, funX, funY)
+            grid = Ec[:,1]
+        elif 'CC' == self.type:
+            grid = call2(funX, self.M.gridCC)
+        elif 'N' == self.type:
+            grid = call2(funX, self.M.gridN)
+
+        comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
+
+        err = np.linalg.norm((comp - anal), np.inf)
+        return err
+
+    def test_orderCC(self):
+        self.type = 'CC'
+        self.name = 'Interpolation 2D CYLMESH: CC'
+        self.orderTest()
+
+    def test_orderN(self):
+        self.type = 'N'
+        self.name = 'Interpolation 2D CYLMESH: N'
+        self.orderTest()
+
+    def test_orderFx(self):
+        self.type = 'Fx'
+        self.name = 'Interpolation 2D CYLMESH: Fx'
+        self.orderTest()
+
+    def test_orderFz(self):
+        self.type = 'Fz'
+        self.name = 'Interpolation 2D CYLMESH: Fz'
+        self.orderTest()
+
+    def test_orderEy(self):
+        self.type = 'Ey'
+        self.name = 'Interpolation 2D CYLMESH: Ey'
+        self.orderTest()
+
 class TestInterpolation3D(OrderTest):
     name = "Interpolation"
     LOCS = np.random.rand(50,3)*0.6+0.2

SimPEG/Tests/test_utils.py

@@ -164,12 +164,34 @@ class TestSequenceFunctions(unittest.TestCase):
             Z = B2*A - sp.identity(M.nC*3)
             self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
 
-    def test_isFloat(self):
+    def test_isScalar(self):
         self.assertTrue(isScalar(1.))
         self.assertTrue(isScalar(1))
         self.assertTrue(isScalar(long(1)))
         self.assertTrue(isScalar(np.r_[1.]))
         self.assertTrue(isScalar(np.r_[1]))
 
+    def test_asArray_N_x_Dim(self):
+
+        true = np.array([[1,2,3]])
+
+        listArray = asArray_N_x_Dim([1,2,3],3)
+        self.assertTrue(np.all(true == listArray))
+        self.assertTrue(true.shape == listArray.shape)
+
+        listArray = asArray_N_x_Dim(np.r_[1,2,3],3)
+        self.assertTrue(np.all(true == listArray))
+        self.assertTrue(true.shape == listArray.shape)
+
+        listArray = asArray_N_x_Dim(np.array([[1,2,3.]]),3)
+        self.assertTrue(np.all(true == listArray))
+        self.assertTrue(true.shape == listArray.shape)
+
+        true = np.array([[1,2],[4,5]])
+
+        listArray = asArray_N_x_Dim([[1,2],[4,5]],2)
+        self.assertTrue(np.all(true == listArray))
+        self.assertTrue(true.shape == listArray.shape)
+
 if __name__ == '__main__':
     unittest.main()

SimPEG/Utils/__init__.py

@@ -144,6 +144,21 @@ def dependentProperty(name, value, children, doc):
         setattr(self, name, val)
     return property(fget=fget, fset=fset, doc=doc)
 
+
+def asArray_N_x_Dim(pts, dim):
+        if type(pts) == list:
+            pts = np.array(pts)
+        assert type(pts) == np.ndarray, "pts must be a numpy array"
+
+        if dim > 1:
+            pts = np.atleast_2d(pts)
+        elif len(pts.shape) == 1:
+            pts = pts[:,np.newaxis]
+
+        assert pts.shape[1] == dim, "pts must be a column vector of shape (nPts, %d) not (%d, %d)" % ((dim,)+pts.shape)
+
+        return pts
+
 def requires(var):
     """
         Use this to wrap a funciton::