# Source code for projection

"""
An implementation of the projection. 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.py]</contents/LIBRARY/ptc/mathischeap_ptc/projection.py>.
Dependence may exist. In case of error, check import and install required

"""

import numpy as np
from mimetic_basis_polynomials import grid

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

:param bf: The basis functions in :math:\\Omega_{\\mathrm{ref}}.
:type bf: MimeticBasisPolynomials
:param ct: The coordinate transformation representing the
mapping :math:\\Phi,

.. 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 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')
>>> rd = Reduction(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)
>>> rd.NP(pressure) # doctest: +ELLIPSIS
array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00,...
>>> rd.EP((velocity_x, velocity_y, velocity_z)) # doctest: +ELLIPSIS
array([ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,...
>>> rd.FP((velocity_x, velocity_y, velocity_z)) # doctest: +ELLIPSIS
array([ 0.03986372, -0.00413826,  0.03344094, -0.03986372,...
>>> rd.VP(pressure) # doctest: +ELLIPSIS
array([0.00445833, 0.00454444, 0.00079032, 0.00454444, 0.01872708,
0.00615386, 0.00079032, 0.00615386, 0.00079032, ...

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

sbn0 = self.bf.nodes[0]
sbn1 = self.bf.nodes[1]
sbn2 = self.bf.nodes[2]
if d_ == 'x':
a = sbn0[1:] - sbn0[:-1]
a = a.ravel('F')
b = (self.bf.degree[1] + 1) * (self.bf.degree[2] + 1)
edge_size_x = np.tile(a, b)
snb_x = b * self.bf.degree[0]
D = quad_nodes[0][:, np.newaxis].repeat(snb_x, axis=1) + 1
assert np.shape(D)[1] == len(edge_size_x)
xi1 = D * edge_size_x / 2
xi2 = np.tile(sbn0[:-1], b)
xi = xi1 + xi2
eta = np.tile(np.tile(
(self.bf.degree[0], 1)).reshape(
(quad_num_nodes[0], self.bf.degree[0] * (self.bf.degree[1] + 1)),
order='F'), (1, self.bf.degree[2] + 1))
sigma = sbn2.repeat(self.bf.degree[0] * (self.bf.degree[1] + 1))[np.newaxis,
return xi, eta, sigma, edge_size_x, quad_weights
elif d_ == 'y':
edge_size_y = np.tile(np.repeat((sbn1[1:] - sbn1[:-1]),
self.bf.degree[0] + 1),
self.bf.degree[2] + 1)
xi = np.tile(sbn0, self.bf.degree[1] * (self.bf.degree[2] + 1))[np.newaxis,
sn_by = (self.bf.degree[0]+1)*self.bf.degree[1]*(self.bf.degree[2]+1)
axis=1) + 1) * edge_size_y / 2
eta2 = np.tile(np.repeat(sbn1[:-1], (self.bf.degree[0] + 1)),
(self.bf.degree[2] + 1))
eta = eta1 + eta2
sigma = sbn2.repeat(self.bf.degree[1] * (self.bf.degree[0] + 1))[np.newaxis,
return xi, eta, sigma, edge_size_y, quad_weights
elif d_ == 'z':
edge_size_z = np.repeat((sbn2[1:] - sbn2[:-1]),
self.bf.degree[0] + 1).repeat(self.bf.degree[1] + 1)
xi = np.tile(sbn0, (self.bf.degree[1] + 1) * (self.bf.degree[2]))[
eta = np.tile(np.repeat(sbn1, (self.bf.degree[0] + 1)), self.bf.degree[2])[
sn_bz = (self.bf.degree[0]+1)*(self.bf.degree[1]+1)*self.bf.degree[2]
axis=1) + 1) * edge_size_z / 2
sigma2 = sbn2[:-1].repeat((self.bf.degree[0] + 1) * (self.bf.degree[1] + 1))
sigma = sigma1 + sigma2
return xi, eta, sigma, edge_size_z, quad_weights
else:
raise Exception()

@staticmethod
def _HELP_F_(uvw, qw1, qw2, area):
return np.einsum('jkl, kl, j -> j',
uvw, np.tensordot(qw1, qw2, axes=0),
area * 0.25, optimize='optimal'
)

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

:param scalar: A function in :math:H^1(\\Omega) to be reduced
to :math:\\text{NP}_{N}(\\Omega).
:type scalar: function
:returns: A 1d np.array representing the local coefficients/dofs
of the discrete scalar.
"""
nodes = self.bf.nodes
nodes = np.meshgrid(*nodes, indexing='ij')
nodes = [nodes[i].ravel('F') for i in range(3)]
xyz = self.ct.mapping(*nodes)
coefficients = scalar(*xyz)
return coefficients

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

:param vector: A list or tuple of three functions which
represents the vector in :math:H(\\mathrm{curl};\\Omega)
to be reduced to :math:\\text{FP}_{N-1}(\\Omega).
:type vector: list, tuple
:returns: A 1d np.array representing the local coefficients/dofs
of the discrete vector.
"""
xi_x, eta_x, sigma_x, edge_size_d_xi, quad_weights = \
xi_y, eta_y, sigma_y, edge_size_d_et, quad_weights = \
xi_z, eta_z, sigma_z, edge_size_d_sg, quad_weights = \
edge_size = (edge_size_d_xi, edge_size_d_et, edge_size_d_sg)

smctm = self.ct.mapping(xi_x, eta_x, sigma_x)
J = self.ct.Jacobian_matrix(xi_x, eta_x, sigma_x)
J = (J[0][0], J[1][0], J[2][0])
u = vector[0](*smctm)
v = vector[1](*smctm)
w = vector[2](*smctm)
local_dx = np.einsum(
'jk, j, k -> k', J[0]*u + J[1]*v + J[2]*w, quad_weights[0],
edge_size[0] * 0.5, optimize='greedy')

smctm = self.ct.mapping(xi_y, eta_y, sigma_y)
J = self.ct.Jacobian_matrix(xi_y, eta_y, sigma_y)
J = (J[0][1], J[1][1], J[2][1])
u = vector[0](*smctm)
v = vector[1](*smctm)
w = vector[2](*smctm)
local_dy = np.einsum(
'jk, j, k -> k', J[0]*u + J[1]*v + J[2]*w, quad_weights[1],
edge_size[1]*0.5, optimize='greedy')

smctm = self.ct.mapping(xi_z, eta_z, sigma_z)
J = self.ct.Jacobian_matrix(xi_z, eta_z, sigma_z)
J = (J[0][2], J[1][2], J[2][2])
u = vector[0](*smctm)
v = vector[1](*smctm)
w = vector[2](*smctm)
local_dz = np.einsum(
'jk, j, k -> k', J[0]*u + J[1]*v + J[2]*w, quad_weights[2],
edge_size[2]*0.5, optimize='greedy')

return np.hstack((local_dx, local_dy, local_dz))

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

:param vector: A list or tuple of three functions which
represents the vector in :math:H(\\mathrm{div};\\Omega)
to be reduced to :math:\\text{FP}_{N-1}(\\Omega).
:type vector: list, tuple
:returns: A 1d np.array representing the local coefficients/dofs
of the discrete vector.
"""
N = self.bf.degree
qn0, qw0 = Gauss(p[0])
qn1, qw1 = Gauss(p[1])
qn2, qw2 = Gauss(p[2])

xi = np.zeros(((N[0]+1)*N[1]*N[2], p[1] + 1, p[2] + 1))
et = np.zeros(((N[0]+1)*N[1]*N[2], p[1] + 1, p[2] + 1))
si = np.zeros(((N[0]+1)*N[1]*N[2], p[1] + 1, p[2] + 1))
area_dy_dz = np.zeros((N[0] + 1) * N[1] * N[2])
for k in range(self.bf.degree[2]):
for j in range(self.bf.degree[1]):
for i in range(self.bf.degree[0] + 1):
m = i + j * (self.bf.degree[0] + 1) + k * (self.bf.degree[0] + 1) * self.bf.degree[1]
xi[m, ...] = np.ones((p[1] + 1, p[2] + 1)) * self.bf.nodes[0][i]
et[m, ...] = (quad_nodes[1][:, np.newaxis].repeat(p[2] + 1, axis=1) + 1) * (
self.bf.nodes[1][j + 1] - self.bf.nodes[1][j]) / 2 + \
self.bf.nodes[1][j]
si[m, ...] = (quad_nodes[2][np.newaxis, :].repeat((p[1] + 1), axis=0) + 1) * (
self.bf.nodes[2][k + 1] - self.bf.nodes[2][k]) / 2 + \
self.bf.nodes[2][k]
area_dy_dz[m] = (self.bf.nodes[2][k + 1] - self.bf.nodes[2][k]) \
* (self.bf.nodes[1][j + 1] - self.bf.nodes[1][j])
xi_x, et_x, si_x = xi, et, si

xi = np.zeros((N[0]*(N[1]+1)*N[2], p[0] + 1, p[2] + 1))
et = np.zeros((N[0]*(N[1]+1)*N[2], p[0] + 1, p[2] + 1))
si = np.zeros((N[0]*(N[1]+1)*N[2], p[0] + 1, p[2] + 1))
area_dz_dx = np.zeros(N[0] * (N[1] + 1) * N[2])
for k in range(self.bf.degree[2]):
for j in range(self.bf.degree[1] + 1):
for i in range(self.bf.degree[0]):
m = i + j * self.bf.degree[0] + k * (self.bf.degree[1] + 1) * self.bf.degree[0]
xi[m, ...] = (quad_nodes[0][:, np.newaxis].repeat(p[2] + 1, axis=1) + 1) * (
self.bf.nodes[0][i + 1] - self.bf.nodes[0][i]) / 2 + \
self.bf.nodes[0][i]
et[m, ...] = np.ones((p[0] + 1, p[2] + 1)) * self.bf.nodes[1][j]
si[m, ...] = (quad_nodes[2][np.newaxis, :].repeat(p[0] + 1, axis=0) + 1) * (
self.bf.nodes[2][k + 1] - self.bf.nodes[2][k]) / 2 + \
self.bf.nodes[2][k]
area_dz_dx[m] = (self.bf.nodes[2][k + 1] - self.bf.nodes[2][k]) \
* (self.bf.nodes[0][i + 1] - self.bf.nodes[0][i])
xi_y, et_y, si_y = xi, et, si

xi = np.zeros((N[0]*N[1]*(N[2]+1), p[0] + 1, p[1] + 1))
et = np.zeros((N[0]*N[1]*(N[2]+1), p[0] + 1, p[1] + 1))
si = np.zeros((N[0]*N[1]*(N[2]+1), p[0] + 1, p[1] + 1))
area_dx_dy = np.zeros(N[0] * N[1] * (N[2] + 1))
for k in range(self.bf.degree[2] + 1):
for j in range(self.bf.degree[1]):
for i in range(self.bf.degree[0]):
m = i + j * self.bf.degree[0] + k * self.bf.degree[1] * self.bf.degree[0]
xi[m, ...] = (quad_nodes[0][:, np.newaxis].repeat(p[1] + 1, axis=1) + 1) * (
self.bf.nodes[0][i + 1] - self.bf.nodes[0][i]) / 2 + \
self.bf.nodes[0][i]
et[m, ...] = (quad_nodes[1][np.newaxis, :].repeat(p[0] + 1, axis=0) + 1) * (
self.bf.nodes[1][j + 1] - self.bf.nodes[1][j]) / 2 + \
self.bf.nodes[1][j]
si[m, ...] = np.ones((p[0] + 1, p[1] + 1)) * self.bf.nodes[2][k]
area_dx_dy[m] = (self.bf.nodes[1][j + 1] - self.bf.nodes[1][j]) \
* (self.bf.nodes[0][i + 1] - self.bf.nodes[0][i])
xi_z, et_z, si_z = xi, et, si

J = self.ct.Jacobian_matrix(xi_x, et_x, si_x)
Jx_0 = J[1][1] * J[2][2] - J[1][2] * J[2][1]
Jx_1 = J[2][1] * J[0][2] - J[2][2] * J[0][1]
Jx_2 = J[0][1] * J[1][2] - J[0][2] * J[1][1]
J = self.ct.Jacobian_matrix(xi_y, et_y, si_y)
Jy_0 = J[1][2] * J[2][0] - J[1][0] * J[2][2]
Jy_1 = J[2][2] * J[0][0] - J[2][0] * J[0][2]
Jy_2 = J[0][2] * J[1][0] - J[0][0] * J[1][2]
J = self.ct.Jacobian_matrix(xi_z, et_z, si_z)
Jz_0 = J[1][0] * J[2][1] - J[1][1] * J[2][0]
Jz_1 = J[2][0] * J[0][1] - J[2][1] * J[0][0]
Jz_2 = J[0][0] * J[1][1] - J[0][1] * J[1][0]

smctm_x = self.ct.mapping(xi_x, et_x, si_x)
smctm_y = self.ct.mapping(xi_y, et_y, si_y)
smctm_z = self.ct.mapping(xi_z, et_z, si_z)

u = vector[0](*smctm_x)
v = vector[1](*smctm_x)
w = vector[2](*smctm_x)
uvw_dy_dz = Jx_0 * u + Jx_1 * v + Jx_2 * w
local_dy_dz = self._HELP_F_(
)

u = vector[0](*smctm_y)
v = vector[1](*smctm_y)
w = vector[2](*smctm_y)
uvw_dz_dx = Jy_0 * u + Jy_1 * v + Jy_2 * w
local_dz_dx = self._HELP_F_(
)

u = vector[0](*smctm_z)
v = vector[1](*smctm_z)
w = vector[2](*smctm_z)
uvw_dx_dy = Jz_0 * u + Jz_1 * v + Jz_2 * w
local_dx_dy = self._HELP_F_(
)

return np.hstack((local_dy_dz, local_dz_dx, local_dx_dy))

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

:param scalar: A function in :math:L^2(\\Omega) to be reduced
to :math:\\text{VP}_{N-1}(\\Omega).
:type scalar: function
:returns: A 1d np.array representing the local coefficients/dofs
of the discrete scalar.
"""
N = self.bf.degree
qn0, qw0 = Gauss(p[0])
qn1, qw1 = Gauss(p[1])
qn2, qw2 = Gauss(p[2])
magic_factor = 0.125
xi = np.zeros((N[0]*N[1]*N[2], p[0]+1, p[1]+1, p[2]+1))
et = np.zeros((N[0]*N[1]*N[2], p[0]+1, p[1]+1, p[2]+1))
si = np.zeros((N[0]*N[1]*N[2], p[0]+1, p[1]+1, p[2]+1))
volume = np.zeros(N[0]*N[1]*N[2])
for k in range(self.bf.degree[2]):
for j in range(self.bf.degree[1]):
for i in range(self.bf.degree[0]):
m = i + j*self.bf.degree[0] + k*self.bf.degree[0]*self.bf.degree[1]
axis=1)[:,:,np.newaxis].repeat(p[2]+1, axis=2) + 1)\
* (self.bf.nodes[0][i+1]-self.bf.nodes[0][i]
)/2 + self.bf.nodes[0][i]
axis=0)[:,:,np.newaxis].repeat(p[2]+1, axis=2) + 1)\
* (self.bf.nodes[1][j+1]-self.bf.nodes[1][j]
)/2 + self.bf.nodes[1][j]
axis=0)[np.newaxis,:,:].repeat(p[0]+1, axis=0) + 1)\
* (self.bf.nodes[2][k+1]-self.bf.nodes[2][k]
)/2 + self.bf.nodes[2][k]
volume[m] = (self.bf.nodes[0][i+1]-self.bf.nodes[0][i]) \
* (self.bf.nodes[1][j+1]-self.bf.nodes[1][j]) \
* (self.bf.nodes[2][k+1]-self.bf.nodes[2][k]) * magic_factor

detJ = self.ct.Jacobian(xi, et, si)
xyz = self.ct.mapping(xi, et, si)
fxyz = scalar(*xyz)
return np.einsum('jklm, k, l, m, j -> j',
volume, optimize='greedy'
)

[docs]
class Reconstruction(object):
"""A wrapper of reconstruction functions in :math:\\Omega.

:param bf: The basis functions in :math:\\Omega_{\\mathrm{ref}}.
:type bf: MimeticBasisPolynomials
:param ct: The coordinate transformation representing the
mapping :math:\\Phi,

.. math::

\\Phi: \\Omega_{\\mathrm{ref}}\\to\\Omega

:type ct: CoordinateTransformation

:example:

>>> import numpy as np
>>> 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-5', 'Lobatto-5', 'Lobatto-5')
>>> 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)
>>> loc_dofs_N = rd.NP(pressure)
>>> loc_dofs_E = rd.EP((velocity_x, velocity_y, velocity_z))
>>> loc_dofs_F = rd.FP((velocity_x, velocity_y, velocity_z))
>>> loc_dofs_V = rd.VP(source)
>>> rc = Reconstruction(bf, ct)
>>> xi = np.linspace(-1, 1, 50)
>>> et = np.linspace(-1, 1, 50)
>>> sg = np.linspace(-1, 1, 50)
>>> xyz_p, p = rc.NP(loc_dofs_N, xi, et, sg)
>>> xyz_w, w = rc.EP(loc_dofs_E, xi, et, sg)
>>> xyz_u, u = rc.FP(loc_dofs_F, xi, et, sg)
>>> xyz_f, f = rc.VP(loc_dofs_V, xi, et, sg)
>>> p - pressure(*xyz_p) # the error # doctest: +ELLIPSIS
array([[[ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, ...,
0.00000000e+00,  0.00000000e+00,  0.00000000e+00],...
>>> w[0] - velocity_x(*xyz_w) # the error # doctest: +ELLIPSIS
array([[[ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, ...,
0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
[ 0.00000000e+00,  2.29166833e-02,  2.98219834e-02, ...,
3.36914832e-03, -5.18668625e-03,  6.18969815e-18],...
>>> u[0] - velocity_x(*xyz_u) # the error # doctest: +ELLIPSIS
array([[[ 1.34796048e-05,  4.22664447e-04,  8.32751318e-04, ...,
8.32751318e-04,  4.22664447e-04,  1.34796048e-05],
[ 4.22664447e-04,  3.56786323e-04,  3.72185374e-04, ...,
3.72185374e-04,  3.56786323e-04,  4.22664447e-04],...
>>> f - source(*xyz_f) # the error # doctest: +ELLIPSIS
array([[[-0.77424798, -1.49814167, -2.09148356, ..., -1.60496992,
-3.0320961 , -4.59503016],
[-1.49814167, -1.29548405, -1.11149995, ..., -0.90717392,
-1.6210195 , -2.40674137],...

"""
def __init__(self, bf, ct):
assert bf.__class__.__name__ == 'MimeticBasisPolynomials'
assert ct.__class__.__name__ == 'CoordinateTransformation'
self.bf = bf
self.ct = 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:\\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).

