"""
Given a :class:`CoordinateTransformation` instance which stands for a
mapping :math:`\\Phi:\\Omega_{\mathrm{ref}}\\to\\Omega`, we compute the
mass matrices, :math:`\\mathsf{M}_{\\text{N}}`,
:math:`\\mathsf{M}_{\\text{E}}`, :math:`\\mathsf{M}_{\\text{F}}` and
:math:`\\mathsf{M}_{\\text{V}}`, of spaces
:math:`\\text{NP}_{N}(\\Omega)`, :math:`\\text{EP}_{N-1}(\\Omega)`,
:math:`\\text{FP}_{N-1}(\\Omega)` and :math:`\\text{VP}_{N-1}(\\Omega)`.
⭕ To access the source code, click on the [source] button at the right
side or click on
:download:`[mass_matrices.py]</contents/LIBRARY/ptc/mathischeap_ptc/mass_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 quadrature import Gauss
from mimetic_basis_polynomials import grid
from scipy import sparse as spspa
[docs]
class MassMatrices(object):
""" The mass matrices.
:param bf:
: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-2', 'Lobatto-2', 'Lobatto-2')
>>> mm = MassMatrices(bf, ct)
>>> mm.quad_degree
[4, 4, 4]
>>> MN = mm.NP
>>> np.shape(MN)
(27, 27)
>>> ME = mm.EP
>>> np.shape(ME)
(54, 54)
>>> MF = mm.FP
>>> np.shape(MF)
(36, 36)
>>> MV = mm.VP
>>> np.shape(MV)
(8, 8)
>>> MN[13, :].toarray() # doctest: +ELLIPSIS
array([[0.00069202, 0.00092145, 0.00016439, 0.00092145, 0.01896296,
0.00381929, 0.00016439, 0.00381929, 0.00016439, ...
>>> ME[25, :].toarray() # doctest: +ELLIPSIS
array([[ 2.84993722e-04, 2.18286425e-02, 2.31867356e-02,
8.84765301e-03, 7.36065685e-04, -6.34982123e-03,...
>>> MF[30,:].toarray() # doctest: +ELLIPSIS
array([[ 9.57306626e-03, 1.40151619e-02, -3.83134271e-03,
-1.56689038e-01, 5.83237762e-02, 7.05662159e-02,...
>>> MV.toarray() # doctest: +ELLIPSIS
array([[12.32813585, -0.99033139, -0.99033139, 0.29781084, -0.99033139,
0.29781084, 0.29781084, -0.04159915],
[-0.99033139, 16.33440736, 0.29781084, -3.27422609, 0.29781084,
-3.27422609, -0.04159915, 0.29781084],...
"""
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_xi, qw_xi = Gauss(quad_degree[0])
qn_et, qw_et = Gauss(quad_degree[1])
qn_sg, qw_sg = Gauss(quad_degree[2])
self.quad_weights = np.kron(qw_sg, np.kron(qw_et, qw_xi))
self.quad_nodes = [qn_xi, qn_et, qn_sg]
self._xi_eta_sigma_ = grid(qn_xi, qn_et, qn_sg)
@staticmethod
def _helper_E_(quad_weights, sqrt_g_g, bfO, bfS):
M = np.einsum('m, im, jm -> ij',
quad_weights*sqrt_g_g, bfO, bfS,
optimize='greedy')
return spspa.csc_matrix(M)
@staticmethod
def _helper_F_(quad_weights, sqrtg, g, bfO, bfS):
M = np.einsum('m, im, jm -> ij',
quad_weights*sqrtg*g, bfO, bfS,
optimize='greedy')
return spspa.csc_matrix(M)
@property
def NP(self):
"""The mass matrix :math:`\\mathsf{M}_{\\text{N}}`.
:return: Return a csc_matrix representing the mass matrix.
"""
bf_N = self.bf.node_polynomials(*self.quad_nodes)
detJ = self.ct.Jacobian(*self._xi_eta_sigma_)
M_N = np.einsum('im, jm, m -> ij',
bf_N, bf_N, detJ*self.quad_weights,
optimize='greedy')
return spspa.csc_matrix(M_N)
@property
def EP(self):
"""The mass matrix :math:`\\mathsf{M}_{\\text{E}}`.
:return: Return a csc_matrix representing the mass matrix.
"""
bf_E = self.bf.edge_polynomials(*self.quad_nodes)
sqrtg = self.ct.Jacobian(*self._xi_eta_sigma_)
g = self.ct.inverse_metric_matrix(*self._xi_eta_sigma_)
M00 = self._helper_E_(self.quad_weights, sqrtg*g[0][0],
bf_E[0], bf_E[0])
M11 = self._helper_E_(self.quad_weights, sqrtg*g[1][1],
bf_E[1], bf_E[1])
M22 = self._helper_E_(self.quad_weights, sqrtg*g[2][2],
bf_E[2], bf_E[2])
M01 = self._helper_E_(self.quad_weights, sqrtg*g[0][1],
bf_E[0], bf_E[1])
M02 = self._helper_E_(self.quad_weights, sqrtg*g[0][2],
bf_E[0], bf_E[2])
M12 = self._helper_E_(self.quad_weights, sqrtg*g[1][2],
bf_E[1], bf_E[2])
M10 = M01.T
M20 = M02.T
M21 = M12.T
M_E = spspa.bmat([(M00, M01, M02),
(M10, M11, M12),
(M20, M21, M22)], format='csc')
return M_E
@property
def FP(self):
"""The mass matrix :math:`\\mathsf{M}_{\\text{F}}`.
:return: Return a csc_matrix representing the mass matrix.
"""
bf_F = self.bf.face_polynomials(*self.quad_nodes)
sqrtg = self.ct.Jacobian(*self._xi_eta_sigma_)
g = self.ct.inverse_metric_matrix(*self._xi_eta_sigma_)
M00 = self._helper_F_(self.quad_weights, sqrtg,
g[1][1]*g[2][2]-g[1][2]*g[2][1], bf_F[0], bf_F[0])
M11 = self._helper_F_(self.quad_weights, sqrtg,
g[2][2]*g[0][0]-g[2][0]*g[0][2], bf_F[1], bf_F[1])
M22 = self._helper_F_(self.quad_weights, sqrtg,
g[0][0]*g[1][1]-g[0][1]*g[1][0], bf_F[2], bf_F[2])
g12_20_g10_22 = g[1][2] * g[2][0] - g[1][0] * g[2][2]
g10_21_g11_20 = g[1][0] * g[2][1] - g[1][1] * g[2][0]
g20_01_g21_00 = g[2][0] * g[0][1] - g[2][1] * g[0][0]
M01 = self._helper_F_(self.quad_weights, sqrtg,
g12_20_g10_22, bf_F[0], bf_F[1])
M02 = self._helper_F_(self.quad_weights, sqrtg,
g10_21_g11_20, bf_F[0], bf_F[2])
M12 = self._helper_F_(self.quad_weights, sqrtg,
g20_01_g21_00, bf_F[1], bf_F[2])
M10 = M01.T
M20 = M02.T
M21 = M12.T
M_F = spspa.bmat([(M00, M01, M02),
(M10, M11, M12),
(M20, M21, M22)], format='csc')
return M_F
@property
def VP(self):
"""The mass matrix :math:`\\mathsf{M}_{\\text{V}}`.
:return: Return a csc_matrix representing the mass matrix.
"""
bf_V = self.bf.volume_polynomials(*self.quad_nodes)
detJ = self.ct.Jacobian(*self._xi_eta_sigma_)
M_V = np.einsum('im, jm, m -> ij',
bf_V, bf_V, np.reciprocal(detJ)*self.quad_weights,
optimize='greedy')
return spspa.csc_matrix(M_V)
if __name__ == '__main__':
import doctest
doctest.testmod()
