test cyl2cly for mesh to mesh interpolation

This commit is contained in:
Lindsey Heagy
2016-03-31 11:16:14 -07:00
parent f972a50fce
commit bf300061bc
2 changed files with 311 additions and 206 deletions
+12 -6
View File
@@ -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')
+299 -200
View File
@@ -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__':