"""
xi_etg_sigma = grid(xi, et, sg)
xyz = self.ct.mapping(*xi_etg_sigma)
basis = self.bf.node_polynomials(xi, et, sg)
v = np.einsum('ij, i -> j', basis, loc_dofs, optimize='optimal')
if ravel:
pass
else:
shape = [len(xi), len(et), len(sg)]
xyz = [xyz[j].reshape(shape, order='F') for j in range(3)]
v = v.reshape(shape, order='F')
return xyz, v

[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:\\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).
"""
N0, N1, N2 = self.bf.degree
xi_etg_sigma = grid(xi, et, sg)
xyz = self.ct.mapping(*xi_etg_sigma)
basis = self.bf.edge_polynomials(xi, et, sg)
num_dofs_x = N0 * (N1+1) * (N2+1)
num_dofs_y = (N0+1) * N1 * (N2+1)
loc_dofs_x = loc_dofs[ : num_dofs_x]
loc_dofs_y = loc_dofs[num_dofs_x : num_dofs_x+num_dofs_y]
loc_dofs_z = loc_dofs[num_dofs_x+num_dofs_y : ]
iJ = self.ct.inverse_Jacobian_matrix(*xi_etg_sigma)
u = np.einsum('ij, i -> j', basis[0], loc_dofs_x, optimize='greedy')
v = np.einsum('ij, i -> j', basis[1], loc_dofs_y, optimize='greedy')
w = np.einsum('ij, i -> j', basis[2], loc_dofs_z, optimize='greedy')
value = [None, None, None]
for j in range(3):
value[j] = u*iJ[0][j] + v*iJ[1][j] + w*iJ[2][j]
if ravel:
pass
else:
shape = [len(xi), len(et), len(sg)]
xyz = [xyz[j].reshape(shape, order='F') for j in range(3)]
# noinspection PyUnresolvedReferences
value = [value[j].reshape(shape, order='F') for j in range(3)]
return xyz, value

