mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-27 19:48:52 +08:00
mesh2mesh interpolation a bit slow and memory heavy right now -> might not want to return the mat, just perform the interpolation
This commit is contained in:
Lindsey Heagy
@@ -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])
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
+44
-35
@@ -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."""
|
||||
|
||||
@@ -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()
|
||||
Reference in New Issue
Block a user