"""
An implementation of the projection for the dual mimetic spaces. The
projection can be divided into two processes, *reduction* and
*reconstruction*.
The *reduction* process will use the given scalar or vector to obtain
the coefficients of its projection in the discrete space.
The *reconstruction* process assembles the coefficients with the correct
basis functions and leads to evaluable function.
⭕ To access the source code, click on the [source] button at the right
side or click on
:download:`[projection_dual.py]</contents/LIBRARY/ptc/mathischeap_ptc/projection_dual.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 projection import Reduction, Reconstruction
from projection_trace import ReductionTrace, ReconstructionTrace
from mass_matrices import MassMatrices
from scipy.sparse import linalg as spspalinalg
[docs]
class ReductionDual(object):
"""The reduction for spaces
:math:`\\widetilde{\\text{NP}}_{N}(\\Omega)`,
:math:`\\widetilde{\\text{EP}}_{N-1}(\\Omega)`,
:math:`\\widetilde{\\text{FP}}_{N-1}(\\Omega)` and
:math:`\\widetilde{\\text{VP}}_{N-1}(\\Omega)`.
:example:
>>> from numpy import sin, cos, pi
>>> 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-3', 'Lobatto-3', 'Lobatto-3')
>>> rdd = ReductionDual(bf, ct)
>>> def pressure(x,y,z):
... return sin(np.pi*x) * sin(pi*y) * sin(pi*z)
>>> def velocity_x(x,y,z):
... return pi * cos(pi*x) * sin(pi*y) * sin(pi*z)
>>> def velocity_y(x,y,z):
... return pi * sin(pi*x) * cos(pi*y) * sin(pi*z)
>>> def velocity_z(x,y,z):
... return pi * sin(pi*x) * sin(pi*y) * cos(pi*z)
>>> rdd.NP(pressure) # doctest: +ELLIPSIS
array([ 4.55293100e-05, 4.23149147e-04, 4.96681526e-05,...
>>> rdd.EP((velocity_x, velocity_y, velocity_z)) # doctest: +ELLIPSIS
array([-4.46767360e-03, 3.73876731e-03, 1.82174958e-03,...
>>> rdd.FP((velocity_x, velocity_y, velocity_z)) # doctest: +ELLIPSIS
array([ 0.04649943, -0.03142164, -0.07862941, 0.03226957,...
>>> rdd.VP(pressure) # doctest: +ELLIPSIS
array([ 0.14087023, 0.17787218, -0.02315247, 0.17787218,...
"""
def __init__(self, bf, ct, quad_degree=None):
self._primal_reduction_ = Reduction(bf, ct, quad_degree=quad_degree)
self._mass_matrices_ = MassMatrices(bf, ct)
[docs]
def NP(self, scalar):
"""Reduce a scalar to :math:`\\widetilde{\\text{NP}}_{N}(\\Omega)`.
:param scalar: A function in :math:`H^1(\\Omega)` to be reduced
to :math:`\\widetilde{\\text{NP}}_{N}(\\Omega)`.
:type scalar: function
:returns: A 1d np.array representing the local coefficients/dofs
of the discrete scalar.
"""
M = self._mass_matrices_.NP
dofs = self._primal_reduction_.NP(scalar)
return M @ dofs
[docs]
def EP(self, vector):
"""Reduce a vector to :math:`\\widetilde{\\text{EP}}_{N-1}(\\Omega)`.
:param vector: A vector in :math:`H(\\mathrm{curl};\\Omega)` to be reduced
to :math:`\\widetilde{\\text{EP}}_{N-1}(\\Omega)`.
:type vector: tuple, list
:returns: A 1d np.array representing the local coefficients/dofs
of the discrete vector.
"""
M = self._mass_matrices_.EP
dofs = self._primal_reduction_.EP(vector)
return M @ dofs
[docs]
def FP(self, vector):
"""Reduce a vector to :math:`\\widetilde{\\text{FP}}_{N-1}(\\Omega)`.
:param vector: A vector in :math:`H(\\mathrm{div};\\Omega)` to be reduced
to :math:`\\widetilde{\\text{FP}}_{N-1}(\\Omega)`.
:type vector: tuple, list
:returns: A 1d np.array representing the local coefficients/dofs
of the discrete vector.
"""
M = self._mass_matrices_.FP
dofs = self._primal_reduction_.FP(vector)
return M @ dofs
[docs]
def VP(self, scalar):
"""Reduce a scalar to :math:`\\widetilde{\\text{VP}}_{N-1}(\\Omega)`.
:param scalar: A function in :math:`L^2(\\Omega)` to be reduced
to :math:`\\widetilde{\\text{VP}}_{N-1}(\\Omega)`.
:type scalar: function
:returns: A 1d np.array representing the local coefficients/dofs
of the discrete scalar.
"""
M = self._mass_matrices_.VP
dofs = self._primal_reduction_.VP(scalar)
return M @ dofs
[docs]
class ReconstructionDual(object):
"""Reconstruction functions for spaces
:math:`\\widetilde{\\text{NP}}_{N}(\\Omega)`,
:math:`\\widetilde{\\text{EP}}_{N-1}(\\Omega)`,
:math:`\\widetilde{\\text{FP}}_{N-1}(\\Omega)` and
:math:`\\widetilde{\\text{VP}}_{N-1}(\\Omega)`.
:example:
>>> from numpy import sin, cos, pi
>>> 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-3', 'Lobatto-3', 'Lobatto-3')
>>> def pressure(x,y,z):
... return sin(np.pi*x) * sin(pi*y) * sin(pi*z)
>>> def velocity_x(x,y,z):
... return pi * cos(pi*x) * sin(pi*y) * sin(pi*z)
>>> def velocity_y(x,y,z):
... return pi * sin(pi*x) * cos(pi*y) * sin(pi*z)
>>> def velocity_z(x,y,z):
... return pi * sin(pi*x) * sin(pi*y) * cos(pi*z)
>>> def source(x,y,z):
... return -3 * pi**2 * sin(np.pi*x) * sin(pi*y) * sin(pi*z)
>>> rd = Reduction(bf, ct)
>>> dofs_N = rd.NP(pressure)
>>> dofs_E = rd.EP((velocity_x, velocity_y, velocity_z))
>>> dofs_F = rd.FP((velocity_x, velocity_y, velocity_z))
>>> dofs_V = rd.VP(source)
>>> rdd = ReductionDual(bf, ct)
>>> d_dofs_N = rdd.NP(pressure)
>>> d_dofs_E = rdd.EP((velocity_x, velocity_y, velocity_z))
>>> d_dofs_F = rdd.FP((velocity_x, velocity_y, velocity_z))
>>> d_dofs_V = rdd.VP(source)
>>> xi = np.linspace(-1, 1, 50)
>>> et = np.linspace(-1, 1, 50)
>>> sg = np.linspace(-1, 1, 50)
>>> rc = Reconstruction(bf, ct)
>>> xyz_p, p = rc.NP(dofs_N, xi, et, sg)
>>> xyz_w, w = rc.EP(dofs_E, xi, et, sg)
>>> xyz_u, u = rc.FP(dofs_F, xi, et, sg)
>>> xyz_f, f = rc.VP(dofs_V, xi, et, sg)
>>> rcd = ReconstructionDual(bf, ct)
>>> XYZ_p, P = rcd.NP(d_dofs_N, xi, et, sg)
>>> XYZ_w, W = rcd.EP(d_dofs_E, xi, et, sg)
>>> XYZ_u, U = rcd.FP(d_dofs_F, xi, et, sg)
>>> XYZ_f, F = rcd.VP(d_dofs_V, xi, et, sg)
>>> np.testing.assert_array_almost_equal(p, P)
>>> np.testing.assert_array_almost_equal(w[0], W[0])
>>> np.testing.assert_array_almost_equal(w[1], W[1])
>>> np.testing.assert_array_almost_equal(w[2], W[2])
>>> np.testing.assert_array_almost_equal(u[0], U[0])
>>> np.testing.assert_array_almost_equal(u[1], U[1])
>>> np.testing.assert_array_almost_equal(u[2], U[2])
>>> np.testing.assert_array_almost_equal(f, F)
>>> np.testing.assert_array_almost_equal(xyz_p[1], XYZ_f[1])
>>> np.testing.assert_array_almost_equal(XYZ_w[0], xyz_u[0])
>>> np.testing.assert_array_almost_equal(XYZ_f[2], XYZ_u[2])
"""
def __init__(self, bf, ct):
self._primal_reconstruction_ = Reconstruction(bf, ct)
self._mass_matrices_ = MassMatrices(bf, ct)
[docs]
def NP(self, loc_dofs, xi, et, sg, ravel=False):
"""Reconstruct a node polynomial evaluated at
:math:`\\Phi\\circ \\text{meshgrid}(\\xi, \\eta, \\varsigma)``.
:param loc_dofs: A 1d np.array containing the coefficients of
the discrete scalar in :math:`\\widetilde{\\text{NP}}_{N}(\\Omega)`.
:type loc_dofs: np.array
:param xi: :math:`\\xi`.
:param et: :math:`\\eta`.
:param sg: :math:`\\varsigma`. The reconstructed scalar will be
evaluated at :math:`\\Phi\\circ \\text{meshgrid}(
\\xi, \\eta, \\varsigma)``.
:type xi: 1d np.array
:type et: 1d np.array
:type sg: 1d np.array
:param ravel: (default: ``False``) If ``ravel`` is ``True``, we
will flat the outputs (as a 1d array) according to local
numbering. Otherwise, you get 3d outputs corresponding to
the indexing.
:returns: A tuple of two outputs:
* :math:`(x,y,z)`: The reconstructed scalar is evaluated at
:math:`(x,y,z):=\\Phi\\circ \\text{meshgrid}(
\\xi, \\eta, \\varsigma)`.
* values: The values of the reconstructed scalar evaluated
at :math:`(x,y,z)`.
"""
M = self._mass_matrices_.NP
invM = spspalinalg.inv(M)
primal_dofs = invM @ loc_dofs
return self._primal_reconstruction_.NP(
primal_dofs, xi, et, sg, ravel=ravel)
[docs]
def EP(self, loc_dofs, xi, et, sg, ravel=False):
"""Reconstruct a vector of edge polynomials evaluated at
:math:`\\Phi\\circ \\text{meshgrid}(\\xi, \\eta, \\varsigma)``.
:param loc_dofs: A 1d np.array containing the coefficients of
the discrete vector in :math:`\\widetilde{\\text{EP}}_{N-1}(\\Omega)`.
:type loc_dofs: np.array
:param xi: :math:`\\xi`.
:param et: :math:`\\eta`.
:param sg: :math:`\\varsigma`. The reconstructed vector will be
evaluated at :math:`\\Phi\\circ \\text{meshgrid}(
\\xi, \\eta, \\varsigma)``.
:type xi: 1d np.array
:type et: 1d np.array
:type sg: 1d np.array
:param ravel: (default: ``False``) If ``ravel`` is ``True``, we
will flat the outputs (as a 1d array) according to local
numbering. Otherwise, you get 3d outputs corresponding to
the indexing.
:returns: A tuple of two outputs:
* :math:`(x,y,z)`: The reconstructed scalar is evaluated at
:math:`(x,y,z):=\\Phi\\circ \\text{meshgrid}(
\\xi, \\eta, \\varsigma)`.
* :math:`(u,v,w)`: Three components of the reconstructed
vector evaluated at :math:`(x,y,z)`.
"""
M = self._mass_matrices_.EP
invM = spspalinalg.inv(M)
primal_dofs = invM @ loc_dofs
return self._primal_reconstruction_.EP(
primal_dofs, xi, et, sg, ravel=ravel)
[docs]
def FP(self, loc_dofs, xi, et, sg, ravel=False):
"""Reconstruct a vector of face polynomials evaluated at
:math:`\\Phi\\circ \\text{meshgrid}(\\xi, \\eta, \\varsigma)``.
:param loc_dofs: A 1d np.array containing the coefficients of
the discrete vector in :math:`\\widetilde{\\text{FP}}_{N-1}(\\Omega)`.
:type loc_dofs: np.array
:param xi: :math:`\\xi`.
:param et: :math:`\\eta`.
:param sg: :math:`\\varsigma`. The reconstructed vector will be
evaluated at :math:`\\Phi\\circ \\text{meshgrid}(
\\xi, \\eta, \\varsigma)``.
:type xi: 1d np.array
:type et: 1d np.array
:type sg: 1d np.array
:param ravel: (default: ``False``) If ``ravel`` is ``True``, we
will flat the outputs (as a 1d array) according to local
numbering. Otherwise, you get 3d outputs corresponding to
the indexing.
:returns: A tuple of two outputs:
* :math:`(x,y,z)`: The reconstructed scalar is evaluated at
:math:`(x,y,z):=\\Phi\\circ \\text{meshgrid}(
\\xi, \\eta, \\varsigma)`.
* :math:`(u,v,w)`: Three components of the reconstructed
vector evaluated at :math:`(x,y,z)`.
"""
M = self._mass_matrices_.FP
invM = spspalinalg.inv(M)
primal_dofs = invM @ loc_dofs
return self._primal_reconstruction_.FP(
primal_dofs, xi, et, sg, ravel=ravel)
[docs]
def VP(self, loc_dofs, xi, et, sg, ravel=False):
"""Reconstruct a volume polynomial evaluated at
:math:`\\Phi\\circ \\text{meshgrid}(\\xi, \\eta, \\varsigma)``.
:param loc_dofs: A 1d np.array containing the coefficients of
the discrete scala in :math:`\\widetilde{\\text{VP}}_{N-1}(\\Omega)`.
:type loc_dofs: np.array
:param xi: :math:`\\xi`.
:param et: :math:`\\eta`.
:param sg: :math:`\\varsigma`. The reconstructed scalar will be
evaluated at :math:`\\Phi\\circ \\text{meshgrid}(
\\xi, \\eta, \\varsigma)``.
:type xi: 1d np.array
:type et: 1d np.array
:type sg: 1d np.array
:param ravel: (default: ``False``) If ``ravel`` is ``True``, we
will flat the outputs (as a 1d array) according to local
numbering. Otherwise, you get 3d outputs corresponding to
the indexing.
:returns: A tuple of two outputs:
* :math:`(x,y,z)`: The reconstructed scalar is evaluated at
:math:`(x,y,z):=\\Phi\\circ \\text{meshgrid}(
\\xi, \\eta, \\varsigma)`.
* values: The values of the reconstructed scalar evaluated
at :math:`(x,y,z)`.
"""
M = self._mass_matrices_.VP
invM = spspalinalg.inv(M)
primal_dofs = invM @ loc_dofs
if 'float' in primal_dofs.__class__.__name__ :
primal_dofs = np.array([primal_dofs,])
return self._primal_reconstruction_.VP(
primal_dofs, xi, et, sg, ravel=ravel)
if __name__ == '__main__':
import doctest
doctest.testmod()
