diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 4fc9451f..4f4e98d3 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -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 diff --git a/SimPEG/Utils/interputils.py b/SimPEG/Utils/interputils.py index 8fff1dcb..17b9b469 100644 --- a/SimPEG/Utils/interputils.py +++ b/SimPEG/Utils/interputils.py @@ -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, diff --git a/SimPEG/Utils/interputils_cython.c b/SimPEG/Utils/interputils_cython.c index ab5fdd3a..d03d682c 100644 --- a/SimPEG/Utils/interputils_cython.c +++ b/SimPEG/Utils/interputils_cython.c @@ -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;} diff --git a/SimPEG/Utils/interputils_cython.pyx b/SimPEG/Utils/interputils_cython.pyx index f7d59bb3..accf8cd6 100644 --- a/SimPEG/Utils/interputils_cython.pyx +++ b/SimPEG/Utils/interputils_cython.pyx @@ -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, diff --git a/tests/mesh/test_TreeInterpolation.py b/tests/mesh/test_TreeInterpolation.py new file mode 100644 index 00000000..90860959 --- /dev/null +++ b/tests/mesh/test_TreeInterpolation.py @@ -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() diff --git a/tests/mesh/test_TreeMesh.py b/tests/mesh/test_TreeMesh.py index 82288f99..61c2f9be 100644 --- a/tests/mesh/test_TreeMesh.py +++ b/tests/mesh/test_TreeMesh.py @@ -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__':