dubiner.icc
Go to the documentation of this file.
#ifndef _RHEOLEF_DUBINER_ICC
#define _RHEOLEF_DUBINER_ICC
//
// evaluate the Proriol-Dubiner basis on a point of the reference element
//
// author: Pierre.Saramito@imag.fr
//
// date: 5 september 2017
//
#include "rheolef/quadrature.h"
#include "basis_option.h"
namespace rheolef {
// ----------------------------------------------
// coefficients in the recurrence formula of the Jacobi polynomials:
// see War-2010, page 5:4
// ----------------------------------------------
namespace details {
namespace legendre {
template<class T>
static inline T a (size_t i) {
return T(int(i))/sqrt(T(2*int(i)+1)*T(2*int(i)-1));
}
} // namespace legendre
namespace jacobi {
template<class T>
static inline T a (T alpha, T beta, size_t degree) {
return (2*degree + 1 + alpha + beta)*(2*degree + 2 + alpha + beta)
/(2*(degree + 1)*(degree + 1 + alpha + beta));
}
template<class T>
static inline T b (T alpha, T beta, size_t degree) {
return (sqr(alpha) - sqr(beta))*(2*degree + 1 + alpha + beta)
/(2*(degree + 1)*(2*degree + alpha + beta)*(degree + 1 + alpha + beta));
}
template<class T>
static inline T c (T alpha, T beta, size_t degree) {
return (degree + alpha)*(degree + beta)*(2*degree + 2 + alpha + beta)
/((degree + 1)*(degree + 1 + alpha + beta)*(2*degree + alpha + beta));
}
} // namespace jacobi
} // namespace details
// ----------------------------------------------
// Legendre/Dubiner basis evaluation:
// - Dubiner basis [Kir-2010] when K=t,T
// - Legendre basis when K=e,q,H
// - product basis when K=P
// ----------------------------------------------
// ortho-normalized Legendre polynoms on a segment: see e.g. HesWar-2008, p. 45
// Note: HesWar-2008 refers to tilde_e=[-1,1] instead of hat_e=[0,1]
template<class T, class Container>
static
void
eval_dubiner_basis_e (
const T& tilde_x, // in [-1,1]
size_t degree,
const std::vector<size_t>& perm,
Container& value,
std::vector<point_basic<size_t> >& power_index)
{
typedef reference_element::size_type size_type;
typedef typename float_traits<T>::type float_type;
value.resize (degree+1);
power_index.resize (value.size());
size_t loc_idof;
loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(0));
power_index[perm[loc_idof]] = point_basic<size_t>(0);
value [perm[loc_idof]] = 1/sqrt(float_type(2));
if (degree == 0) return;
loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(1));
power_index[perm[loc_idof]] = point_basic<size_t>(1);
value [perm[loc_idof]] = sqrt(float_type(3)/float_type(2))*tilde_x;
for (size_type i = 1; i+1 <= degree; i++) {
namespace L = rheolef::details::legendre;
size_type loc_idof_i = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
size_type loc_idof_im1 = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i-1));
loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i+1));
power_index[perm[loc_idof]] = point_basic<size_t>(i+1);
value [perm[loc_idof]] = (tilde_x*value[perm[loc_idof_i]] - L::a<float_type>(i)*value[perm[loc_idof_im1]])/L::a<float_type>(i+1);
}
}
// orthogonal Dubiner polynoms on a triangle: see Kir-2010, p. 5:7
// Notes: polynoms are not normalized => mass matrix is diag but not identity
// Kir-2010 refers to [-1,1]^2 instead of hat_K subset [0,1]^2
template<class T, class Container>
static
void
eval_dubiner_basis_t (
const point_basic<T>& tilde_x, // in [-1,1]^2
size_t degree,
const std::vector<size_t>& perm,
Container& value,
std::vector<point_basic<size_t> >& power_index)
{
typedef reference_element::size_type size_type;
typedef typename float_traits<T>::type float_type;
value.resize ((degree+1)*(degree+2)/2);
power_index.resize (value.size());
// Dubiner, by algorithm Kir-2010, page 5:7
T xi0 = tilde_x[0],
xi1 = tilde_x[1];
size_t loc_idof;
loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(0,0));
power_index[perm[loc_idof]] = point_basic<size_t>(0,0);
value [perm[loc_idof]] = 1;
if (degree == 0) return;
loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(1,0));
power_index[perm[loc_idof]] = point_basic<size_t>(1,0);
value [perm[loc_idof]] = (1 + 2*xi0 + xi1)/float_type(2);
for (size_type i = 1; i < degree; i++) {
size_t loc_idof_i_0 = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i ,0));
size_t loc_idof_im1_0 = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i-1,0));
loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i+1,0));
power_index[perm[loc_idof]] = point_basic<size_t>(i+1,0);
value [perm[loc_idof]] = float_type(2*int(i)+1)/float_type(int(i)+1)*(1 + 2*xi0 + xi1)/float_type(2)*value[perm[loc_idof_i_0]]
- float_type( int(i) )/float_type(int(i)+1)*sqr((1 - xi1)/float_type(2)) *value[perm[loc_idof_im1_0]];
}
for (size_type i = 0; i < degree; i++) {
size_t loc_idof_i_0 = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,0));
loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,1));
power_index[perm[loc_idof]] = point_basic<size_t>(i,1);
value [perm[loc_idof]] = value[perm[loc_idof_i_0]]*(1 + 2*i + (3+2*i)*xi1)/float_type(2);
}
for (size_type i = 0; i < degree; i++) {
for (size_type j = 1; i+j < degree; j++) {
namespace J = rheolef::details::jacobi;
size_t loc_idof_i_j = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j ));
size_t loc_idof_i_jm1 = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j-1));
loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j+1));
power_index[perm[loc_idof]] = point_basic<size_t>(i,j+1);
value [perm[loc_idof]] = (J::a<float_type>(2*int(i)+1,0,j)*xi1 + J::b<float_type>(2*int(i)+1,0,j))*value[perm[loc_idof_i_j]]
- J::c<float_type>(2*int(i)+1,0,j) *value[perm[loc_idof_i_jm1]];
}
}
}
template<class T, class Container>
static
void
eval_dubiner_basis_T (
const point_basic<T>& tilde_x, // in [-1,1]^3
size_t degree,
const std::vector<size_t>& perm,
Container& value,
std::vector<point_basic<size_t> >& power_index)
{
typedef reference_element::size_type size_type;
typedef typename float_traits<T>::type float_type;
value.resize (reference_element::n_node(reference_element::T,degree));
power_index.resize (value.size());
// Dubiner, by algorithm Kir-2010, page 5:7
T xi0 = tilde_x[0],
xi1 = tilde_x[1],
xi2 = tilde_x[2];
T F1 = (2 + 2*xi0 + xi1 + xi2)/2,
F2 = sqr((xi1 + xi2)/2),
F3 = (2 + 3*xi1 + xi2)/2,
F4 = (1 + 2*xi1 + xi2)/2,
F5 = (1 - xi2)/2;
size_t loc_idof;
loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(0,0,0));
power_index[perm[loc_idof]] = point_basic<size_t>(0,0,0);
value [perm[loc_idof]] = 1;
if (degree == 0) return;
loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(1,0,0));
power_index[perm[loc_idof]] = point_basic<size_t>(1,0,0);
value [perm[loc_idof]] = F1;
for (size_type i = 1; i+1 <= degree; i++) {
size_t loc_idof_i_0_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i ,0,0));
size_t loc_idof_im1_0_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i-1,0,0));
loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i+1,0,0));
power_index[perm[loc_idof]] = point_basic<size_t>(i+1,0,0);
value [perm[loc_idof]] = float_type(2*int(i)+1)/float_type(int(i)+1)*F1*value[perm[loc_idof_i_0_0]]
- float_type( int(i) )/float_type(int(i)+1)*F2*value[perm[loc_idof_im1_0_0]];
}
for (size_type i = 0; i+1 <= degree; i++) {
size_t loc_idof_i_0_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,0,0));
loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,1,0));
power_index[perm[loc_idof]] = point_basic<size_t>(i,1,0);
value [perm[loc_idof]] = (i*(1 + xi1) + F3)*value[perm[loc_idof_i_0_0]];
}
for (size_type i = 0; i+2 <= degree; i++) {
for (size_type j = 1; i+j+1 <= degree; j++) {
namespace J = rheolef::details::jacobi;
size_t loc_idof_i_j_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j ,0));
size_t loc_idof_i_jm1_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j-1,0));
loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j+1,0));
power_index[perm[loc_idof]] = point_basic<size_t>(i,j+1,0);
value [perm[loc_idof]] = (J::a<float_type>(2*int(i)+1,0,j)*F4 + J::b<float_type>(2*int(i)+1,0,j)*F5)*value[perm[loc_idof_i_j_0]]
- J::c<float_type>(2*int(i)+1,0,j)*sqr(F5) *value[perm[loc_idof_i_jm1_0]];
}}
for (size_type i = 0; i+1 <= degree; i++) {
for (size_type j = 0; i+j+1 <= degree; j++) {
size_t loc_idof_i_j_0 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,0));
loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,1));
power_index[perm[loc_idof]] = point_basic<size_t>(i,j,1);
value [perm[loc_idof]] = (1+i+j + (2+i+j)*xi2)*value[perm[loc_idof_i_j_0]];
}}
for (size_type i = 0; i+2 <= degree; i++) {
for (size_type j = 0; i+j+2 <= degree; j++) {
for (size_type k = 1; i+j+k+1 <= degree; k++) {
namespace J = rheolef::details::jacobi;
size_t loc_idof_i_j_k = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
size_t loc_idof_i_j_km1 = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,k-1));
loc_idof = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,k+1));
power_index[perm[loc_idof]] = point_basic<size_t>(i,j,k+1);
value [perm[loc_idof]] = (J::a<float_type>(2*int(i)+2*int(j)+2,0,k)*xi2 + J::b<float_type>(2*int(i)+2*int(j)+2,0,k))*value[perm[loc_idof_i_j_k]]
- J::c<float_type>(2*int(i)+2*int(j)+2,0,k) *value[perm[loc_idof_i_j_km1]];
}}}
}
// --------------------------------------------------
// main call: evaluate all the basis at a given point
// --------------------------------------------------
template<class T, class Container>
static
void
eval_dubiner_basis (
const point_basic<T>& hat_x,
reference_element hat_K,
size_t degree,
const std::vector<size_t>& perm,
Container& value,
std::vector<point_basic<size_t> >& power_index,
bool do_tilde = true)
{
typedef reference_element::size_type size_type;
value.resize (reference_element::n_node(hat_K.variant(), degree));
power_index.resize (value.size());
switch (hat_K.variant()) {
case reference_element::p: {
power_index[perm[0]] = point_basic<size_t>();
value [perm[0]] = 1;
break;
}
case reference_element::e: {
// HesWar-2008 refers to [-1,1] instead of hat_K=[0,1]:
T tilde_x = do_tilde ? 2*hat_x[0] - 1: hat_x[0];
eval_dubiner_basis_e (tilde_x, degree, perm, value, power_index);
break;
}
case reference_element::t: {
// Kir-2010 refers to [-1,1]^2 instead of [0,1]^2:
point_basic<T> tilde_x
= do_tilde ? point_basic<T>(2*hat_x[0]-1, 2*hat_x[1]-1) : hat_x;
eval_dubiner_basis_t (tilde_x, degree, perm, value, power_index);
break;
}
case reference_element::T: {
// Dubiner, by algorithm Kir-2010, page 5:7
// Kir-2010 refers to [-1,1]^2 instead of hat_K subset [0,1]^2
point_basic<T> tilde_x
= do_tilde ? point_basic<T>(2*hat_x[0]-1, 2*hat_x[1]-1, 2*hat_x[2]-1) : hat_x;
eval_dubiner_basis_T (tilde_x, degree, perm, value, power_index);
break;
}
case reference_element::q: {
Container value0, value1;
std::vector<size_t> id_e (degree+1);
for (size_type i = 0; i < id_e.size(); i++) { id_e [i] = i; }
std::vector<point_basic<size_t> > power_index_e;
eval_dubiner_basis_e (hat_x[0], degree, id_e, value0, power_index_e);
eval_dubiner_basis_e (hat_x[1], degree, id_e, value1, power_index_e);
for (size_type i = 0; i <= degree; i++) {
for (size_type j = 0; j <= degree; j++) {
size_t loc_idof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
size_t loc_jdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(j));
size_t loc_idof = reference_element_q::ilat2loc_inod (degree, point_basic<size_t>(i,j));
power_index[perm[loc_idof]] = point_basic<size_t>(i,j);
value [perm[loc_idof]] = value0[loc_idof_e]*value1[loc_jdof_e];
}}
break;
}
case reference_element::H: {
Container value0, value1, value2;
std::vector<size_t> id_e (degree+1);
for (size_type i = 0; i < id_e.size(); i++) { id_e [i] = i; }
std::vector<point_basic<size_t> > power_index_e;
eval_dubiner_basis_e (hat_x[0], degree, id_e, value0, power_index_e);
eval_dubiner_basis_e (hat_x[1], degree, id_e, value1, power_index_e);
eval_dubiner_basis_e (hat_x[2], degree, id_e, value2, power_index_e);
for (size_type i = 0; i <= degree; i++) {
for (size_type j = 0; j <= degree; j++) {
for (size_type k = 0; k <= degree; k++) {
size_t loc_idof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
size_t loc_jdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(j));
size_t loc_kdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(k));
size_t loc_idof = reference_element_H::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
power_index[perm[loc_idof]] = point_basic<size_t>(i,j,k);
value [perm[loc_idof]] = value0[loc_idof_e]*value1[loc_jdof_e]*value2[loc_kdof_e];
}}}
break;
}
case reference_element::P: {
Container value01, value2;
std::vector<size_t> id_e (degree+1), id_t ((degree+1)*(degree+2)/2);
for (size_type iloc = 0; iloc < id_e.size(); iloc ++) { id_e [iloc ] = iloc ; }
for (size_type iloc = 0; iloc < id_t.size(); iloc ++) { id_t [iloc ] = iloc ; }
point_basic<T> tilde_x
= do_tilde ? point_basic<T>(2*hat_x[0]-1, 2*hat_x[1]-1)
: point_basic<T>( hat_x[0], hat_x[1]);
std::vector<point_basic<size_t> > power_index_e, power_index_t;
eval_dubiner_basis_t (tilde_x, degree, id_t, value01, power_index_t);
eval_dubiner_basis_e (hat_x[2], degree, id_e, value2, power_index_e);
for (size_type i = 0; i <= degree; i++) {
for (size_type j = 0; i+j <= degree; j++) {
for (size_type k = 0; k <= degree; k++) {
size_t loc_ijdof_t = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j));
size_t loc_kdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(k));
size_t loc_idof = reference_element_P::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
power_index[perm[loc_idof]] = point_basic<size_t>(i,j,k);
value [perm[loc_idof]] = value01[loc_ijdof_t]*value2[loc_kdof_e];
}}}
break;
}
default:
error_macro ("unsupported element: "<<hat_K.name());
}
}
// omit the power_index last arg
template<class T, class Container>
static
void
eval_dubiner_basis (
const point_basic<T>& hat_x,
reference_element hat_K,
size_t degree,
const std::vector<size_t>& perm,
Container& value,
bool do_tilde = true)
{
std::vector<point_basic<size_t> > power_index;
eval_dubiner_basis (hat_x, hat_K, degree, perm, value, power_index, do_tilde);
}
}// namespace rheolef
#endif // _RHEOLEF_DUBINER_ICC
basis_option.h
basis_option - finite element method options
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::jacobi
Definition: jacobi.icc:24
rheolef::point_basic< size_t >
rheolef::reference_element_H::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:790
rheolef::reference_element_P::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:695
rheolef::reference_element_T::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:582
rheolef::reference_element_e::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:350
rheolef::reference_element_q::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:481
rheolef::reference_element_t::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:413
rheolef::reference_element::H
static const variant_type H
Definition: reference_element.h:81
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::reference_element::p
static const variant_type p
Definition: reference_element.h:75
rheolef::reference_element::size_type
std::vector< int >::size_type size_type
Definition: reference_element.h:71
rheolef::reference_element::n_node
static size_type n_node(variant_type variant, size_type order)
Definition: reference_element.cc:201
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
rheolef::reference_element::P
static const variant_type P
Definition: reference_element.h:80
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
rheolef::float_type
typename float_traits< value_type >::type float_type
Definition: field_expr_recursive.h:501
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
rheolef::value
rheolef::std value
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
T
Expr1::float_type T
Definition: field_expr.h:261
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
mkgeo_obstacle.L
string L
Definition: mkgeo_obstacle.sh:133
rheolef::details::jacobi
Definition: dubiner.icc:49
rheolef::details::legendre
Definition: dubiner.icc:42
rheolef::space_constant::vector
@ vector
Definition: space_constant.h:137
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rk::beta
Float beta[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:60
a
Definition: diffusion_isotropic.h:25
rheolef::float_traits::type
T type
Definition: Float.h:94