mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-03 11:45:39 +08:00
more object oriented innerProducts
This commit is contained in:
rowanc1
+199
-240
@@ -93,12 +93,7 @@ class InnerProducts(object):
|
||||
'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)
|
||||
proj = getattr(self, '_getFaceP' + ('x'*d))()
|
||||
|
||||
elif projType == 'E':
|
||||
locs = {
|
||||
@@ -111,12 +106,7 @@ class InnerProducts(object):
|
||||
'011': [ None , None , ('eX3', 'eY2', 'eZ2')],
|
||||
'111': [ None , None , ('eX3', 'eY3', 'eZ3')]
|
||||
}
|
||||
if d == 1:
|
||||
proj = _getEdgePx(self)
|
||||
elif d == 2:
|
||||
proj = _getEdgePxx(self)
|
||||
elif d == 3:
|
||||
proj = _getEdgePxxx(self)
|
||||
proj = getattr(self, '_getEdgeP' + ('x'*d))()
|
||||
|
||||
return [V*proj(*locs[node][d-1]) for node in nodes]
|
||||
|
||||
@@ -260,283 +250,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 _getEdgePx(M):
|
||||
assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh'
|
||||
return _getEdgePx_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 _getEdgePx_Rectangular(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
|
||||
# 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
|
||||
|
||||
def _getEdgePxx_Rectangular(M):
|
||||
i, j = np.int64(range(M.nCx)), np.int64(range(M.nCy))
|
||||
posX = 0 if xFace == 'fXm' else 1
|
||||
posY = 0 if yFace == 'fYm' else 1
|
||||
posZ = 0 if zFace == 'fZm' else 1
|
||||
|
||||
iijj = ndgrid(i, j)
|
||||
ii, jj = iijj[:, 0], iijj[:, 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
|
||||
|
||||
if M._meshType == 'LRM':
|
||||
eT1 = M.r(M.tangents, 'E', 'Ex', 'M')
|
||||
eT2 = M.r(M.tangents, 'E', 'Ey', 'M')
|
||||
IND = np.r_[ind1, ind2, ind3].flatten()
|
||||
|
||||
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
|
||||
PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, M.nF)).tocsr()
|
||||
|
||||
ind1 = sub2ind(M.vnEx, np.c_[ii, jj + posX])
|
||||
ind2 = sub2ind(M.vnEy, np.c_[ii + posY, jj]) + M.nEx
|
||||
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
|
||||
|
||||
IND = np.r_[ind1, ind2].flatten()
|
||||
return PXXX
|
||||
return Pxxx
|
||||
|
||||
PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, M.nE)).tocsr()
|
||||
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
|
||||
|
||||
Reference in New Issue
Block a user