Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
bernstein.icc
Go to the documentation of this file.
1#ifndef _RHEOLEF_BERNSTEIN_ICC
2#define _RHEOLEF_BERNSTEIN_ICC
23//
24// Bernstein basis, as initial local FEM basis
25//
26// author: Pierre.Saramito@imag.fr
27//
28// date: 2 september 2017
29//
30#include "rheolef/reference_element.h"
31
32namespace rheolef {
33using namespace std;
34
35// pre-compute factorial[i] = i! for i=0..degree
36// thus, avoid some inner loops
37template<class T>
38static
39void
40precompute_factorial (
41 size_t degree,
42 std::vector<T>& factorial)
43{
45 factorial.resize (degree+1);
46 factorial[0] = 1;
47 for (size_type i = 1; i <= degree; i++) {
48 factorial[i] = factorial[i-1]*i;
49 }
50}
51// lambda[1][mu] = i-th barycentric coordinate, mu=0..d-1
52// lambda[i][mu] = (lambda[1][mu])^i, i=0..degree, mu=0..d-1
53// manage also variants for tensor-product elements (qPH)
54template<class T>
55static
56void
57precompute_power_bernstein (
58 reference_element hat_K,
59 size_t d,
60 const point_basic<T>& x,
61 size_t degree,
62 std::vector<std::array<T,6> >& lambda)
63{
64 lambda.resize (degree+1);
65 size_t n_comp = 0;
66 switch (hat_K.variant()) {
67 case reference_element::p: n_comp = 0; break;
70 case reference_element::T: n_comp = d+1; break;
72 case reference_element::H: n_comp = 2*d; break;
73 case reference_element::P: n_comp = 5; break;
74 default: error_macro ("unsupported element: "<<hat_K.name());
75 }
76 for (size_t mu = 0; mu < n_comp; mu++) {
77 lambda[0][mu] = 1;
78 }
79 if (degree == 0) return;
80 switch (hat_K.variant()) {
82 break;
84 lambda[1][0] = x[0];
85 lambda[1][1] = 1-x[0];
86 break;
88 lambda[1][0] = x[0];
89 lambda[1][1] = x[1];
90 lambda[1][2] = 1-x[0]-x[1];
91 break;
93 lambda[1][0] = x[0];
94 lambda[1][1] = x[1];
95 lambda[1][2] = x[2];
96 lambda[1][3] = 1-x[0]-x[1]-x[2];
97 break;
99 T x0 = (1+x[0])/2,
100 x1 = (1+x[1])/2;
101 lambda[1][0] = x0;
102 lambda[1][1] = 1-x0;
103 lambda[1][2] = x1;
104 lambda[1][3] = 1-x1;
105 break;
106 }
108 T x_2 = (1+x[2])/2;
109 lambda[1][0] = x[0];
110 lambda[1][1] = x[1];
111 lambda[1][2] = 1-x[0]-x[1];
112 lambda[1][3] = x_2;
113 lambda[1][4] = 1-x_2;
114 break;
115 }
117 T x0 = (1+x[0])/2,
118 x1 = (1+x[1])/2,
119 x2 = (1+x[2])/2;
120 lambda[1][0] = x0;
121 lambda[1][1] = 1-x0;
122 lambda[1][2] = x1;
123 lambda[1][3] = 1-x1;
124 lambda[1][4] = x2;
125 lambda[1][5] = 1-x2;
126 break;
127 }
128 }
129 for (size_t i = 2; i <= degree; i++) {
130 for (size_t mu = 0; mu < n_comp; mu++) {
131 lambda[i][mu] = lambda[i-1][mu]*lambda[1][mu];
132 }
133 }
134}
135// ----------------------------------------------
136// Bernstein basis evaluation
137// result = c * prod_{i=0}^{d} lambda[i]^m[i]
138// with C = s!/(prod m[i]!)
139// s = sum_{i=0}^d m[i]!
140// m[0] = sum_{i=1}^d m[i]!
141// use precomputed i! and x[i]^m for efficiency
142// Bernstein basis functions on quadrilaterals
143// or hexahedra are defined using
144// a tensor product construction [ChaWar-2016]
145// ----------------------------------------------
146template<class T, class U>
147static
148U
149eval_bernstein_internal (
150 reference_element hat_K,
151 size_t d,
152 const std::vector<std::array<U,6> >& lambda_m,
153 const std::vector<T>& factorial,
154 const point_basic<size_t>& m,
155 size_t degree)
156{
157 switch (hat_K.variant()) {
159 return 1;
160 }
162 T deno = factorial[ m[0]]
163 *factorial[degree-m[0]];
164 T c = factorial[degree]/deno;
165 return c*lambda_m[ m[0]][0]
166 *lambda_m[degree-m[0]][1];
167 }
169 T deno = factorial[ m[0] ]
170 *factorial[ m[1]]
171 *factorial[degree-m[0]-m[1]];
172 T c = factorial[degree]/deno;
173 return c*lambda_m[ m[0] ][0]
174 *lambda_m[ m[1]][1]
175 *lambda_m[degree-m[0]-m[1]][2];
176 }
178 T deno = factorial[ m[0] ]
179 *factorial[ m[1] ]
180 *factorial[ m[2]]
181 *factorial[degree-m[0]-m[1]-m[2]];
182 T c = factorial[degree]/deno;
183 return c*lambda_m[ m[0] ][0]
184 *lambda_m[ m[1] ][1]
185 *lambda_m[ m[2]][2]
186 *lambda_m[degree-m[0]-m[1]-m[2]][3];
187 }
189 T deno = factorial[ m[0]]
190 *factorial[degree-m[0]]
191 *factorial[ m[1]]
192 *factorial[degree-m[1]];
193 T c = sqr(factorial[degree])/deno;
194 return c*lambda_m[ m[0]][0]
195 *lambda_m[degree-m[0]][1]
196 *lambda_m[ m[1]][2]
197 *lambda_m[degree-m[1]][3];
198 }
200 T deno = factorial[ m[0]]
201 *factorial[degree-m[0]]
202 *factorial[ m[1]]
203 *factorial[degree-m[1]]
204 *factorial[ m[2]]
205 *factorial[degree-m[2]];
206 T c = pow(factorial[degree],3)/deno;
207 return c*lambda_m[ m[0]][0]
208 *lambda_m[degree-m[0]][1]
209 *lambda_m[ m[1]][2]
210 *lambda_m[degree-m[1]][3]
211 *lambda_m[ m[2]][4]
212 *lambda_m[degree-m[2]][5];
213 }
215 T deno = factorial[ m[0] ]
216 *factorial[ m[1]]
217 *factorial[degree-m[0]-m[1]]
218 *factorial[ m[2]]
219 *factorial[degree-m[2]];
220 T c = sqr(factorial[degree])/deno;
221 return c*lambda_m[ m[0] ][0]
222 *lambda_m[ m[1]][1]
223 *lambda_m[degree-m[0]-m[1]][2]
224 *lambda_m[ m[2] ][3]
225 *lambda_m[degree-m[2] ][4];
226 }
227 default:
228 error_macro ("unsupported element: "<<hat_K.name());
229 return 0;
230 }
231}
232
233}// namespace rheolef
234#endif // _RHEOLEF_BERNSTEIN_ICC
field::size_type size_type
Definition branch.cc:430
static const variant_type H
static const variant_type q
static const variant_type e
static const variant_type p
std::vector< int >::size_type size_type
static const variant_type T
static const variant_type P
static const variant_type t
#define error_macro(message)
Definition dis_macros.h:49
Expr1::float_type T
Definition field_expr.h:230
This file is part of Rheolef.
STL namespace.