Source code for mass_matrices

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