# Source code for projection_trace

"""
An implementation of the projection for the trace spaces on a 3D surface
:math:\\Gamma. 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 first assembles the coefficients with the
correct basis functions and evaluates the discrete variable at given
points.

⭕ To access the source code, click on the [source] button at the right
side or click on
:download:[projection_trace.py]</contents/LIBRARY/ptc/mathischeap_ptc/projection_trace.py>.
Dependence may exist. In case of error, check import and install required

"""

import numpy as np
from mimetic_basis_polynomials_2d import grid2d

[docs]
class ReductionTrace(object):
"""A wrapper of reduction functions.

:param bf: The basis functions in :math:\\Pi_{\\mathrm{ref}}.
:type bf: MimeticBasisPolynomials2D
:param ct: The coordinate transformation representing the
mapping :math:\\Psi,

.. math::

\\Psi: \\Pi_{\\mathrm{ref}}\\to\\Gamma

:type ct: CoordinateTransformationSurface
:param quad_degree: (default: None) The degree used for the
numerical integral. It should be a list or tuple of two
positive integers. If it is None, a suitable degree will be
obtained from bf.

:example:

>>> def p(x,y,z): return np.sin(np.pi*x) * np.sin(np.pi*y) * np.sin(np.pi*z)
>>> from mimetic_basis_polynomials_2d import MimeticBasisPolynomials2D
>>> from coordinate_transformation_surface import CoordinateTransformationSurface
>>> from coordinate_transformation_surface import Psi, d_Psi
>>> bf = MimeticBasisPolynomials2D('Lobatto-3', 'Lobatto-3')
>>> ct = CoordinateTransformationSurface(Psi, d_Psi)
>>> Rd = ReductionTrace(bf, ct)
>>> Rd.TF(p)
array([0.03873144, 0.02717954, 0.02752846, 0.02717954, 0.04952738,
0.04752636, 0.02752846, 0.04752636, 0.01838736])

"""
assert bf.__class__.__name__ == 'MimeticBasisPolynomials2D'
assert ct.__class__.__name__ == 'CoordinateTransformationSurface'
self.bf = bf
self.ct = ct
bf_N = bf.degree
quad_degree = [bf_N[i]+2 for i in range(2)]
else:
pass

[docs]
def TN(self, scalar):
"""Reduce a scalar on :math:\\Gamma to
:math:\\text{TN}_{N}(\\Gamma).

:param scalar:
:return: A 1d np.array representing the local coefficients/dofs
of the discrete scalar.

"""
raise NotImplementedError("Could you code it?")

[docs]
def TE(self, vector):
"""Reduce a vector on :math:\\Gamma to
:math:\\text{TE}_{N-1}(\\Gamma).

:param vector:
:return: A 1d np.array representing the local coefficients/dofs
of the discrete vector.

"""
raise NotImplementedError("Could you code it?")

[docs]
def TF(self, scalar):
"""Reduce a scalar on :math:\\Gamma to
:math:\\text{TF}_{N-1}(\\Gamma).

:param scalar:
:return: A 1d np.array representing the local coefficients/dofs
of the discrete scalar.

"""
N = self.bf.degree
NUM_basis = N[0]*N[1]

nodes = self.bf.nodes

qn0, qw0 = Gauss(p[0])
qn1, qw1 = Gauss(p[1])
magic_factor = 0.25

rho = np.zeros((NUM_basis, p[0] + 1, p[1] + 1))
tau = np.zeros((NUM_basis, p[0] + 1, p[1] + 1))
volume = np.zeros(NUM_basis)
for j in range(N[1]):
for i in range(N[0]):
m = i + j*N[0]
rho[m,...] = (quad_nodes[0][:,np.newaxis].repeat(p[1]+1, axis=1) + 1)\
* (nodes[0][i+1]-nodes[0][i]
)/2 + nodes[0][i]
tau[m,...] = (quad_nodes[1][np.newaxis,:].repeat(p[0]+1, axis=0) + 1)\
* (nodes[1][j+1]-nodes[1][j]
)/2 + nodes[1][j]
volume[m] = (nodes[0][i+1]-nodes[0][i]) \
* (nodes[1][j+1]-nodes[1][j])  * magic_factor

g = self.ct.metric(rho, tau)
xyz = self.ct.mapping(rho, tau)
fxyz = scalar(*xyz)
return np.einsum('jkl, k, l, j -> j',
volume, optimize='greedy'
)

[docs]
class ReconstructionTrace(object):
"""A wrapper of reconstruction functions.

:param bf: The basis functions in :math:\\Pi_{\\mathrm{ref}}.
:type bf: MimeticBasisPolynomials2D
:param ct: The coordinate transformation representing the
mapping :math:\\Psi,

.. math::

\\Psi: \\Pi_{\\mathrm{ref}}\\to\\Gamma

:type ct: CoordinateTransformationSurface

:example:

>>> def p(x,y,z): return np.sin(np.pi*x) * np.sin(np.pi*y) * np.sin(np.pi*z)
>>> from mimetic_basis_polynomials_2d import MimeticBasisPolynomials2D
>>> from coordinate_transformation_surface import CoordinateTransformationSurface
>>> from coordinate_transformation_surface import Psi, d_Psi
>>> bf = MimeticBasisPolynomials2D('Lobatto-20', 'Lobatto-20')
>>> ct = CoordinateTransformationSurface(Psi, d_Psi)
>>> Rd = ReductionTrace(bf, ct)
>>> dofs = Rd.TF(p)
>>> Rc = ReconstructionTrace(bf, ct)
>>> rho = np.linspace(-1,1,5)
>>> tau = np.linspace(-1,1,5)
>>> xyz, v = Rc.TF(dofs, rho, tau)
>>> np.max(np.abs(p(*xyz) - v)) # the max error from the exact function # doctest: +ELLIPSIS
0.0002989...

"""
def __init__(self, bf, ct):
assert bf.__class__.__name__ == 'MimeticBasisPolynomials2D'
assert ct.__class__.__name__ == 'CoordinateTransformationSurface'
self.bf = bf
self.ct = ct

