mirror of
https://github.com/wassname/simpeg.git
synced 2026-09-10 12:37:30 +08:00
InnerProducts working as an operator. Simplifications and generalizations in inner product code.
This commit is contained in:
rowanc1
+119
-138
@@ -10,112 +10,43 @@ class InnerProducts(object):
|
||||
def __init__(self):
|
||||
raise Exception('InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class.')
|
||||
|
||||
def getFaceInnerProduct(self, prop=None, returnP=False,
|
||||
invProp=False, invMat=False, doFast=True):
|
||||
def getFaceInnerProduct(self, prop=None, invProp=False, invMat=False, doFast=True):
|
||||
"""
|
||||
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param bool returnP: returns the projection matrices
|
||||
:param bool invProp: inverts the material property
|
||||
:param bool invMat: inverts the matrix
|
||||
:param bool doFast: do a faster implementation if available.
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
fast = None
|
||||
return self._getInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat, doFast=True)
|
||||
|
||||
if returnP is False and hasattr(self, '_fastFaceInnerProduct') and doFast:
|
||||
fast = self._fastFaceInnerProduct(prop=prop, invProp=invProp, invMat=invMat)
|
||||
|
||||
if fast is not None:
|
||||
return fast
|
||||
|
||||
if invProp:
|
||||
prop = invPropertyTensor(self, prop)
|
||||
|
||||
Mu = makePropertyTensor(self, prop)
|
||||
|
||||
d = self.dim
|
||||
# We will multiply by sqrt on each side to keep symmetry
|
||||
V = sp.kron(sp.identity(d), sdiag(np.sqrt((2**(-d))*self.vol)))
|
||||
|
||||
if d == 1:
|
||||
fP = _getFacePx(self)
|
||||
P000 = V*fP('fXm')
|
||||
P100 = V*fP('fXp')
|
||||
elif d == 2:
|
||||
fP = _getFacePxx(self)
|
||||
P000 = V*fP('fXm', 'fYm')
|
||||
P100 = V*fP('fXp', 'fYm')
|
||||
P010 = V*fP('fXm', 'fYp')
|
||||
P110 = V*fP('fXp', 'fYp')
|
||||
elif d == 3:
|
||||
fP = _getFacePxxx(self)
|
||||
P000 = V*fP('fXm', 'fYm', 'fZm')
|
||||
P100 = V*fP('fXp', 'fYm', 'fZm')
|
||||
P010 = V*fP('fXm', 'fYp', 'fZm')
|
||||
P110 = V*fP('fXp', 'fYp', 'fZm')
|
||||
P001 = V*fP('fXm', 'fYm', 'fZp')
|
||||
P101 = V*fP('fXp', 'fYm', 'fZp')
|
||||
P011 = V*fP('fXm', 'fYp', 'fZp')
|
||||
P111 = V*fP('fXp', 'fYp', 'fZp')
|
||||
|
||||
A = P000.T*Mu*P000 + P100.T*Mu*P100
|
||||
P = [P000, P100]
|
||||
|
||||
if d > 1:
|
||||
A = A + P010.T*Mu*P010 + P110.T*Mu*P110
|
||||
P += [P010, P110]
|
||||
if d > 2:
|
||||
A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111
|
||||
P += [P001, P101, P011, P111]
|
||||
|
||||
if invMat and tensorType(self, prop) < 3:
|
||||
A = sdInv(A)
|
||||
elif invMat and tensorType(self, prop) == 3:
|
||||
raise Exception('Solver needed to invert A.')
|
||||
|
||||
if returnP:
|
||||
return A, P
|
||||
else:
|
||||
return A
|
||||
|
||||
def getFaceInnerProductDeriv(self, prop=None, v=None, P=None, doFast=True):
|
||||
def getEdgeInnerProduct(self, prop=None, invProp=False, invMat=False, doFast=True):
|
||||
"""
|
||||
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param numpy.array v: vector to multiply (required in the general implementation)
|
||||
:param list P: list of projection matrices
|
||||
:param bool doFast: do a faster implementation if available.
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: dMdm, the derivative of the inner product matrix (nF, nC*nA)
|
||||
"""
|
||||
fast = None
|
||||
|
||||
if hasattr(self, '_fastFaceInnerProductDeriv') and doFast:
|
||||
fast = self._fastFaceInnerProductDeriv(prop=prop, v=v)
|
||||
|
||||
if fast is not None:
|
||||
return fast
|
||||
|
||||
if P is None:
|
||||
M, P = self.getFaceInnerProduct(prop=prop, returnP=True)
|
||||
|
||||
return self._getInnerProductDeriv(prop, v, P, self.nF)
|
||||
|
||||
def getEdgeInnerProduct(self, prop=None, returnP=False,
|
||||
invProp=False, invMat=False, doFast=True):
|
||||
"""
|
||||
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param bool returnP: returns the projection matrices
|
||||
:param bool invProp: inverts the material property
|
||||
:param bool invMat: inverts the matrix
|
||||
:param bool doFast: do a faster implementation if available.
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nE, nE)
|
||||
"""
|
||||
return self._getInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat)
|
||||
|
||||
def _getInnerProduct(self, projType, prop=None, invProp=False, invMat=False, doFast=True):
|
||||
"""
|
||||
:param str projType: 'F' for faces 'E' for edges
|
||||
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param bool invProp: inverts the material property
|
||||
:param bool invMat: inverts the matrix
|
||||
:param bool doFast: do a faster implementation if available.
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nE, nE)
|
||||
"""
|
||||
assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
|
||||
fast = None
|
||||
|
||||
if returnP is False and hasattr(self, '_fastEdgeInnerProduct') and doFast:
|
||||
fast = self._fastEdgeInnerProduct(prop=prop, invProp=invProp, invMat=invMat)
|
||||
if hasattr(self, '_fastInnerProduct') and doFast:
|
||||
fast = self._fastInnerProduct(projType, prop=prop, invProp=invProp, invMat=invMat)
|
||||
|
||||
if fast is not None:
|
||||
return fast
|
||||
@@ -123,72 +54,122 @@ class InnerProducts(object):
|
||||
if invProp:
|
||||
prop = invPropertyTensor(self, prop)
|
||||
|
||||
tensorType = TensorType(self, prop)
|
||||
Mu = makePropertyTensor(self, prop)
|
||||
|
||||
Ps = self._getInnerProductProjectionMatrices(projType, tensorType)
|
||||
|
||||
A = np.sum([P.T * Mu * P for P in Ps])
|
||||
|
||||
if invMat and tensorType < 3:
|
||||
A = sdInv(A)
|
||||
elif invMat and tensorType == 3:
|
||||
raise Exception('Solver needed to invert A.')
|
||||
|
||||
return A
|
||||
|
||||
def _getInnerProductProjectionMatrices(self, projType, tensorType):
|
||||
"""
|
||||
:param str projType: 'F' for faces 'E' for edges
|
||||
:param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
|
||||
"""
|
||||
assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.'
|
||||
assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
|
||||
|
||||
d = self.dim
|
||||
# We will multiply by sqrt on each side to keep symmetry
|
||||
V = sp.kron(sp.identity(d), sdiag(np.sqrt((2**(-d))*self.vol)))
|
||||
|
||||
if d == 1:
|
||||
raise NotImplementedError('getEdgeInnerProduct not implemented for 1D')
|
||||
elif d == 2:
|
||||
eP = _getEdgePxx(self)
|
||||
P000 = V*eP('eX0', 'eY0')
|
||||
P100 = V*eP('eX0', 'eY1')
|
||||
P010 = V*eP('eX1', 'eY0')
|
||||
P110 = V*eP('eX1', 'eY1')
|
||||
elif d == 3:
|
||||
eP = _getEdgePxxx(self)
|
||||
P000 = V*eP('eX0', 'eY0', 'eZ0')
|
||||
P100 = V*eP('eX0', 'eY1', 'eZ1')
|
||||
P010 = V*eP('eX1', 'eY0', 'eZ2')
|
||||
P110 = V*eP('eX1', 'eY1', 'eZ3')
|
||||
P001 = V*eP('eX2', 'eY2', 'eZ0')
|
||||
P101 = V*eP('eX2', 'eY3', 'eZ1')
|
||||
P011 = V*eP('eX3', 'eY2', 'eZ2')
|
||||
P111 = V*eP('eX3', 'eY3', 'eZ3')
|
||||
nodes = ['000', '100', '010', '110', '001', '101', '011', '111'][:2**d]
|
||||
|
||||
Mu = makePropertyTensor(self, prop)
|
||||
A = P000.T*Mu*P000 + P100.T*Mu*P100 + P010.T*Mu*P010 + P110.T*Mu*P110
|
||||
P = [P000, P100, P010, P110]
|
||||
if d == 3:
|
||||
A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111
|
||||
P += [P001, P101, P011, P111]
|
||||
if projType == 'F':
|
||||
locs = {
|
||||
'000': [('fXm',), ('fXm', 'fYm'), ('fXm', 'fYm', 'fZm')],
|
||||
'100': [('fXp',), ('fXp', 'fYm'), ('fXp', 'fYm', 'fZm')],
|
||||
'010': [ None , ('fXm', 'fYp'), ('fXm', 'fYp', 'fZm')],
|
||||
'110': [ None , ('fXp', 'fYp'), ('fXp', 'fYp', 'fZm')],
|
||||
'001': [ None , None , ('fXm', 'fYm', 'fZp')],
|
||||
'101': [ None , None , ('fXp', 'fYm', 'fZp')],
|
||||
'011': [ None , None , ('fXm', 'fYp', 'fZp')],
|
||||
'111': [ None , None , ('fXp', 'fYp', 'fZp')]
|
||||
}
|
||||
if d == 1:
|
||||
proj = _getFacePx(self)
|
||||
elif d == 2:
|
||||
proj = _getFacePxx(self)
|
||||
elif d == 3:
|
||||
proj = _getFacePxxx(self)
|
||||
|
||||
if invMat and tensorType(self, prop) < 3:
|
||||
A = sdInv(A)
|
||||
elif invMat and tensorType(self, prop) == 3:
|
||||
raise Exception('Solver needed to invert A.')
|
||||
elif projType == 'E':
|
||||
locs = {
|
||||
'000': [ None , ('eX0', 'eY0'), ('eX0', 'eY0', 'eZ0')],
|
||||
'100': [ None , ('eX0', 'eY1'), ('eX0', 'eY1', 'eZ1')],
|
||||
'010': [ None , ('eX1', 'eY0'), ('eX1', 'eY0', 'eZ2')],
|
||||
'110': [ None , ('eX1', 'eY1'), ('eX1', 'eY1', 'eZ3')],
|
||||
'001': [ None , None , ('eX2', 'eY2', 'eZ0')],
|
||||
'101': [ None , None , ('eX2', 'eY3', 'eZ1')],
|
||||
'011': [ None , None , ('eX3', 'eY2', 'eZ2')],
|
||||
'111': [ None , None , ('eX3', 'eY3', 'eZ3')]
|
||||
}
|
||||
if d == 1:
|
||||
raise NotImplementedError('getEdgeInnerProduct not implemented for 1D')
|
||||
elif d == 2:
|
||||
proj = _getEdgePxx(self)
|
||||
elif d == 3:
|
||||
proj = _getEdgePxxx(self)
|
||||
|
||||
if returnP:
|
||||
return A, P
|
||||
else:
|
||||
return A
|
||||
return [V*proj(*locs[node][d-1]) for node in nodes]
|
||||
|
||||
def getEdgeInnerProductDeriv(self, prop=None, v=None, P=None, doFast=True):
|
||||
|
||||
def getFaceInnerProductDeriv(self, tensorType, P=None, doFast=True):
|
||||
"""
|
||||
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param numpy.array v: vector to multiply (required in the general implementation)
|
||||
:param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
|
||||
:param list P: list of projection matrices
|
||||
:param bool doFast: do a faster implementation if available.
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
|
||||
:return: dMdm, the derivative of the inner product matrix (nF, nC*nA)
|
||||
"""
|
||||
|
||||
assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.'
|
||||
fast = None
|
||||
|
||||
if hasattr(self, '_fastEdgeInnerProductDeriv') and doFast:
|
||||
fast = self._fastEdgeInnerProductDeriv(prop=prop, v=v)
|
||||
if hasattr(self, '_fastInnerProductDeriv') and doFast:
|
||||
fast = self._fastInnerProductDeriv('F', tensorType)
|
||||
|
||||
if fast is not None:
|
||||
return fast
|
||||
|
||||
if P is None:
|
||||
M, P = self.getEdgeInnerProduct(prop=prop, returnP=True)
|
||||
P = self._getInnerProductProjectionMatrices('F', tensorType=tensorType)
|
||||
|
||||
return self._getInnerProductDeriv(prop, v, P, self.nE)
|
||||
def innerProductDeriv(v):
|
||||
return self._getInnerProductDeriv(tensorType, P, self.nF, v)
|
||||
return DerivOperator(innerProductDeriv)
|
||||
|
||||
def _getInnerProductDeriv(self, prop, v, P, n):
|
||||
|
||||
def getEdgeInnerProductDeriv(self, tensorType, P=None, doFast=True):
|
||||
"""
|
||||
:param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
|
||||
:param list P: list of projection matrices
|
||||
:param bool doFast: do a faster implementation if available.
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
|
||||
"""
|
||||
assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.'
|
||||
fast = None
|
||||
|
||||
if hasattr(self, '_fastInnerProductDeriv') and doFast:
|
||||
fast = self._fastInnerProductDeriv('E', tensorType)
|
||||
|
||||
if fast is not None:
|
||||
return fast
|
||||
|
||||
if P is None:
|
||||
P = self._getInnerProductProjectionMatrices('E', tensorType=tensorType)
|
||||
def innerProductDeriv(v):
|
||||
return self._getInnerProductDeriv(tensorType, P, self.nE, v)
|
||||
return DerivOperator(innerProductDeriv)
|
||||
|
||||
def _getInnerProductDeriv(self, tensorType, P, n, v):
|
||||
"""
|
||||
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param numpy.array v: vector to multiply (required in the general implementation)
|
||||
@@ -197,7 +178,7 @@ class InnerProducts(object):
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: dMdm, the derivative of the inner product matrix (n, nC*nA)
|
||||
"""
|
||||
if prop is None:
|
||||
if tensorType == -1:
|
||||
return None
|
||||
|
||||
if v is None:
|
||||
@@ -206,24 +187,24 @@ class InnerProducts(object):
|
||||
d = self.dim
|
||||
Z = spzeros(self.nC, self.nC)
|
||||
|
||||
if isScalar(prop):
|
||||
if tensorType == 0:
|
||||
dMdm = spzeros(n, 1)
|
||||
for i, p in enumerate(P):
|
||||
dMdm = dMdm + sp.csr_matrix((p.T * (p * v), (range(n), np.zeros(n))), shape=(n,1))
|
||||
if d == 1:
|
||||
if prop.size == self.nC:
|
||||
if tensorType == 1:
|
||||
dMdm = spzeros(n, self.nC)
|
||||
for i, p in enumerate(P):
|
||||
dMdm = dMdm + p.T * sdiag( p * v )
|
||||
elif d == 2:
|
||||
if prop.size == self.nC:
|
||||
if tensorType == 1:
|
||||
dMdm = spzeros(n, self.nC)
|
||||
for i, p in enumerate(P):
|
||||
Y = p * v
|
||||
y1 = Y[:self.nC]
|
||||
y2 = Y[self.nC:]
|
||||
dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 )))
|
||||
elif prop.size == self.nC*2:
|
||||
elif tensorType == 2:
|
||||
dMdms = [spzeros(n, self.nC) for _ in range(2)]
|
||||
for i, p in enumerate(P):
|
||||
Y = p * v
|
||||
@@ -232,7 +213,7 @@ class InnerProducts(object):
|
||||
dMdms[0] = dMdms[0] + p.T * sp.vstack(( sdiag( y1 ), Z))
|
||||
dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 )))
|
||||
dMdm = sp.hstack(dMdms)
|
||||
elif prop.size == self.nC*3:
|
||||
elif tensorType == 3:
|
||||
dMdms = [spzeros(n, self.nC) for _ in range(3)]
|
||||
for i, p in enumerate(P):
|
||||
Y = p * v
|
||||
@@ -243,7 +224,7 @@ class InnerProducts(object):
|
||||
dMdms[2] = dMdms[2] + p.T * sp.vstack(( sdiag( y2 ), sdiag( y1 )))
|
||||
dMdm = sp.hstack(dMdms)
|
||||
elif d == 3:
|
||||
if prop.size == self.nC:
|
||||
if tensorType == 1:
|
||||
dMdm = spzeros(n, self.nC)
|
||||
for i, p in enumerate(P):
|
||||
Y = p * v
|
||||
@@ -251,7 +232,7 @@ class InnerProducts(object):
|
||||
y2 = Y[self.nC:self.nC*2]
|
||||
y3 = Y[self.nC*2:]
|
||||
dMdm = dMdm + p.T * sp.vstack((sdiag( y1 ), sdiag( y2 ), sdiag( y3 )))
|
||||
elif prop.size == self.nC*3:
|
||||
elif tensorType == 2:
|
||||
dMdms = [spzeros(n, self.nC) for _ in range(3)]
|
||||
for i, p in enumerate(P):
|
||||
Y = p * v
|
||||
@@ -262,7 +243,7 @@ class InnerProducts(object):
|
||||
dMdms[1] = dMdms[1] + p.T * sp.vstack(( Z, sdiag( y2 ), Z))
|
||||
dMdms[2] = dMdms[2] + p.T * sp.vstack(( Z, Z, sdiag( y3 )))
|
||||
dMdm = sp.hstack(dMdms)
|
||||
elif prop.size == self.nC*6:
|
||||
elif tensorType == 3:
|
||||
dMdms = [spzeros(n, self.nC) for _ in range(6)]
|
||||
for i, p in enumerate(P):
|
||||
Y = p * v
|
||||
|
||||
+27
-70
@@ -260,49 +260,21 @@ class BaseTensorMesh(BaseRectangularMesh):
|
||||
return Q.tocsr()
|
||||
|
||||
|
||||
def _fastFaceInnerProduct(self, prop=None, invProp=False, invMat=False):
|
||||
def _fastInnerProduct(self, projType, prop=None, invProp=False, invMat=False):
|
||||
"""
|
||||
Fast version of getFaceInnerProduct.
|
||||
This does not handle the case of a full tensor prop.
|
||||
|
||||
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param str projType: 'E' or 'F'
|
||||
:param bool returnP: returns the projection matrices
|
||||
:param bool invProp: inverts the material property
|
||||
:param bool invMat: inverts the matrix
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
return self._fastInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat)
|
||||
assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
|
||||
|
||||
|
||||
def _fastEdgeInnerProduct(self, prop=None, invProp=False, invMat=False):
|
||||
"""
|
||||
Fast version of getEdgeInnerProduct.
|
||||
This does not handle the case of a full tensor prop.
|
||||
|
||||
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param bool returnP: returns the projection matrices
|
||||
:param bool invProp: inverts the material property
|
||||
:param bool invMat: inverts the matrix
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nE, nE)
|
||||
"""
|
||||
return self._fastInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat)
|
||||
|
||||
|
||||
def _fastInnerProduct(self, AvType, prop=None, invProp=False, invMat=False):
|
||||
"""
|
||||
Fast version of getFaceInnerProduct.
|
||||
This does not handle the case of a full tensor prop.
|
||||
|
||||
:param numpy.array prop: 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 invProp: inverts the material property
|
||||
:param bool invMat: inverts the matrix
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
if prop is None:
|
||||
prop = np.ones(self.nC)
|
||||
|
||||
@@ -313,11 +285,11 @@ class BaseTensorMesh(BaseRectangularMesh):
|
||||
prop = prop*np.ones(self.nC)
|
||||
|
||||
if prop.size == self.nC:
|
||||
Av = getattr(self, 'ave'+AvType+'2CC')
|
||||
Av = getattr(self, 'ave'+projType+'2CC')
|
||||
Vprop = self.vol * Utils.mkvc(prop)
|
||||
M = self.dim * Utils.sdiag(Av.T * Vprop)
|
||||
elif prop.size == self.nC*self.dim:
|
||||
Av = getattr(self, 'ave'+AvType+'2CCV')
|
||||
Av = getattr(self, 'ave'+projType+'2CCV')
|
||||
V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol))
|
||||
M = Utils.sdiag(Av.T * V * Utils.mkvc(prop))
|
||||
else:
|
||||
@@ -328,55 +300,40 @@ class BaseTensorMesh(BaseRectangularMesh):
|
||||
else:
|
||||
return M
|
||||
|
||||
def _fastFaceInnerProductDeriv(self, prop=None, v=None):
|
||||
def _fastInnerProductDeriv(self, projType, tensorType):
|
||||
"""
|
||||
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param str projType: 'E' or 'F'
|
||||
:param TensorType tensorType: type of the tensor
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: M, the inner product matrix (nF, nF)
|
||||
"""
|
||||
return self._fastInnerProductDeriv('F', prop=prop, v=v)
|
||||
|
||||
|
||||
def _fastEdgeInnerProductDeriv(self, prop=None, v=None):
|
||||
"""
|
||||
:param numpy.array prop: 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', prop=prop, v=v)
|
||||
|
||||
|
||||
def _fastInnerProductDeriv(self, AvType, prop=None, v=None):
|
||||
"""
|
||||
:param str AvType: 'E' or 'F'
|
||||
:param numpy.array prop: 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 prop is None:
|
||||
assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
|
||||
if tensorType == -1:
|
||||
return None
|
||||
|
||||
if Utils.isScalar(prop):
|
||||
Av = getattr(self, 'ave'+AvType+'2CC')
|
||||
if tensorType == 0:
|
||||
Av = getattr(self, 'ave'+projType+'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 v is None:
|
||||
# return self.dim * Av.T * V * ones
|
||||
def scalarInnerProductDeriv(v):
|
||||
return Utils.sdiag(v) * self.dim * Av.T * V * ones
|
||||
return Utils.DerivOperator(scalarInnerProductDeriv)
|
||||
|
||||
if prop.size == self.nC:
|
||||
Av = getattr(self, 'ave'+AvType+'2CC')
|
||||
if tensorType == 1:
|
||||
Av = getattr(self, 'ave'+projType+'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
|
||||
def isotropicInnerProductDeriv(v):
|
||||
return Utils.sdiag(v) * self.dim * Av.T * V
|
||||
return Utils.DerivOperator(isotropicInnerProductDeriv)
|
||||
|
||||
if prop.size == self.nC*self.dim: # anisotropic
|
||||
Av = getattr(self, 'ave'+AvType+'2CCV')
|
||||
if tensorType == 2: # anisotropic
|
||||
Av = getattr(self, 'ave'+projType+'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
|
||||
def anisotropicInnerProductDeriv(v):
|
||||
return Utils.sdiag(v) * Av.T * V
|
||||
return Utils.DerivOperator(anisotropicInnerProductDeriv)
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -6,130 +6,93 @@ from TestUtils import checkDerivative
|
||||
|
||||
class TestInnerProductsDerivs(unittest.TestCase):
|
||||
|
||||
def doTestFace(self, h, rep, vec, fast):
|
||||
def doTestFace(self, h, rep, fast):
|
||||
mesh = Mesh.TensorMesh(h)
|
||||
v = np.random.rand(mesh.nF)
|
||||
sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
|
||||
Md = mesh.getFaceInnerProductDeriv(Utils.TensorType(mesh, sig), doFast=fast)
|
||||
def fun(sig):
|
||||
M = mesh.getFaceInnerProduct(sig)
|
||||
if vec:
|
||||
Md = mesh.getFaceInnerProductDeriv(sig, v=v, doFast=fast)
|
||||
return M*v, Md
|
||||
Md = mesh.getFaceInnerProductDeriv(sig, doFast=fast)
|
||||
return M*v, Utils.sdiag(v)*Md
|
||||
sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
|
||||
return M*v, Md*v
|
||||
return checkDerivative(fun, sig, num=5, plotIt=False)
|
||||
|
||||
def doTestEdge(self, h, rep, vec, fast):
|
||||
def doTestEdge(self, h, rep, fast):
|
||||
mesh = Mesh.TensorMesh(h)
|
||||
v = np.random.rand(mesh.nE)
|
||||
sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
|
||||
Md = mesh.getEdgeInnerProductDeriv(Utils.TensorType(mesh, sig), doFast=fast)
|
||||
def fun(sig):
|
||||
M = mesh.getEdgeInnerProduct(sig)
|
||||
if vec:
|
||||
Md = mesh.getEdgeInnerProductDeriv(sig, v=v, doFast=fast)
|
||||
return M*v, Md
|
||||
Md = mesh.getEdgeInnerProductDeriv(sig, doFast=fast)
|
||||
return M*v, Utils.sdiag(v)*Md
|
||||
sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep)
|
||||
return M*v, Md*v
|
||||
return checkDerivative(fun, sig, num=5, plotIt=False)
|
||||
|
||||
def test_FaceIP_1D_float(self):
|
||||
self.assertTrue(self.doTestFace([10],0,True, False))
|
||||
self.assertTrue(self.doTestFace([10],0, False))
|
||||
def test_FaceIP_2D_float(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0,True, False))
|
||||
self.assertTrue(self.doTestFace([10, 4],0, False))
|
||||
def test_FaceIP_3D_float(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0,True, False))
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, False))
|
||||
def test_FaceIP_1D_isotropic(self):
|
||||
self.assertTrue(self.doTestFace([10],1,True, False))
|
||||
self.assertTrue(self.doTestFace([10],1, False))
|
||||
def test_FaceIP_2D_isotropic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1,True, False))
|
||||
self.assertTrue(self.doTestFace([10, 4],1, False))
|
||||
def test_FaceIP_3D_isotropic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1,True, False))
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, False))
|
||||
def test_FaceIP_2D_anisotropic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2,True, False))
|
||||
self.assertTrue(self.doTestFace([10, 4],2, False))
|
||||
def test_FaceIP_3D_anisotropic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3,True, False))
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, False))
|
||||
def test_FaceIP_2D_tensor(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],3,True, False))
|
||||
self.assertTrue(self.doTestFace([10, 4],3, False))
|
||||
def test_FaceIP_3D_tensor(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],6,True, False))
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],6, False))
|
||||
|
||||
def test_FaceIP_1D_float_fast(self):
|
||||
self.assertTrue(self.doTestFace([10],0, False, True))
|
||||
self.assertTrue(self.doTestFace([10],0, True))
|
||||
def test_FaceIP_2D_float_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, False, True))
|
||||
self.assertTrue(self.doTestFace([10, 4],0, True))
|
||||
def test_FaceIP_3D_float_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, False, True))
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, True))
|
||||
def test_FaceIP_1D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10],1, False, True))
|
||||
self.assertTrue(self.doTestFace([10],1, True))
|
||||
def test_FaceIP_2D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, False, True))
|
||||
self.assertTrue(self.doTestFace([10, 4],1, True))
|
||||
def test_FaceIP_3D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, False, True))
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, True))
|
||||
def test_FaceIP_2D_anisotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, False, True))
|
||||
self.assertTrue(self.doTestFace([10, 4],2, True))
|
||||
def test_FaceIP_3D_anisotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, False, True))
|
||||
|
||||
def test_FaceIP_1D_float_fast_vec(self):
|
||||
self.assertTrue(self.doTestFace([10],0, True, True))
|
||||
def test_FaceIP_2D_float_fast_vec(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, True, True))
|
||||
def test_FaceIP_3D_float_fast_vec(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, True, True))
|
||||
def test_FaceIP_1D_isotropic_fast_vec(self):
|
||||
self.assertTrue(self.doTestFace([10],1, True, True))
|
||||
def test_FaceIP_2D_isotropic_fast_vec(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, True, True))
|
||||
def test_FaceIP_3D_isotropic_fast_vec(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, True, True))
|
||||
def test_FaceIP_2D_anisotropic_fast_vec(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, True, True))
|
||||
def test_FaceIP_3D_anisotropic_fast_vec(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, True, True))
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, True))
|
||||
|
||||
def test_EdgeIP_2D_float(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0,True, False))
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, False))
|
||||
def test_EdgeIP_3D_float(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0,True, False))
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0, False))
|
||||
def test_EdgeIP_2D_isotropic(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1,True, False))
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, False))
|
||||
def test_EdgeIP_3D_isotropic(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1,True, False))
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1, False))
|
||||
def test_EdgeIP_2D_anisotropic(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2,True, False))
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, False))
|
||||
def test_EdgeIP_3D_anisotropic(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3,True, False))
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, False))
|
||||
def test_EdgeIP_2D_tensor(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],3,True, False))
|
||||
self.assertTrue(self.doTestEdge([10, 4],3, False))
|
||||
def test_EdgeIP_3D_tensor(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],6,True, False))
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],6, False))
|
||||
|
||||
def test_EdgeIP_2D_float_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, False, True))
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, True))
|
||||
def test_EdgeIP_3D_float_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0, False, True))
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0, True))
|
||||
def test_EdgeIP_2D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, False, True))
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, True))
|
||||
def test_EdgeIP_3D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1, False, True))
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1, True))
|
||||
def test_EdgeIP_2D_anisotropic_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, False, True))
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, True))
|
||||
def test_EdgeIP_3D_anisotropic_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, False, True))
|
||||
|
||||
def test_EdgeIP_2D_float_fast_vec(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, True, True))
|
||||
def test_EdgeIP_3D_float_fast_vec(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0, True, True))
|
||||
def test_EdgeIP_2D_isotropic_fast_vec(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, True, True))
|
||||
def test_EdgeIP_3D_isotropic_fast_vec(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1, True, True))
|
||||
def test_EdgeIP_2D_anisotropic_fast_vec(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, True, True))
|
||||
def test_EdgeIP_3D_anisotropic_fast_vec(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, True, True))
|
||||
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, True))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
@@ -160,7 +160,7 @@ class TestSequenceFunctions(unittest.TestCase):
|
||||
Z = B2*A - sp.identity(M.nC*2)
|
||||
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
|
||||
|
||||
def test_tensorType2D(self):
|
||||
def test_TensorType2D(self):
|
||||
M = Mesh.TensorMesh([6, 6])
|
||||
a1 = np.random.rand(M.nC)
|
||||
a2 = np.random.rand(M.nC)
|
||||
@@ -170,12 +170,12 @@ class TestSequenceFunctions(unittest.TestCase):
|
||||
prop3 = np.c_[a1, a2, a3]
|
||||
|
||||
for ii, prop in enumerate([4, prop1, prop2, prop3]):
|
||||
self.assertTrue(tensorType(M, prop) == ii)
|
||||
self.assertTrue(TensorType(M, prop) == ii)
|
||||
|
||||
self.assertRaises(Exception, tensorType, M, np.c_[a1, a2, a3, a3])
|
||||
self.assertTrue(tensorType(M, None) == -1)
|
||||
self.assertRaises(Exception, TensorType, M, np.c_[a1, a2, a3, a3])
|
||||
self.assertTrue(TensorType(M, None) == -1)
|
||||
|
||||
def test_tensorType3D(self):
|
||||
def test_TensorType3D(self):
|
||||
M = Mesh.TensorMesh([6, 6, 7])
|
||||
a1 = np.random.rand(M.nC)
|
||||
a2 = np.random.rand(M.nC)
|
||||
@@ -188,10 +188,10 @@ class TestSequenceFunctions(unittest.TestCase):
|
||||
prop3 = np.c_[a1, a2, a3, a4, a5, a6]
|
||||
|
||||
for ii, prop in enumerate([4, prop1, prop2, prop3]):
|
||||
self.assertTrue(tensorType(M, prop) == ii)
|
||||
self.assertTrue(TensorType(M, prop) == ii)
|
||||
|
||||
self.assertRaises(Exception, tensorType, M, np.c_[a1, a2, a3, a3])
|
||||
self.assertTrue(tensorType(M, None) == -1)
|
||||
self.assertRaises(Exception, TensorType, M, np.c_[a1, a2, a3, a3])
|
||||
self.assertTrue(TensorType(M, None) == -1)
|
||||
|
||||
|
||||
def test_invPropertyTensor3D(self):
|
||||
|
||||
+36
-21
@@ -251,25 +251,34 @@ def inv2X2BlockDiagonal(a11, a12, a21, a22, returnMatrix=True):
|
||||
return sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))),
|
||||
sp.hstack((sdiag(b21), sdiag(b22)))))
|
||||
|
||||
def tensorType(M, tensor):
|
||||
if tensor is None: # default is ones
|
||||
return -1
|
||||
|
||||
if isScalar(tensor):
|
||||
return 0
|
||||
|
||||
if tensor.size == M.nC:
|
||||
return 1
|
||||
|
||||
if ((M.dim == 2 and tensor.size == M.nC*2) or
|
||||
(M.dim == 3 and tensor.size == M.nC*3)):
|
||||
return 2
|
||||
|
||||
if ((M.dim == 2 and tensor.size == M.nC*3) or
|
||||
(M.dim == 3 and tensor.size == M.nC*6)):
|
||||
return 3
|
||||
|
||||
raise Exception('Unexpected shape of tensor')
|
||||
class TensorType(object):
|
||||
def __init__(self, M, tensor):
|
||||
if tensor is None: # default is ones
|
||||
self._tt = -1
|
||||
self._tts = 'none'
|
||||
elif isScalar(tensor):
|
||||
self._tt = 0
|
||||
self._tts = 'scalar'
|
||||
elif tensor.size == M.nC:
|
||||
self._tt = 1
|
||||
self._tts = 'isotropic'
|
||||
elif ((M.dim == 2 and tensor.size == M.nC*2) or
|
||||
(M.dim == 3 and tensor.size == M.nC*3)):
|
||||
self._tt = 2
|
||||
self._tts = 'anisotropic'
|
||||
elif ((M.dim == 2 and tensor.size == M.nC*3) or
|
||||
(M.dim == 3 and tensor.size == M.nC*6)):
|
||||
self._tt = 3
|
||||
self._tts = 'tensor'
|
||||
else:
|
||||
raise Exception('Unexpected shape of tensor')
|
||||
def __str__(self):
|
||||
return 'TensorType[%i]: %s' % (self._tt, self._tts)
|
||||
def __eq__(self, v): return self._tt == v
|
||||
def __le__(self, v): return self._tt <= v
|
||||
def __ge__(self, v): return self._tt >= v
|
||||
def __lt__(self, v): return self._tt < v
|
||||
def __gt__(self, v): return self._tt > v
|
||||
|
||||
def makePropertyTensor(M, tensor):
|
||||
if tensor is None: # default is ones
|
||||
@@ -278,7 +287,7 @@ def makePropertyTensor(M, tensor):
|
||||
if isScalar(tensor):
|
||||
tensor = tensor * np.ones(M.nC)
|
||||
|
||||
propType = tensorType(M, tensor)
|
||||
propType = TensorType(M, tensor)
|
||||
if propType == 1: # Isotropic!
|
||||
Sigma = sp.kron(sp.identity(M.dim), sdiag(mkvc(tensor)))
|
||||
elif propType == 2: # Diagonal tensor
|
||||
@@ -302,7 +311,7 @@ def makePropertyTensor(M, tensor):
|
||||
|
||||
def invPropertyTensor(M, tensor, returnMatrix=False):
|
||||
|
||||
propType = tensorType(M, tensor)
|
||||
propType = TensorType(M, tensor)
|
||||
|
||||
if isScalar(tensor):
|
||||
T = 1./tensor
|
||||
@@ -341,3 +350,9 @@ class SimPEGLinearOperator(LinearOperator):
|
||||
def T(self):
|
||||
return self.__class__((self.shape[1],self.shape[0]),self.rmatvec,rmatvec=self.matvec,matmat=self.matmat)
|
||||
|
||||
|
||||
class DerivOperator(object):
|
||||
def __init__(self, f):
|
||||
self.f = f
|
||||
def __mul__(self, v):
|
||||
return self.f(v)
|
||||
|
||||
Reference in New Issue
Block a user