"""
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
packages or download required scripts. © mathischeap.com
"""
import numpy as np
from quadrature import Gauss
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``.
:type quad_degree: list, tuple
: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, ...
"""
def __init__(self, bf, ct, quad_degree=None):
assert bf.__class__.__name__ == 'MimeticBasisPolynomials'
assert ct.__class__.__name__ == 'CoordinateTransformation'
self.bf = bf
self.ct = ct
if quad_degree is None:
bf_N = bf.degree
quad_degree = [bf_N[i]+2 for i in range(3)]
else:
pass
self.quad_degree = quad_degree
def ___HELPER_E___(self, d_, quad_degree):
qn0, qw0 = Gauss(quad_degree[0])
qn1, qw1 = Gauss(quad_degree[1])
qn2, qw2 = Gauss(quad_degree[2])
quad_nodes, quad_weights = [qn0, qn1, qn2], [qw0, qw1, qw2]
quad_num_nodes = [len(quad_nodes_i) for quad_nodes_i in
quad_nodes]
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(
sbn1[:, np.newaxis].repeat(quad_num_nodes[0], axis=1).T,
(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,
:].repeat(quad_num_nodes[0], axis=0)
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,
:].repeat(quad_num_nodes[1], axis=0)
sn_by = (self.bf.degree[0]+1)*self.bf.degree[1]*(self.bf.degree[2]+1)
eta1 = (quad_nodes[1][:, np.newaxis].repeat(sn_by,
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,
:].repeat(quad_num_nodes[1], axis=0)
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]))[
np.newaxis, :].repeat(quad_num_nodes[2], axis=0)
eta = np.tile(np.repeat(sbn1, (self.bf.degree[0] + 1)), self.bf.degree[2])[
np.newaxis, :].repeat(quad_num_nodes[2], axis=0)
sn_bz = (self.bf.degree[0]+1)*(self.bf.degree[1]+1)*self.bf.degree[2]
sigma1 = (quad_nodes[2][:, np.newaxis].repeat(sn_bz,
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 = \
self.___HELPER_E___('x', self.quad_degree)
xi_y, eta_y, sigma_y, edge_size_d_et, quad_weights = \
self.___HELPER_E___('y', self.quad_degree)
xi_z, eta_z, sigma_z, edge_size_d_sg, quad_weights = \
self.___HELPER_E___('z', self.quad_degree)
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
p = self.quad_degree
qn0, qw0 = Gauss(p[0])
qn1, qw1 = Gauss(p[1])
qn2, qw2 = Gauss(p[2])
quad_nodes, quad_weights = [qn0, qn1, qn2], [qw0, qw1, qw2]
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_(
uvw_dy_dz, quad_weights[1], quad_weights[2], area_dy_dz
)
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_(
uvw_dz_dx, quad_weights[0], quad_weights[2], area_dz_dx
)
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_(
uvw_dx_dy, quad_weights[0], quad_weights[1], area_dx_dy
)
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.
"""
p = self.quad_degree
N = self.bf.degree
qn0, qw0 = Gauss(p[0])
qn1, qw1 = Gauss(p[1])
qn2, qw2 = Gauss(p[2])
quad_nodes, quad_weights = [qn0, qn1, qn2], [qw0, qw1, qw2]
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]
xi[m,...] = (quad_nodes[0][:,np.newaxis].repeat(p[1]+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]
et[m,...] = (quad_nodes[1][np.newaxis,:].repeat(p[0]+1,
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]
si[m,...] = (quad_nodes[2][np.newaxis,:].repeat(p[1]+1,
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',
fxyz * detJ, quad_weights[0], quad_weights[1],
quad_weights[2],
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()