[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:\\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).
"""
N0, N1, N2 = self.bf.degree
xi_etg_sigma = grid(xi, et, sg)
xyz = self.ct.mapping(*xi_etg_sigma)
basis = self.bf.face_polynomials(xi, et, sg)
num_dofs_x = (N0+1) * N1 * N2
num_dofs_y = N0 * (N1+1) * N2
loc_dofs_x = loc_dofs[ : num_dofs_x]
loc_dofs_y = loc_dofs[num_dofs_x : num_dofs_x+num_dofs_y]
loc_dofs_z = loc_dofs[num_dofs_x+num_dofs_y : ]
iJ = self.ct.inverse_Jacobian_matrix(*xi_etg_sigma)
u = np.einsum('ij, i -> j', basis[0], loc_dofs_x, optimize='greedy')
v = np.einsum('ij, i -> j', basis[1], loc_dofs_y, optimize='greedy')
w = np.einsum('ij, i -> j', basis[2], loc_dofs_z, optimize='greedy')
_0u = iJ[1][1] * iJ[2][2] - iJ[1][2] * iJ[2][1]
_0v = iJ[2][1] * iJ[0][2] - iJ[2][2] * iJ[0][1]
_0w = iJ[0][1] * iJ[1][2] - iJ[0][2] * iJ[1][1]
_1u = iJ[1][2] * iJ[2][0] - iJ[1][0] * iJ[2][2]
_1v = iJ[2][2] * iJ[0][0] - iJ[2][0] * iJ[0][2]
_1w = iJ[0][2] * iJ[1][0] - iJ[0][0] * iJ[1][2]
_2u = iJ[1][0] * iJ[2][1] - iJ[1][1] * iJ[2][0]
_2v = iJ[2][0] * iJ[0][1] - iJ[2][1] * iJ[0][0]
_2w = iJ[0][0] * iJ[1][1] - iJ[0][1] * iJ[1][0]
value_x = u * _0u +  v * _0v +  w * _0w
value_y = u * _1u +  v * _1v +  w * _1w
value_z = u * _2u +  v * _2v +  w * _2w
if ravel:
pass
else:
shape = [len(xi), len(et), len(sg)]
xyz = [xyz[j].reshape(shape, order='F') for j in range(3)]
value_x = value_x.reshape(shape, order='F')
value_y = value_y.reshape(shape, order='F')
value_z = value_z.reshape(shape, order='F')
return xyz, (value_x, value_y, value_z)

[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:\\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).
"""
xi_etg_sigma = grid(xi, et, sg)
xyz = self.ct.mapping(*xi_etg_sigma)
basis = self.bf.volume_polynomials(xi, et, sg)
det_iJ = self.ct.inverse_Jacobian(*xi_etg_sigma)
v = np.einsum('ij, i -> j', basis, loc_dofs,
optimize='greedy') * det_iJ
if ravel:
pass
else:
shape = [len(xi), len(et), len(sg)]
xyz = [xyz[j].reshape(shape, order='F') for j in range(3)]
v = v.reshape(shape, order='F')
return xyz, v

if __name__ == '__main__':
import doctest
doctest.testmod()