mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-01 10:19:17 +08:00
Merge pull request #88 from simpeg/feat/CylMeshInterpolation
cyl mesh interpolation
This commit is contained in:
Rowan Cockett
+61
-2
@@ -1,7 +1,7 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
from scipy.constants import pi
|
||||
from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, ddx, av, avExtrap
|
||||
from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, spzeros, ddx, av, avExtrap
|
||||
from TensorMesh import BaseTensorMesh
|
||||
from InnerProducts import InnerProducts
|
||||
from View import CylView
|
||||
@@ -27,11 +27,14 @@ class CylMesh(BaseTensorMesh, InnerProducts, CylView):
|
||||
|
||||
_unitDimensions = [1, 2*np.pi, 1]
|
||||
|
||||
def __init__(self, h, x0=None):
|
||||
def __init__(self, h, x0=None, cartesianOrigin=None):
|
||||
BaseTensorMesh.__init__(self, h, x0)
|
||||
assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi"
|
||||
if self.dim == 2:
|
||||
print 'Warning, a disk mesh has not been tested thoroughly.'
|
||||
cartesianOrigin = np.zeros(self.dim) if cartesianOrigin is None else cartesianOrigin
|
||||
assert len(cartesianOrigin) == self.dim, "cartesianOrigin must be the same length as the dimension of the mesh."
|
||||
self.cartesianOrigin = np.array(cartesianOrigin, dtype=float)
|
||||
|
||||
|
||||
@property
|
||||
@@ -327,6 +330,62 @@ class CylMesh(BaseTensorMesh, InnerProducts, CylView):
|
||||
raise NotImplementedError('wrapping in the averaging is not yet implemented')
|
||||
return self._aveF2CCV
|
||||
|
||||
def getInterpolationMatCartMesh(self, Mrect, locType='CC'):
|
||||
"""
|
||||
Takes a cartesian mesh and returns a projection to translate onto the cartesian grid.
|
||||
"""
|
||||
|
||||
assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes"
|
||||
|
||||
|
||||
if locType == 'F':
|
||||
# do this three times for each component
|
||||
X = self.getInterpolationMatCartMesh(Mrect, locType='Fx')
|
||||
Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy')
|
||||
Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz')
|
||||
return sp.vstack((X,Y,Z))
|
||||
if locType == 'E':
|
||||
X = self.getInterpolationMatCartMesh(Mrect, locType='Ex')
|
||||
Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey')
|
||||
Z = spzeros(Mrect.nEz, self.nE)
|
||||
return sp.vstack((X,Y,Z))
|
||||
|
||||
grid = getattr(Mrect, 'grid' + locType)
|
||||
# This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice
|
||||
theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2
|
||||
theta[theta < 0] += np.pi*2.0
|
||||
r = ((grid[:,0] - self.cartesianOrigin[0])**2 + (grid[:,1] - self.cartesianOrigin[1])**2)**0.5
|
||||
|
||||
if locType in ['CC', 'N', 'Fz', 'Ez']:
|
||||
G, proj = np.c_[r, theta, grid[:,2]], np.ones(r.size)
|
||||
else:
|
||||
dotMe = {
|
||||
'Fx': Mrect.normals[:Mrect.nFx,:],
|
||||
'Fy': Mrect.normals[Mrect.nFx:(Mrect.nFx+Mrect.nFy),:],
|
||||
'Fz': Mrect.normals[-Mrect.nFz:,:],
|
||||
'Ex': Mrect.tangents[:Mrect.nEx,:],
|
||||
'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:],
|
||||
'Ez': Mrect.tangents[-Mrect.nEz:,:],
|
||||
}[locType]
|
||||
if 'F' in locType:
|
||||
normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)]
|
||||
proj = ( normals * dotMe ).sum(axis=1)
|
||||
if 'E' in locType:
|
||||
tangents = np.c_[-np.sin(theta), np.cos(theta), np.zeros(theta.size)]
|
||||
proj = ( tangents * dotMe ).sum(axis=1)
|
||||
G = np.c_[r, theta, grid[:,2]]
|
||||
|
||||
interpType = locType
|
||||
if interpType == 'Fy':
|
||||
interpType = 'Fx'
|
||||
elif interpType == 'Ex':
|
||||
interpType = 'Ey'
|
||||
|
||||
Pc2r = self.getInterpolationMat(G, interpType)
|
||||
Proj = sdiag(proj)
|
||||
return Proj * Pc2r
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
|
||||
@@ -571,35 +571,58 @@ class DiffOperators(object):
|
||||
@property
|
||||
def aveF2CC(self):
|
||||
"Construct the averaging operator on cell faces to cell centers."
|
||||
if getattr(self, '_aveF2CC', None) is None:
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
self._aveF2CC = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])),
|
||||
sp.kron(av(n[1]), speye(n[0]))), format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
|
||||
kron3(speye(n[2]), av(n[1]), speye(n[0])),
|
||||
kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
|
||||
return self._aveF2CC
|
||||
|
||||
if(self.dim == 1):
|
||||
return self.aveFx2CC
|
||||
elif(self.dim == 2):
|
||||
return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC), format="csr")
|
||||
elif(self.dim == 3):
|
||||
return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr")
|
||||
|
||||
@property
|
||||
def aveF2CCV(self):
|
||||
"Construct the averaging operator on cell faces to cell centers."
|
||||
if getattr(self, '_aveF2CCV', None) is None:
|
||||
if(self.dim == 1):
|
||||
return self.aveFx2CC
|
||||
elif(self.dim == 2):
|
||||
return sp.block_diag((self.aveFx2CC, self.aveFy2CC), format="csr")
|
||||
elif(self.dim == 3):
|
||||
return sp.block_diag((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr")
|
||||
|
||||
@property
|
||||
def aveFx2CC(self):
|
||||
"Construct the averaging operator on cell faces in the x direction to cell centers."
|
||||
if getattr(self, '_aveFx2CC', None) is None:
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
self._aveF2CCV = av(n[0])
|
||||
self._aveFx2CC = av(n[0])
|
||||
elif(self.dim == 2):
|
||||
self._aveF2CCV = sp.block_diag((sp.kron(speye(n[1]), av(n[0])),
|
||||
sp.kron(av(n[1]), speye(n[0]))), format="csr")
|
||||
self._aveFx2CC = sp.kron(speye(n[1]), av(n[0]))
|
||||
elif(self.dim == 3):
|
||||
self._aveF2CCV = sp.block_diag((kron3(speye(n[2]), speye(n[1]), av(n[0])),
|
||||
kron3(speye(n[2]), av(n[1]), speye(n[0])),
|
||||
kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
|
||||
return self._aveF2CCV
|
||||
self._aveFx2CC = kron3(speye(n[2]), speye(n[1]), av(n[0]))
|
||||
return self._aveFx2CC
|
||||
|
||||
@property
|
||||
def aveFy2CC(self):
|
||||
"Construct the averaging operator on cell faces in the y direction to cell centers."
|
||||
if self.dim < 2: return None
|
||||
if getattr(self, '_aveFy2CC', None) is None:
|
||||
n = self.vnC
|
||||
if(self.dim == 2):
|
||||
self._aveFy2CC = sp.kron(av(n[1]), speye(n[0]))
|
||||
elif(self.dim == 3):
|
||||
self._aveFy2CC = kron3(speye(n[2]), av(n[1]), speye(n[0]))
|
||||
return self._aveFy2CC
|
||||
|
||||
@property
|
||||
def aveFz2CC(self):
|
||||
"Construct the averaging operator on cell faces in the z direction to cell centers."
|
||||
if self.dim < 3: return None
|
||||
if getattr(self, '_aveFz2CC', None) is None:
|
||||
n = self.vnC
|
||||
if(self.dim == 3):
|
||||
self._aveFz2CC = kron3(av(n[2]), speye(n[1]), speye(n[0]))
|
||||
return self._aveFz2CC
|
||||
|
||||
|
||||
@property
|
||||
def aveCC2F(self):
|
||||
@@ -620,36 +643,60 @@ class DiffOperators(object):
|
||||
@property
|
||||
def aveE2CC(self):
|
||||
"Construct the averaging operator on cell edges to cell centers."
|
||||
if getattr(self, '_aveE2CC', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
self._aveE2CC = speye(n[0])
|
||||
elif(self.dim == 2):
|
||||
self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])),
|
||||
sp.kron(speye(n[1]), av(n[0]))), format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
|
||||
kron3(av(n[2]), speye(n[1]), av(n[0])),
|
||||
kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
|
||||
return self._aveE2CC
|
||||
if(self.dim == 1):
|
||||
return self.aveEx2CC
|
||||
elif(self.dim == 2):
|
||||
return 0.5*sp.hstack((self.aveEx2CC, self.aveEy2CC), format="csr")
|
||||
elif(self.dim == 3):
|
||||
return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr")
|
||||
|
||||
@property
|
||||
def aveE2CCV(self):
|
||||
"Construct the averaging operator on cell edges to cell centers."
|
||||
if getattr(self, '_aveE2CCV', None) is None:
|
||||
if(self.dim == 1):
|
||||
return self.aveEx2CC
|
||||
elif(self.dim == 2):
|
||||
return sp.block_diag((self.aveEx2CC, self.aveEy2CC), format="csr")
|
||||
elif(self.dim == 3):
|
||||
return sp.block_diag((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr")
|
||||
|
||||
@property
|
||||
def aveEx2CC(self):
|
||||
"Construct the averaging operator on cell edges in the x direction to cell centers."
|
||||
if getattr(self, '_aveEx2CC', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
|
||||
self._aveEx2CC = speye(n[0])
|
||||
elif(self.dim == 2):
|
||||
self._aveE2CCV = sp.block_diag((sp.kron(av(n[1]), speye(n[0])),
|
||||
sp.kron(speye(n[1]), av(n[0]))), format="csr")
|
||||
self._aveEx2CC = sp.kron(av(n[1]), speye(n[0]))
|
||||
elif(self.dim == 3):
|
||||
self._aveE2CCV = sp.block_diag((kron3(av(n[2]), av(n[1]), speye(n[0])),
|
||||
kron3(av(n[2]), speye(n[1]), av(n[0])),
|
||||
kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
|
||||
return self._aveE2CCV
|
||||
self._aveEx2CC = kron3(av(n[2]), av(n[1]), speye(n[0]))
|
||||
return self._aveEx2CC
|
||||
|
||||
@property
|
||||
def aveEy2CC(self):
|
||||
"Construct the averaging operator on cell edges in the y direction to cell centers."
|
||||
if self.dim < 2: return None
|
||||
if getattr(self, '_aveEy2CC', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
if(self.dim == 2):
|
||||
self._aveEy2CC = sp.kron(speye(n[1]), av(n[0]))
|
||||
elif(self.dim == 3):
|
||||
self._aveEy2CC = kron3(av(n[2]), speye(n[1]), av(n[0]))
|
||||
return self._aveEy2CC
|
||||
|
||||
@property
|
||||
def aveEz2CC(self):
|
||||
"Construct the averaging operator on cell edges in the z direction to cell centers."
|
||||
if self.dim < 3: return None
|
||||
if getattr(self, '_aveEz2CC', None) is None:
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
if(self.dim == 3):
|
||||
self._aveEz2CC = kron3(speye(n[2]), av(n[1]), av(n[0]))
|
||||
return self._aveEz2CC
|
||||
|
||||
@property
|
||||
def aveN2CC(self):
|
||||
|
||||
@@ -168,21 +168,21 @@ class BaseTensorMesh(BaseRectangularMesh):
|
||||
|
||||
"""
|
||||
|
||||
if key is 'Fx':
|
||||
if key == 'Fx':
|
||||
ten = [self.vectorNx , self.vectorCCy, self.vectorCCz]
|
||||
elif key is 'Fy':
|
||||
elif key == 'Fy':
|
||||
ten = [self.vectorCCx, self.vectorNy , self.vectorCCz]
|
||||
elif key is 'Fz':
|
||||
elif key == 'Fz':
|
||||
ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ]
|
||||
elif key is 'Ex':
|
||||
elif key == 'Ex':
|
||||
ten = [self.vectorCCx, self.vectorNy , self.vectorNz ]
|
||||
elif key is 'Ey':
|
||||
elif key == 'Ey':
|
||||
ten = [self.vectorNx , self.vectorCCy, self.vectorNz ]
|
||||
elif key is 'Ez':
|
||||
elif key == 'Ez':
|
||||
ten = [self.vectorNx , self.vectorNy , self.vectorCCz]
|
||||
elif key is 'CC':
|
||||
elif key == 'CC':
|
||||
ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz]
|
||||
elif key is 'N':
|
||||
elif key == 'N':
|
||||
ten = [self.vectorNx , self.vectorNy , self.vectorNz ]
|
||||
|
||||
return [t for t in ten if t is not None]
|
||||
@@ -204,10 +204,12 @@ class BaseTensorMesh(BaseRectangularMesh):
|
||||
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]]
|
||||
tensors[1] = np.r_[tensors[1], 2.0*np.pi]
|
||||
|
||||
inside = np.ones(pts.shape[0],dtype=bool)
|
||||
for i, tensor in enumerate(tensors):
|
||||
inside = inside & (pts[:,i] >= tensor.min()) & (pts[:,i] <= tensor.max())
|
||||
TOL = np.diff(tensor).min() * 1.0e-10
|
||||
inside = inside & (pts[:,i] >= tensor.min()-TOL) & (pts[:,i] <= tensor.max()+TOL)
|
||||
return inside
|
||||
|
||||
def getInterpolationMat(self, loc, locType, zerosOutside=False):
|
||||
|
||||
+20
-3
@@ -109,6 +109,10 @@ class TensorView(object):
|
||||
vc = (self.aveN2CC*v).reshape(self.vnC, order='F')
|
||||
elif vType in ['Fx', 'Fy', 'Fz', 'Ex', 'Ey', 'Ez']:
|
||||
aveOp = 'ave' + vType[0] + '2CCV'
|
||||
# n = getattr(self,'vn'+vType[0])
|
||||
# if 'x' in vType: v = np.r_[v,np.zeros(n[1]),np.zeros(n[2])]
|
||||
# if 'y' in vType: v = np.r_[np.zeros(n[0]),v,np.zeros(n[2])]
|
||||
# if 'z' in vType: v = np.r_[np.zeros(n[0]),np.zeros(n[1]),v]
|
||||
v = getattr(self,aveOp)*v # average to cell centers
|
||||
ind_xyz = {'x':0,'y':1,'z':2}[vType[1]]
|
||||
vc = self.r(v.reshape((self.nC,-1),order='F'), 'CC','CC','M')[ind_xyz]
|
||||
@@ -190,9 +194,16 @@ class TensorView(object):
|
||||
M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
|
||||
|
||||
"""
|
||||
if type(vType) in [list, tuple]:
|
||||
assert ax is None, "cannot specify an axis to plot on with this function."
|
||||
fig, axs = plt.subplots(1,len(vType))
|
||||
out = []
|
||||
for vTypeI, ax in zip(vType, axs):
|
||||
out += [self.plotSlice(v,vType=vTypeI, normal=normal, ind=ind, grid=grid, view=view, ax=ax, clim=clim, showIt=False, pcolorOpts=pcolorOpts, streamOpts=streamOpts, gridOpts=gridOpts)]
|
||||
return out
|
||||
viewOpts = ['real','imag','abs','vec']
|
||||
normalOpts = ['X', 'Y', 'Z']
|
||||
vTypeOpts = ['CC', 'CCv','F','E']
|
||||
vTypeOpts = ['CC', 'CCv','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez']
|
||||
|
||||
# Some user error checking
|
||||
assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts)
|
||||
@@ -219,9 +230,15 @@ class TensorView(object):
|
||||
elif vType == 'CCv':
|
||||
assert view == 'vec', 'Other types for CCv not supported'
|
||||
else:
|
||||
# Now just deal with 'F' and 'E'
|
||||
# Now just deal with 'F' and 'E' (x,y,z, maybe...)
|
||||
aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC')
|
||||
v = getattr(self,aveOp)*v # average to cell centers (might be a vector)
|
||||
Av = getattr(self,aveOp)
|
||||
if v.size == Av.shape[1]:
|
||||
v = Av * v
|
||||
else:
|
||||
v = self.r(v,vType[0],vType) # get specific component
|
||||
v = Av * v
|
||||
# we should now be averaged to cell centers (might be a vector)
|
||||
v = self.r(v.reshape((self.nC,-1),order='F'),'CC','CC','M')
|
||||
if view == 'vec':
|
||||
outSlice = []
|
||||
|
||||
@@ -133,6 +133,80 @@ class TestCyl2DMesh(unittest.TestCase):
|
||||
def test_lightOperators(self):
|
||||
self.assertTrue(self.mesh.nodalGrad is None)
|
||||
|
||||
def test_getInterpMatCartMesh_Cells(self):
|
||||
|
||||
Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
|
||||
Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
|
||||
|
||||
mc = np.arange(Mc.nC)
|
||||
xr = np.linspace(0,0.4,50)
|
||||
xc = np.linspace(0,0.4,50) + 0.2
|
||||
Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'CC')
|
||||
Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC')
|
||||
Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC')
|
||||
|
||||
assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
|
||||
|
||||
def test_getInterpMatCartMesh_Faces(self):
|
||||
|
||||
Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
|
||||
Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
|
||||
|
||||
Pf = Mc.getInterpolationMatCartMesh(Mr, 'F')
|
||||
mf = np.ones(Mc.nF)
|
||||
|
||||
frect = Pf * mf
|
||||
|
||||
fxcc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx')
|
||||
fycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy')
|
||||
fzcc = Mr.r(frect, 'F', 'Fz')
|
||||
|
||||
indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
|
||||
indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
|
||||
|
||||
TOL = 1e-2
|
||||
assert np.abs(float(fxcc[indX]) - 1) < TOL
|
||||
assert np.abs(float(fxcc[indY]) - 0) < TOL
|
||||
assert np.abs(float(fycc[indX]) - 0) < TOL
|
||||
assert np.abs(float(fycc[indY]) - 1) < TOL
|
||||
assert np.abs((fzcc - 1).sum()) < TOL
|
||||
|
||||
mag = (fxcc**2 + fycc**2)**0.5
|
||||
dist = ((Mr.gridCC[:,0] + 0.2)**2 + (Mr.gridCC[:,1] + 0.2)**2)**0.5
|
||||
|
||||
assert np.abs(mag[dist > 0.1].max() - 1) < TOL
|
||||
assert np.abs(mag[dist > 0.1].min() - 1) < TOL
|
||||
|
||||
|
||||
def test_getInterpMatCartMesh_Edges(self):
|
||||
|
||||
Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
|
||||
Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
|
||||
|
||||
Pe = Mc.getInterpolationMatCartMesh(Mr, 'E')
|
||||
me = np.ones(Mc.nE)
|
||||
|
||||
erect = Pe * me
|
||||
|
||||
excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex')
|
||||
eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey')
|
||||
ezcc = Mr.r(erect, 'E', 'Ez')
|
||||
|
||||
indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
|
||||
indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
|
||||
|
||||
TOL = 1e-2
|
||||
assert np.abs(float(excc[indX]) - 0) < TOL
|
||||
assert np.abs(float(excc[indY]) + 1) < TOL
|
||||
assert np.abs(float(eycc[indX]) - 1) < TOL
|
||||
assert np.abs(float(eycc[indY]) - 0) < TOL
|
||||
assert np.abs(ezcc.sum()) < TOL
|
||||
|
||||
mag = (excc**2 + eycc**2)**0.5
|
||||
dist = ((Mr.gridCC[:,0] + 0.2)**2 + (Mr.gridCC[:,1] + 0.2)**2)**0.5
|
||||
|
||||
assert np.abs(mag[dist > 0.1].max() - 1) < TOL
|
||||
assert np.abs(mag[dist > 0.1].min() - 1) < TOL
|
||||
|
||||
|
||||
MESHTYPES = ['uniformCylMesh']
|
||||
|
||||
Reference in New Issue
Block a user