Loading...
Searching...
No Matches
bernstein.icc
Go to the documentation of this file.
#ifndef _RHEOLEF_BERNSTEIN_ICC
#define _RHEOLEF_BERNSTEIN_ICC
//
// Bernstein basis, as initial local FEM basis
//
// author: Pierre.Saramito@imag.fr
//
// date: 2 september 2017
//
#include "rheolef/reference_element.h"
namespace rheolef {
using namespace std;
// pre-compute factorial[i] = i! for i=0..degree
// thus, avoid some inner loops
template<class T>
static
void
precompute_factorial (
size_t degree,
std::vector<T>& factorial)
{
typedef reference_element::size_type size_type;
factorial.resize (degree+1);
factorial[0] = 1;
for (size_type i = 1; i <= degree; i++) {
factorial[i] = factorial[i-1]*i;
}
}
// lambda[1][mu] = i-th barycentric coordinate, mu=0..d-1
// lambda[i][mu] = (lambda[1][mu])^i, i=0..degree, mu=0..d-1
// manage also variants for tensor-product elements (qPH)
template<class T>
static
void
precompute_power_bernstein (
reference_element hat_K,
size_t d,
const point_basic<T>& x,
size_t degree,
std::vector<std::array<T,6> >& lambda)
{
lambda.resize (degree+1);
size_t n_comp = 0;
switch (hat_K.variant()) {
case reference_element::p: n_comp = 0; break;
case reference_element::e:
case reference_element::t:
case reference_element::T: n_comp = d+1; break;
case reference_element::q:
case reference_element::H: n_comp = 2*d; break;
case reference_element::P: n_comp = 5; break;
default: error_macro ("unsupported element: "<<hat_K.name());
}
for (size_t mu = 0; mu < n_comp; mu++) {
lambda[0][mu] = 1;
}
if (degree == 0) return;
switch (hat_K.variant()) {
case reference_element::p:
break;
case reference_element::e:
lambda[1][0] = x[0];
lambda[1][1] = 1-x[0];
break;
case reference_element::t:
lambda[1][0] = x[0];
lambda[1][1] = x[1];
lambda[1][2] = 1-x[0]-x[1];
break;
case reference_element::T:
lambda[1][0] = x[0];
lambda[1][1] = x[1];
lambda[1][2] = x[2];
lambda[1][3] = 1-x[0]-x[1]-x[2];
break;
case reference_element::q: {
T x0 = (1+x[0])/2,
x1 = (1+x[1])/2;
lambda[1][0] = x0;
lambda[1][1] = 1-x0;
lambda[1][2] = x1;
lambda[1][3] = 1-x1;
break;
}
case reference_element::P: {
T x_2 = (1+x[2])/2;
lambda[1][0] = x[0];
lambda[1][1] = x[1];
lambda[1][2] = 1-x[0]-x[1];
lambda[1][3] = x_2;
lambda[1][4] = 1-x_2;
break;
}
case reference_element::H: {
T x0 = (1+x[0])/2,
x1 = (1+x[1])/2,
x2 = (1+x[2])/2;
lambda[1][0] = x0;
lambda[1][1] = 1-x0;
lambda[1][2] = x1;
lambda[1][3] = 1-x1;
lambda[1][4] = x2;
lambda[1][5] = 1-x2;
break;
}
}
for (size_t i = 2; i <= degree; i++) {
for (size_t mu = 0; mu < n_comp; mu++) {
lambda[i][mu] = lambda[i-1][mu]*lambda[1][mu];
}
}
}
// ----------------------------------------------
// Bernstein basis evaluation
// result = c * prod_{i=0}^{d} lambda[i]^m[i]
// with C = s!/(prod m[i]!)
// s = sum_{i=0}^d m[i]!
// m[0] = sum_{i=1}^d m[i]!
// use precomputed i! and x[i]^m for efficiency
// Bernstein basis functions on quadrilaterals
// or hexahedra are defined using
// a tensor product construction [ChaWar-2016]
// ----------------------------------------------
template<class T, class U>
static
U
eval_bernstein_internal (
reference_element hat_K,
size_t d,
const std::vector<std::array<U,6> >& lambda_m,
const std::vector<T>& factorial,
const point_basic<size_t>& m,
size_t degree)
{
switch (hat_K.variant()) {
case reference_element::p: {
return 1;
}
case reference_element::e: {
T deno = factorial[ m[0]]
*factorial[degree-m[0]];
T c = factorial[degree]/deno;
return c*lambda_m[ m[0]][0]
*lambda_m[degree-m[0]][1];
}
case reference_element::t: {
T deno = factorial[ m[0] ]
*factorial[ m[1]]
*factorial[degree-m[0]-m[1]];
T c = factorial[degree]/deno;
return c*lambda_m[ m[0] ][0]
*lambda_m[ m[1]][1]
*lambda_m[degree-m[0]-m[1]][2];
}
case reference_element::T: {
T deno = factorial[ m[0] ]
*factorial[ m[1] ]
*factorial[ m[2]]
*factorial[degree-m[0]-m[1]-m[2]];
T c = factorial[degree]/deno;
return c*lambda_m[ m[0] ][0]
*lambda_m[ m[1] ][1]
*lambda_m[ m[2]][2]
*lambda_m[degree-m[0]-m[1]-m[2]][3];
}
case reference_element::q: {
T deno = factorial[ m[0]]
*factorial[degree-m[0]]
*factorial[ m[1]]
*factorial[degree-m[1]];
T c = sqr(factorial[degree])/deno;
return c*lambda_m[ m[0]][0]
*lambda_m[degree-m[0]][1]
*lambda_m[ m[1]][2]
*lambda_m[degree-m[1]][3];
}
case reference_element::H: {
T deno = factorial[ m[0]]
*factorial[degree-m[0]]
*factorial[ m[1]]
*factorial[degree-m[1]]
*factorial[ m[2]]
*factorial[degree-m[2]];
T c = pow(factorial[degree],3)/deno;
return c*lambda_m[ m[0]][0]
*lambda_m[degree-m[0]][1]
*lambda_m[ m[1]][2]
*lambda_m[degree-m[1]][3]
*lambda_m[ m[2]][4]
*lambda_m[degree-m[2]][5];
}
case reference_element::P: {
T deno = factorial[ m[0] ]
*factorial[ m[1]]
*factorial[degree-m[0]-m[1]]
*factorial[ m[2]]
*factorial[degree-m[2]];
T c = sqr(factorial[degree])/deno;
return c*lambda_m[ m[0] ][0]
*lambda_m[ m[1]][1]
*lambda_m[degree-m[0]-m[1]][2]
*lambda_m[ m[2] ][3]
*lambda_m[degree-m[2] ][4];
}
default:
error_macro ("unsupported element: "<<hat_K.name());
return 0;
}
}
}// namespace rheolef
#endif // _RHEOLEF_BERNSTEIN_ICC
size_type
field::size_type size_type
Definition branch.cc:430
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::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
error_macro
#define error_macro(message)
Definition dis_macros.h:49
T
Expr1::float_type T
Definition field_expr.h:230
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
std
STL namespace.
d
size_t d
Definition neumann-laplace-lambda.cc:28
lambda
Definition yield_slip_circle.h:34