sherwin.icc
Go to the documentation of this file.
#ifndef _RHEOLEF_SHERWIN_ICC
#define _RHEOLEF_SHERWIN_ICC
//
// evaluate the Dubiner-Sherwin basis on a point of the reference element
// see SheKar-2005, pages 52 (1d) & 114 (2d) & app. D
//
// author: Pierre.Saramito@imag.fr
//
// date: 5 september 2017
//
#include "rheolef/quadrature.h"
#include "basis_option.h"
#include "reference_element_aux.icc"
#ifdef SHERWIN_USE_GINAC
#include <ginac/ginac.h>
#endif // SHERWIN_USE_GINAC
namespace rheolef {
// --------------------------------------------------------------------------
// compute all [0:order] Jacobi polynomials
// --------------------------------------------------------------------------
namespace details {
template<class T0, class T, class Container>
static
void
eval_jacobi (
const T0& tilde_x, // in [-1,1]
size_t order,
const T& alpha,
const T& beta,
Container& value)
{
typedef reference_element::size_type size_type;
value.resize (order+1);
value[0] = 1;
if (order == 0) return;
value[1] = ((alpha+beta+2)*tilde_x + alpha - beta)/2;
for (size_t r = 1; r < order; r++) {
T a = (beta*beta - alpha*alpha)/((alpha+beta+T(2.*r))*(alpha+beta+T(2.*r+2)));
T b = 2*(alpha+T(1.*r))*(beta+T(1.*r))/((alpha+beta+T(2.*r))*(alpha+beta+T(2.*r+1)));
T c = 2*T(r+1.)*(alpha+beta+T(r+1.))/((alpha+beta+T(2.*r+1))*(alpha+beta+T(2.*r+2)));
value[r+1] = ((tilde_x-a)*value[r] - b*value[r-1])/c;
}
}
} // namespace details
// --------------------------------------------------------------------------
// Note: SheKar-2005 refers to tilde_e=[-1,1] instead of hat_e=[0,1]
// --------------------------------------------------------------------------
template<class T0, class T, class Container>
static
void
eval_sherwin_basis_e (
const T0& tilde_x, // in [-1,1]
size_t degree,
const T& alpha,
const T& beta,
Container& jacobi, // working array
Container& value)
{
typedef reference_element::size_type size_type;
typedef typename Container::value_type value_type;
typedef typename float_traits<value_type>::type float_type;
value.resize (degree+1);
size_t loc_idof;
if (degree == 0) {
loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(0));
value [loc_idof] = 1;
return;
}
loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(0));
value [loc_idof] = (1-tilde_x)/2;
loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(degree));
value [loc_idof] = (1+tilde_x)/2;
if (degree == 1) return;
details::eval_jacobi (tilde_x, degree-2, alpha, beta, jacobi);
for (size_type i = 1; i+1 <= degree; i++) {
loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
value [loc_idof] = 0.25*(1-tilde_x)*(1+tilde_x)*jacobi [i-1];
}
}
template<class T, class Container>
inline
bool
eval_sherwin_basis_is_singular_point (
const point_basic<T>& tilde_x, // in [-1,1]^2
reference_element hat_K,
size_t degree,
Container& value)
{
typedef typename Container::value_type value_type;
typedef typename float_traits<value_type>::type float_type;
float_type eps = std::numeric_limits<float_type>::epsilon();
size_t d = hat_K.dimension();
if (d == 0) return false;
T distance_to_corner = fabs(tilde_x[d-1]-1);
if (distance_to_corner >= eps) {
return false;
}
// we are at top singular vertex:
for (size_t i = 0; i < size_t(value.size()); ++i) {
value[i] = value_type(0);
}
point_basic<size_t> ilat (0,0,0);
ilat [d] = degree;
size_t loc_idof = ilat2loc_inod (hat_K, degree, ilat);
value[loc_idof] = 1;
return true;
}
#ifdef SHERWIN_USE_GINAC
// is_singular is not relevant by ginac (symbolic computations):
template<>
inline
bool
eval_sherwin_basis_is_singular_point<GiNaC::ex,std::vector<GiNaC::ex> > (
const point_basic<GiNaC::ex>& tilde_x, // in [-1,1]^2
reference_element hat_K,
size_t degree,
std::vector<GiNaC::ex>& value)
{
return false;
}
#endif // SHERWIN_USE_GINAC
// Dubiner polynoms on a triangle with C0 junctions
// See KarShe-2005 p. 585
template<class T0, class T, class Container>
static
void
eval_sherwin_basis_t (
const point_basic<T0>& tilde_x, // in [-1,1]^2
size_t degree,
const T& alpha,
const T& beta,
Container& jacobi_x, // working array
Container& jacobi_y, // working array
Container& value)
{
typedef reference_element::size_type size_type;
typedef typename Container::value_type value_type;
typedef typename float_traits<value_type>::type float_type;
value.resize ((degree+1)*(degree+2)/2);
if (degree == 0) {
value[0] = 1;
return;
}
if (eval_sherwin_basis_is_singular_point (tilde_x, reference_element(reference_element::t), degree, value)) return;
value_type eta0 = 2*(1+tilde_x[0])/(1-tilde_x[1])-1,
eta1 = tilde_x[1];
size_t loc_idof;
// 1. verticies polynomials
loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(0,0));
value [loc_idof] = 0.25*(1-eta0)*(1-eta1);
loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(degree,0));
value [loc_idof] = 0.25*(1+eta0)*(1-eta1);
loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(0,degree));
value [loc_idof] = 0.5*(1+eta1);
// 2. edge polynomials
if (degree <= 1) return;
// 2.1. 1st edge: (-1,-1) --> (1,-1)
details::eval_jacobi (eta0, degree-2, T(1), T(1), jacobi_x);
value_type coef = 0.5*(1-eta1);
for (size_type i = 1; i <= degree-1; ++i) {
loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,0));
coef *= 0.5*(1-eta1);
value [loc_idof] = 0.25*(1-eta0)*(1+eta0)*coef*jacobi_x[i-1];
}
// 2.2. 2nd edge: (1,-1) --> (-1,1)
details::eval_jacobi (eta1, degree-2, T(1), T(1), jacobi_y);
coef = 0.5*(1+eta1);
for (size_type j = 1; j <= degree-1; ++j) {
loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(degree-j,j));
value [loc_idof] = 0.25*(1+eta0)*(1-eta1)*coef*jacobi_y[j-1];
}
// 2.3. 3rd edge: (-1,1) --> (-1,-1)
for (size_type j = 1; j <= degree-1; ++j) {
loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(0,j));
value [loc_idof] = 0.25*(1-eta0)*(1-eta1)*coef*jacobi_y[j-1];
}
// 3. interior polynomials
coef = 0.5*(1-eta1);
for (size_type i = 1; i <= degree-1; ++i) {
coef *= 0.5*(1-eta1);
details::eval_jacobi (eta1, degree-2, T(2*int(i)+1), T(1), jacobi_y);
for (size_type j = 1; i+j <= degree-1; ++j) {
loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j));
value [loc_idof] = 0.25*(1-eta0)*(1+eta0)*jacobi_x[i-1]*coef
*0.5*(1+eta1)*jacobi_y[j-1];
}
}
}
#ifdef TODO
template<class T>
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,
Eigen::Matrix<T,Eigen::Dynamic,1>& value,
std::vector<point_basic<size_t> >& power_index)
{
}
#endif // TODO
// --------------------------------------------------
// main call: evaluate all the basis at a given point
// --------------------------------------------------
template<class T0, class T, class Container>
static
void
eval_sherwin_basis (
const point_basic<T0>& hat_x,
reference_element hat_K,
size_t degree,
const T& alpha,
const T& beta,
Container& jacobi_x, // working array
Container& jacobi_y, // working array
Container& jacobi_z, // working array
Container& value,
bool do_tilde = true)
{
typedef reference_element::size_type size_type;
typedef typename Container::value_type value_type;
typedef typename float_traits<value_type>::type float_type;
value.resize (reference_element::n_node(hat_K.variant(), degree));
switch (hat_K.variant()) {
case reference_element::p: {
value [0] = 1;
break;
}
case reference_element::e: {
// HesWar-2008 refers to [-1,1] instead of hat_K=[0,1]:
value_type tilde_x = do_tilde ? 2*hat_x[0] - 1: hat_x[0];
eval_sherwin_basis_e (tilde_x, degree, alpha, beta, jacobi_x, value);
break;
}
case reference_element::t: {
// Kir-2010 refers to [-1,1]^2 instead of [0,1]^2:
point_basic<value_type> tilde_x
= do_tilde ? point_basic<value_type>(2*hat_x[0]-1, 2*hat_x[1]-1)
: point_basic<value_type>( hat_x[0], hat_x[1]);
eval_sherwin_basis_t (tilde_x, degree, alpha, beta, jacobi_x, jacobi_y, value);
break;
}
#ifdef TODO
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: {
Eigen::Matrix<T,Eigen::Dynamic,1> 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: {
Eigen::Matrix<T,Eigen::Dynamic,1> 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: {
Eigen::Matrix<T,Eigen::Dynamic,1> 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;
}
#endif // TODO
default:
error_macro ("unsupported element: "<<hat_K.name());
}
}
}// namespace rheolef
#endif // _RHEOLEF_SHERWIN_ICC
basis_option.h
basis_option - finite element method options
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::point_basic
Definition: point.h:87
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_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
see the reference_element page for the full documentation
Definition: reference_element.h:66
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::dimension
size_type dimension() const
Definition: reference_element.h:101
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
rheolef::value_type
result_type value_type
Definition: field_expr_recursive.h:499
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
T
Expr1::float_type T
Definition: field_expr.h:261
reference_element_aux.icc
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
mkgeo_ball.a
int a
Definition: mkgeo_ball.sh:151
rheolef::space_constant::vector
@ vector
Definition: space_constant.h:137
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::eval_sherwin_basis_is_singular_point
bool eval_sherwin_basis_is_singular_point(const point_basic< T > &tilde_x, reference_element hat_K, size_t degree, Container &value)
Definition: sherwin.icc:109
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
epsilon
Float epsilon
Definition: transmission_error.cc:25