"""
Building incidence matrices, :math:`\\mathsf{E}_{(\\nabla)}`,
:math:`\\mathsf{E}_{(\\nabla\\times)}` and
:math:`\\mathsf{E}_{(\\nabla\\cdot)}` in the 3d reference domain
:math:`\\Omega_{\\mathrm{ref}}=\\left[-1,1\\right]^3`. We build them
through the local numberings which are also implemented here in this
script.
⭕ To access the source code, click on the [source] button at the right
side or click on
:download:`[incidence_matrices.py]</contents/LIBRARY/ptc/mathischeap_ptc/incidence_matrices.py>`.
Dependence may exist. In case of error, check import and install required
packages or download required scripts. © mathischeap.com
"""
import numpy as np
from scipy.sparse import csc_matrix
[docs]
def local_numbering_NP(N_xi, N_et, N_sg):
"""Generating the local numbering for a node polynomials
in :math:`\\text{NP}_{N}(\\Omega_{\\mathrm{ref}})`
constructed on three sets of nodes,
:math:`\\left\\lbrace\\xi_0,\\xi_1,\\cdots, \\xi_{N0}
\\right\\rbrace`, :math:`\\left\\lbrace\\eta_0,\\eta_1,\\cdots,
\\eta_{N1}\\right\\rbrace` and :math:`\\left\\lbrace\\varsigma_0,
\\varsigma_1,\\cdots, \\varsigma_{N2}\\right\\rbrace`.
:param N_xi: ``N_xi + 1`` nodes in set :math:`\\left\\lbrace\\xi_0,
\\xi_1,\\cdots, \\xi_{N_\\xi} \\right\\rbrace`.
:param N_et: ``N_et + 1`` nodes in set :math:`\\left\\lbrace\\eta_0,
\\eta_1,\\cdots, \\eta_{N_\\eta} \\right\\rbrace`.
:param N_sg: ``N_sg + 1`` nodes in set :math:`\\left\\lbrace
\\varsigma_0, \\varsigma_1,\\cdots, \\varsigma_{N_\\varsigma}
\\right\\rbrace`.
:type N_xi: positive integer
:type N_et: positive integer
:type N_sg: positive integer
:return: A 3d np.array contain the local numbering. Its three
dimensions refer to three axes :math:`\\left(\\xi,\\eta,
\\varsigma\\right)`.
:example:
>>> ln = local_numbering_NP(3,3,3)
>>> ln[:,:,0]
array([[ 0, 4, 8, 12],
[ 1, 5, 9, 13],
[ 2, 6, 10, 14],
[ 3, 7, 11, 15]])
>>> ln[:,:,1]
array([[16, 20, 24, 28],
[17, 21, 25, 29],
[18, 22, 26, 30],
[19, 23, 27, 31]])
>>> ln[:,:,2]
array([[32, 36, 40, 44],
[33, 37, 41, 45],
[34, 38, 42, 46],
[35, 39, 43, 47]])
>>> ln[:,:,3]
array([[48, 52, 56, 60],
[49, 53, 57, 61],
[50, 54, 58, 62],
[51, 55, 59, 63]])
"""
N = [N_xi, N_et, N_sg]
p = [N[i] + 1 for i in range(3)]
ln = np.arange((N_xi + 1) * (N_et + 1) * (N_sg + 1)).reshape(
p, order='F')
return ln
[docs]
def local_numbering_EP(N_xi, N_et, N_sg):
"""Generating the local numbering for a vector of 3D edge
polynomials in :math:`\\text{EP}_{N-1}(\\Omega_{\\mathrm{ref}})`
constructed on three sets of nodes,
:math:`\\left\\lbrace\\xi_0,\\xi_1,\\cdots, \\xi_{N0}
\\right\\rbrace`, :math:`\\left\\lbrace\\eta_0,\\eta_1,\\cdots,
\\eta_{N1}\\right\\rbrace` and :math:`\\left\\lbrace\\varsigma_0,
\\varsigma_1,\\cdots, \\varsigma_{N2}\\right\\rbrace`.
:param N_xi: ``N_xi + 1`` nodes in set :math:`\\left\\lbrace\\xi_0,
\\xi_1,\\cdots, \\xi_{N_\\xi} \\right\\rbrace`.
:param N_et: ``N_et + 1`` nodes in set :math:`\\left\\lbrace\\eta_0,
\\eta_1,\\cdots, \\eta_{N_\\eta} \\right\\rbrace`.
:param N_sg: ``N_sg + 1`` nodes in set :math:`\\left\\lbrace
\\varsigma_0, \\varsigma_1,\\cdots, \\varsigma_{N_\\varsigma}
\\right\\rbrace`.
:type N_xi: positive integer
:type N_et: positive integer
:type N_sg: positive integer
:return: A tuple of three outputs:
#. A 3d np.array contain the local numbering for the first
componement of the vector. Its three dimensions refer to
three axes :math:`\\left(\\xi,\\eta, \\varsigma\\right)`.
#. A 3d np.array contain the local numbering for the second
componement of the vector. Its three dimensions refer to
three axes :math:`\\left(\\xi,\\eta, \\varsigma\\right)`.
#. A 3d np.array contain the local numbering for the third
componement of the vector. Its three dimensions refer to
three axes :math:`\\left(\\xi,\\eta, \\varsigma\\right)`.
:example:
>>> ln0, ln1, ln2 = local_numbering_EP(2,2,2)
>>> ln0[:,:,0]
array([[0, 2, 4],
[1, 3, 5]])
>>> ln0[:,:,1]
array([[ 6, 8, 10],
[ 7, 9, 11]])
>>> ln0[:,:,2]
array([[12, 14, 16],
[13, 15, 17]])
>>> ln1[:,:,0]
array([[18, 21],
[19, 22],
[20, 23]])
>>> ln1[:,:,1]
array([[24, 27],
[25, 28],
[26, 29]])
>>> ln1[:,:,2]
array([[30, 33],
[31, 34],
[32, 35]])
>>> ln2[:,:,0]
array([[36, 39, 42],
[37, 40, 43],
[38, 41, 44]])
>>> ln2[:,:,1]
array([[45, 48, 51],
[46, 49, 52],
[47, 50, 53]])
"""
ln0 = np.arange(N_xi * (N_et + 1) * (N_sg + 1)).reshape(
[N_xi, N_et + 1, N_sg + 1], order='F')
ln1 = np.arange((N_xi + 1) * N_et * (N_sg + 1)).reshape(
[N_xi + 1, N_et, N_sg + 1], order='F') \
+ N_xi * (N_et + 1) * (N_sg + 1)
ln2 = np.arange((N_xi + 1) * (N_et + 1) * N_sg).reshape(
[N_xi + 1, N_et + 1, N_sg], order='F') \
+ N_xi * (N_et + 1) * (N_sg + 1) + \
(N_xi + 1) * N_et * (N_sg + 1)
return ln0, ln1, ln2
[docs]
def local_numbering_FP(N_xi, N_et, N_sg):
"""Generating the local numbering for a vector of 3D face
polynomials in :math:`\\text{FP}_{N-1}(\\Omega_{\\mathrm{ref}})`
constructed on three sets of nodes,
:math:`\\left\\lbrace\\xi_0,\\xi_1,\\cdots, \\xi_{N0}
\\right\\rbrace`, :math:`\\left\\lbrace\\eta_0,\\eta_1,\\cdots,
\\eta_{N1}\\right\\rbrace` and :math:`\\left\\lbrace\\varsigma_0,
\\varsigma_1,\\cdots, \\varsigma_{N2}\\right\\rbrace`.
:param N_xi: ``N_xi + 1`` nodes in set :math:`\\left\\lbrace\\xi_0,
\\xi_1,\\cdots, \\xi_{N_\\xi} \\right\\rbrace`.
:param N_et: ``N_et + 1`` nodes in set :math:`\\left\\lbrace\\eta_0,
\\eta_1,\\cdots, \\eta_{N_\\eta} \\right\\rbrace`.
:param N_sg: ``N_sg + 1`` nodes in set :math:`\\left\\lbrace
\\varsigma_0, \\varsigma_1,\\cdots, \\varsigma_{N_\\varsigma}
\\right\\rbrace`.
:type N_xi: positive integer
:type N_et: positive integer
:type N_sg: positive integer
:return: A tuple of three outputs:
#. A 3d np.array contain the local numbering for the first
componement of the vector. Its three dimensions refer to
three axes :math:`\\left(\\xi,\\eta, \\varsigma\\right)`.
#. A 3d np.array contain the local numbering for the second
componement of the vector. Its three dimensions refer to
three axes :math:`\\left(\\xi,\\eta, \\varsigma\\right)`.
#. A 3d np.array contain the local numbering for the third
componement of the vector. Its three dimensions refer to
three axes :math:`\\left(\\xi,\\eta, \\varsigma\\right)`.
:example:
>>> ln0, ln1, ln2 = local_numbering_FP(2,2,2)
>>> ln0[:,:,0]
array([[0, 3],
[1, 4],
[2, 5]])
>>> ln0[:,:,1]
array([[ 6, 9],
[ 7, 10],
[ 8, 11]])
>>> ln1[:,:,0]
array([[12, 14, 16],
[13, 15, 17]])
>>> ln1[:,:,1]
array([[18, 20, 22],
[19, 21, 23]])
>>> ln2[:,:,0]
array([[24, 26],
[25, 27]])
>>> ln2[:,:,1]
array([[28, 30],
[29, 31]])
>>> ln2[:,:,2]
array([[32, 34],
[33, 35]])
"""
ln0 = np.arange((N_xi + 1) * N_et * N_sg).reshape(
[N_xi + 1, N_et, N_sg], order='F')
ln1 = np.arange(N_xi * (N_et + 1) * N_sg).reshape(
[N_xi, N_et + 1, N_sg], order='F') +\
(N_xi + 1) * N_et * N_sg
ln2 = np.arange(N_xi * N_et * (N_sg + 1)).reshape(
[N_xi, N_et, N_sg + 1], order='F') +\
(N_xi + 1) * N_et * N_sg +\
N_xi * (N_et + 1) * N_sg
return ln0, ln1, ln2
[docs]
def local_numbering_VP(N_xi, N_et, N_sg):
"""Generating the local numbering for a volume polynomials
in :math:`\\text{VP}_{N-1}(\\Omega_{\\mathrm{ref}})`
constructed on three sets of nodes,
:math:`\\left\\lbrace\\xi_0,\\xi_1,\\cdots, \\xi_{N0}
\\right\\rbrace`, :math:`\\left\\lbrace\\eta_0,\\eta_1,\\cdots,
\\eta_{N1}\\right\\rbrace` and :math:`\\left\\lbrace\\varsigma_0,
\\varsigma_1,\\cdots, \\varsigma_{N2}\\right\\rbrace`.
:param N_xi: ``N_xi + 1`` nodes in set :math:`\\left\\lbrace\\xi_0,
\\xi_1,\\cdots, \\xi_{N_\\xi} \\right\\rbrace`.
:param N_et: ``N_et + 1`` nodes in set :math:`\\left\\lbrace\\eta_0,
\\eta_1,\\cdots, \\eta_{N_\\eta} \\right\\rbrace`.
:param N_sg: ``N_sg + 1`` nodes in set :math:`\\left\\lbrace
\\varsigma_0, \\varsigma_1,\\cdots, \\varsigma_{N_\\varsigma}
\\right\\rbrace`.
:type N_xi: positive integer
:type N_et: positive integer
:type N_sg: positive integer
:return: A 3d np.array contain the local numbering. Its three
dimensions refer to three axes :math:`\\left(\\xi,\\eta,
\\varsigma\\right)`.
:example:
>>> ln = local_numbering_VP(3,3,3)
>>> ln[:,:,0]
array([[0, 3, 6],
[1, 4, 7],
[2, 5, 8]])
>>> ln[:,:,1]
array([[ 9, 12, 15],
[10, 13, 16],
[11, 14, 17]])
>>> ln[:,:,2]
array([[18, 21, 24],
[19, 22, 25],
[20, 23, 26]])
"""
N = [N_xi, N_et, N_sg]
_ln_ = np.arange(N_xi * N_et * N_sg).reshape(N, order='F')
return _ln_
[docs]
def E_grad(N_xi, N_et, N_sg):
"""The incidence matrix :math:`\\mathsf{E}_{(\\nabla)}` for mimetic
polynomials constructed on three sets of nodes,
:math:`\\left\\lbrace\\xi_0,\\xi_1,\\cdots, \\xi_{N0}
\\right\\rbrace`, :math:`\\left\\lbrace\\eta_0,\\eta_1,\\cdots,
\\eta_{N1}\\right\\rbrace` and :math:`\\left\\lbrace\\varsigma_0,
\\varsigma_1,\\cdots, \\varsigma_{N2}\\right\\rbrace`.
:param N_xi: ``N_xi + 1`` nodes in set :math:`\\left\\lbrace\\xi_0,
\\xi_1,\\cdots, \\xi_{N_\\xi} \\right\\rbrace`.
:param N_et: ``N_et + 1`` nodes in set :math:`\\left\\lbrace\\eta_0,
\\eta_1,\\cdots, \\eta_{N_\\eta} \\right\\rbrace`.
:param N_sg: ``N_sg + 1`` nodes in set :math:`\\left\\lbrace
\\varsigma_0, \\varsigma_1,\\cdots, \\varsigma_{N_\\varsigma}
\\right\\rbrace`.
:type N_xi: positive integer
:type N_et: positive integer
:type N_sg: positive integer
:return: A csc_matrix :math:`\\mathsf{E}_{(\\nabla)}`.
:example:
>>> E = E_grad(1,1,1)
>>> E.toarray()
array([[-1, 1, 0, 0, 0, 0, 0, 0],
[ 0, 0, -1, 1, 0, 0, 0, 0],
[ 0, 0, 0, 0, -1, 1, 0, 0],
[ 0, 0, 0, 0, 0, 0, -1, 1],
[-1, 0, 1, 0, 0, 0, 0, 0],
[ 0, -1, 0, 1, 0, 0, 0, 0],
[ 0, 0, 0, 0, -1, 0, 1, 0],
[ 0, 0, 0, 0, 0, -1, 0, 1],
[-1, 0, 0, 0, 1, 0, 0, 0],
[ 0, -1, 0, 0, 0, 1, 0, 0],
[ 0, 0, -1, 0, 0, 0, 1, 0],
[ 0, 0, 0, -1, 0, 0, 0, 1]], dtype=int32)
"""
sn = local_numbering_NP(N_xi, N_et, N_sg)
dn = local_numbering_EP(N_xi, N_et, N_sg)
E = np.zeros((N_xi * (N_et + 1) * (N_sg + 1) +
(N_xi + 1) * N_et * (N_sg + 1) +
(N_xi + 1) * (N_et + 1) * N_sg,
(N_xi + 1) * (N_et + 1) * (N_sg + 1)), dtype=int)
I, J, K = np.shape(dn[0])
for k in range(K):
for j in range(J):
for i in range(I):
E[dn[0][i,j,k], sn[i,j,k]] = -1 # North
E[dn[0][i,j,k], sn[i+1,j,k]] = +1 # South
I, J, K = np.shape(dn[1])
for k in range(K):
for j in range(J):
for i in range(I):
E[dn[1][i,j,k], sn[i,j,k]] = -1 # West
E[dn[1][i,j,k], sn[i,j+1,k]] = +1 # East
I, J, K = np.shape(dn[2])
for k in range(K):
for j in range(J):
for i in range(I):
E[dn[2][i,j,k], sn[i,j,k]] = -1 # Back
E[dn[2][i,j,k], sn[i,j,k+1]] = +1 # Front
return csc_matrix(E)
[docs]
def E_curl(N_xi, N_et, N_sg):
"""The incidence matrix :math:`\\mathsf{E}_{(\\nabla\\times)}` for mimetic
polynomials constructed on three sets of nodes,
:math:`\\left\\lbrace\\xi_0,\\xi_1,\\cdots, \\xi_{N0}
\\right\\rbrace`, :math:`\\left\\lbrace\\eta_0,\\eta_1,\\cdots,
\\eta_{N1}\\right\\rbrace` and :math:`\\left\\lbrace\\varsigma_0,
\\varsigma_1,\\cdots, \\varsigma_{N2}\\right\\rbrace`.
:param N_xi: ``N_xi + 1`` nodes in set :math:`\\left\\lbrace\\xi_0,
\\xi_1,\\cdots, \\xi_{N_\\xi} \\right\\rbrace`.
:param N_et: ``N_et + 1`` nodes in set :math:`\\left\\lbrace\\eta_0,
\\eta_1,\\cdots, \\eta_{N_\\eta} \\right\\rbrace`.
:param N_sg: ``N_sg + 1`` nodes in set :math:`\\left\\lbrace
\\varsigma_0, \\varsigma_1,\\cdots, \\varsigma_{N_\\varsigma}
\\right\\rbrace`.
:type N_xi: positive integer
:type N_et: positive integer
:type N_sg: positive integer
:return: A csc_matrix :math:`\\mathsf{E}_{(\\nabla\\times)}`.
:example:
>>> E = E_curl(1,1,1)
>>> E.toarray()
array([[ 0, 0, 0, 0, 1, 0, -1, 0, -1, 0, 1, 0],
[ 0, 0, 0, 0, 0, 1, 0, -1, 0, -1, 0, 1],
[-1, 0, 1, 0, 0, 0, 0, 0, 1, -1, 0, 0],
[ 0, -1, 0, 1, 0, 0, 0, 0, 0, 0, 1, -1],
[ 1, -1, 0, 0, -1, 1, 0, 0, 0, 0, 0, 0],
[ 0, 0, 1, -1, 0, 0, -1, 1, 0, 0, 0, 0]], dtype=int32)
"""
sn = local_numbering_EP(N_xi, N_et, N_sg)
dn = local_numbering_FP(N_xi, N_et, N_sg)
E = np.zeros(((N_xi + 1) * N_et * N_sg +
N_xi * (N_et + 1) * N_sg +
N_xi * N_et * (N_sg + 1),
N_xi * (N_et + 1) * (N_sg + 1) +
(N_xi + 1) * N_et * (N_sg + 1) +
(N_xi + 1) * (N_et + 1) * N_sg), dtype=int)
I, J, K = np.shape(dn[0])
for k in range(K):
for j in range(J):
for i in range(I):
E[dn[0][i,j,k], sn[1][i,j,k ]] = +1 # Back
E[dn[0][i,j,k], sn[1][i,j,k+1]] = -1 # Front
E[dn[0][i,j,k], sn[2][i,j ,k]] = -1 # West
E[dn[0][i,j,k], sn[2][i,j+1,k]] = +1 # East
I, J, K = np.shape(dn[1])
for k in range(K):
for j in range(J):
for i in range(I):
E[dn[1][i,j,k], sn[0][i,j,k ]] = -1 # Back
E[dn[1][i,j,k], sn[0][i,j,k+1]] = +1 # Front
E[dn[1][i,j,k], sn[2][i ,j,k]] = +1 # North
E[dn[1][i,j,k], sn[2][i+1,j,k]] = -1 # South
I, J, K = np.shape(dn[2])
for k in range(K):
for j in range(J):
for i in range(I):
E[dn[2][i,j,k], sn[0][i,j ,k]] = +1 # West
E[dn[2][i,j,k], sn[0][i,j+1,k]] = -1 # East
E[dn[2][i,j,k], sn[1][i ,j,k]] = -1 # North
E[dn[2][i,j,k], sn[1][i+1,j,k]] = +1 # South
return csc_matrix(E)
[docs]
def E_div(N_xi, N_et, N_sg):
"""The incidence matrix :math:`\\mathsf{E}_{(\\nabla\\cdot)}` for mimetic
polynomials constructed on three sets of nodes,
:math:`\\left\\lbrace\\xi_0,\\xi_1,\\cdots, \\xi_{N0}
\\right\\rbrace`, :math:`\\left\\lbrace\\eta_0,\\eta_1,\\cdots,
\\eta_{N1}\\right\\rbrace` and :math:`\\left\\lbrace\\varsigma_0,
\\varsigma_1,\\cdots, \\varsigma_{N2}\\right\\rbrace`.
:param N_xi: ``N_xi + 1`` nodes in set :math:`\\left\\lbrace\\xi_0,
\\xi_1,\\cdots, \\xi_{N_\\xi} \\right\\rbrace`.
:param N_et: ``N_et + 1`` nodes in set :math:`\\left\\lbrace\\eta_0,
\\eta_1,\\cdots, \\eta_{N_\\eta} \\right\\rbrace`.
:param N_sg: ``N_sg + 1`` nodes in set :math:`\\left\\lbrace
\\varsigma_0, \\varsigma_1,\\cdots, \\varsigma_{N_\\varsigma}
\\right\\rbrace`.
:type N_xi: positive integer
:type N_et: positive integer
:type N_sg: positive integer
:return: A csc_matrix :math:`\\mathsf{E}_{(\\nabla\\cdot)}`.
:example:
>>> E = E_div(1,1,1)
>>> E.toarray()
array([[-1, 1, -1, 1, -1, 1]], dtype=int32)
"""
sn = local_numbering_FP(N_xi, N_et, N_sg)
dn = local_numbering_VP(N_xi, N_et, N_sg)
E = np.zeros((N_xi * N_et * N_sg,
(N_xi + 1) * N_et * N_sg +
N_xi * (N_et + 1) * N_sg +
N_xi * N_et * (N_sg + 1)), dtype=int)
I, J, K = np.shape(dn)
for k in range(K):
for j in range(J):
for i in range(I):
E[dn[i,j,k], sn[0][i ,j,k]] = -1 # North
E[dn[i,j,k], sn[0][i+1,j,k]] = +1 # South
E[dn[i,j,k], sn[1][i,j ,k]] = -1 # West
E[dn[i,j,k], sn[1][i,j+1,k]] = +1 # East
E[dn[i,j,k], sn[2][i,j,k ]] = -1 # Back
E[dn[i,j,k], sn[2][i,j,k+1]] = +1 # Front
return csc_matrix(E)
if __name__ == '__main__':
import doctest
doctest.testmod()
import random
N0 = random.randint(1,5)
N1 = random.randint(1,5)
N2 = random.randint(1,5)
Eg = E_grad(N0, N1, N2)
Ec = E_curl(N0, N1, N2)
Ed = E_div(N0, N1, N2)
E_curl_E_grad = Ec @ Eg
E_div_E_curl = Ed @ Ec
print(E_curl_E_grad.nnz) # must be 0 because curl(grad()) = 0
print(E_div_E_curl.nnz) # must be 0 because div(curl()) = 0
assert E_curl_E_grad.nnz == 0
assert E_div_E_curl.nnz == 0
# print(E_grad(1, 1, 1).toarray())
# print(E_curl(1, 1, 1).toarray())
print(E_div(1, 1, 1).toarray())
