# Source code for incidence_matrices

"""
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

"""

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]
"""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.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)

Ec = E_curl(N0, N1, N2)
Ed = E_div(N0, N1, N2)

E_div_E_curl = Ed @ Ec

print(E_div_E_curl.nnz) # must be 0 because div(curl()) = 0

assert E_div_E_curl.nnz == 0