[docs]
def TN(self, loc_dofs, rho, tau, ravel=False):
"""Reconstruct a discrete trace polynomial in
:math:\\text{TN}_{N}(\\Gamma)
evaluated at
:math:\\Psi\\circ \\text{grid2d}(\\varrho, \\tau).

:param loc_dofs: A 1d np.array containing the coefficients of
the discrete variable in :math:\\text{TN}_{N}(\\Gamma).
:type loc_dofs: np.array
:param rho: :math:\\varrho.
:param tau: :math:\\tau. The reconstruction will be
evaluated at :math:\\Psi\\circ \\text{grid2d}(
\\varrho, \\tau).
:type rho: 1d np.array
:type tau: 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 2d outputs corresponding to
the indexing.
:returns: A tuple of two outputs:

* :math:(x,y,z): The reconstruction is evaluated at
:math:(x,y,z):=\\Psi\\circ \\text{grid2d}(\\varrho, \\tau).
* values: The values of the reconstructed variable at :math:(x,y,z).

"""
raise NotImplementedError("Could you code it?")

[docs]
def TE(self, loc_dofs, rho, tau, ravel=False):
"""Reconstruct a discrete trace polynomial in
:math:\\text{TE}_{N-1}(\\Gamma)
evaluated at
:math:\\Psi\\circ \\text{grid2d}(\\varrho, \\tau).

:param loc_dofs: A 1d np.array containing the coefficients of
the discrete variable in :math:\\text{TE}_{N-1}(\\Gamma).
:type loc_dofs: np.array
:param rho: :math:\\varrho.
:param tau: :math:\\tau. The reconstruction will be
evaluated at :math:\\Psi\\circ \\text{grid2d}(
\\varrho, \\tau).
:type rho: 1d np.array
:type tau: 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 2d outputs corresponding to
the indexing.
:returns: A tuple of two outputs:

* :math:(x,y,z): The reconstruction is evaluated at
:math:(x,y,z):=\\Psi\\circ \\text{grid2d}(\\varrho, \\tau).
* values: The values of the reconstructed variable at :math:(x,y,z).

"""
raise NotImplementedError("Could you code it?")

[docs]
def TF(self, loc_dofs, rho, tau, ravel=False):
"""Reconstruct a discrete trace polynomial in
:math:\\text{TF}_{N-1}(\\Gamma)
evaluated at
:math:\\Psi\\circ \\text{grid2d}(\\varrho, \\tau).

:param loc_dofs: A 1d np.array containing the coefficients of
the discrete variable in :math:\\text{TF}_{N-1}(\\Gamma).
:type loc_dofs: np.array
:param rho: :math:\\varrho.
:param tau: :math:\\tau. The reconstruction will be
evaluated at :math:\\Psi\\circ \\text{grid2d}(
\\varrho, \\tau).
:type rho: 1d np.array
:type tau: 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 2d outputs corresponding to
the indexing.
:returns: A tuple of two outputs:

* :math:(x,y,z): The reconstruction is evaluated at
:math:(x,y,z):=\\Psi\\circ \\text{grid2d}(\\varrho, \\tau).
* values: The values of the reconstructed variable at :math:(x,y,z).

"""
basis = self.bf.face_polynomials(rho, tau)
rho_tau = grid2d(rho, tau)
xyz = self.ct.mapping(*rho_tau)
g = self.ct.metric(*rho_tau)
v = np.einsum('ij, i -> j', basis, loc_dofs,
optimize='greedy') * np.reciprocal(np.sqrt(g))

if ravel:
pass
else:
shape = [len(rho), len(tau)]
xyz = [xyz[j].reshape(shape, order='F') for j in range(3)]
v = v.reshape(shape, order='F')

return xyz, v

if __name__ == '__main__':
def p(x,y,z): return np.sin(np.pi*x) + 0 * y * z # * np.sin(np.pi*y) * np.sin(np.pi*z)

from mimetic_basis_polynomials_2d import MimeticBasisPolynomials2D
from coordinate_transformation_surface import CoordinateTransformationSurface
from coordinate_transformation_surface import Psi, d_Psi

bf = MimeticBasisPolynomials2D('Lobatto-10', 'Lobatto-10')
ct = CoordinateTransformationSurface(Psi, d_Psi)

Rd = ReductionTrace(bf, ct)

dofs = Rd.TF(p)

Rc = ReconstructionTrace(bf, ct)
rho = np.linspace(-1,1,100)
tau = np.linspace(-1,1,100)

xyz, v = Rc.TF(dofs, rho, tau)

import matplotlib.pyplot as plt
from matplotlib import cm

fig = plt.figure(figsize=(8, 6))

V = (v + 1) / 2
cmap = cm.bwr

ax.plot_surface(*xyz, facecolors=cmap(V))

mappable = cm.ScalarMappable(cmap=cmap)
mappable.set_array(V)
cbar = plt.colorbar(mappable, ax=ax, ticks=[0, 0.5, 1], shrink=1, aspect=20,
extend='both',
orientation='vertical', label='Some Units')

cbar.ax.set_yticklabels(
[r"-1", "0", "1"])  # vertically oriented colorbar

ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$z$')
plt.show()

# print(v[:,0]- v[:,1])