diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index a8e29768..69c4cd49 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/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') diff --git a/SimPEG/Tests/test_interpolation.py b/SimPEG/Tests/test_interpolation.py index e2634cb7..4ea0ef6d 100644 --- a/SimPEG/Tests/test_interpolation.py +++ b/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 diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py index ad737489..0ac03ef3 100644 --- a/SimPEG/Tests/test_utils.py +++ b/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() diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 05796ef4..0d6a8dda 100644 --- a/SimPEG/Utils/__init__.py +++ b/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::