mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-06 05:16:51 +08:00
rearranging of interpolation code.
This commit is contained in:
rowanc1
+105
-102
@@ -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
|
||||
|
||||
+81
-78
@@ -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):
|
||||
"""
|
||||
|
||||
@@ -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)
|
||||
|
||||
Reference in New Issue
Block a user