diff --git a/SimPEG/Mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py index 14b3aecf..0ca82ce1 100644 --- a/SimPEG/Mesh/BaseMesh.py +++ b/SimPEG/Mesh/BaseMesh.py @@ -1,4 +1,5 @@ import numpy as np +import scipy.sparse as sp from SimPEG import Utils @@ -594,3 +595,54 @@ class BaseRectangularMesh(BaseMesh): return out else: return switchKernal(x) + + def getInterpolationMatMesh2Mesh(self, mesh2, locType='CC'): + """ + Interpolates variables from the current mesh to a new mesh (mesh2) + + :param Mesh mesh2: SimPEG mesh which we interpolate values to + :param string locType: location of variables 'CC', 'E', 'F', 'N' + :rtype: scipy.sparse.csr_matrix + :return P: interpolation matrix + """ + + # Error Checking + if self._meshType == 'Cyl': + assert self.isSymmetric, "Currently, we do not support non-symmetric cyl meshes" + if mesh2._meshType == 'Cyl': + assert self._meshType == 'Cyl', "Interpolation from 3D mesh to Cyl mesh is not supported" + + # if Cyl to cart call + if self._meshType == 'Cyl' and mesh2._meshType != 'Cyl': + return self.getInterpolationMatCartMesh(mesh2, locType) + + # Scalars + if locType in ['CC', 'N', 'Fx', 'Fy', 'Fz', 'Ex', 'Ey', 'Ez']: + grid = getattr(mesh2, 'grid%s'%locType) + return self.getInterpolationMat(grid, locType) + + # Vectors + else: + if self._meshType == 'Cyl': + if locType == 'F': + X = self.getInterpolationMatMesh2Mesh(mesh2, locType='Fx') + Z = self.getInterpolationMatMesh2Mesh(mesh2, locType='Fz') + return sp.block_diag([X, Z]) + elif locType == 'E': + return self.getInterpolationMatMesh2Mesh(mesh2, locType='Ey') + + if self.dim == 1: + return self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType) + elif self.dim == 2: + X = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType) + Y = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sy'%locType) + return sp.block_diag([X, Y]) + elif self.dim == 3: + X = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sx'%locType) + Y = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sy'%locType) + Z = self.getInterpolationMatMesh2Mesh(mesh2, locType='%sz'%locType) + return sp.block_diag([X, Y, Z]) + + + + diff --git a/SimPEG/Tests.py b/SimPEG/Tests.py index 804ca6b8..7c395ef3 100644 --- a/SimPEG/Tests.py +++ b/SimPEG/Tests.py @@ -82,14 +82,14 @@ class OrderTest(unittest.TestCase): _meshType = meshTypes[0] meshDimension = 3 - def setupMesh(self, nc): + def makeMesh(self, nc, meshType=_meshType, meshDimension=meshDimension): """ For a given number of cells nc, generate a TensorMesh with uniform cells with edge length h=1/nc. """ - if 'TensorMesh' in self._meshType: - if 'uniform' in self._meshType: + if 'TensorMesh' in meshType: + if 'uniform' in meshType: h = [nc, nc, nc] - elif 'random' in self._meshType: + elif 'random' in meshType: h1 = np.random.rand(nc)*nc*0.5 + nc*0.5 h2 = np.random.rand(nc)*nc*0.5 + nc*0.5 h3 = np.random.rand(nc)*nc*0.5 + nc*0.5 @@ -97,46 +97,46 @@ class OrderTest(unittest.TestCase): else: raise Exception('Unexpected meshType') - self.M = TensorMesh(h[:self.meshDimension]) - max_h = max([np.max(hi) for hi in self.M.h]) - return max_h + M = TensorMesh(h[:meshDimension]) + max_h = max([np.max(hi) for hi in M.h]) + return M, max_h - elif 'CylMesh' in self._meshType: - if 'uniform' in self._meshType: + elif 'CylMesh' in meshType: + if 'uniform' in meshType: h = [nc, nc, nc] else: raise Exception('Unexpected meshType') - if self.meshDimension == 2: - self.M = CylMesh([h[0], 1, h[2]]) - max_h = max([np.max(hi) for hi in [self.M.hx, self.M.hz]]) - elif self.meshDimension == 3: - self.M = CylMesh(h) - max_h = max([np.max(hi) for hi in self.M.h]) - return max_h + if meshDimension == 2: + M = CylMesh([h[0], 1, h[2]]) + max_h = max([np.max(hi) for hi in [M.hx, M.hz]]) + elif meshDimension == 3: + M = CylMesh(h) + max_h = max([np.max(hi) for hi in M.h]) + return M, max_h - elif 'Curv' in self._meshType: - if 'uniform' in self._meshType: + elif 'Curv' in meshType: + if 'uniform' in meshType: kwrd = 'rect' - elif 'rotate' in self._meshType: + elif 'rotate' in meshType: kwrd = 'rotate' else: raise Exception('Unexpected meshType') - if self.meshDimension == 1: + if meshDimension == 1: raise Exception('Lom not supported for 1D') - elif self.meshDimension == 2: + elif meshDimension == 2: X, Y = Utils.exampleLrmGrid([nc, nc], kwrd) - self.M = CurvilinearMesh([X, Y]) - elif self.meshDimension == 3: + M = CurvilinearMesh([X, Y]) + elif meshDimension == 3: X, Y, Z = Utils.exampleLrmGrid([nc, nc, nc], kwrd) - self.M = CurvilinearMesh([X, Y, Z]) - return 1./nc + M = CurvilinearMesh([X, Y, Z]) + return M, 1./nc - elif 'Tree' in self._meshType: + elif 'Tree' in meshType: nc *= 2 - if 'uniform' in self._meshType or 'notatree' in self._meshType: + if 'uniform' in meshType or 'notatree' in meshType: h = [nc, nc, nc] - elif 'random' in self._meshType: + elif 'random' in meshType: h1 = np.random.rand(nc)*nc*0.5 + nc*0.5 h2 = np.random.rand(nc)*nc*0.5 + nc*0.5 h3 = np.random.rand(nc)*nc*0.5 + nc*0.5 @@ -145,20 +145,29 @@ class OrderTest(unittest.TestCase): raise Exception('Unexpected meshType') levels = int(np.log(nc)/np.log(2)) - self.M = Tree(h[:self.meshDimension], levels=levels) + M = Tree(h[:meshDimension], levels=levels) def function(cell): - if 'notatree' in self._meshType: + if 'notatree' in meshType: return levels - 1 r = cell.center - np.array([0.5]*len(cell.center)) dist = np.sqrt(r.dot(r)) if dist < 0.2: return levels return levels - 1 - self.M.refine(function,balance=False) - self.M.number(balance=False) - # self.M.plotGrid(showIt=True) - max_h = max([np.max(hi) for hi in self.M.h]) - return max_h + M.refine(function,balance=False) + M.number(balance=False) + # M.plotGrid(showIt=True) + max_h = max([np.max(hi) for hi in M.h]) + return M, max_h + + + def setupMesh(self, nc): + """ + For a given number of cells nc, generate a TensorMesh with uniform cells with edge length h=1/nc. + """ + M, h = self.makeMesh(nc, meshType=self._meshType, meshDimension=self.meshDimension) + self.M = M + return h def getError(self): """For given h, generate A[h], f and A(f) and return norm of error.""" diff --git a/tests/mesh/test_InterpolationMesh2Mesh.py b/tests/mesh/test_InterpolationMesh2Mesh.py new file mode 100644 index 00000000..e0b6a817 --- /dev/null +++ b/tests/mesh/test_InterpolationMesh2Mesh.py @@ -0,0 +1,255 @@ +import numpy as np +import unittest +from SimPEG.Utils import mkvc +from SimPEG import Mesh, Tests +import unittest + +call1 = lambda fun, xyz: fun(xyz) +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))) + +TOL = 1e-7 + +class TestInterpolationMesh2Mesh_Tensor1D(Tests.OrderTest): + + name = 'Mesh2Mesh Tensor1D' + meshSizes = [8, 16, 32] + meshTypes = ['uniformTensorMesh'] + meshDimension = 1 + + def getError(self): + funX = lambda x: np.cos(2*np.pi*x) + + mesh2, _ = self.makeMesh(self.M.nC-1, meshType=self._meshType, meshDimension=self.meshDimension ) + ana = call1(funX, getattr(mesh2, 'grid%s'%self.type)) + + v = call1(funX, getattr(self.M, 'grid%s'%self.type)) + P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type) + num = P*v + + return np.linalg.norm((num - ana), np.inf) + + def test_orderCC_1D(self): + self.type = 'CC' + self.name = 'Mesh2Mesh Tensor1D: CC' + self.orderTest() + + def test_orderN_1D(self): + self.type = 'N' + self.name = 'Mesh2Mesh Tensor1D: N' + self.orderTest() + + def test_orderEx_1D(self): + self.type = 'Ex' + self.name = 'Mesh2Mesh Tensor1D: Ex' + self.orderTest() + + def test_orderFx_1D(self): + self.type = 'Fx' + self.name = 'Mesh2Mesh Tensor1D: Fx' + self.orderTest() + +class TestInterpolationMesh2Mesh_Tensor2D(Tests.OrderTest): + + name = 'Mesh2Mesh Tensor2D' + meshSizes = [4, 8, 16] + meshTypes = ['uniformTensorMesh'] + meshDimension = 2 + + def getError(self): + funX = lambda x, y: np.cos(2*np.pi*y) + funY = lambda x, y: np.cos(2*np.pi*x) + + mesh2, _ = self.makeMesh(self.M.nC-1, meshType=self._meshType, meshDimension=self.meshDimension ) + + if 'x' in self.type: + ana = call2(funX, getattr(mesh2, 'grid%s'%self.type)) + elif 'y' in self.type: + ana = call2(funY, getattr(mesh2, 'grid%s'%self.type)) + elif 'F' in self.type: + ana = cartF2(mesh2, funX, funY) + ana = mesh2.projectFaceVector(ana) + elif 'E' in self.type: + ana = cartE2(mesh2, funX, funY) + ana = mesh2.projectFaceVector(ana) + else: + ana = call2(funX, getattr(mesh2, 'grid%s'%self.type)) + + + if 'F' in self.type: + v = cartF2(self.M, funX, funY) + if 'x' in self.type or 'y' in self.type: + v = self.M.projectFaceVector(v) + else: + v = mkvc(v) + elif 'E' in self.type: + v = cartE2(self.M, funX, funY) + if 'x' in self.type or 'y' in self.type: + v = self.M.projectFaceVector(v) + else: + v = mkvc(v) + elif 'CC' == self.type: + v = call2(funX, self.M.gridCC) + elif 'N' == self.type: + v = call2(funX, self.M.gridN) + + P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type) + # print P.shape, v.shape + num = P*v + + return np.linalg.norm((num - ana), np.inf) + + def test_orderCC_2D(self): + self.type = 'CC' + self.name = 'Mesh2Mesh Tensor2D: CC' + self.orderTest() + + def test_orderN_2D(self): + self.type = 'N' + self.name = 'Mesh2Mesh Tensor2D: N' + self.orderTest() + + def test_orderE_2D(self): + self.type = 'E' + self.name = 'Mesh2Mesh Tensor2D: E' + self.orderTest() + + def test_orderEx_2D(self): + self.type = 'Ex' + self.name = 'Mesh2Mesh Tensor2D: Ex' + self.orderTest() + + def test_orderEy_2D(self): + self.type = 'Ey' + self.name = 'Mesh2Mesh Tensor2D: Ey' + self.orderTest() + + def test_orderF_2D(self): + self.type = 'F' + self.name = 'Mesh2Mesh Tensor2D: F' + self.orderTest() + + def test_orderFx_2D(self): + self.type = 'Fx' + self.name = 'Mesh2Mesh Tensor2D: Fx' + self.orderTest() + + def test_orderFy_2D(self): + self.type = 'Fy' + self.name = 'Mesh2Mesh Tensor2D: Fy' + self.orderTest() + +class TestInterpolationMesh2Mesh_Tensor3D(Tests.OrderTest): + + name = 'Mesh2Mesh Tensor3D' + meshSizes = [4, 8, 16] + meshTypes = ['uniformTensorMesh'] + meshDimension = 3 + + def getError(self): + funX = lambda x, y, z: np.cos(2*np.pi*y) + funY = lambda x, y, z: np.cos(2*np.pi*z) + funZ = lambda x, y, z: np.cos(2*np.pi*x) + + mesh2, _ = self.makeMesh(self.M.nC-1, meshType=self._meshType, meshDimension=self.meshDimension ) + + if 'x' in self.type: + ana = call3(funX, getattr(mesh2, 'grid%s'%self.type)) + elif 'y' in self.type: + ana = call3(funY, getattr(mesh2, 'grid%s'%self.type)) + elif 'z' in self.type: + ana = call3(funZ, getattr(mesh2, 'grid%s'%self.type)) + elif 'F' in self.type: + ana = cartF3(mesh2, funX, funY, funZ) + ana = mesh2.projectFaceVector(ana) + elif 'E' in self.type: + ana = cartE3(mesh2, funX, funY, funZ) + ana = mesh2.projectFaceVector(ana) + else: + ana = call3(funX, getattr(mesh2, 'grid%s'%self.type)) + + + if 'F' in self.type: + v = cartF3(self.M, funX, funY, funZ) + if 'x' in self.type or 'y' in self.type or 'z' in self.type: + v = self.M.projectFaceVector(v) + else: + v = mkvc(v) + elif 'E' in self.type: + v = cartE3(self.M, funX, funY, funZ) + if 'x' in self.type or 'y' in self.type or 'z' in self.type: + v = self.M.projectFaceVector(v) + else: + v = mkvc(v) + elif 'CC' == self.type: + v = call3(funX, self.M.gridCC) + elif 'N' == self.type: + v = call3(funX, self.M.gridN) + + P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type) + # print P.shape, v.shape + num = P*v + + return np.linalg.norm((num - ana), np.inf) + + def test_orderCC_3D(self): + self.type = 'CC' + self.name = 'Mesh2Mesh Tensor3D: CC' + self.orderTest() + + def test_orderN_3D(self): + self.type = 'N' + self.name = 'Mesh2Mesh Tensor3D: N' + self.orderTest() + + def test_orderE_3D(self): + self.type = 'E' + self.name = 'Mesh2Mesh Tensor3D: E' + self.orderTest() + + def test_orderEx_3D(self): + self.type = 'Ex' + self.name = 'Mesh2Mesh Tensor3D: Ex' + self.orderTest() + + def test_orderEy_3D(self): + self.type = 'Ey' + self.name = 'Mesh2Mesh Tensor3D: Ey' + self.orderTest() + + def test_orderEz_3D(self): + self.type = 'Ez' + self.name = 'Mesh2Mesh Tensor3D: Ez' + self.orderTest() + + def test_orderF_3D(self): + self.type = 'F' + self.name = 'Mesh2Mesh Tensor3D: F' + self.orderTest() + + def test_orderFx_3D(self): + self.type = 'Fx' + self.name = 'Mesh2Mesh Tensor3D: Fx' + self.orderTest() + + def test_orderFy_3D(self): + self.type = 'Fy' + self.name = 'Mesh2Mesh Tensor3D: Fy' + self.orderTest() + + def test_orderFz_3D(self): + self.type = 'Fz' + self.name = 'Mesh2Mesh Tensor3D: Fz' + self.orderTest() + + +if __name__ == '__main__': + unittest.main()