# Source code for projection_dual

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

"""

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,...

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