From bf300061bc712afc51cf94877c1d4d83f0bbc3da Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 31 Mar 2016 11:16:14 -0700 Subject: [PATCH] test cyl2cly for mesh to mesh interpolation --- SimPEG/Mesh/BaseMesh.py | 18 +- tests/mesh/test_InterpolationMesh2Mesh.py | 499 +++++++++++++--------- 2 files changed, 311 insertions(+), 206 deletions(-) diff --git a/SimPEG/Mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py index 0ca82ce1..bccbd4cc 100644 --- a/SimPEG/Mesh/BaseMesh.py +++ b/SimPEG/Mesh/BaseMesh.py @@ -596,6 +596,7 @@ class BaseRectangularMesh(BaseMesh): else: return switchKernal(x) + def getInterpolationMatMesh2Mesh(self, mesh2, locType='CC'): """ Interpolates variables from the current mesh to a new mesh (mesh2) @@ -606,15 +607,20 @@ class BaseRectangularMesh(BaseMesh): :return P: interpolation matrix """ + # import warnings + # warnings.warn( + # "`getInterpolationMatMesh2Mesh` will be slow. If you want to interpolate a vector from one mesh to another, use `InterpolateVecMesh2Mesh`", + # RuntimeWarning) + # Error Checking - if self._meshType == 'Cyl': + 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 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) + 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']: @@ -623,7 +629,7 @@ class BaseRectangularMesh(BaseMesh): # Vectors else: - if self._meshType == 'Cyl': + if self._meshType == 'CYL': if locType == 'F': X = self.getInterpolationMatMesh2Mesh(mesh2, locType='Fx') Z = self.getInterpolationMatMesh2Mesh(mesh2, locType='Fz') diff --git a/tests/mesh/test_InterpolationMesh2Mesh.py b/tests/mesh/test_InterpolationMesh2Mesh.py index e0b6a817..79d3e665 100644 --- a/tests/mesh/test_InterpolationMesh2Mesh.py +++ b/tests/mesh/test_InterpolationMesh2Mesh.py @@ -4,6 +4,10 @@ from SimPEG.Utils import mkvc from SimPEG import Mesh, Tests import unittest +test1D = False +test2D = True +test3D = False + 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]) @@ -18,237 +22,332 @@ cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_ TOL = 1e-7 -class TestInterpolationMesh2Mesh_Tensor1D(Tests.OrderTest): +if test1D: + class TestInterpolationMesh2Mesh_Tensor1D(Tests.OrderTest): - name = 'Mesh2Mesh Tensor1D' - meshSizes = [8, 16, 32] - meshTypes = ['uniformTensorMesh'] - meshDimension = 1 + name = 'Mesh2Mesh Tensor1D' + meshSizes = [8, 16, 32] + meshTypes = ['uniformTensorMesh'] + meshDimension = 1 - def getError(self): - funX = lambda x: np.cos(2*np.pi*x) + 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)) + 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 + 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) + 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_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_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_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() + def test_orderFx_1D(self): + self.type = 'Fx' + self.name = 'Mesh2Mesh Tensor1D: Fx' + self.orderTest() -class TestInterpolationMesh2Mesh_Tensor2D(Tests.OrderTest): +if test2D: + class TestInterpolationMesh2Mesh_Tensor2D(Tests.OrderTest): - name = 'Mesh2Mesh Tensor2D' - meshSizes = [4, 8, 16] - meshTypes = ['uniformTensorMesh'] - meshDimension = 2 + 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) + 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 ) + 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) + 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.projectEdgeVector(ana) 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) + 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.projectEdgeVector(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_Cyl(Tests.OrderTest): + + name = 'Mesh2Mesh Cyl' + meshSizes = [4, 8, 16] + meshTypes = ['uniformCylMesh'] + 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 'z' in self.type: + ana = call2(funY, getattr(mesh2, 'grid%s'%self.type)) + elif 'F' in self.type: + ana = cartF2Cyl(mesh2, funX, funY) + ana = np.c_[ana[:,0], np.zeros_like(ana[:,0]), ana[:,1]] + ana = mesh2.projectFaceVector(ana) + elif 'E' in self.type: + ana = cartE2Cyl(mesh2, funX, funY) + ana = np.c_[np.zeros_like(ana[:,1]),ana[:,1],np.zeros_like(ana[:,1])] + ana = mesh2.projectEdgeVector(ana) 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)) + ana = call2(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) + if 'F' in self.type: + v = cartF2Cyl(self.M, funX, funY) + v = np.c_[v[:,0], np.zeros_like(v[:,0]),v[:,1]] + if 'x' in self.type or 'z' in self.type: + v = self.M.projectFaceVector(v) + else: + v = np.c_[v[:,0], v[:,2]] + v = mkvc(v) + elif 'E' in self.type: + v = cartE2Cyl(self.M, funX, funY) + v = np.c_[np.zeros_like(v[:,1]), v[:,1],np.zeros_like(v[:,1])] + v = self.M.projectEdgeVector(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_Cyl(self): + self.type = 'CC' + self.name = 'Mesh2Mesh Tensor2D: CC' + self.orderTest() + + def test_orderN_Cyl(self): + self.type = 'N' + self.name = 'Mesh2Mesh Tensor2D: N' + self.orderTest() + + def test_orderE_Cyl(self): + self.type = 'E' + self.name = 'Mesh2Mesh Tensor2D: E' + self.orderTest() + + def test_orderEy_Cyl(self): + self.type = 'Ey' + self.name = 'Mesh2Mesh Tensor2D: Ey' + self.orderTest() + + def test_orderF_Cyl(self): + self.type = 'F' + self.name = 'Mesh2Mesh Tensor2D: F' + self.orderTest() + + def test_orderFx_Cyl(self): + self.type = 'Fx' + self.name = 'Mesh2Mesh Tensor2D: Fx' + self.orderTest() + + def test_orderFz_Cyl(self): + self.type = 'Fz' + self.name = 'Mesh2Mesh Tensor2D: Fz' + self.orderTest() + +if test3D: + 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: - 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) + ana = call3(funX, getattr(mesh2, 'grid%s'%self.type)) - P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type) - # print P.shape, v.shape - num = P*v - return np.linalg.norm((num - ana), np.inf) + 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) - def test_orderCC_3D(self): - self.type = 'CC' - self.name = 'Mesh2Mesh Tensor3D: CC' - self.orderTest() + P = self.M.getInterpolationMatMesh2Mesh(mesh2, locType=self.type) + # print P.shape, v.shape + num = P*v - def test_orderN_3D(self): - self.type = 'N' - self.name = 'Mesh2Mesh Tensor3D: N' - self.orderTest() + return np.linalg.norm((num - ana), np.inf) - def test_orderE_3D(self): - self.type = 'E' - self.name = 'Mesh2Mesh Tensor3D: E' - self.orderTest() + def test_orderCC_3D(self): + self.type = 'CC' + self.name = 'Mesh2Mesh Tensor3D: CC' + self.orderTest() - def test_orderEx_3D(self): - self.type = 'Ex' - self.name = 'Mesh2Mesh Tensor3D: Ex' - self.orderTest() + def test_orderN_3D(self): + self.type = 'N' + self.name = 'Mesh2Mesh Tensor3D: N' + self.orderTest() - def test_orderEy_3D(self): - self.type = 'Ey' - self.name = 'Mesh2Mesh Tensor3D: Ey' - self.orderTest() + def test_orderE_3D(self): + self.type = 'E' + self.name = 'Mesh2Mesh Tensor3D: E' + self.orderTest() - def test_orderEz_3D(self): - self.type = 'Ez' - self.name = 'Mesh2Mesh Tensor3D: Ez' - self.orderTest() + def test_orderEx_3D(self): + self.type = 'Ex' + self.name = 'Mesh2Mesh Tensor3D: Ex' + self.orderTest() - def test_orderF_3D(self): - self.type = 'F' - self.name = 'Mesh2Mesh Tensor3D: F' - self.orderTest() + def test_orderEy_3D(self): + self.type = 'Ey' + self.name = 'Mesh2Mesh Tensor3D: Ey' + self.orderTest() - def test_orderFx_3D(self): - self.type = 'Fx' - self.name = 'Mesh2Mesh Tensor3D: Fx' - 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() - 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__':