From 5bcdc5d46acb32f90b3454b7399d9195233c703d Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 6 Mar 2014 19:05:58 -0800 Subject: [PATCH] rearranging of interpolation code. --- SimPEG/Mesh/CylMesh.py | 207 ++++++++++++++++++----------------- SimPEG/Mesh/TensorMesh.py | 159 ++++++++++++++------------- SimPEG/Tests/test_cylMesh.py | 1 - 3 files changed, 186 insertions(+), 181 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 267533a6..c44afebe 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -4,6 +4,37 @@ from scipy.constants import pi from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, ddx, av, avExtrap from TensorMesh import BaseTensorMesh + + # def getNearest(self, loc, locType): + # """ Returns the index of the closest face or edge to a given location + + # :param numpy.ndarray loc: Test point + # :param str locType: Type of location desired (see below) + # :rtype: int + # :return: ind: + + # locType can be:: + + # 'fz' -> location of nearest z-face + # 'fr' -> location of nearest r-face + # 'et' -> location of nearest edge + # """ + + # if locType=='et': + # dr = self.gridN[:,0] - loc[0] + # dz = self.gridN[:,1] - loc[1] + # elif locType=='fz': + # dr = self.gridFz[:,0] - loc[0] + # dz = self.gridFz[:,1] - loc[1] + # elif locType=='fr': + # dr = self.gridFr[:,0] - loc[0] + # dz = self.gridFr[:,1] - loc[1] + # else: + # raise ValueError('Invalid locType') + # R = np.sqrt(dr**2 + dz**2) + # ind = np.argmin(R) + # return ind + class CylMesh(BaseTensorMesh): """ CylMesh is a mesh class for cylindrical problems @@ -313,125 +344,97 @@ class CylMesh(BaseTensorMesh): """mass matrix for products of face functions w'*M(materialProp)*f""" return self.getMass(loc='f', materialProp=materialProp) - def getInterpolationMat(self, loc, locType='fz'): - """ Produces intrpolation matrix + # def getInterpolationMat(self, loc, locType='fz'): + # """ Produces intrpolation matrix - :param numpy.ndarray loc: Location of points to interpolate to - :param str locType: What to interpolate (see below) - :rtype: scipy.sparse.csr.csr_matrix - :return: M, the intrpolation matrix + # :param numpy.ndarray loc: Location of points to interpolate to + # :param str locType: What to interpolate (see below) + # :rtype: scipy.sparse.csr.csr_matrix + # :return: M, the intrpolation matrix - locType can be:: + # locType can be:: - 'fz' -> z-component of field defined on faces - 'fr' -> r-component of field defined on faces - 'et' -> theta-component of field defined on edges - """ + # 'fz' -> z-component of field defined on faces + # 'fr' -> r-component of field defined on faces + # 'et' -> theta-component of field defined on edges + # """ - loc = np.atleast_2d(loc) + # loc = np.atleast_2d(loc) - assert np.all(loc[:,0]<=self.vectorNx.max()) & \ - np.all(loc[:,1]>=self.vectorNz.min()) & \ - np.all(loc[:,1]<=self.vectorNz.max()), \ - "Points outside of mesh" + # assert np.all(loc[:,0]<=self.vectorNx.max()) & \ + # np.all(loc[:,1]>=self.vectorNz.min()) & \ + # np.all(loc[:,1]<=self.vectorNz.max()), \ + # "Points outside of mesh" - if locType=='fz': - Q = sp.lil_matrix((loc.shape[0], self.nF), dtype=float) + # if locType=='fz': + # Q = sp.lil_matrix((loc.shape[0], self.nF), dtype=float) - for i, iloc in enumerate(loc): - # Point is on a z-interface - if np.any(np.abs(self.vectorNz-iloc[1])<0.001): - dFz = self.gridFz-iloc #Distance to z faces - dFz[dFz[:,0]>0,:] = np.inf #Looking for next face to the left... - indL = np.argmin(np.sum(dFz**2, axis=1)) #Closest one - if self.gridFz[indL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation) - zFL = self.gridFz[indL,:] - zFLL = self.gridFz[indL-1,:] - Q[i, indL+self.nFr] = (iloc[0] - zFLL[0])/(zFL[0] - zFLL[0]) - Q[i, indL+self.nFr-1] = -(iloc[0] - zFL[0])/(zFL[0] - zFLL[0]) - else: - zFL = self.gridFz[indL,:] - zFR = self.gridFz[indL+1,:] - Q[i,indL+self.nFr] = (zFR[0] - iloc[0])/(zFR[0] - zFL[0]) - Q[i,indL+self.nFr+1] = (iloc[0] - zFL[0])/(zFR[0] - zFL[0]) - # Point is in a cell - else: - dFz = self.gridFz-iloc - dFz[dFz>0] = np.inf - dFz = np.sum(dFz**2, axis=1) + # for i, iloc in enumerate(loc): + # # Point is on a z-interface + # if np.any(np.abs(self.vectorNz-iloc[1])<0.001): + # dFz = self.gridFz-iloc #Distance to z faces + # dFz[dFz[:,0]>0,:] = np.inf #Looking for next face to the left... + # indL = np.argmin(np.sum(dFz**2, axis=1)) #Closest one + # if self.gridFz[indL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation) + # zFL = self.gridFz[indL,:] + # zFLL = self.gridFz[indL-1,:] + # Q[i, indL+self.nFr] = (iloc[0] - zFLL[0])/(zFL[0] - zFLL[0]) + # Q[i, indL+self.nFr-1] = -(iloc[0] - zFL[0])/(zFL[0] - zFLL[0]) + # else: + # zFL = self.gridFz[indL,:] + # zFR = self.gridFz[indL+1,:] + # Q[i,indL+self.nFr] = (zFR[0] - iloc[0])/(zFR[0] - zFL[0]) + # Q[i,indL+self.nFr+1] = (iloc[0] - zFL[0])/(zFR[0] - zFL[0]) + # # Point is in a cell + # else: + # dFz = self.gridFz-iloc + # dFz[dFz>0] = np.inf + # dFz = np.sum(dFz**2, axis=1) - indBL = np.argmin(dFz) # Face below and to the left - indAL = indBL + self.nCx # Face above and to the left + # indBL = np.argmin(dFz) # Face below and to the left + # indAL = indBL + self.nCx # Face above and to the left - zF_BL = self.gridFz[indBL,:] - zF_AL = self.gridFz[indAL,:] + # zF_BL = self.gridFz[indBL,:] + # zF_AL = self.gridFz[indAL,:] - dzB = iloc[1] - zF_BL[1] # z-distance to face below - dzA = zF_AL[1] - iloc[1] # z-distance to face above + # dzB = iloc[1] - zF_BL[1] # z-distance to face below + # dzA = zF_AL[1] - iloc[1] # z-distance to face above - if self.gridFz[indBL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation) - zF_BLL = self.gridFz[indBL-1,:] - zF_ALL = self.gridFz[indAL-1,:] + # if self.gridFz[indBL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation) + # zF_BLL = self.gridFz[indBL-1,:] + # zF_ALL = self.gridFz[indAL-1,:] - DZ = zF_AL[1] - zF_BL[1] - DR = zF_AL[0] - zF_ALL[0] + # DZ = zF_AL[1] - zF_BL[1] + # DR = zF_AL[0] - zF_ALL[0] - drL = iloc[0] - zF_AL[0] - drLL = iloc[0] - zF_ALL[0] + # drL = iloc[0] - zF_AL[0] + # drLL = iloc[0] - zF_ALL[0] - Q[i, indBL+self.nFr-1] = -(1 - dzB/DZ)*(drL/DR) - Q[i, indBL+self.nFr] = (1 - dzB/DZ)*(drLL/DR) - Q[i, indAL+self.nFr-1] = -(dzB/DZ)*(drL/DR) - Q[i, indAL+self.nFr] = (dzB/DZ)*(drLL/DR) - else: - indBR = indBL+1 # Face below and to the right - indAR = indAL + 1 # Face above and to the right - zF_BR = self.gridFz[indBR,:] + # Q[i, indBL+self.nFr-1] = -(1 - dzB/DZ)*(drL/DR) + # Q[i, indBL+self.nFr] = (1 - dzB/DZ)*(drLL/DR) + # Q[i, indAL+self.nFr-1] = -(dzB/DZ)*(drL/DR) + # Q[i, indAL+self.nFr] = (dzB/DZ)*(drLL/DR) + # else: + # indBR = indBL+1 # Face below and to the right + # indAR = indAL + 1 # Face above and to the right + # zF_BR = self.gridFz[indBR,:] - drL = iloc[0] - zF_BL[0] # r-distance to face on left - drR = zF_BR[0] - iloc[0] # r-distance to face on right + # drL = iloc[0] - zF_BL[0] # r-distance to face on left + # drR = zF_BR[0] - iloc[0] # r-distance to face on right - drz = (drL + drR)*(dzB + dzA) - Q[i,indBL+self.nFr] = drR*dzA/drz - Q[i,indBR+self.nFr] = drL*dzA/drz - Q[i,indAL+self.nFr] = drR*dzB/drz - Q[i,indAR+self.nFr] = drL*dzB/drz + # drz = (drL + drR)*(dzB + dzA) + # Q[i,indBL+self.nFr] = drR*dzA/drz + # Q[i,indBR+self.nFr] = drL*dzA/drz + # Q[i,indAL+self.nFr] = drR*dzB/drz + # Q[i,indAR+self.nFr] = drL*dzB/drz - elif locType=='fr': - raise NotImplementedError('locType==fr') - elif locType=='et': - raise NotImplementedError('locType==et') - else: - raise ValueError('Invalid locType') - return Q.tocsr() + # elif locType=='fr': + # raise NotImplementedError('locType==fr') + # elif locType=='et': + # raise NotImplementedError('locType==et') + # else: + # raise ValueError('Invalid locType') + # return Q.tocsr() - def getNearest(self, loc, locType): - """ Returns the index of the closest face or edge to a given location - :param numpy.ndarray loc: Test point - :param str locType: Type of location desired (see below) - :rtype: int - :return: ind: - - locType can be:: - - 'fz' -> location of nearest z-face - 'fr' -> location of nearest r-face - 'et' -> location of nearest edge - """ - - if locType=='et': - dr = self.gridN[:,0] - loc[0] - dz = self.gridN[:,1] - loc[1] - elif locType=='fz': - dr = self.gridFz[:,0] - loc[0] - dz = self.gridFz[:,1] - loc[1] - elif locType=='fr': - dr = self.gridFr[:,0] - loc[0] - dz = self.gridFr[:,1] - loc[1] - else: - raise ValueError('Invalid locType') - R = np.sqrt(dr**2 + dz**2) - ind = np.argmin(R) - return ind diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 377c948b..cc3bcbb1 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -170,6 +170,87 @@ class BaseTensorMesh(BaseRectangularMesh): return [t for t in ten if t is not None] + # --------------- Methods --------------------- + + 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 + """ + + tensors = self.getTensor(locType) + if type(pts) == list: + pts = np.array(pts) + assert type(pts) == np.ndarray, "must be a numpy array" + if self.dim > 1: + assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" + elif len(pts.shape) == 1: + pts = pts[:,np.newaxis] + else: + assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" + + inside = np.ones(pts.shape[0],dtype=bool) + for i, tensor in enumerate(tensors): + inside = inside & (pts[:,i] >= tensor.min()) & (pts[:,i] <= tensor.max()) + return inside + + def getInterpolationMat(self, loc, locType, zerosOutside=False): + """ Produces interpolation matrix + + :param numpy.ndarray loc: 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 type(loc) == list: + loc = np.array(loc) + assert type(loc) == np.ndarray, "must be a numpy array" + if self.dim > 1: + assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" + elif len(loc.shape) == 1: + loc = loc[:,np.newaxis] + else: + assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" + + if zerosOutside is False: + assert np.all(self.isInside(loc)), "Points outside of mesh" + else: + indZeros = np.logical_not(self.isInside(loc)) + loc[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.' + nF_nE = self.vnF if 'F' in locType else self.vnE + components = [Utils.spzeros(loc.shape[0], n) for n in nF_nE] + components[ind] = Utils.interpmat(loc, *self.getTensor(locType)) + Q = sp.hstack(components) + elif locType in ['CC', 'N']: + Q = Utils.interpmat(loc, *self.getTensor(locType)) + else: + raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) + + if zerosOutside: + Q[indZeros, :] = 0 + + return Q.tocsr() + + class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts): """ @@ -320,84 +401,6 @@ class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts): self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2), Utils.mkvc(l3)] return self._edge - # --------------- Methods --------------------- - - 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 - """ - - tensors = self.getTensor(locType) - if type(pts) == list: - pts = np.array(pts) - assert type(pts) == np.ndarray, "must be a numpy array" - if self.dim > 1: - assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" - elif len(pts.shape) == 1: - pts = pts[:,np.newaxis] - else: - assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" - - inside = np.ones(pts.shape[0],dtype=bool) - for i, tensor in enumerate(tensors): - inside = inside & (pts[:,i] >= tensor.min()) & (pts[:,i] <= tensor.max()) - return inside - - def getInterpolationMat(self, loc, locType, zerosOutside=False): - """ Produces interpolation matrix - - :param numpy.ndarray loc: 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 type(loc) == list: - loc = np.array(loc) - assert type(loc) == np.ndarray, "must be a numpy array" - if self.dim > 1: - assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" - elif len(loc.shape) == 1: - loc = loc[:,np.newaxis] - else: - assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" - - if zerosOutside is False: - assert np.all(self.isInside(loc)), "Points outside of mesh" - else: - indZeros = np.logical_not(self.isInside(loc)) - loc[indZeros, :] = np.array([v.mean() for v in self.getTensor('CC')]) - - ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1 - if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind: - nF_nE = self.vnF if 'F' in locType else self.vnE - components = [Utils.spzeros(loc.shape[0], n) for n in nF_nE] - components[ind] = Utils.interpmat(loc, *self.getTensor(locType)) - Q = sp.hstack(components) - elif locType in ['CC', 'N']: - Q = Utils.interpmat(loc, *self.getTensor(locType)) - else: - raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) - if zerosOutside: - Q[indZeros, :] = 0 - return Q.tocsr() - - @property def faceBoundaryInd(self): """ diff --git a/SimPEG/Tests/test_cylMesh.py b/SimPEG/Tests/test_cylMesh.py index 89311bac..ed316775 100644 --- a/SimPEG/Tests/test_cylMesh.py +++ b/SimPEG/Tests/test_cylMesh.py @@ -163,7 +163,6 @@ class TestFaceDiv2D(OrderTest): F = self.M.projectFaceVector(Fc) divF = self.M.faceDiv.dot(F) - tm = Mesh.TensorMesh([self.M.nCx, self.M.nCz]) divF_anal = call3(sol, self.M.gridCC) err = np.linalg.norm((divF-divF_anal), np.inf)