mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-27 19:48:52 +08:00
Merge pull request #75 from simpeg/innerProducts
Inner products derivative
This commit is contained in:
Rowan Cockett
@@ -620,7 +620,7 @@ class DiffOperators(object):
|
||||
# The number of cell centers in each direction
|
||||
n = self.vnC
|
||||
if(self.dim == 1):
|
||||
raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
|
||||
self._aveE2CC = speye(n[0])
|
||||
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")
|
||||
|
||||
+318
-360
@@ -10,194 +10,172 @@ 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=doFast)
|
||||
|
||||
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, doFast=doFast)
|
||||
|
||||
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
|
||||
|
||||
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')]
|
||||
}
|
||||
proj = getattr(self, '_getFaceP' + ('x'*d))()
|
||||
|
||||
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': [('eX0',), ('eX0', 'eY0'), ('eX0', 'eY0', 'eZ0')],
|
||||
'100': [('eX0',), ('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')]
|
||||
}
|
||||
proj = getattr(self, '_getEdgeP' + ('x'*d))()
|
||||
|
||||
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, prop, doFast=True, invProp=False, invMat=False):
|
||||
"""
|
||||
: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.
|
||||
:param bool invProp: inverts the material property
|
||||
:param bool invMat: inverts the matrix
|
||||
:rtype: function
|
||||
:return: dMdmu(u), the derivative of the inner product matrix (u)
|
||||
|
||||
Given u, dMdmu returns (nF, nC*nA)
|
||||
|
||||
:param np.ndarray u: vector that multiplies dMdmu
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: dMdmu, the derivative of the inner product matrix for a certain u
|
||||
"""
|
||||
return self._getInnerProductDeriv(prop, 'F', doFast=doFast, invProp=invProp, invMat=invMat)
|
||||
|
||||
|
||||
def getEdgeInnerProductDeriv(self, prop, doFast=True, invProp=False, invMat=False):
|
||||
"""
|
||||
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param bool doFast: do a faster implementation if available.
|
||||
:param bool invProp: inverts the material property
|
||||
:param bool invMat: inverts the matrix
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
|
||||
"""
|
||||
return self._getInnerProductDeriv(prop, 'E', doFast=doFast, invProp=invProp, invMat=invMat)
|
||||
|
||||
def _getInnerProductDeriv(self, prop, projType, doFast=True, invProp=False, invMat=False):
|
||||
"""
|
||||
:param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
|
||||
:param str projType: 'F' for faces 'E' for edges
|
||||
:param bool doFast: do a faster implementation if available.
|
||||
:param bool invProp: inverts the material property
|
||||
:param bool invMat: inverts the matrix
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
|
||||
"""
|
||||
fast = None
|
||||
|
||||
if hasattr(self, '_fastEdgeInnerProductDeriv') and doFast:
|
||||
fast = self._fastEdgeInnerProductDeriv(prop=prop, v=v)
|
||||
|
||||
if hasattr(self, '_fastInnerProductDeriv') and doFast:
|
||||
fast = self._fastInnerProductDeriv(projType, prop, invProp=invProp, invMat=invMat)
|
||||
if fast is not None:
|
||||
return fast
|
||||
|
||||
if P is None:
|
||||
M, P = self.getEdgeInnerProduct(prop=prop, returnP=True)
|
||||
if invProp or invMat:
|
||||
raise NotImplementedError('inverting the property or the matrix is not yet implemented for this mesh/tensorType. You should write it!')
|
||||
|
||||
return self._getInnerProductDeriv(prop, v, P, self.nE)
|
||||
tensorType = TensorType(self, prop)
|
||||
P = self._getInnerProductProjectionMatrices(projType, tensorType=tensorType)
|
||||
def innerProductDeriv(v):
|
||||
return self._getInnerProductDerivFunction(tensorType, P, projType, v)
|
||||
return innerProductDeriv
|
||||
|
||||
def _getInnerProductDeriv(self, prop, v, P, n):
|
||||
def _getInnerProductDerivFunction(self, tensorType, P, projType, 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)
|
||||
:param list P: list of projection matrices
|
||||
:param int n: nF or nE
|
||||
:param str projType: 'F' for faces 'E' for edges
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: dMdm, the derivative of the inner product matrix (n, nC*nA)
|
||||
"""
|
||||
if prop is None:
|
||||
assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
|
||||
n = getattr(self,'n'+projType)
|
||||
|
||||
if tensorType == -1:
|
||||
return None
|
||||
|
||||
if v is None:
|
||||
@@ -206,24 +184,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 +210,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 +221,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 +229,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 +240,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
|
||||
@@ -279,272 +257,252 @@ class InnerProducts(object):
|
||||
|
||||
return dMdm
|
||||
|
||||
# ------------------------ Geometries ------------------------------
|
||||
#
|
||||
#
|
||||
# node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1)
|
||||
# / /
|
||||
# / / |
|
||||
# edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k)
|
||||
# / / |
|
||||
# / / |
|
||||
# node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k)
|
||||
# | | |
|
||||
# | | node(i+1,j+1,k+1)
|
||||
# | | /
|
||||
# edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k)
|
||||
# | | /
|
||||
# | | /
|
||||
# | |/
|
||||
# node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
|
||||
# ------------------------ Geometries ------------------------------
|
||||
#
|
||||
#
|
||||
# node(i,j,k+1) ------ edge2(i,j,k+1) ----- node(i,j+1,k+1)
|
||||
# / /
|
||||
# / / |
|
||||
# edge3(i,j,k) face1(i,j,k) edge3(i,j+1,k)
|
||||
# / / |
|
||||
# / / |
|
||||
# node(i,j,k) ------ edge2(i,j,k) ----- node(i,j+1,k)
|
||||
# | | |
|
||||
# | | node(i+1,j+1,k+1)
|
||||
# | | /
|
||||
# edge1(i,j,k) face3(i,j,k) edge1(i,j+1,k)
|
||||
# | | /
|
||||
# | | /
|
||||
# | |/
|
||||
# node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k)
|
||||
|
||||
def _getFacePx(M):
|
||||
assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh'
|
||||
return _getFacePx_Rectangular(M)
|
||||
|
||||
def _getFacePxx(M):
|
||||
if M._meshType == 'TREE':
|
||||
return M._getFacePxx
|
||||
|
||||
return _getFacePxx_Rectangular(M)
|
||||
|
||||
def _getFacePxxx(M):
|
||||
if M._meshType == 'TREE':
|
||||
return M._getFacePxxx
|
||||
|
||||
return _getFacePxxx_Rectangular(M)
|
||||
|
||||
def _getEdgePxx(M):
|
||||
if M._meshType == 'TREE':
|
||||
return M._getEdgePxx
|
||||
|
||||
return _getEdgePxx_Rectangular(M)
|
||||
|
||||
def _getEdgePxxx(M):
|
||||
if M._meshType == 'TREE':
|
||||
return M._getEdgePxxx
|
||||
|
||||
return _getEdgePxxx_Rectangular(M)
|
||||
|
||||
def _getFacePx_Rectangular(M):
|
||||
"""Returns a function for creating projection matrices
|
||||
|
||||
"""
|
||||
ii = np.int64(range(M.nCx))
|
||||
|
||||
def Px(xFace):
|
||||
"""
|
||||
xFace is 'fXp' or 'fXm'
|
||||
"""
|
||||
posFx = 0 if xFace == 'fXm' else 1
|
||||
IND = ii + posFx
|
||||
PX = sp.csr_matrix((np.ones(M.nC), (range(M.nC), IND)), shape=(M.nC, M.nF))
|
||||
return PX
|
||||
|
||||
return Px
|
||||
|
||||
def _getFacePxx_Rectangular(M):
|
||||
"""returns a function for creating projection matrices
|
||||
|
||||
Mats takes you from faces a subset of all faces on only the
|
||||
faces that you ask for.
|
||||
|
||||
These are centered around a single nodes.
|
||||
|
||||
For example, if this was your entire mesh:
|
||||
|
||||
f3(Yp)
|
||||
2_______________3
|
||||
| |
|
||||
| |
|
||||
| |
|
||||
f0(Xm) | x | f1(Xp)
|
||||
| |
|
||||
| |
|
||||
|_______________|
|
||||
0 1
|
||||
f2(Ym)
|
||||
|
||||
Pxx('fXm','fYm') = | 1, 0, 0, 0 |
|
||||
| 0, 0, 1, 0 |
|
||||
|
||||
Pxx('fXp','fYm') = | 0, 1, 0, 0 |
|
||||
| 0, 0, 1, 0 |
|
||||
def _getFacePx(M):
|
||||
"""Returns a function for creating projection matrices
|
||||
|
||||
"""
|
||||
i, j = np.int64(range(M.nCx)), np.int64(range(M.nCy))
|
||||
ii = np.arange(M.nCx)
|
||||
|
||||
iijj = ndgrid(i, j)
|
||||
ii, jj = iijj[:, 0], iijj[:, 1]
|
||||
def Px(xFace):
|
||||
"""
|
||||
xFace is 'fXp' or 'fXm'
|
||||
"""
|
||||
posFx = 0 if xFace == 'fXm' else 1
|
||||
IND = ii + posFx
|
||||
PX = sp.csr_matrix((np.ones(M.nC), (range(M.nC), IND)), shape=(M.nC, M.nF))
|
||||
return PX
|
||||
|
||||
if M._meshType == 'LRM':
|
||||
fN1 = M.r(M.normals, 'F', 'Fx', 'M')
|
||||
fN2 = M.r(M.normals, 'F', 'Fy', 'M')
|
||||
return Px
|
||||
|
||||
def Pxx(xFace, yFace):
|
||||
"""
|
||||
xFace is 'fXp' or 'fXm'
|
||||
yFace is 'fYp' or 'fYm'
|
||||
"""
|
||||
# no | node | f1 | f2
|
||||
# 00 | i ,j | i , j | i, j
|
||||
# 10 | i+1,j | i+1, j | i, j
|
||||
# 01 | i ,j+1 | i , j | i, j+1
|
||||
# 11 | i+1,j+1 | i+1, j | i, j+1
|
||||
def _getFacePxx(M):
|
||||
"""returns a function for creating projection matrices
|
||||
|
||||
posFx = 0 if xFace == 'fXm' else 1
|
||||
posFy = 0 if yFace == 'fYm' else 1
|
||||
Mats takes you from faces a subset of all faces on only the
|
||||
faces that you ask for.
|
||||
|
||||
ind1 = sub2ind(M.vnFx, np.c_[ii + posFx, jj])
|
||||
ind2 = sub2ind(M.vnFy, np.c_[ii, jj + posFy]) + M.nFx
|
||||
These are centered around a single nodes.
|
||||
|
||||
IND = np.r_[ind1, ind2].flatten()
|
||||
For example, if this was your entire mesh:
|
||||
|
||||
PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF))
|
||||
f3(Yp)
|
||||
2_______________3
|
||||
| |
|
||||
| |
|
||||
| |
|
||||
f0(Xm) | x | f1(Xp)
|
||||
| |
|
||||
| |
|
||||
|_______________|
|
||||
0 1
|
||||
f2(Ym)
|
||||
|
||||
Pxx('fXm','fYm') = | 1, 0, 0, 0 |
|
||||
| 0, 0, 1, 0 |
|
||||
|
||||
Pxx('fXp','fYm') = | 0, 1, 0, 0 |
|
||||
| 0, 0, 1, 0 |
|
||||
|
||||
"""
|
||||
i, j = np.arange(M.nCx), np.arange(M.nCy)
|
||||
|
||||
iijj = ndgrid(i, j)
|
||||
ii, jj = iijj[:, 0], iijj[:, 1]
|
||||
|
||||
if M._meshType == 'LRM':
|
||||
I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]),
|
||||
getSubArray(fN2[0], [i, j + posFy]), getSubArray(fN2[1], [i, j + posFy]))
|
||||
PXX = I2x2 * PXX
|
||||
fN1 = M.r(M.normals, 'F', 'Fx', 'M')
|
||||
fN2 = M.r(M.normals, 'F', 'Fy', 'M')
|
||||
|
||||
return PXX
|
||||
def Pxx(xFace, yFace):
|
||||
"""
|
||||
xFace is 'fXp' or 'fXm'
|
||||
yFace is 'fYp' or 'fYm'
|
||||
"""
|
||||
# no | node | f1 | f2
|
||||
# 00 | i ,j | i , j | i, j
|
||||
# 10 | i+1,j | i+1, j | i, j
|
||||
# 01 | i ,j+1 | i , j | i, j+1
|
||||
# 11 | i+1,j+1 | i+1, j | i, j+1
|
||||
|
||||
return Pxx
|
||||
posFx = 0 if xFace == 'fXm' else 1
|
||||
posFy = 0 if yFace == 'fYm' else 1
|
||||
|
||||
def _getFacePxxx_Rectangular(M):
|
||||
"""returns a function for creating projection matrices
|
||||
ind1 = sub2ind(M.vnFx, np.c_[ii + posFx, jj])
|
||||
ind2 = sub2ind(M.vnFy, np.c_[ii, jj + posFy]) + M.nFx
|
||||
|
||||
Mats takes you from faces a subset of all faces on only the
|
||||
faces that you ask for.
|
||||
IND = np.r_[ind1, ind2].flatten()
|
||||
|
||||
These are centered around a single nodes.
|
||||
"""
|
||||
PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nF))
|
||||
|
||||
i, j, k = np.int64(range(M.nCx)), np.int64(range(M.nCy)), np.int64(range(M.nCz))
|
||||
if M._meshType == 'LRM':
|
||||
I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]),
|
||||
getSubArray(fN2[0], [i, j + posFy]), getSubArray(fN2[1], [i, j + posFy]))
|
||||
PXX = I2x2 * PXX
|
||||
|
||||
iijjkk = ndgrid(i, j, k)
|
||||
ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
|
||||
return PXX
|
||||
|
||||
if M._meshType == 'LRM':
|
||||
fN1 = M.r(M.normals, 'F', 'Fx', 'M')
|
||||
fN2 = M.r(M.normals, 'F', 'Fy', 'M')
|
||||
fN3 = M.r(M.normals, 'F', 'Fz', 'M')
|
||||
return Pxx
|
||||
|
||||
def Pxxx(xFace, yFace, zFace):
|
||||
"""
|
||||
xFace is 'fXp' or 'fXm'
|
||||
yFace is 'fYp' or 'fYm'
|
||||
zFace is 'fZp' or 'fZm'
|
||||
def _getFacePxxx(M):
|
||||
"""returns a function for creating projection matrices
|
||||
|
||||
Mats takes you from faces a subset of all faces on only the
|
||||
faces that you ask for.
|
||||
|
||||
These are centered around a single nodes.
|
||||
"""
|
||||
|
||||
# no | node | f1 | f2 | f3
|
||||
# 000 | i ,j ,k | i , j, k | i, j , k | i, j, k
|
||||
# 100 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k
|
||||
# 010 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k
|
||||
# 110 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k
|
||||
# 001 | i ,j ,k+1 | i , j, k | i, j , k | i, j, k+1
|
||||
# 101 | i+1,j ,k+1 | i+1, j, k | i, j , k | i, j, k+1
|
||||
# 011 | i ,j+1,k+1 | i , j, k | i, j+1, k | i, j, k+1
|
||||
# 111 | i+1,j+1,k+1 | i+1, j, k | i, j+1, k | i, j, k+1
|
||||
i, j, k = np.arange(M.nCx), np.arange(M.nCy), np.arange(M.nCz)
|
||||
|
||||
posX = 0 if xFace == 'fXm' else 1
|
||||
posY = 0 if yFace == 'fYm' else 1
|
||||
posZ = 0 if zFace == 'fZm' else 1
|
||||
|
||||
ind1 = sub2ind(M.vnFx, np.c_[ii + posX, jj, kk])
|
||||
ind2 = sub2ind(M.vnFy, np.c_[ii, jj + posY, kk]) + M.nFx
|
||||
ind3 = sub2ind(M.vnFz, np.c_[ii, jj, kk + posZ]) + M.nFx + M.nFy
|
||||
|
||||
IND = np.r_[ind1, ind2, ind3].flatten()
|
||||
|
||||
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr()
|
||||
iijjkk = ndgrid(i, j, k)
|
||||
ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
|
||||
|
||||
if M._meshType == 'LRM':
|
||||
I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]),
|
||||
getSubArray(fN2[0], [i, j + posY, k]), getSubArray(fN2[1], [i, j + posY, k]), getSubArray(fN2[2], [i, j + posY, k]),
|
||||
getSubArray(fN3[0], [i, j, k + posZ]), getSubArray(fN3[1], [i, j, k + posZ]), getSubArray(fN3[2], [i, j, k + posZ]))
|
||||
PXXX = I3x3 * PXXX
|
||||
fN1 = M.r(M.normals, 'F', 'Fx', 'M')
|
||||
fN2 = M.r(M.normals, 'F', 'Fy', 'M')
|
||||
fN3 = M.r(M.normals, 'F', 'Fz', 'M')
|
||||
|
||||
return PXXX
|
||||
return Pxxx
|
||||
def Pxxx(xFace, yFace, zFace):
|
||||
"""
|
||||
xFace is 'fXp' or 'fXm'
|
||||
yFace is 'fYp' or 'fYm'
|
||||
zFace is 'fZp' or 'fZm'
|
||||
"""
|
||||
|
||||
def _getEdgePxx_Rectangular(M):
|
||||
i, j = np.int64(range(M.nCx)), np.int64(range(M.nCy))
|
||||
# no | node | f1 | f2 | f3
|
||||
# 000 | i ,j ,k | i , j, k | i, j , k | i, j, k
|
||||
# 100 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k
|
||||
# 010 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k
|
||||
# 110 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k
|
||||
# 001 | i ,j ,k+1 | i , j, k | i, j , k | i, j, k+1
|
||||
# 101 | i+1,j ,k+1 | i+1, j, k | i, j , k | i, j, k+1
|
||||
# 011 | i ,j+1,k+1 | i , j, k | i, j+1, k | i, j, k+1
|
||||
# 111 | i+1,j+1,k+1 | i+1, j, k | i, j+1, k | i, j, k+1
|
||||
|
||||
iijj = ndgrid(i, j)
|
||||
ii, jj = iijj[:, 0], iijj[:, 1]
|
||||
posX = 0 if xFace == 'fXm' else 1
|
||||
posY = 0 if yFace == 'fYm' else 1
|
||||
posZ = 0 if zFace == 'fZm' else 1
|
||||
|
||||
if M._meshType == 'LRM':
|
||||
eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
|
||||
eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
|
||||
ind1 = sub2ind(M.vnFx, np.c_[ii + posX, jj, kk])
|
||||
ind2 = sub2ind(M.vnFy, np.c_[ii, jj + posY, kk]) + M.nFx
|
||||
ind3 = sub2ind(M.vnFz, np.c_[ii, jj, kk + posZ]) + M.nFx + M.nFy
|
||||
|
||||
def Pxx(xEdge, yEdge):
|
||||
# no | node | e1 | e2
|
||||
# 00 | i ,j | i ,j | i ,j
|
||||
# 10 | i+1,j | i ,j | i+1,j
|
||||
# 01 | i ,j+1 | i ,j+1 | i ,j
|
||||
# 11 | i+1,j+1 | i ,j+1 | i+1,j
|
||||
posX = 0 if xEdge == 'eX0' else 1
|
||||
posY = 0 if yEdge == 'eY0' else 1
|
||||
IND = np.r_[ind1, ind2, ind3].flatten()
|
||||
|
||||
ind1 = sub2ind(M.vnEx, np.c_[ii, jj + posX])
|
||||
ind2 = sub2ind(M.vnEy, np.c_[ii + posY, jj]) + M.nEx
|
||||
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr()
|
||||
|
||||
IND = np.r_[ind1, ind2].flatten()
|
||||
if M._meshType == 'LRM':
|
||||
I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]),
|
||||
getSubArray(fN2[0], [i, j + posY, k]), getSubArray(fN2[1], [i, j + posY, k]), getSubArray(fN2[2], [i, j + posY, k]),
|
||||
getSubArray(fN3[0], [i, j, k + posZ]), getSubArray(fN3[1], [i, j, k + posZ]), getSubArray(fN3[2], [i, j, k + posZ]))
|
||||
PXXX = I3x3 * PXXX
|
||||
|
||||
PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr()
|
||||
return PXXX
|
||||
return Pxxx
|
||||
|
||||
def _getEdgePx(M):
|
||||
"""Returns a function for creating projection matrices"""
|
||||
def Px(xEdge):
|
||||
assert xEdge == 'eX0', 'xEdge = %s, not eX0' % xEdge
|
||||
return sp.identity(M.nC)
|
||||
return Px
|
||||
|
||||
def _getEdgePxx(M):
|
||||
i, j = np.arange(M.nCx), np.arange(M.nCy)
|
||||
|
||||
iijj = ndgrid(i, j)
|
||||
ii, jj = iijj[:, 0], iijj[:, 1]
|
||||
|
||||
if M._meshType == 'LRM':
|
||||
I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]),
|
||||
getSubArray(eT2[0], [i + posY, j]), getSubArray(eT2[1], [i + posY, j]))
|
||||
PXX = I2x2 * PXX
|
||||
eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
|
||||
eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
|
||||
|
||||
return PXX
|
||||
return Pxx
|
||||
def Pxx(xEdge, yEdge):
|
||||
# no | node | e1 | e2
|
||||
# 00 | i ,j | i ,j | i ,j
|
||||
# 10 | i+1,j | i ,j | i+1,j
|
||||
# 01 | i ,j+1 | i ,j+1 | i ,j
|
||||
# 11 | i+1,j+1 | i ,j+1 | i+1,j
|
||||
posX = 0 if xEdge == 'eX0' else 1
|
||||
posY = 0 if yEdge == 'eY0' else 1
|
||||
|
||||
def _getEdgePxxx_Rectangular(M):
|
||||
i, j, k = np.int64(range(M.nCx)), np.int64(range(M.nCy)), np.int64(range(M.nCz))
|
||||
ind1 = sub2ind(M.vnEx, np.c_[ii, jj + posX])
|
||||
ind2 = sub2ind(M.vnEy, np.c_[ii + posY, jj]) + M.nEx
|
||||
|
||||
iijjkk = ndgrid(i, j, k)
|
||||
ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
|
||||
IND = np.r_[ind1, ind2].flatten()
|
||||
|
||||
if M._meshType == 'LRM':
|
||||
eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
|
||||
eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
|
||||
eT3 = M.r(M.tangents, 'E', 'Ez', 'M')
|
||||
PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr()
|
||||
|
||||
def Pxxx(xEdge, yEdge, zEdge):
|
||||
if M._meshType == 'LRM':
|
||||
I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]),
|
||||
getSubArray(eT2[0], [i + posY, j]), getSubArray(eT2[1], [i + posY, j]))
|
||||
PXX = I2x2 * PXX
|
||||
|
||||
# no | node | e1 | e2 | e3
|
||||
# 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k
|
||||
# 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k
|
||||
# 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k
|
||||
# 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k
|
||||
# 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k
|
||||
# 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k
|
||||
# 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k
|
||||
# 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k
|
||||
return PXX
|
||||
return Pxx
|
||||
|
||||
posX = [0,0] if xEdge == 'eX0' else [1, 0] if xEdge == 'eX1' else [0,1] if xEdge == 'eX2' else [1,1]
|
||||
posY = [0,0] if yEdge == 'eY0' else [1, 0] if yEdge == 'eY1' else [0,1] if yEdge == 'eY2' else [1,1]
|
||||
posZ = [0,0] if zEdge == 'eZ0' else [1, 0] if zEdge == 'eZ1' else [0,1] if zEdge == 'eZ2' else [1,1]
|
||||
def _getEdgePxxx(M):
|
||||
i, j, k = np.arange(M.nCx), np.arange(M.nCy), np.arange(M.nCz)
|
||||
|
||||
ind1 = sub2ind(M.vnEx, np.c_[ii, jj + posX[0], kk + posX[1]])
|
||||
ind2 = sub2ind(M.vnEy, np.c_[ii + posY[0], jj, kk + posY[1]]) + M.nEx
|
||||
ind3 = sub2ind(M.vnEz, np.c_[ii + posZ[0], jj + posZ[1], kk]) + M.nEx + M.nEy
|
||||
|
||||
IND = np.r_[ind1, ind2, ind3].flatten()
|
||||
|
||||
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr()
|
||||
iijjkk = ndgrid(i, j, k)
|
||||
ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2]
|
||||
|
||||
if M._meshType == 'LRM':
|
||||
I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]),
|
||||
getSubArray(eT2[0], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[1], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[2], [i + posY[0], j, k + posY[1]]),
|
||||
getSubArray(eT3[0], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[1], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[2], [i + posZ[0], j + posZ[1], k]))
|
||||
PXXX = I3x3 * PXXX
|
||||
eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
|
||||
eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
|
||||
eT3 = M.r(M.tangents, 'E', 'Ez', 'M')
|
||||
|
||||
return PXXX
|
||||
return Pxxx
|
||||
def Pxxx(xEdge, yEdge, zEdge):
|
||||
|
||||
# no | node | e1 | e2 | e3
|
||||
# 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k
|
||||
# 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k
|
||||
# 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k
|
||||
# 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k
|
||||
# 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k
|
||||
# 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k
|
||||
# 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k
|
||||
# 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k
|
||||
|
||||
posX = [0,0] if xEdge == 'eX0' else [1, 0] if xEdge == 'eX1' else [0,1] if xEdge == 'eX2' else [1,1]
|
||||
posY = [0,0] if yEdge == 'eY0' else [1, 0] if yEdge == 'eY1' else [0,1] if yEdge == 'eY2' else [1,1]
|
||||
posZ = [0,0] if zEdge == 'eZ0' else [1, 0] if zEdge == 'eZ1' else [0,1] if zEdge == 'eZ2' else [1,1]
|
||||
|
||||
ind1 = sub2ind(M.vnEx, np.c_[ii, jj + posX[0], kk + posX[1]])
|
||||
ind2 = sub2ind(M.vnEy, np.c_[ii + posY[0], jj, kk + posY[1]]) + M.nEx
|
||||
ind3 = sub2ind(M.vnEz, np.c_[ii + posZ[0], jj + posZ[1], kk]) + M.nEx + M.nEy
|
||||
|
||||
IND = np.r_[ind1, ind2, ind3].flatten()
|
||||
|
||||
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nE)).tocsr()
|
||||
|
||||
if M._meshType == 'LRM':
|
||||
I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]),
|
||||
getSubArray(eT2[0], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[1], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[2], [i + posY[0], j, k + posY[1]]),
|
||||
getSubArray(eT3[0], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[1], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[2], [i + posZ[0], j + posZ[1], k]))
|
||||
PXXX = I3x3 * PXXX
|
||||
|
||||
return PXXX
|
||||
return Pxxx
|
||||
|
||||
if __name__ == '__main__':
|
||||
from TensorMesh import TensorMesh
|
||||
|
||||
+45
-71
@@ -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,57 @@ class BaseTensorMesh(BaseRectangularMesh):
|
||||
else:
|
||||
return M
|
||||
|
||||
def _fastFaceInnerProductDeriv(self, prop=None, v=None):
|
||||
def _fastInnerProductDeriv(self, projType, prop, invProp=False, invMat=False):
|
||||
"""
|
||||
: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)
|
||||
:param str projType: 'E' or 'F'
|
||||
:param TensorType tensorType: type of the tensor
|
||||
:param bool invProp: inverts the material property
|
||||
:param bool invMat: inverts the matrix
|
||||
:rtype: function
|
||||
:return: dMdmu, the derivative of the inner product matrix
|
||||
"""
|
||||
return self._fastInnerProductDeriv('F', prop=prop, v=v)
|
||||
assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
|
||||
tensorType = Utils.TensorType(self, prop)
|
||||
|
||||
dMdprop = None
|
||||
|
||||
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)
|
||||
if invMat:
|
||||
MI = self._fastInnerProduct(projType, prop, invProp=invProp, invMat=invMat)
|
||||
|
||||
|
||||
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:
|
||||
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 not invMat and not invProp:
|
||||
dMdprop = self.dim * Av.T * V * ones
|
||||
elif invMat and invProp:
|
||||
dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * ones * Utils.sdiag(1./prop**2)
|
||||
|
||||
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
|
||||
if not invMat and not invProp:
|
||||
dMdprop = self.dim * Av.T * V
|
||||
elif invMat and invProp:
|
||||
dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2)
|
||||
|
||||
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
|
||||
if not invMat and not invProp:
|
||||
dMdprop = Av.T * V
|
||||
elif invMat and invProp:
|
||||
dMdprop = Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2)
|
||||
|
||||
if dMdprop is not None:
|
||||
def innerProductDeriv(v=None):
|
||||
if v is None:
|
||||
print 'Depreciation Warning: TensorMesh.innerProductDeriv. You should be supplying a vector. Use: sdiag(u)*dMdprop'
|
||||
return dMdprop
|
||||
return Utils.sdiag(v) * dMdprop
|
||||
return innerProductDeriv
|
||||
else:
|
||||
return None
|
||||
|
||||
|
||||
|
||||
|
||||
+30
-22
@@ -1022,31 +1022,39 @@ class TreeMesh(InnerProducts, BaseMesh):
|
||||
V += [1./lenj]*lenj
|
||||
return sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nE))
|
||||
|
||||
def _getFacePxx(self, xFace, yFace):
|
||||
self.number()
|
||||
xP = self._getFaceP(xFace, yFace, None)
|
||||
yP = self._getFaceP(yFace, xFace, None)
|
||||
return sp.vstack((xP, yP))
|
||||
def _getFacePxx(self):
|
||||
def Pxx(xFace, yFace):
|
||||
self.number()
|
||||
xP = self._getFaceP(xFace, yFace, None)
|
||||
yP = self._getFaceP(yFace, xFace, None)
|
||||
return sp.vstack((xP, yP))
|
||||
return Pxx
|
||||
|
||||
def _getEdgePxx(self, xEdge, yEdge):
|
||||
self.number()
|
||||
xP = self._getEdgeP(xEdge, yEdge, None)
|
||||
yP = self._getEdgeP(yEdge, xEdge, None)
|
||||
return sp.vstack((xP, yP))
|
||||
def _getEdgePxx(self):
|
||||
def Pxx(xEdge, yEdge):
|
||||
self.number()
|
||||
xP = self._getEdgeP(xEdge, yEdge, None)
|
||||
yP = self._getEdgeP(yEdge, xEdge, None)
|
||||
return sp.vstack((xP, yP))
|
||||
return Pxx
|
||||
|
||||
def _getFacePxxx(self, xFace, yFace, zFace):
|
||||
self.number()
|
||||
xP = self._getFaceP(xFace, yFace, zFace)
|
||||
yP = self._getFaceP(yFace, xFace, zFace)
|
||||
zP = self._getFaceP(zFace, xFace, yFace)
|
||||
return sp.vstack((xP, yP, zP))
|
||||
def _getFacePxxx(self):
|
||||
def Pxxx(xFace, yFace, zFace):
|
||||
self.number()
|
||||
xP = self._getFaceP(xFace, yFace, zFace)
|
||||
yP = self._getFaceP(yFace, xFace, zFace)
|
||||
zP = self._getFaceP(zFace, xFace, yFace)
|
||||
return sp.vstack((xP, yP, zP))
|
||||
return Pxxx
|
||||
|
||||
def _getEdgePxxx(self, xEdge, yEdge, zEdge):
|
||||
self.number()
|
||||
xP = self._getEdgeP(xEdge, yEdge, zEdge)
|
||||
yP = self._getEdgeP(yEdge, xEdge, zEdge)
|
||||
zP = self._getEdgeP(zEdge, xEdge, yEdge)
|
||||
return sp.vstack((xP, yP, zP))
|
||||
def _getEdgePxxx(self):
|
||||
def Pxxx(xEdge, yEdge, zEdge):
|
||||
self.number()
|
||||
xP = self._getEdgeP(xEdge, yEdge, zEdge)
|
||||
yP = self._getEdgeP(yEdge, xEdge, zEdge)
|
||||
zP = self._getEdgeP(zEdge, xEdge, yEdge)
|
||||
return sp.vstack((xP, yP, zP))
|
||||
return Pxxx
|
||||
|
||||
def plotGrid(self, ax=None, text=False, centers=False, faces=False, edges=False, lines=True, nodes=False, showIt=False):
|
||||
self.number()
|
||||
|
||||
@@ -260,56 +260,60 @@ class TestInnerProducts1D(OrderTest):
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
if __name__ == '__main__' and False:
|
||||
import sympy
|
||||
###################################################
|
||||
#### Uncomment to Reevaluate the InnerProducts ####
|
||||
###################################################
|
||||
|
||||
x,y,z = sympy.symbols(['x','y','z'])
|
||||
ex = x**2+y*z
|
||||
ey = (z**2)*x+y*z
|
||||
ez = y**2+x*z
|
||||
e = sympy.Matrix([ex,ey,ez])
|
||||
# if __name__ == '__main__':
|
||||
# import sympy
|
||||
|
||||
sigma1 = x*y+1
|
||||
sigma2 = x*z+2
|
||||
sigma3 = 3+z*y
|
||||
sigma4 = 0.1*x*y*z
|
||||
sigma5 = 0.2*x*y
|
||||
sigma6 = 0.1*z
|
||||
# x,y,z = sympy.symbols(['x','y','z'])
|
||||
# ex = x**2+y*z
|
||||
# ey = (z**2)*x+y*z
|
||||
# ez = y**2+x*z
|
||||
# e = sympy.Matrix([ex,ey,ez])
|
||||
|
||||
S1 = sympy.Matrix([[sigma1,0,0],[0,sigma1,0],[0,0,sigma1]])
|
||||
S2 = sympy.Matrix([[sigma1,0,0],[0,sigma2,0],[0,0,sigma3]])
|
||||
S3 = sympy.Matrix([[sigma1,sigma4,sigma5],[sigma4,sigma2,sigma6],[sigma5,sigma6,sigma3]])
|
||||
# sigma1 = x*y+1
|
||||
# sigma2 = x*z+2
|
||||
# sigma3 = 3+z*y
|
||||
# sigma4 = 0.1*x*y*z
|
||||
# sigma5 = 0.2*x*y
|
||||
# sigma6 = 0.1*z
|
||||
|
||||
print '3D'
|
||||
print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1)), (z,0,1))
|
||||
print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1)), (z,0,1))
|
||||
print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1)), (z,0,1))
|
||||
# S1 = sympy.Matrix([[sigma1,0,0],[0,sigma1,0],[0,0,sigma1]])
|
||||
# S2 = sympy.Matrix([[sigma1,0,0],[0,sigma2,0],[0,0,sigma3]])
|
||||
# S3 = sympy.Matrix([[sigma1,sigma4,sigma5],[sigma4,sigma2,sigma6],[sigma5,sigma6,sigma3]])
|
||||
|
||||
# print '3D'
|
||||
# print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1)), (z,0,1))
|
||||
# print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1)), (z,0,1))
|
||||
# print sympy.integrate(sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1)), (z,0,1))
|
||||
|
||||
|
||||
z = 5
|
||||
ex = x**2+y*z
|
||||
ey = (z**2)*x+y*z
|
||||
e = sympy.Matrix([ex,ey])
|
||||
# z = 5
|
||||
# ex = x**2+y*z
|
||||
# ey = (z**2)*x+y*z
|
||||
# e = sympy.Matrix([ex,ey])
|
||||
|
||||
sigma1 = x*y+1
|
||||
sigma2 = x*z+2
|
||||
sigma3 = 3+z*y
|
||||
# sigma1 = x*y+1
|
||||
# sigma2 = x*z+2
|
||||
# sigma3 = 3+z*y
|
||||
|
||||
S1 = sympy.Matrix([[sigma1,0],[0,sigma1]])
|
||||
S2 = sympy.Matrix([[sigma1,0],[0,sigma2]])
|
||||
S3 = sympy.Matrix([[sigma1,sigma3],[sigma3,sigma2]])
|
||||
# S1 = sympy.Matrix([[sigma1,0],[0,sigma1]])
|
||||
# S2 = sympy.Matrix([[sigma1,0],[0,sigma2]])
|
||||
# S3 = sympy.Matrix([[sigma1,sigma3],[sigma3,sigma2]])
|
||||
|
||||
print '2D'
|
||||
print sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1))
|
||||
print sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1))
|
||||
print sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1))
|
||||
# print '2D'
|
||||
# print sympy.integrate(sympy.integrate(e.T*S1*e, (x,0,1)), (y,0,1))
|
||||
# print sympy.integrate(sympy.integrate(e.T*S2*e, (x,0,1)), (y,0,1))
|
||||
# print sympy.integrate(sympy.integrate(e.T*S3*e, (x,0,1)), (y,0,1))
|
||||
|
||||
y = 12
|
||||
z = 5
|
||||
ex = x**2+y*z
|
||||
e = ex
|
||||
# y = 12
|
||||
# z = 5
|
||||
# ex = x**2+y*z
|
||||
# e = ex
|
||||
|
||||
sigma1 = x*y+1
|
||||
# sigma1 = x*y+1
|
||||
|
||||
print '1D'
|
||||
print sympy.integrate(e*sigma1*e, (x,0,1))
|
||||
# print '1D'
|
||||
# print sympy.integrate(e*sigma1*e, (x,0,1))
|
||||
|
||||
@@ -6,131 +6,259 @@ from TestUtils import checkDerivative
|
||||
|
||||
class TestInnerProductsDerivs(unittest.TestCase):
|
||||
|
||||
def doTestFace(self, h, rep, vec, fast):
|
||||
mesh = Mesh.TensorMesh(h)
|
||||
def doTestFace(self, h, rep, fast, meshType, invProp=False, invMat=False):
|
||||
if meshType == 'LRM':
|
||||
hRect = Utils.exampleLrmGrid(h,'rotate')
|
||||
mesh = Mesh.LogicallyRectMesh(hRect)
|
||||
elif meshType == 'Tree':
|
||||
mesh = Mesh.TreeMesh(h)
|
||||
elif meshType == 'Tensor':
|
||||
mesh = Mesh.TensorMesh(h)
|
||||
v = np.random.rand(mesh.nF)
|
||||
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)
|
||||
def fun(sig):
|
||||
M = mesh.getFaceInnerProduct(sig, invProp=invProp, invMat=invMat)
|
||||
Md = mesh.getFaceInnerProductDeriv(sig, invProp=invProp, invMat=invMat, doFast=fast)
|
||||
return M*v, Md(v)
|
||||
print meshType, 'Face', h, rep, fast, ('harmonic' if invProp and invMat else 'standard')
|
||||
return checkDerivative(fun, sig, num=5, plotIt=False)
|
||||
|
||||
def doTestEdge(self, h, rep, vec, fast):
|
||||
mesh = Mesh.TensorMesh(h)
|
||||
def doTestEdge(self, h, rep, fast, meshType, invProp=False, invMat=False):
|
||||
if meshType == 'LRM':
|
||||
hRect = Utils.exampleLrmGrid(h,'rotate')
|
||||
mesh = Mesh.LogicallyRectMesh(hRect)
|
||||
elif meshType == 'Tree':
|
||||
mesh = Mesh.TreeMesh(h)
|
||||
elif meshType == 'Tensor':
|
||||
mesh = Mesh.TensorMesh(h)
|
||||
v = np.random.rand(mesh.nE)
|
||||
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)
|
||||
def fun(sig):
|
||||
M = mesh.getEdgeInnerProduct(sig, invProp=invProp, invMat=invMat)
|
||||
Md = mesh.getEdgeInnerProductDeriv(sig, invProp=invProp, invMat=invMat, doFast=fast)
|
||||
return M*v, Md(v)
|
||||
print meshType, 'Edge', h, rep, fast, ('harmonic' if invProp and invMat else 'standard')
|
||||
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, 'Tensor'))
|
||||
def test_FaceIP_2D_float(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0,True, False))
|
||||
self.assertTrue(self.doTestFace([10, 4],0, False, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
def test_FaceIP_1D_isotropic(self):
|
||||
self.assertTrue(self.doTestFace([10],1,True, False))
|
||||
self.assertTrue(self.doTestFace([10],1, False, 'Tensor'))
|
||||
def test_FaceIP_2D_isotropic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1,True, False))
|
||||
self.assertTrue(self.doTestFace([10, 4],1, False, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
def test_FaceIP_2D_anisotropic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2,True, False))
|
||||
self.assertTrue(self.doTestFace([10, 4],2, False, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
def test_FaceIP_2D_tensor(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],3,True, False))
|
||||
self.assertTrue(self.doTestFace([10, 4],3, False, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
|
||||
def test_FaceIP_1D_float_fast(self):
|
||||
self.assertTrue(self.doTestFace([10],0, False, True))
|
||||
self.assertTrue(self.doTestFace([10],0, True, 'Tensor'))
|
||||
def test_FaceIP_2D_float_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, False, True))
|
||||
self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
def test_FaceIP_1D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10],1, False, True))
|
||||
self.assertTrue(self.doTestFace([10],1, True, 'Tensor'))
|
||||
def test_FaceIP_2D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, False, True))
|
||||
self.assertTrue(self.doTestFace([10, 4],1, True, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
def test_FaceIP_2D_anisotropic_fast(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, False, True))
|
||||
self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
|
||||
def test_EdgeIP_1D_float(self):
|
||||
self.assertTrue(self.doTestEdge([10],0, False, 'Tensor'))
|
||||
def test_EdgeIP_2D_float(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0,True, False))
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, False, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
def test_EdgeIP_1D_isotropic(self):
|
||||
self.assertTrue(self.doTestEdge([10],1, False, 'Tensor'))
|
||||
def test_EdgeIP_2D_isotropic(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1,True, False))
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, False, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
def test_EdgeIP_2D_anisotropic(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2,True, False))
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, False, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
def test_EdgeIP_2D_tensor(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],3,True, False))
|
||||
self.assertTrue(self.doTestEdge([10, 4],3, False, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
|
||||
def test_EdgeIP_1D_float_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10],0, True, 'Tensor'))
|
||||
def test_EdgeIP_2D_float_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, False, True))
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, True, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
def test_EdgeIP_1D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10],1, True, 'Tensor'))
|
||||
def test_EdgeIP_2D_isotropic_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, False, True))
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, True, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
def test_EdgeIP_2D_anisotropic_fast(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, False, True))
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, True, 'Tensor'))
|
||||
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, 'Tensor'))
|
||||
|
||||
|
||||
|
||||
def test_FaceIP_1D_float_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10],0, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_2D_float_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_3D_float_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_1D_isotropic_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10],1, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_2D_isotropic_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_3D_isotropic_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_2D_anisotropic_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor', invProp=True, invMat=True))
|
||||
def test_FaceIP_3D_anisotropic_fast_harmonic(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tensor', invProp=True, invMat=True))
|
||||
|
||||
|
||||
|
||||
def test_FaceIP_2D_float_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, False, 'LRM'))
|
||||
def test_FaceIP_3D_float_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'LRM'))
|
||||
def test_FaceIP_2D_isotropic_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, False, 'LRM'))
|
||||
def test_FaceIP_3D_isotropic_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'LRM'))
|
||||
def test_FaceIP_2D_anisotropic_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, False, 'LRM'))
|
||||
def test_FaceIP_3D_anisotropic_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'LRM'))
|
||||
def test_FaceIP_2D_tensor_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],3, False, 'LRM'))
|
||||
def test_FaceIP_3D_tensor_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'LRM'))
|
||||
|
||||
def test_FaceIP_2D_float_fast_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, True, 'LRM'))
|
||||
def test_FaceIP_3D_float_fast_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'LRM'))
|
||||
def test_FaceIP_2D_isotropic_fast_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, True, 'LRM'))
|
||||
def test_FaceIP_3D_isotropic_fast_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'LRM'))
|
||||
def test_FaceIP_2D_anisotropic_fast_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, True, 'LRM'))
|
||||
def test_FaceIP_3D_anisotropic_fast_LRM(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'LRM'))
|
||||
|
||||
def test_EdgeIP_2D_float_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, False, 'LRM'))
|
||||
def test_EdgeIP_3D_float_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'LRM'))
|
||||
def test_EdgeIP_2D_isotropic_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, False, 'LRM'))
|
||||
def test_EdgeIP_3D_isotropic_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'LRM'))
|
||||
def test_EdgeIP_2D_anisotropic_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, False, 'LRM'))
|
||||
def test_EdgeIP_3D_anisotropic_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'LRM'))
|
||||
def test_EdgeIP_2D_tensor_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],3, False, 'LRM'))
|
||||
def test_EdgeIP_3D_tensor_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'LRM'))
|
||||
|
||||
def test_EdgeIP_2D_float_fast_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, True, 'LRM'))
|
||||
def test_EdgeIP_3D_float_fast_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'LRM'))
|
||||
def test_EdgeIP_2D_isotropic_fast_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, True, 'LRM'))
|
||||
def test_EdgeIP_3D_isotropic_fast_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'LRM'))
|
||||
def test_EdgeIP_2D_anisotropic_fast_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, True, 'LRM'))
|
||||
def test_EdgeIP_3D_anisotropic_fast_LRM(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'LRM'))
|
||||
|
||||
|
||||
|
||||
|
||||
def test_FaceIP_2D_float_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, False, 'Tree'))
|
||||
def test_FaceIP_3D_float_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Tree'))
|
||||
def test_FaceIP_2D_isotropic_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, False, 'Tree'))
|
||||
def test_FaceIP_3D_isotropic_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Tree'))
|
||||
def test_FaceIP_2D_anisotropic_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, False, 'Tree'))
|
||||
def test_FaceIP_3D_anisotropic_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Tree'))
|
||||
def test_FaceIP_2D_tensor_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],3, False, 'Tree'))
|
||||
def test_FaceIP_3D_tensor_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Tree'))
|
||||
|
||||
def test_FaceIP_2D_float_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],0, True, 'Tree'))
|
||||
def test_FaceIP_3D_float_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tree'))
|
||||
def test_FaceIP_2D_isotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],1, True, 'Tree'))
|
||||
def test_FaceIP_3D_isotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tree'))
|
||||
def test_FaceIP_2D_anisotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4],2, True, 'Tree'))
|
||||
def test_FaceIP_3D_anisotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tree'))
|
||||
|
||||
def test_EdgeIP_2D_float_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, False, 'Tree'))
|
||||
def test_EdgeIP_3D_float_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'Tree'))
|
||||
def test_EdgeIP_2D_isotropic_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, False, 'Tree'))
|
||||
def test_EdgeIP_3D_isotropic_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'Tree'))
|
||||
def test_EdgeIP_2D_anisotropic_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, False, 'Tree'))
|
||||
def test_EdgeIP_3D_anisotropic_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'Tree'))
|
||||
def test_EdgeIP_2D_tensor_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],3, False, 'Tree'))
|
||||
def test_EdgeIP_3D_tensor_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Tree'))
|
||||
|
||||
def test_EdgeIP_2D_float_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],0, True, 'Tree'))
|
||||
def test_EdgeIP_3D_float_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'Tree'))
|
||||
def test_EdgeIP_2D_isotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],1, True, 'Tree'))
|
||||
def test_EdgeIP_3D_isotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'Tree'))
|
||||
def test_EdgeIP_2D_anisotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4],2, True, 'Tree'))
|
||||
def test_EdgeIP_3D_anisotropic_fast_Tree(self):
|
||||
self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Tree'))
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.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):
|
||||
|
||||
+30
-22
@@ -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
|
||||
@@ -340,4 +349,3 @@ class SimPEGLinearOperator(LinearOperator):
|
||||
@property
|
||||
def T(self):
|
||||
return self.__class__((self.shape[1],self.shape[0]),self.rmatvec,rmatvec=self.matvec,matmat=self.matmat)
|
||||
|
||||
|
||||
@@ -94,7 +94,10 @@ def closestPoints(mesh, pts, gridLoc='CC'):
|
||||
nodeInds = np.empty(pts.shape[0], dtype=int)
|
||||
|
||||
for i, pt in enumerate(pts):
|
||||
nodeInds[i] = ((np.tile(pt, (grid.shape[0],1)) - grid)**2).sum(axis=1).argmin()
|
||||
if mesh.dim == 1:
|
||||
nodeInds[i] = ((pt - grid)**2).argmin()
|
||||
else:
|
||||
nodeInds[i] = ((np.tile(pt, (grid.shape[0],1)) - grid)**2).sum(axis=1).argmin()
|
||||
|
||||
return nodeInds
|
||||
|
||||
@@ -107,7 +110,7 @@ def readUBCTensorMesh(fileName):
|
||||
|
||||
Output:
|
||||
:param SimPEG TensorMesh object
|
||||
:return
|
||||
:return
|
||||
"""
|
||||
|
||||
# Interal function to read cell size lines for the UBC mesh files.
|
||||
@@ -119,8 +122,8 @@ def readUBCTensorMesh(fileName):
|
||||
re = np.array(sp[0],dtype=int)*(' ' + sp[1])
|
||||
line = line.replace(st,re.strip())
|
||||
return np.array(line.split(),dtype=float)
|
||||
|
||||
# Read the file as line strings, remove lines with comment = !
|
||||
|
||||
# Read the file as line strings, remove lines with comment = !
|
||||
msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!')
|
||||
|
||||
# Fist line is the size of the model
|
||||
@@ -134,7 +137,7 @@ def readUBCTensorMesh(fileName):
|
||||
h3 = h3temp[::-1] # Invert the indexing of the vector to start from the bottom.
|
||||
# Adjust the reference point to the bottom south west corner
|
||||
x0[2] = x0[2] - np.sum(h3)
|
||||
# Make the mesh
|
||||
# Make the mesh
|
||||
from SimPEG import Mesh
|
||||
tensMsh = Mesh.TensorMesh([h1,h2,h3],x0)
|
||||
return tensMsh
|
||||
|
||||
Reference in New Issue
Block a user