# 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

"""

import numpy as np
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.

: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)
[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],...

"""
assert bf.__class__.__name__ == 'MimeticBasisPolynomials'
assert ct.__class__.__name__ == 'CoordinateTransformation'
self.bf = bf
self.ct = ct
bf_N = bf.degree
quad_degree = [bf_N[i]+2 for i in range(3)]
else:
pass
self._xi_eta_sigma_ = grid(qn_xi, qn_et, qn_sg)

@staticmethod
M = np.einsum('m, im, jm -> ij',
optimize='greedy')
return spspa.csc_matrix(M)

@staticmethod
def _helper_F_(quad_weights, sqrtg, g, bfO, bfS):
M = np.einsum('m, im, jm -> ij',
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.
"""
detJ = self.ct.Jacobian(*self._xi_eta_sigma_)
M_N = np.einsum('im, jm, m -> ij',
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.
"""
sqrtg = self.ct.Jacobian(*self._xi_eta_sigma_)
g = self.ct.inverse_metric_matrix(*self._xi_eta_sigma_)
bf_E[0], bf_E[0])
bf_E[1], bf_E[1])
bf_E[2], bf_E[2])
bf_E[0], bf_E[1])
bf_E[0], bf_E[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.
"""
sqrtg = self.ct.Jacobian(*self._xi_eta_sigma_)
g = self.ct.inverse_metric_matrix(*self._xi_eta_sigma_)
g[1][1]*g[2][2]-g[1][2]*g[2][1], bf_F[0], bf_F[0])
g[2][2]*g[0][0]-g[2][0]*g[0][2], bf_F[1], bf_F[1])
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]
g12_20_g10_22, bf_F[0], bf_F[1])
g10_21_g11_20, bf_F[0], bf_F[2])
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.
"""