wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-04 13:13:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

Interpolation

Browse Source
This commit is contained in:
Rowan Cockett
2015-11-17 19:47:10 -08:00
parent 5afaa54ad9
commit 1718e3506d
6 changed files with 629 additions and 84 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Mesh/TreeMesh.py
+292 -18
View File
@@ -104,9 +104,12 @@ import time
MAX_BITS = 20
class TreeMesh(BaseMesh, InnerProducts):
_meshType = 'TREE'
def __init__(self, h_in, x0_in=None, levels=3):
assert type(h_in) is list, 'h_in must be a list'
assert len(h_in) > 1, "len(h_in) must be greater than 1"
assert len(h_in) in [2,3], "There is only support for TreeMesh in 2D or 3D."
h = range(len(h_in))
for i, h_i in enumerate(h_in):
@@ -119,7 +122,7 @@ class TreeMesh(BaseMesh, InnerProducts):
assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i)
assert len(h_i) == 2**levels, "must make h and levels match"
h[i] = h_i[:] # make a copy.
self.h = h
self._h = h
x0 = np.zeros(len(h))
if x0_in is not None:
@@ -149,7 +152,12 @@ class TreeMesh(BaseMesh, InnerProducts):
@property
def __dirty__(self):
return self.__dirtyFaces__ or self.__dirtyEdges__ or self.__dirtyNodes__ or self.__dirtyHanging__ or self.__dirtySets__
return (self.__dirtyFaces__ or
self.__dirtyEdges__ or
self.__dirtyNodes__ or
self.__dirtyCells__ or
self.__dirtyHanging__ or
self.__dirtySets__)
@__dirty__.setter
def __dirty__(self, val):
@@ -157,6 +165,7 @@ class TreeMesh(BaseMesh, InnerProducts):
self.__dirtyFaces__ = True
self.__dirtyEdges__ = True
self.__dirtyNodes__ = True
self.__dirtyCells__ = True
self.__dirtyHanging__ = True
self.__dirtySets__ = True
@@ -175,6 +184,95 @@ class TreeMesh(BaseMesh, InnerProducts):
@property
def levels(self): return self._levels
@property
def h(self):
"""h is a list containing the cell widths of the tensor mesh in each dimension."""
return self._h
@property
def hx(self):
"Width of cells in the x direction"
return self._h[0]
@property
def hy(self):
"Width of cells in the y direction"
return self._h[1]
@property
def hz(self):
"Width of cells in the z direction"
return None if self.dim < 3 else self._h[2]
@property
def vectorNx(self):
"""Nodal grid vector (1D) in the x direction."""
return np.r_[0., self.hx.cumsum()] + self.x0[0]
@property
def vectorNy(self):
"""Nodal grid vector (1D) in the y direction."""
return np.r_[0., self.hy.cumsum()] + self.x0[1]
@property
def vectorNz(self):
"""Nodal grid vector (1D) in the z direction."""
return None if self.dim < 3 else np.r_[0., self.hz.cumsum()] + self.x0[2]
@property
def vectorCCx(self):
"""Cell-centered grid vector (1D) in the x direction."""
return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + self.x0[0]
@property
def vectorCCy(self):
"""Cell-centered grid vector (1D) in the y direction."""
return np.r_[0, self.hy[:-1].cumsum()] + self.hy*0.5 + self.x0[1]
@property
def vectorCCz(self):
"""Cell-centered grid vector (1D) in the z direction."""
return None if self.dim < 3 else np.r_[0, self.hz[:-1].cumsum()] + self.hz*0.5 + self.x0[2]
def getTensor(self, key):
""" Returns a tensor list.
:param str key: What tensor (see below)
:rtype: list
:return: list of the tensors that make up the mesh.
key can be::
'CC' -> scalar field defined on cell centers
'N' -> scalar field defined on nodes
'Fx' -> x-component of field defined on faces
'Fy' -> y-component of field defined on faces
'Fz' -> z-component of field defined on faces
'Ex' -> x-component of field defined on edges
'Ey' -> y-component of field defined on edges
'Ez' -> z-component of field defined on edges
"""
if key == 'Fx':
ten = [self.vectorNx , self.vectorCCy, self.vectorCCz]
elif key == 'Fy':
ten = [self.vectorCCx, self.vectorNy , self.vectorCCz]
elif key == 'Fz':
ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ]
elif key == 'Ex':
ten = [self.vectorCCx, self.vectorNy , self.vectorNz ]
elif key == 'Ey':
ten = [self.vectorNx , self.vectorCCy, self.vectorNz ]
elif key == 'Ez':
ten = [self.vectorNx , self.vectorNy , self.vectorCCz]
elif key == 'CC':
ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz]
elif key == 'N':
ten = [self.vectorNx , self.vectorNy , self.vectorNz ]
return [t for t in ten if t is not None]
@property
def nC(self): return len(self._cells)
@@ -282,6 +380,10 @@ class TreeMesh(BaseMesh, InnerProducts):
def ntF(self):
return self.ntFx + self.ntFy + (0 if self.dim == 2 else self.ntFz)
@property
def vntF(self):
return [self.ntFx, self.ntFy] + ([] if self.dim == 2 else [self.ntFz])
@property
def ntFx(self):
self.number()
@@ -302,6 +404,10 @@ class TreeMesh(BaseMesh, InnerProducts):
def ntE(self):
return self.ntEx + self.ntEy + (0 if self.dim == 2 else self.ntEz)
@property
def vntE(self):
return [self.ntEx, self.ntEy] + ([] if self.dim == 2 else [self.ntEz])
@property
def ntEx(self):
if self.dim == 2:return self.ntFy
@@ -491,6 +597,7 @@ class TreeMesh(BaseMesh, InnerProducts):
def _parentPointer(self, pointer):
if pointer[-1] == 0: return None
mod = self._levelWidth(pointer[-1] - 1)
return [p - (p % mod) for p in pointer[:-1]] + [pointer[-1]-1]
@@ -760,6 +867,14 @@ class TreeMesh(BaseMesh, InnerProducts):
self.__dirtySets__ = False
def _numberCells(self, force=False):
if not self.__dirtyCells__ and not force: return
self._cc2i = dict()
for ii, c in enumerate(sorted(self._cells)):
self._cc2i[c] = ii
self.__dirtyCells__ = False
def _numberNodes(self, force=False):
if not self.__dirtyNodes__ and not force: return
self._createNumberingSets(force=force)
@@ -878,6 +993,7 @@ class TreeMesh(BaseMesh, InnerProducts):
def _hanging(self, force=False):
if not self.__dirtyHanging__ and not force: return
self._numberCells(force=force)
self._numberNodes(force=force)
self._numberFaces(force=force)
self._numberEdges(force=force)
@@ -1649,7 +1765,7 @@ class TreeMesh(BaseMesh, InnerProducts):
self._aveN2CC = Av*Re
return self._aveN2CC
def _getFaceP(self, xFace, yFace, zFace):
ind1, ind2, ind3 = [], [], []
@@ -1723,12 +1839,165 @@ class TreeMesh(BaseMesh, InnerProducts):
return Pxxx
def isInside(self, pts, locType='N'):
"""
Determines if a set of points are inside a mesh.
:param numpy.ndarray pts: Location of points to test
:rtype numpy.ndarray
:return inside, numpy array of booleans
"""
pts = Utils.asArray_N_x_Dim(pts, self.dim)
tensors = self.getTensor(locType)
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):
TOL = np.diff(tensor).min() * 1.0e-10
inside = inside & (pts[:,i] >= tensor.min()-TOL) & (pts[:,i] <= tensor.max()+TOL)
return inside
def point2index(self, locs):
locs = Utils.asArray_N_x_Dim(locs, self.dim)
TOL = 1e-10
Nx = self.vectorNx
Ny = self.vectorNy
Nz = self.vectorNz
pointers = range(self.dim)
Nx = np.r_[Nx[0] - TOL, Nx[1:-1], Nx[-1] + TOL]
pointers[0] = np.searchsorted(Nx, locs[:,0])
Ny = np.r_[Ny[0] - TOL, Ny[1:-1], Ny[-1] + TOL]
pointers[1] = np.searchsorted(Ny, locs[:,1])
if self.dim == 3:
Nz = np.r_[Nz[0] - TOL, Nz[1:-1], Nz[-1] + TOL]
pointers[2] = np.searchsorted(Nz, locs[:,2])
if np.any([np.any(P == len(N)) or np.any(P == 0) for P,N in zip(pointers,[Nx,Ny,Nz])]):
raise Exception('There are points outside of the mesh.')
out = []
for pointer in zip(*pointers):
for level in range(self.levels+1):
width = self._levelWidth(level)
testPointer = [((p-1)//width)*width for p in pointer] + [level]
test = self._index(testPointer)
if test in self:
out += [test]
break
return out
def getInterpolationMat(self, locs, locType, zerosOutside=False):
""" Produces interpolation matrix
:param numpy.ndarray locs: Location of points to interpolate to
:param str locType: What to interpolate (see below)
:rtype: scipy.sparse.csr.csr_matrix
:return: M, the interpolation matrix
locType can be::
'Ex' -> x-component of field defined on edges
'Ey' -> y-component of field defined on edges
'Ez' -> z-component of field defined on edges
'Fx' -> x-component of field defined on faces
'Fy' -> y-component of field defined on faces
'Fz' -> z-component of field defined on faces
'N' -> scalar field defined on nodes
'CC' -> scalar field defined on cell centers
"""
if 'E' in locType and self.dim == 2: raise Exception('Interpolation for edges is not supported in 2D.')
locs = Utils.asArray_N_x_Dim(locs, self.dim)
TOL = 1e-10
self.number()
cells = self.point2index(locs)
I,J,V=[],[],[]
numberer = getattr(self, '_'+locType.lower()+'2i')
if zerosOutside is False:
assert np.all(self.isInside(locs)), "Points outside of mesh"
else:
indZeros = np.logical_not(self.isInside(locs))
locs[indZeros, :] = np.array([v.mean() for v in self.getTensor('CC')])
if locType in ['Fx','Fy','Fz','Ex','Ey','Ez']:
ind = {'x':0, 'y':1, 'z':2}[locType[1]]
assert self.dim >= ind, 'mesh is not high enough dimension.'
antiInd = {'x':[1,2], 'y':[0,2], 'z':[0,1]}[locType[1]][:self.dim-1]
nF_nE = self.vntF if 'F' in locType else self.vntE
components = [Utils.spzeros(locs.shape[0], n) for n in nF_nE]
for ii, cell in enumerate(cells):
loc = locs[ii,:]
p = self._asPointer(cell)
h, n = self._cellH(p), self._cellN(p)
w = self._levelWidth(p[-1])
if 'E' in locType:
iLocs, weights = Utils.interputils._interpmat2D(np.array([(loc-n-self.x0)[antiInd]]),np.r_[0.,h[antiInd[0]]+TOL],np.r_[0.,h[antiInd[1]]+TOL])
newJ = [numberer[self._index([__+w*iLocs[IND][0] if _ == antiInd[0] else __+w*iLocs[IND][1] if _ == antiInd[1] else __ for _, __ in enumerate(p[:-1])] + [p[-1]])] for IND in range(4)] #sorry
elif 'F' in locType:
_, weights = Utils.interputils._interpmat1D(np.r_[(loc-n-self.x0)[ind]],np.r_[0.,h[ind]+TOL])
plusFace = self._index([__+w if _ == ind else __ for _, __ in enumerate(p[:-1])] + [p[-1]])
newJ = [numberer[cell], numberer[plusFace]]
I += [ii]*len(newJ)
J += newJ
V += weights
components[ind] = sp.csr_matrix((V,(I,J)), shape=(locs.shape[0], nF_nE[ind]))
# remove any zero blocks (hstack complains)
components = [comp for comp in components if comp.shape[1] > 0]
Q = sp.hstack(components).tocsr()
if 'E' in locType:
R = self._deflationMatrix(locType[0],asOnes=False,withHanging=True)
else: # faces
R = self._deflationMatrix(locType[0],asOnes=True,withHanging=True)
elif locType == 'N':
for ii, cell in enumerate(cells):
loc = locs[ii,:]
p = self._asPointer(cell)
h, n = self._cellH(p), self._cellN(p)
w = self._levelWidth(p[-1])
iLocs, weights = Utils.interputils._interpmat3D(np.array([(loc-n-self.x0)]),*[np.r_[0.,h[_]+TOL] for _ in range(3)])
newJ = [numberer[self._index([__+w*iLocs[IND][_] for _, __ in enumerate(p[:-1])] + [p[-1]])] for IND in range(8)] #sorry
I += [ii]*len(newJ)
J += newJ
V += weights
Q = sp.csr_matrix((V,(I,J)), shape=(locs.shape[0], self.ntN))
R = self._deflationMatrix('N',withHanging=True)
elif locType == 'CC':
for ii, cell in enumerate(cells):
I += [ii]
J += [numberer[cell]]
V += [1.0]
Q = sp.csr_matrix((V,(I,J)), shape=(locs.shape[0], self.nC))
R = Utils.Identity()
else:
raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))
if zerosOutside:
Q[indZeros, :] = 0
return Q * R
def plotGrid(self, ax=None, showIt=False,
grid=True,
cells=True, cellLine=False,
nodes=False,
facesX=False, facesY=False, facesZ=False,
edgesX=False, edgesY=False, edgesZ=False):
grid=True,
cells=True, cellLine=False,
nodes=False,
facesX=False, facesY=False, facesZ=False,
edgesX=False, edgesY=False, edgesZ=False):
# self.number()
@@ -1840,7 +2109,6 @@ class TreeMesh(BaseMesh, InnerProducts):
if showIt:plt.show()
def plotImage(self, I, ax=None, showIt=True):
if self.dim == 3: raise Exception()
@@ -1914,7 +2182,6 @@ def SortGrid(grid, offset=0):
return sorted(range(offset,grid.shape[0]+offset), key=K)
class NotBalancedException(Exception):
pass
@@ -1944,19 +2211,26 @@ if __name__ == '__main__':
# T = TreeMesh([[(1,128)],[(1,128)],[(1,128)]],levels=7)
# T = TreeMesh([128,128,128],levels=7)
T = TreeMesh([64,64],levels=6)
# T = TreeMesh([64,64],levels=6)
T = TreeMesh([4,4,4],levels=2)
# T = TreeMesh([[(1,128)],[(1,128)]],levels=7)
# T.refine(lambda xc:1, balance=False)
# T.refine(lambda xc:2, balance=False)
# T._index([0,0,0])
# T._pointer(0)
tic = time.time()
T.refine(function)#, balance=False)
print time.time() - tic
# tic = time.time()
# T.refine(function)#, balance=False)
# print time.time() - tic
# print T.nC
print T.nC
T.plotImage(T.vol,showIt=True)
P = T.getInterpolationMat([0.2,0,0], 'Ex')
print P.todense()
blah
# T.plotImage(np.arange(len(T.vol)),showIt=True)
# print T.getFaceInnerProduct()
# print T.gridFz
@@ -1972,7 +2246,7 @@ if __name__ == '__main__':
# T.__dirty__ = True
print T.gridFx.shape[0], T.nFx
# print T.gridFx.shape[0], T.nFx
SimPEG/Utils/interputils.py
+8 -8
View File
@@ -124,13 +124,13 @@ if not _interpCython:
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
inds += [( ind_x1, ind_y2),
( ind_x1, ind_y1),
inds += [( ind_x1, ind_y1),
( ind_x1, ind_y2),
( ind_x2, ind_y1),
( ind_x2, ind_y2)]
vals += [wx1*wy2,
wx1*wy1,
vals += [wx1*wy1,
wx1*wy2,
wx2*wy1,
wx2*wy2]
@@ -152,8 +152,8 @@ if not _interpCython:
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
ind_z1, ind_z2, wz1, wz2 = _interp_point_1D(z, locs[i, 2])
inds += [( ind_x1, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z1),
inds += [( ind_x1, ind_y1, ind_z1),
( ind_x1, ind_y2, ind_z1),
( ind_x2, ind_y1, ind_z1),
( ind_x2, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z2),
@@ -161,8 +161,8 @@ if not _interpCython:
( ind_x2, ind_y1, ind_z2),
( ind_x2, ind_y2, ind_z2)]
vals += [wx1*wy2*wz1,
wx1*wy1*wz1,
vals += [wx1*wy1*wz1,
wx1*wy2*wz1,
wx2*wy1*wz1,
wx2*wy2*wz1,
wx1*wy1*wz2,
SimPEG/Utils/interputils_cython.c
+49 -49
View File
@@ -2546,7 +2546,7 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_4_interpmat2D(CYTH
* ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
* ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1]) # <<<<<<<<<<<<<<
*
* inds += [( ind_x1, ind_y2),
* inds += [( ind_x1, ind_y1),
*/
__pyx_t_9 = __Pyx_GetModuleGlobalName(__pyx_n_s_interp_point_1D); if (unlikely(!__pyx_t_9)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 72; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_9);
@@ -2669,8 +2669,8 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_4_interpmat2D(CYTH
/* "SimPEG/Utils/interputils_cython.pyx":74
* ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
*
* inds += [( ind_x1, ind_y2), # <<<<<<<<<<<<<<
* ( ind_x1, ind_y1),
* inds += [( ind_x1, ind_y1), # <<<<<<<<<<<<<<
* ( ind_x1, ind_y2),
* ( ind_x2, ind_y1),
*/
__pyx_t_3 = PyTuple_New(2); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 74; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
@@ -2678,14 +2678,14 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_4_interpmat2D(CYTH
__Pyx_INCREF(__pyx_v_ind_x1);
__Pyx_GIVEREF(__pyx_v_ind_x1);
PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_v_ind_x1);
__Pyx_INCREF(__pyx_v_ind_y2);
__Pyx_GIVEREF(__pyx_v_ind_y2);
PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_v_ind_y2);
__Pyx_INCREF(__pyx_v_ind_y1);
__Pyx_GIVEREF(__pyx_v_ind_y1);
PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_v_ind_y1);
/* "SimPEG/Utils/interputils_cython.pyx":75
*
* inds += [( ind_x1, ind_y2),
* ( ind_x1, ind_y1), # <<<<<<<<<<<<<<
* inds += [( ind_x1, ind_y1),
* ( ind_x1, ind_y2), # <<<<<<<<<<<<<<
* ( ind_x2, ind_y1),
* ( ind_x2, ind_y2)]
*/
@@ -2694,13 +2694,13 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_4_interpmat2D(CYTH
__Pyx_INCREF(__pyx_v_ind_x1);
__Pyx_GIVEREF(__pyx_v_ind_x1);
PyTuple_SET_ITEM(__pyx_t_11, 0, __pyx_v_ind_x1);
__Pyx_INCREF(__pyx_v_ind_y1);
__Pyx_GIVEREF(__pyx_v_ind_y1);
PyTuple_SET_ITEM(__pyx_t_11, 1, __pyx_v_ind_y1);
__Pyx_INCREF(__pyx_v_ind_y2);
__Pyx_GIVEREF(__pyx_v_ind_y2);
PyTuple_SET_ITEM(__pyx_t_11, 1, __pyx_v_ind_y2);
/* "SimPEG/Utils/interputils_cython.pyx":76
* inds += [( ind_x1, ind_y2),
* ( ind_x1, ind_y1),
* inds += [( ind_x1, ind_y1),
* ( ind_x1, ind_y2),
* ( ind_x2, ind_y1), # <<<<<<<<<<<<<<
* ( ind_x2, ind_y2)]
*
@@ -2715,11 +2715,11 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_4_interpmat2D(CYTH
PyTuple_SET_ITEM(__pyx_t_8, 1, __pyx_v_ind_y1);
/* "SimPEG/Utils/interputils_cython.pyx":77
* ( ind_x1, ind_y1),
* ( ind_x1, ind_y2),
* ( ind_x2, ind_y1),
* ( ind_x2, ind_y2)] # <<<<<<<<<<<<<<
*
* vals += [wx1*wy2, wx1*wy1, wx2*wy1, wx2*wy2]
* vals += [wx1*wy1, wx1*wy2, wx2*wy1, wx2*wy2]
*/
__pyx_t_1 = PyTuple_New(2); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 77; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
@@ -2733,8 +2733,8 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_4_interpmat2D(CYTH
/* "SimPEG/Utils/interputils_cython.pyx":74
* ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
*
* inds += [( ind_x1, ind_y2), # <<<<<<<<<<<<<<
* ( ind_x1, ind_y1),
* inds += [( ind_x1, ind_y1), # <<<<<<<<<<<<<<
* ( ind_x1, ind_y2),
* ( ind_x2, ind_y1),
*/
__pyx_t_9 = PyList_New(4); if (unlikely(!__pyx_t_9)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 74; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
@@ -2760,13 +2760,13 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_4_interpmat2D(CYTH
/* "SimPEG/Utils/interputils_cython.pyx":79
* ( ind_x2, ind_y2)]
*
* vals += [wx1*wy2, wx1*wy1, wx2*wy1, wx2*wy2] # <<<<<<<<<<<<<<
* vals += [wx1*wy1, wx1*wy2, wx2*wy1, wx2*wy2] # <<<<<<<<<<<<<<
*
* return inds, vals
*/
__pyx_t_1 = PyNumber_Multiply(__pyx_v_wx1, __pyx_v_wy2); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 79; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__pyx_t_1 = PyNumber_Multiply(__pyx_v_wx1, __pyx_v_wy1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 79; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_1);
__pyx_t_9 = PyNumber_Multiply(__pyx_v_wx1, __pyx_v_wy1); if (unlikely(!__pyx_t_9)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 79; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__pyx_t_9 = PyNumber_Multiply(__pyx_v_wx1, __pyx_v_wy2); if (unlikely(!__pyx_t_9)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 79; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_9);
__pyx_t_8 = PyNumber_Multiply(__pyx_v_wx2, __pyx_v_wy1); if (unlikely(!__pyx_t_8)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 79; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_8);
@@ -2794,7 +2794,7 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_4_interpmat2D(CYTH
}
/* "SimPEG/Utils/interputils_cython.pyx":81
* vals += [wx1*wy2, wx1*wy1, wx2*wy1, wx2*wy2]
* vals += [wx1*wy1, wx1*wy2, wx2*wy1, wx2*wy2]
*
* return inds, vals # <<<<<<<<<<<<<<
*
@@ -3376,7 +3376,7 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_6_interpmat3D(CYTH
* ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
* ind_z1, ind_z2, wz1, wz2 = _interp_point_1D(z, locs[i, 2]) # <<<<<<<<<<<<<<
*
* inds += [( ind_x1, ind_y2, ind_z1),
* inds += [( ind_x1, ind_y1, ind_z1),
*/
__pyx_t_11 = __Pyx_GetModuleGlobalName(__pyx_n_s_interp_point_1D); if (unlikely(!__pyx_t_11)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 99; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_11);
@@ -3499,8 +3499,8 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_6_interpmat3D(CYTH
/* "SimPEG/Utils/interputils_cython.pyx":101
* ind_z1, ind_z2, wz1, wz2 = _interp_point_1D(z, locs[i, 2])
*
* inds += [( ind_x1, ind_y2, ind_z1), # <<<<<<<<<<<<<<
* ( ind_x1, ind_y1, ind_z1),
* inds += [( ind_x1, ind_y1, ind_z1), # <<<<<<<<<<<<<<
* ( ind_x1, ind_y2, ind_z1),
* ( ind_x2, ind_y1, ind_z1),
*/
__pyx_t_3 = PyTuple_New(3); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 101; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
@@ -3508,17 +3508,17 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_6_interpmat3D(CYTH
__Pyx_INCREF(__pyx_v_ind_x1);
__Pyx_GIVEREF(__pyx_v_ind_x1);
PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_v_ind_x1);
__Pyx_INCREF(__pyx_v_ind_y2);
__Pyx_GIVEREF(__pyx_v_ind_y2);
PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_v_ind_y2);
__Pyx_INCREF(__pyx_v_ind_y1);
__Pyx_GIVEREF(__pyx_v_ind_y1);
PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_v_ind_y1);
__Pyx_INCREF(__pyx_v_ind_z1);
__Pyx_GIVEREF(__pyx_v_ind_z1);
PyTuple_SET_ITEM(__pyx_t_3, 2, __pyx_v_ind_z1);
/* "SimPEG/Utils/interputils_cython.pyx":102
*
* inds += [( ind_x1, ind_y2, ind_z1),
* ( ind_x1, ind_y1, ind_z1), # <<<<<<<<<<<<<<
* inds += [( ind_x1, ind_y1, ind_z1),
* ( ind_x1, ind_y2, ind_z1), # <<<<<<<<<<<<<<
* ( ind_x2, ind_y1, ind_z1),
* ( ind_x2, ind_y2, ind_z1),
*/
@@ -3527,16 +3527,16 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_6_interpmat3D(CYTH
__Pyx_INCREF(__pyx_v_ind_x1);
__Pyx_GIVEREF(__pyx_v_ind_x1);
PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_v_ind_x1);
__Pyx_INCREF(__pyx_v_ind_y1);
__Pyx_GIVEREF(__pyx_v_ind_y1);
PyTuple_SET_ITEM(__pyx_t_1, 1, __pyx_v_ind_y1);
__Pyx_INCREF(__pyx_v_ind_y2);
__Pyx_GIVEREF(__pyx_v_ind_y2);
PyTuple_SET_ITEM(__pyx_t_1, 1, __pyx_v_ind_y2);
__Pyx_INCREF(__pyx_v_ind_z1);
__Pyx_GIVEREF(__pyx_v_ind_z1);
PyTuple_SET_ITEM(__pyx_t_1, 2, __pyx_v_ind_z1);
/* "SimPEG/Utils/interputils_cython.pyx":103
* inds += [( ind_x1, ind_y2, ind_z1),
* ( ind_x1, ind_y1, ind_z1),
* inds += [( ind_x1, ind_y1, ind_z1),
* ( ind_x1, ind_y2, ind_z1),
* ( ind_x2, ind_y1, ind_z1), # <<<<<<<<<<<<<<
* ( ind_x2, ind_y2, ind_z1),
* ( ind_x1, ind_y1, ind_z2),
@@ -3554,7 +3554,7 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_6_interpmat3D(CYTH
PyTuple_SET_ITEM(__pyx_t_8, 2, __pyx_v_ind_z1);
/* "SimPEG/Utils/interputils_cython.pyx":104
* ( ind_x1, ind_y1, ind_z1),
* ( ind_x1, ind_y2, ind_z1),
* ( ind_x2, ind_y1, ind_z1),
* ( ind_x2, ind_y2, ind_z1), # <<<<<<<<<<<<<<
* ( ind_x1, ind_y1, ind_z2),
@@ -3634,7 +3634,7 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_6_interpmat3D(CYTH
* ( ind_x2, ind_y1, ind_z2),
* ( ind_x2, ind_y2, ind_z2)] # <<<<<<<<<<<<<<
*
* vals += [wx1*wy2*wz1,
* vals += [wx1*wy1*wz1,
*/
__pyx_t_19 = PyTuple_New(3); if (unlikely(!__pyx_t_19)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 108; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_19);
@@ -3651,8 +3651,8 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_6_interpmat3D(CYTH
/* "SimPEG/Utils/interputils_cython.pyx":101
* ind_z1, ind_z2, wz1, wz2 = _interp_point_1D(z, locs[i, 2])
*
* inds += [( ind_x1, ind_y2, ind_z1), # <<<<<<<<<<<<<<
* ( ind_x1, ind_y1, ind_z1),
* inds += [( ind_x1, ind_y1, ind_z1), # <<<<<<<<<<<<<<
* ( ind_x1, ind_y2, ind_z1),
* ( ind_x2, ind_y1, ind_z1),
*/
__pyx_t_20 = PyList_New(8); if (unlikely(!__pyx_t_20)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 101; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
@@ -3690,11 +3690,11 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_6_interpmat3D(CYTH
/* "SimPEG/Utils/interputils_cython.pyx":110
* ( ind_x2, ind_y2, ind_z2)]
*
* vals += [wx1*wy2*wz1, # <<<<<<<<<<<<<<
* wx1*wy1*wz1,
* vals += [wx1*wy1*wz1, # <<<<<<<<<<<<<<
* wx1*wy2*wz1,
* wx2*wy1*wz1,
*/
__pyx_t_19 = PyNumber_Multiply(__pyx_v_wx1, __pyx_v_wy2); if (unlikely(!__pyx_t_19)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 110; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__pyx_t_19 = PyNumber_Multiply(__pyx_v_wx1, __pyx_v_wy1); if (unlikely(!__pyx_t_19)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 110; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_19);
__pyx_t_20 = PyNumber_Multiply(__pyx_t_19, __pyx_v_wz1); if (unlikely(!__pyx_t_20)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 110; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_20);
@@ -3702,20 +3702,20 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_6_interpmat3D(CYTH
/* "SimPEG/Utils/interputils_cython.pyx":111
*
* vals += [wx1*wy2*wz1,
* wx1*wy1*wz1, # <<<<<<<<<<<<<<
* vals += [wx1*wy1*wz1,
* wx1*wy2*wz1, # <<<<<<<<<<<<<<
* wx2*wy1*wz1,
* wx2*wy2*wz1,
*/
__pyx_t_19 = PyNumber_Multiply(__pyx_v_wx1, __pyx_v_wy1); if (unlikely(!__pyx_t_19)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 111; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__pyx_t_19 = PyNumber_Multiply(__pyx_v_wx1, __pyx_v_wy2); if (unlikely(!__pyx_t_19)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 111; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_19);
__pyx_t_18 = PyNumber_Multiply(__pyx_t_19, __pyx_v_wz1); if (unlikely(!__pyx_t_18)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 111; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_18);
__Pyx_DECREF(__pyx_t_19); __pyx_t_19 = 0;
/* "SimPEG/Utils/interputils_cython.pyx":112
* vals += [wx1*wy2*wz1,
* wx1*wy1*wz1,
* vals += [wx1*wy1*wz1,
* wx1*wy2*wz1,
* wx2*wy1*wz1, # <<<<<<<<<<<<<<
* wx2*wy2*wz1,
* wx1*wy1*wz2,
@@ -3727,7 +3727,7 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_6_interpmat3D(CYTH
__Pyx_DECREF(__pyx_t_19); __pyx_t_19 = 0;
/* "SimPEG/Utils/interputils_cython.pyx":113
* wx1*wy1*wz1,
* wx1*wy2*wz1,
* wx2*wy1*wz1,
* wx2*wy2*wz1, # <<<<<<<<<<<<<<
* wx1*wy1*wz2,
@@ -3794,8 +3794,8 @@ static PyObject *__pyx_pf_6SimPEG_5Utils_18interputils_cython_6_interpmat3D(CYTH
/* "SimPEG/Utils/interputils_cython.pyx":110
* ( ind_x2, ind_y2, ind_z2)]
*
* vals += [wx1*wy2*wz1, # <<<<<<<<<<<<<<
* wx1*wy1*wz1,
* vals += [wx1*wy1*wz1, # <<<<<<<<<<<<<<
* wx1*wy2*wz1,
* wx2*wy1*wz1,
*/
__pyx_t_19 = PyList_New(8); if (unlikely(!__pyx_t_19)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 110; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
SimPEG/Utils/interputils_cython.pyx
+7 -7
View File
@@ -71,12 +71,12 @@ def _interpmat2D(np.ndarray[np.float64_t, ndim=2] locs,
ind_x1, ind_x2, wx1, wx2 = _interp_point_1D(x, locs[i, 0])
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
inds += [( ind_x1, ind_y2),
( ind_x1, ind_y1),
inds += [( ind_x1, ind_y1),
( ind_x1, ind_y2),
( ind_x2, ind_y1),
( ind_x2, ind_y2)]
vals += [wx1*wy2, wx1*wy1, wx2*wy1, wx2*wy2]
vals += [wx1*wy1, wx1*wy2, wx2*wy1, wx2*wy2]
return inds, vals
@@ -98,8 +98,8 @@ def _interpmat3D(np.ndarray[np.float64_t, ndim=2] locs,
ind_y1, ind_y2, wy1, wy2 = _interp_point_1D(y, locs[i, 1])
ind_z1, ind_z2, wz1, wz2 = _interp_point_1D(z, locs[i, 2])
inds += [( ind_x1, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z1),
inds += [( ind_x1, ind_y1, ind_z1),
( ind_x1, ind_y2, ind_z1),
( ind_x2, ind_y1, ind_z1),
( ind_x2, ind_y2, ind_z1),
( ind_x1, ind_y1, ind_z2),
@@ -107,8 +107,8 @@ def _interpmat3D(np.ndarray[np.float64_t, ndim=2] locs,
( ind_x2, ind_y1, ind_z2),
( ind_x2, ind_y2, ind_z2)]
vals += [wx1*wy2*wz1,
wx1*wy1*wz1,
vals += [wx1*wy1*wz1,
wx1*wy2*wz1,
wx2*wy1*wz1,
wx2*wy2*wz1,
wx1*wy1*wz2,
tests/mesh/test_TreeInterpolation.py
+182
View File
@@ -0,0 +1,182 @@
import numpy as np
import unittest
from SimPEG import Utils, Tests
MESHTYPES = ['uniformTree'] #['randomTree', 'uniformTree']
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)))
cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, 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)))
plotIt = False
MESHTYPES = ['uniformTree','notatreeTree']
"""
Face interpolation is O(h)
Edge interpolation is O(h^2)
"""
class TestInterpolation2d(Tests.OrderTest):
name = "Interpolation 2D"
np.random.seed(1)
LOCS = np.random.rand(50,2)*0.6+0.2
# LOCS = np.c_[np.ones(100)*0.51, np.linspace(0.3,0.7,100)]
meshTypes = MESHTYPES
# tolerance = TOLERANCES
meshDimension = 2
meshSizes = [8, 16, 32]
expectedOrders = 1
def getError(self):
funX = lambda x, y: np.cos(2.*np.pi*y)*np.cos(2.*np.pi*x) + x
funY = lambda x, y: np.cos(2.*np.pi*x)*np.cos(2.*np.pi*y) + y
# self.LOCS = self.M.gridCC
if 'x' in self.type:
ana = call2(funX, self.LOCS)
elif 'y' in self.type:
ana = call2(funY, self.LOCS)
else:
ana = call2(funX, self.LOCS)
if 'F' in self.type:
Fc = cartF2(self.M, funX, funY)
grid = self.M.projectFaceVector(Fc)
elif 'E' in self.type:
Ec = cartE2(self.M, funX, funY)
grid = self.M.projectEdgeVector(Ec)
elif 'CC' == self.type:
grid = call2(funX, self.M.gridCC)
elif 'N' == self.type:
grid = call2(funX, self.M.gridN)
P = self.M.getInterpolationMat(self.LOCS, self.type)
# print P
comp = P*grid
err = np.linalg.norm((comp - ana), np.inf)
if plotIt:
import matplotlib.pyplot as plt
ax = plt.subplot(211)
self.M.plotGrid(ax=ax)
plt.plot(self.LOCS[:,0],self.LOCS[:,1], 'mx')
# ax = plt.subplot(111)
# self.M.plotImage(call2(funX, self.M.gridCC),ax=ax)
ax = plt.subplot(212)
plt.plot(self.LOCS[:,1],comp, 'bx')
plt.plot(self.LOCS[:,1],ana, 'ro')
plt.show()
return err
def test_orderFx(self):
self.type = 'Fx'
self.name = 'TreeMesh Interpolation 2D: Fx'
self.orderTest()
def test_orderFy(self):
self.type = 'Fy'
self.name = 'TreeMesh Interpolation 2D: Fy'
self.orderTest()
class TestInterpolation3D(Tests.OrderTest):
name = "Interpolation"
LOCS = np.random.rand(50,3)*0.6+0.2
meshTypes = MESHTYPES
# tolerance = TOLERANCES
meshDimension = 3
meshSizes = [8, 16]
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)
if 'x' in self.type:
ana = call3(funX, self.LOCS)
elif 'y' in self.type:
ana = call3(funY, self.LOCS)
elif 'z' in self.type:
ana = call3(funZ, self.LOCS)
else:
ana = call3(funX, self.LOCS)
if 'F' in self.type:
Fc = cartF3(self.M, funX, funY, funZ)
grid = self.M.projectFaceVector(Fc)
elif 'E' in self.type:
Ec = cartE3(self.M, funX, funY, funZ)
grid = self.M.projectEdgeVector(Ec)
elif 'CC' == self.type:
grid = call3(funX, self.M.gridCC)
elif 'N' == self.type:
grid = call3(funX, self.M.gridN)
comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid
err = np.linalg.norm((comp - ana), np.inf)
return err
def test_orderCC(self):
self.type = 'CC'
self.name = 'Interpolation 3D: CC'
self.expectedOrders = 1
self.orderTest()
self.expectedOrders = 2
def test_orderN(self):
self.type = 'N'
self.name = 'Interpolation 3D: N'
self.orderTest()
def test_orderFx(self):
self.type = 'Fx'
self.name = 'Interpolation 3D: Fx'
self.expectedOrders = 1
self.orderTest()
self.expectedOrders = 2
def test_orderFy(self):
self.type = 'Fy'
self.name = 'Interpolation 3D: Fy'
self.expectedOrders = 1
self.orderTest()
self.expectedOrders = 2
def test_orderFz(self):
self.type = 'Fz'
self.name = 'Interpolation 3D: Fz'
self.expectedOrders = 1
self.orderTest()
self.expectedOrders = 2
def test_orderEx(self):
self.type = 'Ex'
self.name = 'Interpolation 3D: Ex'
self.orderTest()
def test_orderEy(self):
self.type = 'Ey'
self.name = 'Interpolation 3D: Ey'
self.orderTest()
def test_orderEz(self):
self.type = 'Ez'
self.name = 'Interpolation 3D: Ez'
self.orderTest()
if __name__ == '__main__':
unittest.main()
tests/mesh/test_TreeMesh.py
+91 -2
View File
@@ -1,9 +1,9 @@
from SimPEG import Mesh
from SimPEG import Mesh, Tests
import numpy as np
import matplotlib.pyplot as plt
import unittest
TOL = 1e-10
TOL = 1e-8
@@ -153,6 +153,95 @@ class TestOcTree(unittest.TestCase):
assert np.max(np.abs((M.faceDiv * M.edgeCurl).todense().flatten())) < TOL
assert np.max(np.abs((Mr.faceDiv * Mr.edgeCurl).todense().flatten())) < TOL
class Test2DInterpolation(unittest.TestCase):
def setUp(self):
def topo(x):
return np.sin(x*(2.*np.pi))*0.3 + 0.5
def function(cell):
r = cell.center - np.array([0.5]*len(cell.center))
dist1 = np.sqrt(r.dot(r)) - 0.08
dist2 = np.abs(cell.center[-1] - topo(cell.center[0]))
dist = min([dist1,dist2])
# if dist < 0.05:
# return 5
if dist < 0.05:
return 6
if dist < 0.2:
return 5
if dist < 0.3:
return 4
if dist < 1.0:
return 3
else:
return 0
M = Mesh.TreeMesh([64,64],levels=6)
M.refine(function)
self.M = M
def test_fx(self):
r = np.random.rand(self.M.nFx)
P = self.M.getInterpolationMat(self.M.gridFx, 'Fx')
assert np.abs(P[:,:self.M.nFx]*r - r).max() < TOL
def test_fy(self):
r = np.random.rand(self.M.nFy)
P = self.M.getInterpolationMat(self.M.gridFy, 'Fy')
assert np.abs(P[:,self.M.nFx:]*r - r).max() < TOL
class Test3DInterpolation(unittest.TestCase):
def setUp(self):
def function(cell):
r = cell.center - np.array([0.5]*len(cell.center))
dist = np.sqrt(r.dot(r))
if dist < 0.2:
return 4
if dist < 0.3:
return 3
if dist < 1.0:
return 2
else:
return 0
M = Mesh.TreeMesh([16,16,16],levels=4)
M.refine(function)
# M.plotGrid(showIt=True)
self.M = M
def test_Fx(self):
r = np.random.rand(self.M.nFx)
P = self.M.getInterpolationMat(self.M.gridFx, 'Fx')
assert np.abs(P[:,:self.M.nFx]*r - r).max() < TOL
def test_Fy(self):
r = np.random.rand(self.M.nFy)
P = self.M.getInterpolationMat(self.M.gridFy, 'Fy')
assert np.abs(P[:,self.M.nFx:(self.M.nFx+self.M.nFy)]*r - r).max() < TOL
def test_Fz(self):
r = np.random.rand(self.M.nFz)
P = self.M.getInterpolationMat(self.M.gridFz, 'Fz')
assert np.abs(P[:,(self.M.nFx+self.M.nFy):]*r - r).max() < TOL
def test_Ex(self):
r = np.random.rand(self.M.nEx)
P = self.M.getInterpolationMat(self.M.gridEx, 'Ex')
assert np.abs(P[:,:self.M.nEx]*r - r).max() < TOL
def test_Ey(self):
r = np.random.rand(self.M.nEy)
P = self.M.getInterpolationMat(self.M.gridEy, 'Ey')
assert np.abs(P[:,self.M.nEx:(self.M.nEx+self.M.nEy)]*r - r).max() < TOL
def test_Ez(self):
r = np.random.rand(self.M.nEz)
P = self.M.getInterpolationMat(self.M.gridEz, 'Ez')
assert np.abs(P[:,(self.M.nEx+self.M.nEy):]*r - r).max() < TOL
if __name__ == '__main__':