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