mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-01 11:10:27 +08:00
Updates to CylMesh. Innerproducts. Remove getMass.
This commit is contained in:
rowanc1
+53
-207
@@ -3,39 +3,10 @@ import scipy.sparse as sp
|
||||
from scipy.constants import pi
|
||||
from SimPEG.Utils import mkvc, ndgrid, sdiag, kron3, speye, ddx, av, avExtrap
|
||||
from TensorMesh import BaseTensorMesh
|
||||
from InnerProducts import InnerProducts
|
||||
|
||||
|
||||
# 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):
|
||||
class CylMesh(BaseTensorMesh, InnerProducts):
|
||||
"""
|
||||
CylMesh is a mesh class for cylindrical problems
|
||||
"""
|
||||
@@ -49,6 +20,10 @@ class CylMesh(BaseTensorMesh):
|
||||
assert self.dim == 3, "dim of mesh must equal 3, for a cylindrically symmetric mesh use [hx, 1, hz]"
|
||||
assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi"
|
||||
|
||||
@property
|
||||
def isSymmetric(self):
|
||||
return self.nCy == 1
|
||||
|
||||
@property
|
||||
def nNx(self):
|
||||
"""
|
||||
@@ -57,7 +32,7 @@ class CylMesh(BaseTensorMesh):
|
||||
:rtype: int
|
||||
:return: nNx
|
||||
"""
|
||||
if self.nCy == 1: return self.nCx
|
||||
if self.isSymmetric: return self.nCx
|
||||
return self.nCx + 1
|
||||
|
||||
@property
|
||||
@@ -68,7 +43,7 @@ class CylMesh(BaseTensorMesh):
|
||||
:rtype: int
|
||||
:return: nNy
|
||||
"""
|
||||
if self.nCy == 1: return 0
|
||||
if self.isSymmetric: return 0
|
||||
return self.nCy
|
||||
|
||||
@property
|
||||
@@ -89,7 +64,7 @@ class CylMesh(BaseTensorMesh):
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: vnEy or None if dim < 2
|
||||
"""
|
||||
nNx = self.nNx if self.nCy == 1 else self.nNx - 1
|
||||
nNx = self.nNx if self.isSymmetric else self.nNx - 1
|
||||
return np.r_[nNx, self.nCy, self.nNz]
|
||||
|
||||
@property
|
||||
@@ -100,7 +75,7 @@ class CylMesh(BaseTensorMesh):
|
||||
:rtype: numpy.array (dim, )
|
||||
:return: vnEz or None if nCy > 1
|
||||
"""
|
||||
if self.nCy == 1:
|
||||
if self.isSymmetric:
|
||||
return np.r_[self.nNx, self.nNy, self.nCz]
|
||||
else:
|
||||
return None
|
||||
@@ -113,7 +88,7 @@ class CylMesh(BaseTensorMesh):
|
||||
:rtype: int
|
||||
:return: nEz
|
||||
"""
|
||||
if self.nCy == 1:
|
||||
if self.isSymmetric:
|
||||
return self.vnEz.prod()
|
||||
return (np.r_[self.nNx-1, self.nNy, self.nCz]).prod() + self.nCz
|
||||
|
||||
@@ -130,14 +105,14 @@ class CylMesh(BaseTensorMesh):
|
||||
@property
|
||||
def vectorNx(self):
|
||||
"""Nodal grid vector (1D) in the x direction."""
|
||||
if self.nCy == 1:
|
||||
if self.isSymmetric:
|
||||
return self.hx.cumsum()
|
||||
return np.r_[0, self.hx].cumsum()
|
||||
|
||||
@property
|
||||
def vectorNy(self):
|
||||
"""Nodal grid vector (1D) in the y direction."""
|
||||
if self.nCy == 1:
|
||||
if self.isSymmetric:
|
||||
# There aren't really any nodes, but all the grids need
|
||||
# somewhere to live, why not zero?!
|
||||
return np.r_[0]
|
||||
@@ -147,7 +122,7 @@ class CylMesh(BaseTensorMesh):
|
||||
def edge(self):
|
||||
"""Edge lengths"""
|
||||
if getattr(self, '_edge', None) is None:
|
||||
if self.nCy == 1:
|
||||
if self.isSymmetric:
|
||||
self._edge = 2*pi*self.gridN[:,0]
|
||||
else:
|
||||
raise NotImplementedError('edges not yet implemented for 3D cyl mesh')
|
||||
@@ -186,7 +161,7 @@ class CylMesh(BaseTensorMesh):
|
||||
# Compute faceDivergence operator on faces
|
||||
D1 = self.faceDivx
|
||||
D3 = self.faceDivz
|
||||
if self.nCy == 1:
|
||||
if self.isSymmetric:
|
||||
D = sp.hstack((D1, D3), format="csr")
|
||||
elif self.nCy > 1:
|
||||
D2 = self.faceDivy
|
||||
@@ -237,7 +212,7 @@ class CylMesh(BaseTensorMesh):
|
||||
def nodalGrad(self):
|
||||
"""Construct gradient operator (nodes to edges)."""
|
||||
# Nodal grad does not make sense for cylindrically symmetric mesh.
|
||||
if self.nCy == 1: return None
|
||||
if self.isSymmetric: return None
|
||||
raise NotImplementedError('nodalGrad not yet implemented')
|
||||
|
||||
@property
|
||||
@@ -263,178 +238,49 @@ class CylMesh(BaseTensorMesh):
|
||||
self._edgeCurl = sp.diags(1/self.area,0)*sp.vstack((Dz, Dr))*sp.diags(self.edge,0)
|
||||
return self._edgeCurl
|
||||
|
||||
# @property
|
||||
# def aveE2CC(self):
|
||||
# """Averaging operator from cell edges to cell centers"""
|
||||
# if getattr(self, '_aveE2CC', None) is None:
|
||||
# if self.isSymmetric:
|
||||
# az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr')
|
||||
# ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr')
|
||||
# ar[0,0] = 1
|
||||
# self._aveE2CC = (0.5)*sp.kron(az, ar).T
|
||||
# else:
|
||||
# raise NotImplementedError('wrapping in the averaging is not yet implemented')
|
||||
# return self._aveE2CC
|
||||
|
||||
@property
|
||||
def aveE2CC(self):
|
||||
"""Averaging operator from cell edges to cell centres"""
|
||||
"Construct the averaging operator on cell edges to cell centers."
|
||||
if getattr(self, '_aveE2CC', None) is None:
|
||||
az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr')
|
||||
ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr')
|
||||
ar[0,0] = 1
|
||||
self._aveE2CC = sp.kron(az, ar).T
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
if self.isSymmetric:
|
||||
avR = av(n[0])[:,1:]
|
||||
avR[0,0] = 1.
|
||||
self._aveE2CC = (0.5)*sp.kron(av(n[2]), avR, format="csr")
|
||||
else:
|
||||
raise NotImplementedError('wrapping in the averaging is not yet implemented')
|
||||
# self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
|
||||
# kron3(av(n[2]), speye(n[1]), av(n[0])),
|
||||
# kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
|
||||
return self._aveE2CC
|
||||
|
||||
|
||||
@property
|
||||
def aveF2CC(self):
|
||||
"""Averaging operator from cell faces to cell centres"""
|
||||
"Construct the averaging operator on cell faces to cell centers."
|
||||
if getattr(self, '_aveF2CC', None) is None:
|
||||
az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr')
|
||||
ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr')
|
||||
ar[0,0] = 1
|
||||
Afr = sp.kron(sp.eye(self.nCz),ar)
|
||||
Afz = sp.kron(az,sp.eye(self.nCx))
|
||||
self._aveF2CC = sp.vstack((Afr,Afz)).T
|
||||
n = self.vnC
|
||||
if self.isSymmetric:
|
||||
avR = av(n[0])[:,1:]
|
||||
avR[0,0] = 1.
|
||||
self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[2]), avR), sp.kron(av(n[2]), speye(n[0]))), format="csr")
|
||||
else:
|
||||
raise NotImplementedError('wrapping in the averaging is not yet implemented')
|
||||
# self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
|
||||
# kron3(speye(n[2]), av(n[1]), speye(n[0])),
|
||||
# kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
|
||||
return self._aveF2CC
|
||||
|
||||
def getFaceMassDeriv(self):
|
||||
Av = self.aveF2CC
|
||||
return Av.T * sdiag(self.vol)
|
||||
|
||||
def getEdgeMassDeriv(self):
|
||||
Av = self.aveE2CC
|
||||
return Av.T * sdiag(self.vol)
|
||||
|
||||
|
||||
####################################################
|
||||
# Methods
|
||||
####################################################
|
||||
|
||||
|
||||
def getMass(self, materialProp=None, loc='e'):
|
||||
""" Produces mass matricies.
|
||||
|
||||
:param None,float,numpy.ndarray materialProp: property to be averaged (see below)
|
||||
:param str loc: Average to location: 'e'-edges, 'f'-faces
|
||||
:rtype: scipy.sparse.csr.csr_matrix
|
||||
:return: M, the mass matrix
|
||||
|
||||
materialProp can be::
|
||||
|
||||
None -> takes materialProp = 1 (default)
|
||||
float -> a constant value for entire domain
|
||||
numpy.ndarray -> if materialProp.size == self.nC
|
||||
3D property model
|
||||
if materialProp.size = self.nCz
|
||||
1D (layered eath) property model
|
||||
"""
|
||||
if materialProp is None:
|
||||
materialProp = np.ones(self.nC)
|
||||
elif type(materialProp) is float:
|
||||
materialProp = np.ones(self.nC)*materialProp
|
||||
elif materialProp.shape == (self.nCz,):
|
||||
materialProp = materialProp.repeat(self.nCx)
|
||||
materialProp = mkvc(materialProp)
|
||||
assert materialProp.shape == (self.nC,), "materialProp incorrect shape"
|
||||
|
||||
if loc=='e':
|
||||
Av = self.aveE2CC
|
||||
elif loc=='f':
|
||||
Av = self.aveF2CC
|
||||
else:
|
||||
raise ValueError('Invalid loc')
|
||||
|
||||
diag = Av.T * (self.vol * mkvc(materialProp))
|
||||
|
||||
return sdiag(diag)
|
||||
|
||||
def getEdgeMass(self, materialProp=None):
|
||||
"""mass matrix for products of edge functions w'*M(materialProp)*e"""
|
||||
return self.getMass(loc='e', materialProp=materialProp)
|
||||
|
||||
def getFaceMass(self, materialProp=None):
|
||||
"""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
|
||||
|
||||
# :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::
|
||||
|
||||
# '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)
|
||||
|
||||
# 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)
|
||||
|
||||
# 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
|
||||
|
||||
# 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
|
||||
|
||||
# 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]
|
||||
|
||||
# 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,:]
|
||||
|
||||
# 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
|
||||
|
||||
# elif locType=='fr':
|
||||
# raise NotImplementedError('locType==fr')
|
||||
# elif locType=='et':
|
||||
# raise NotImplementedError('locType==et')
|
||||
# else:
|
||||
# raise ValueError('Invalid locType')
|
||||
# return Q.tocsr()
|
||||
|
||||
|
||||
|
||||
@@ -623,11 +623,11 @@ class DiffOperators(object):
|
||||
raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
|
||||
elif(self.dim == 2):
|
||||
self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])),
|
||||
sp.kron(speye(n[1]), av(n[0]))), format="csr")
|
||||
sp.kron(speye(n[1]), av(n[0]))), format="csr")
|
||||
elif(self.dim == 3):
|
||||
self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
|
||||
kron3(av(n[2]), speye(n[1]), av(n[0])),
|
||||
kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
|
||||
kron3(av(n[2]), speye(n[1]), av(n[0])),
|
||||
kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
|
||||
return self._aveE2CC
|
||||
|
||||
@property
|
||||
@@ -696,57 +696,3 @@ class DiffOperators(object):
|
||||
kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr")
|
||||
return self._aveN2F
|
||||
|
||||
# --------------- Methods ---------------------
|
||||
|
||||
def getMass(self, materialProp=None, loc='e'):
|
||||
""" Produces mass matricies.
|
||||
|
||||
:param str loc: Average to location: 'e'-edges, 'f'-faces
|
||||
:param None,float,numpy.ndarray materialProp: property to be averaged (see below)
|
||||
:rtype: scipy.sparse.csr.csr_matrix
|
||||
:return: M, the mass matrix
|
||||
|
||||
materialProp can be::
|
||||
|
||||
None -> takes materialProp = 1 (default)
|
||||
float -> a constant value for entire domain
|
||||
numpy.ndarray -> if materialProp.size == self.nC
|
||||
3D property model
|
||||
if materialProp.size = self.nCz
|
||||
1D (layered eath) property model
|
||||
"""
|
||||
if materialProp is None:
|
||||
materialProp = np.ones(self.nC)
|
||||
elif type(materialProp) is float:
|
||||
materialProp = np.ones(self.nC)*materialProp
|
||||
elif materialProp.shape == (self.nCz,):
|
||||
materialProp = materialProp.repeat(self.nCx*self.nCy)
|
||||
materialProp = mkvc(materialProp)
|
||||
assert materialProp.shape == (self.nC,), "materialProp incorrect shape"
|
||||
|
||||
if loc=='e':
|
||||
Av = self.aveE2CC
|
||||
elif loc=='f':
|
||||
Av = self.aveF2CC
|
||||
else:
|
||||
raise ValueError('Invalid loc')
|
||||
|
||||
diag = Av.T * (self.vol * mkvc(materialProp))
|
||||
|
||||
return sdiag(diag)
|
||||
|
||||
def getEdgeMass(self, materialProp=None):
|
||||
"""mass matrix for products of edge functions w'*M(materialProp)*e"""
|
||||
return self.getMass(loc='e', materialProp=materialProp)
|
||||
|
||||
def getFaceMass(self, materialProp=None):
|
||||
"""mass matrix for products of face functions w'*M(materialProp)*f"""
|
||||
return self.getMass(loc='f', materialProp=materialProp)
|
||||
|
||||
def getFaceMassDeriv(self):
|
||||
Av = self.aveF2CC
|
||||
return Av.T * sdiag(self.vol)
|
||||
|
||||
def getEdgeMassDeriv(self):
|
||||
Av = self.aveE2CC
|
||||
return Av.T * sdiag(self.vol)
|
||||
|
||||
+112
-104
@@ -239,6 +239,8 @@ class BaseTensorMesh(BaseRectangularMesh):
|
||||
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))
|
||||
# remove any zero blocks (hstack complains)
|
||||
components = [comp for comp in components if comp.shape[1] > 0]
|
||||
Q = sp.hstack(components)
|
||||
elif locType in ['CC', 'N']:
|
||||
Q = Utils.interpmat(loc, *self.getTensor(locType))
|
||||
@@ -251,6 +253,116 @@ class BaseTensorMesh(BaseRectangularMesh):
|
||||
return Q.tocsr()
|
||||
|
||||
|
||||
def _fastFaceInnerProduct(self, materialProperty=None, invertProperty=False):
|
||||
"""
|
||||
Fast version of getFaceInnerProduct.
|
||||
This does not handle the case of a full tensor materialProperty.
|
||||
|
||||
:param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param bool returnP: returns the projection matrices
|
||||
:param bool invertProperty: inverts the material property
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
return self._fastInnerProduct('F', materialProperty=materialProperty, invertProperty=invertProperty)
|
||||
|
||||
|
||||
def _fastEdgeInnerProduct(self, materialProperty=None, invertProperty=False):
|
||||
"""
|
||||
Fast version of getEdgeInnerProduct.
|
||||
This does not handle the case of a full tensor materialProperty.
|
||||
|
||||
:param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param bool returnP: returns the projection matrices
|
||||
:param bool invertProperty: inverts the material property
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nE, nE)
|
||||
"""
|
||||
return self._fastInnerProduct('E', materialProperty=materialProperty, invertProperty=invertProperty)
|
||||
|
||||
|
||||
def _fastInnerProduct(self, AvType, materialProperty=None, invertProperty=False):
|
||||
"""
|
||||
Fast version of getFaceInnerProduct.
|
||||
This does not handle the case of a full tensor materialProperty.
|
||||
|
||||
:param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param str AvType: 'E' or 'F'
|
||||
:param bool returnP: returns the projection matrices
|
||||
:param bool invertProperty: inverts the material property
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
if materialProperty is None:
|
||||
materialProperty = np.ones(self.nC)
|
||||
|
||||
if invertProperty:
|
||||
materialProperty = 1./materialProperty
|
||||
|
||||
if Utils.isScalar(materialProperty):
|
||||
materialProperty = materialProperty*np.ones(self.nC)
|
||||
|
||||
if materialProperty.size == self.nC:
|
||||
Av = getattr(self, 'ave'+AvType+'2CC')
|
||||
Vprop = self.vol * Utils.mkvc(materialProperty)
|
||||
return self.dim * Utils.sdiag(Av.T * Vprop)
|
||||
if materialProperty.size == self.nC*self.dim:
|
||||
Av = getattr(self, 'ave'+AvType+'2CCV')
|
||||
V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
|
||||
return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty))
|
||||
|
||||
|
||||
def _fastFaceInnerProductDeriv(self, materialProperty=None, v=None):
|
||||
"""
|
||||
:param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
return self._fastInnerProductDeriv('F', materialProperty=materialProperty, v=v)
|
||||
|
||||
|
||||
def _fastEdgeInnerProductDeriv(self, materialProperty=None, v=None):
|
||||
"""
|
||||
:param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nE, nE)
|
||||
"""
|
||||
return self._fastInnerProductDeriv('E', materialProperty=materialProperty, v=v)
|
||||
|
||||
|
||||
def _fastInnerProductDeriv(self, AvType, materialProperty=None, v=None):
|
||||
"""
|
||||
:param str AvType: 'E' or 'F'
|
||||
:param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
if materialProperty is None:
|
||||
return 0.0
|
||||
|
||||
if Utils.isScalar(materialProperty):
|
||||
Av = getattr(self, 'ave'+AvType+'2CC')
|
||||
V = Utils.sdiag(self.vol)
|
||||
ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1))
|
||||
if v is None:
|
||||
return self.dim * Av.T * V * ones
|
||||
return Utils.sdiag(v) * self.dim * Av.T * V * ones
|
||||
|
||||
if materialProperty.size == self.nC:
|
||||
Av = getattr(self, 'ave'+AvType+'2CC')
|
||||
V = Utils.sdiag(self.vol)
|
||||
if v is None:
|
||||
return self.dim * Av.T * V
|
||||
return Utils.sdiag(v) * self.dim * Av.T * V
|
||||
|
||||
if materialProperty.size == self.nC*self.dim: # anisotropic
|
||||
Av = getattr(self, 'ave'+AvType+'2CCV')
|
||||
V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
|
||||
if v is None:
|
||||
return Av.T * V
|
||||
return Utils.sdiag(v) * Av.T * V
|
||||
|
||||
|
||||
|
||||
class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts):
|
||||
"""
|
||||
@@ -449,111 +561,7 @@ class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts):
|
||||
indzu = (self.gridCC[:,2]==max(self.gridCC[:,2]))
|
||||
return indxd, indxu, indyd, indyu, indzd, indzu
|
||||
|
||||
def _fastFaceInnerProduct(self, materialProperty=None, invertProperty=False):
|
||||
"""
|
||||
Fast version of getFaceInnerProduct.
|
||||
This does not handle the case of a full tensor materialProperty.
|
||||
|
||||
:param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param bool returnP: returns the projection matrices
|
||||
:param bool invertProperty: inverts the material property
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
return self._fastInnerProduct('F', materialProperty=materialProperty, invertProperty=invertProperty)
|
||||
|
||||
|
||||
def _fastEdgeInnerProduct(self, materialProperty=None, invertProperty=False):
|
||||
"""
|
||||
Fast version of getEdgeInnerProduct.
|
||||
This does not handle the case of a full tensor materialProperty.
|
||||
|
||||
:param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param bool returnP: returns the projection matrices
|
||||
:param bool invertProperty: inverts the material property
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nE, nE)
|
||||
"""
|
||||
return self._fastInnerProduct('E', materialProperty=materialProperty, invertProperty=invertProperty)
|
||||
|
||||
|
||||
def _fastInnerProduct(self, AvType, materialProperty=None, invertProperty=False):
|
||||
"""
|
||||
Fast version of getFaceInnerProduct.
|
||||
This does not handle the case of a full tensor materialProperty.
|
||||
|
||||
:param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param str AvType: 'E' or 'F'
|
||||
:param bool returnP: returns the projection matrices
|
||||
:param bool invertProperty: inverts the material property
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
if materialProperty is None:
|
||||
materialProperty = np.ones(self.nC)
|
||||
|
||||
if invertProperty:
|
||||
materialProperty = 1./materialProperty
|
||||
|
||||
if Utils.isScalar(materialProperty):
|
||||
materialProperty = materialProperty*np.ones(self.nC)
|
||||
|
||||
if materialProperty.size == self.nC:
|
||||
Av = getattr(self, 'ave'+AvType+'2CC')
|
||||
Vprop = self.vol * Utils.mkvc(materialProperty)
|
||||
return self.dim * Utils.sdiag(Av.T * Vprop)
|
||||
if materialProperty.size == self.nC*self.dim:
|
||||
Av = getattr(self, 'ave'+AvType+'2CCV')
|
||||
V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
|
||||
return Utils.sdiag(Av.T * V * Utils.mkvc(materialProperty))
|
||||
|
||||
|
||||
def _fastFaceInnerProductDeriv(self, materialProperty=None, v=None):
|
||||
"""
|
||||
:param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
return self._fastInnerProductDeriv('F', materialProperty=materialProperty, v=v)
|
||||
|
||||
|
||||
def _fastEdgeInnerProductDeriv(self, materialProperty=None, v=None):
|
||||
"""
|
||||
:param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nE, nE)
|
||||
"""
|
||||
return self._fastInnerProductDeriv('E', materialProperty=materialProperty, v=v)
|
||||
|
||||
|
||||
def _fastInnerProductDeriv(self, AvType, materialProperty=None, v=None):
|
||||
"""
|
||||
:param str AvType: 'E' or 'F'
|
||||
:param numpy.array materialProperty: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
if materialProperty is None:
|
||||
return None
|
||||
if Utils.isScalar(materialProperty):
|
||||
Av = getattr(self, 'ave'+AvType+'2CC')
|
||||
V = Utils.sdiag(self.vol)
|
||||
ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1))
|
||||
if v is None:
|
||||
return self.dim * Av.T * V * ones
|
||||
return Utils.sdiag(v) * self.dim * Av.T * V * ones
|
||||
if materialProperty.size == self.nC:
|
||||
Av = getattr(self, 'ave'+AvType+'2CC')
|
||||
V = Utils.sdiag(self.vol)
|
||||
if v is None:
|
||||
return self.dim * Av.T * V
|
||||
return Utils.sdiag(v) * self.dim * Av.T * V
|
||||
if materialProperty.size == self.nC*self.dim: # anisotropic
|
||||
Av = getattr(self, 'ave'+AvType+'2CCV')
|
||||
V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
|
||||
if v is None:
|
||||
return Av.T * V
|
||||
return Utils.sdiag(v) * Av.T * V
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
+1
-1
@@ -337,7 +337,7 @@ if __name__ == '__main__':
|
||||
|
||||
D = M.faceDiv
|
||||
G = M.cellGrad
|
||||
Msig = M.getFaceMass()
|
||||
Msig = M.getFaceInnerProduct()
|
||||
A = D*Msig*G
|
||||
A[0,0] *= 10 # remove the constant null space from the matrix
|
||||
|
||||
|
||||
@@ -22,7 +22,7 @@ class TestSolver(unittest.TestCase):
|
||||
|
||||
D = M.faceDiv
|
||||
G = M.cellGrad
|
||||
Msig = M.getFaceMass()
|
||||
Msig = M.getFaceInnerProduct()
|
||||
A = D*Msig*G
|
||||
A[0,0] *= 10 # remove the constant null space from the matrix
|
||||
|
||||
|
||||
@@ -206,6 +206,62 @@ class TestEdgeCurl2D(OrderTest):
|
||||
def test_order(self):
|
||||
self.orderTest()
|
||||
|
||||
|
||||
# class TestInnerProducts2D(OrderTest):
|
||||
# """Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts."""
|
||||
|
||||
# meshTypes = MESHTYPES
|
||||
# meshDimension = 2
|
||||
# meshSizes = [4, 8, 16, 32, 64, 128]
|
||||
|
||||
# def getError(self):
|
||||
|
||||
# funR = lambda r, t, z: np.cos(2.0*np.pi*z)
|
||||
# funT = lambda r, t, z: 0*t
|
||||
# funZ = lambda r, t, z: np.sin(2.0*np.pi*r)
|
||||
|
||||
# call = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])
|
||||
|
||||
# sigma1 = lambda r, t, z: z+1
|
||||
# sigma2 = lambda r, t, z: r*z+50
|
||||
# sigma3 = lambda r, t, z: 3+t*r
|
||||
# sigma4 = lambda r, t, z: 0.1*r*t*z
|
||||
# sigma5 = lambda r, t, z: 0.2*z*r*t
|
||||
# sigma6 = lambda r, t, z: 0.1*t
|
||||
|
||||
# Gc = self.M.gridCC
|
||||
# if self.sigmaTest == 1:
|
||||
# sigma = np.c_[call(sigma1, Gc)]
|
||||
# analytic = 144877./360 # Found using sympy. z=5
|
||||
# elif self.sigmaTest == 2:
|
||||
# sigma = np.c_[call(sigma1, Gc), call(sigma2, Gc)]
|
||||
# analytic = 189959./120 # Found using sympy. z=5
|
||||
# elif self.sigmaTest == 3:
|
||||
# sigma = np.r_[call(sigma1, Gc), call(sigma2, Gc), call(sigma3, Gc)]
|
||||
# analytic = 781427./360 # Found using sympy. z=5
|
||||
|
||||
# if self.location == 'edges':
|
||||
# E = call(funT, self.M.gridEy)
|
||||
# A = self.M.getEdgeInnerProduct(sigma)
|
||||
# numeric = E.T.dot(A.dot(E))
|
||||
# elif self.location == 'faces':
|
||||
# Fr = call(funR, self.M.gridFx)
|
||||
# Fz = call(funZ, self.M.gridFz)
|
||||
# A = self.M.getFaceInnerProduct(sigma)
|
||||
# F = np.r_[Fr,Fz]
|
||||
# numeric = F.T.dot(A.dot(F))
|
||||
|
||||
# print numeric
|
||||
# err = np.abs(numeric - analytic)
|
||||
# return err
|
||||
|
||||
# def test_order1_faces(self):
|
||||
# self.name = "2D Face Inner Product - Isotropic"
|
||||
# self.location = 'faces'
|
||||
# self.sigmaTest = 1
|
||||
# self.orderTest()
|
||||
|
||||
|
||||
class TestCyl3DMesh(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
Reference in New Issue
Block a user