"""
⭕ To access the source code, click on the [source] button at the right
side or click on
:download:`[mass_matrices_trace.py]</contents/LIBRARY/ptc/mathischeap_ptc/mass_matrices_trace.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 quadrature import Gauss
from mimetic_basis_polynomials_2d import grid2d, MimeticBasisPolynomials2D
from scipy import sparse as spspa
from coordinate_transformation_surface import extract_surface_coordinate_transformations_of
[docs]
class MassMatricesTrace(object):
""" The mass matrices of trace spaces
:math:`\\text{TN}_{N}(\\partial\\Omega)`,
:math:`\\text{TE}_{N-1}(\\partial\\Omega)` and
:math:`\\text{TF}_{N-1}(\\partial\\Omega)`.
:param bf: A :class:`MimeticBasisPolynomials`
(not :class:`MimeticBasisPolynomials2D`) instance.
:type bf: MimeticBasisPolynomials
:param ct: A CoordinateTransformation instance that represents the
mapping :math:`\\Phi:\\Omega_{\mathrm{ref}}\\to\\Omega`.
:type ct: CoordinateTransformation
:param quad_degree: (default: ``None``) The degree used for the
numerical integral. It should be a list or tuple of three
positive integers. If it is ``None``, a suitable degree will be
obtained from ``bf``.
:type quad_degree: list, tuple
:example:
>>> from coordinate_transformation import CoordinateTransformation
>>> from coordinate_transformation import Phi, d_Phi
>>> from mimetic_basis_polynomials import MimeticBasisPolynomials
>>> ct = CoordinateTransformation(Phi, d_Phi)
>>> bf = MimeticBasisPolynomials('Lobatto-4', 'Lobatto-3', 'Lobatto-2')
>>> MMT = MassMatricesTrace(bf, ct)
>>> MMT.TF # doctest: +ELLIPSIS
<52x52 sparse matrix of type '<class 'numpy.float64'>'...
"""
def __init__(self, bf, ct, quad_degree=None):
assert bf.__class__.__name__ == 'MimeticBasisPolynomials'
assert ct.__class__.__name__ == 'CoordinateTransformation'
self.bf = bf
self.ct = ct
if quad_degree is None:
bf_N = bf.degree
quad_degree = [bf_N[i]+2 for i in range(3)]
else:
pass
self.quad_degree = quad_degree
qn_0, qw_0 = Gauss(quad_degree[0])
qn_1, qw_1 = Gauss(quad_degree[1])
qn_2, qw_2 = Gauss(quad_degree[2])
self.quad_weights = {'NS': np.kron(qw_2, qw_1),
'WE': np.kron(qw_2, qw_0),
'BF': np.kron(qw_1, qw_0)}
self.quad_nodes = {'NS': [qn_1, qn_2],
'WE': [qn_0, qn_2],
'BF': [qn_0, qn_1],}
self._rho_tau_ = {'NS':grid2d(qn_1, qn_2),
'WE':grid2d(qn_0, qn_2),
'BF':grid2d(qn_0, qn_1)}
nodes = self.bf.nodes
self.bft = {'NS':MimeticBasisPolynomials2D(nodes[1], nodes[2]),
'WE':MimeticBasisPolynomials2D(nodes[0], nodes[2]),
'BF':MimeticBasisPolynomials2D(nodes[0], nodes[1])}
@property
def TN(self):
"""The mass matrix :math:`\\mathbb{M}_{\\text{N}}`.
:return: Return a csc_matrix representing the mass matrix.
"""
raise NotImplementedError("Could you code it?")
@property
def TE(self):
"""The mass matrix :math:`\\mathbb{M}_{\\text{E}}`.
:return: Return a csc_matrix representing the mass matrix.
"""
raise NotImplementedError("Could you code it?")
@property
def TF(self):
"""The mass matrix :math:`\\mathbb{M}_{\\text{F}}`.
:return: Return a csc_matrix representing the mass matrix.
"""
N, S, W, E, B, F = extract_surface_coordinate_transformations_of(self.ct)
rho, tau = self.quad_nodes['NS']
bf = self.bft['NS'].face_polynomials(rho, tau)
# North side
g = N.metric(*self._rho_tau_['NS'])
M_N = np.einsum('im, jm, m -> ij',
bf, bf, np.reciprocal(np.sqrt(g)) * self.quad_weights['NS'],
optimize='greedy')
# South side
g = S.metric(*self._rho_tau_['NS'])
M_S = np.einsum('im, jm, m -> ij',
bf, bf, np.reciprocal(np.sqrt(g)) * self.quad_weights['NS'],
optimize='greedy')
rho, tau = self.quad_nodes['WE']
bf = self.bft['WE'].face_polynomials(rho, tau)
# West side
g = W.metric(*self._rho_tau_['WE'])
M_W = np.einsum('im, jm, m -> ij',
bf, bf, np.reciprocal(np.sqrt(g)) * self.quad_weights['WE'],
optimize='greedy')
# East side
g = E.metric(*self._rho_tau_['WE'])
M_E = np.einsum('im, jm, m -> ij',
bf, bf, np.reciprocal(np.sqrt(g)) * self.quad_weights['WE'],
optimize='greedy')
rho, tau = self.quad_nodes['BF']
bf = self.bft['BF'].face_polynomials(rho, tau)
# Back side
g = B.metric(*self._rho_tau_['BF'])
M_B = np.einsum('im, jm, m -> ij',
bf, bf, np.reciprocal(np.sqrt(g)) * self.quad_weights['BF'],
optimize='greedy')
# Front side
g = F.metric(*self._rho_tau_['BF'])
M_F = np.einsum('im, jm, m -> ij',
bf, bf, np.reciprocal(np.sqrt(g)) * self.quad_weights['BF'],
optimize='greedy')
M_N = spspa.csc_matrix(M_N)
M_S = spspa.csc_matrix(M_S)
M_W = spspa.csc_matrix(M_W)
M_E = spspa.csc_matrix(M_E)
M_B = spspa.csc_matrix(M_B)
M_F = spspa.csc_matrix(M_F)
M = spspa.bmat([(M_N , None, None, None, None, None),
(None, M_S , None, None, None, None),
(None, None, M_W , None, None, None),
(None, None, None, M_E , None, None),
(None, None, None, None, M_B , None),
(None, None, None, None, None, M_F )], format='csc')
return M
if __name__ == '__main__':
import doctest
doctest.testmod()
from coordinate_transformation import Phi, d_Phi
from coordinate_transformation import CoordinateTransformation
from mimetic_basis_polynomials import MimeticBasisPolynomials
ct = CoordinateTransformation(Phi, d_Phi)
bf = MimeticBasisPolynomials('Lobatto-4', 'Lobatto-3', 'Lobatto-2')
MMT = MassMatricesTrace(bf, ct)
M2 = MMT.TF
print(M2)
