diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 1a1de3ed..6ffb20b9 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -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