Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
sherwin.icc
Go to the documentation of this file.
1#ifndef _RHEOLEF_SHERWIN_ICC
2#define _RHEOLEF_SHERWIN_ICC
23//
24// evaluate the Dubiner-Sherwin basis on a point of the reference element
25// see SheKar-2005, pages 52 (1d) & 114 (2d) & app. D
26//
27// author: Pierre.Saramito@imag.fr
28//
29// date: 5 september 2017
30//
31#include "rheolef/quadrature.h"
32
33#include "basis_option.h"
35
36#ifdef SHERWIN_USE_GINAC
37#include <ginac/ginac.h>
38#endif // SHERWIN_USE_GINAC
39
40namespace rheolef {
41
42// --------------------------------------------------------------------------
43// compute all [0:order] Jacobi polynomials
44// --------------------------------------------------------------------------
45namespace details {
46
47template<class T0, class T, class Container>
48static
49void
50eval_jacobi (
51 const T0& tilde_x, // in [-1,1]
52 size_t order,
53 const T& alpha,
54 const T& beta,
55 Container& value)
56{
58 value.resize (order+1);
59 value[0] = 1;
60 if (order == 0) return;
61 value[1] = ((alpha+beta+2)*tilde_x + alpha - beta)/2;
62 for (size_t r = 1; r < order; r++) {
63 T a = (beta*beta - alpha*alpha)/((alpha+beta+T(2.*r))*(alpha+beta+T(2.*r+2)));
64 T b = 2*(alpha+T(1.*r))*(beta+T(1.*r))/((alpha+beta+T(2.*r))*(alpha+beta+T(2.*r+1)));
65 T c = 2*T(r+1.)*(alpha+beta+T(r+1.))/((alpha+beta+T(2.*r+1))*(alpha+beta+T(2.*r+2)));
66 value[r+1] = ((tilde_x-a)*value[r] - b*value[r-1])/c;
67 }
68}
69
70} // namespace details
71// --------------------------------------------------------------------------
72// Note: SheKar-2005 refers to tilde_e=[-1,1] instead of hat_e=[0,1]
73// --------------------------------------------------------------------------
74template<class T0, class T, class Container>
75static
76void
77eval_sherwin_basis_e (
78 const T0& tilde_x, // in [-1,1]
79 size_t degree,
80 const T& alpha,
81 const T& beta,
82 Container& jacobi, // working array
83 Container& value)
84{
86 typedef typename Container::value_type value_type;
87 typedef typename float_traits<value_type>::type float_type;
88 value.resize (degree+1);
89 size_t loc_idof;
90 if (degree == 0) {
91 loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(0));
92 value [loc_idof] = 1;
93 return;
94 }
95 loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(0));
96 value [loc_idof] = (1-tilde_x)/2;
97 loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(degree));
98 value [loc_idof] = (1+tilde_x)/2;
99 if (degree == 1) return;
100 details::eval_jacobi (tilde_x, degree-2, alpha, beta, jacobi);
101 for (size_type i = 1; i+1 <= degree; i++) {
102 loc_idof = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
103 value [loc_idof] = 0.25*(1-tilde_x)*(1+tilde_x)*jacobi [i-1];
104 }
105}
106template<class T, class Container>
107inline
108bool
110 const point_basic<T>& tilde_x, // in [-1,1]^2
111 reference_element hat_K,
112 size_t degree,
113 Container& value)
114{
115 typedef typename Container::value_type value_type;
116 typedef typename float_traits<value_type>::type float_type;
117 float_type eps = std::numeric_limits<float_type>::epsilon();
118 size_t d = hat_K.dimension();
119 if (d == 0) return false;
120 T distance_to_corner = fabs(tilde_x[d-1]-1);
121 if (distance_to_corner >= eps) {
122 return false;
123 }
124 // we are at top singular vertex:
125 for (size_t i = 0; i < size_t(value.size()); ++i) {
126 value[i] = value_type(0);
127 }
128 point_basic<size_t> ilat (0,0,0);
129 ilat [d] = degree;
130 size_t loc_idof = ilat2loc_inod (hat_K, degree, ilat);
131 value[loc_idof] = 1;
132 return true;
133}
134#ifdef SHERWIN_USE_GINAC
135// is_singular is not relevant by ginac (symbolic computations):
136template<>
137inline
138bool
139eval_sherwin_basis_is_singular_point<GiNaC::ex,std::vector<GiNaC::ex> > (
140 const point_basic<GiNaC::ex>& tilde_x, // in [-1,1]^2
141 reference_element hat_K,
142 size_t degree,
143 std::vector<GiNaC::ex>& value)
144{
145 return false;
146}
147#endif // SHERWIN_USE_GINAC
148// Dubiner polynoms on a triangle with C0 junctions
149// See KarShe-2005 p. 585
150template<class T0, class T, class Container>
151static
152void
153eval_sherwin_basis_t (
154 const point_basic<T0>& tilde_x, // in [-1,1]^2
155 size_t degree,
156 const T& alpha,
157 const T& beta,
158 Container& jacobi_x, // working array
159 Container& jacobi_y, // working array
160 Container& value)
161{
163 typedef typename Container::value_type value_type;
165 value.resize ((degree+1)*(degree+2)/2);
166 if (degree == 0) {
167 value[0] = 1;
168 return;
169 }
170 if (eval_sherwin_basis_is_singular_point (tilde_x, reference_element(reference_element::t), degree, value)) return;
171 value_type eta0 = 2*(1+tilde_x[0])/(1-tilde_x[1])-1,
172 eta1 = tilde_x[1];
173 size_t loc_idof;
174
175 // 1. verticies polynomials
176 loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(0,0));
177 value [loc_idof] = 0.25*(1-eta0)*(1-eta1);
178 loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(degree,0));
179 value [loc_idof] = 0.25*(1+eta0)*(1-eta1);
180 loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(0,degree));
181 value [loc_idof] = 0.5*(1+eta1);
182
183 // 2. edge polynomials
184 if (degree <= 1) return;
185 // 2.1. 1st edge: (-1,-1) --> (1,-1)
186 details::eval_jacobi (eta0, degree-2, T(1), T(1), jacobi_x);
187 value_type coef = 0.5*(1-eta1);
188 for (size_type i = 1; i <= degree-1; ++i) {
189 loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,0));
190 coef *= 0.5*(1-eta1);
191 value [loc_idof] = 0.25*(1-eta0)*(1+eta0)*coef*jacobi_x[i-1];
192 }
193 // 2.2. 2nd edge: (1,-1) --> (-1,1)
194 details::eval_jacobi (eta1, degree-2, T(1), T(1), jacobi_y);
195 coef = 0.5*(1+eta1);
196 for (size_type j = 1; j <= degree-1; ++j) {
197 loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(degree-j,j));
198 value [loc_idof] = 0.25*(1+eta0)*(1-eta1)*coef*jacobi_y[j-1];
199 }
200 // 2.3. 3rd edge: (-1,1) --> (-1,-1)
201 for (size_type j = 1; j <= degree-1; ++j) {
202 loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(0,j));
203 value [loc_idof] = 0.25*(1-eta0)*(1-eta1)*coef*jacobi_y[j-1];
204 }
205 // 3. interior polynomials
206 coef = 0.5*(1-eta1);
207 for (size_type i = 1; i <= degree-1; ++i) {
208 coef *= 0.5*(1-eta1);
209 details::eval_jacobi (eta1, degree-2, T(2*int(i)+1), T(1), jacobi_y);
210 for (size_type j = 1; i+j <= degree-1; ++j) {
211 loc_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j));
212 value [loc_idof] = 0.25*(1-eta0)*(1+eta0)*jacobi_x[i-1]*coef
213 *0.5*(1+eta1)*jacobi_y[j-1];
214 }
215 }
216}
217#ifdef TODO
218template<class T>
219static
220void
221eval_dubiner_basis_T (
222 const point_basic<T>& tilde_x, // in [-1,1]^3
223 size_t degree,
224 const std::vector<size_t>& perm,
225 Eigen::Matrix<T,Eigen::Dynamic,1>& value,
226 std::vector<point_basic<size_t> >& power_index)
227{
228}
229#endif // TODO
230// --------------------------------------------------
231// main call: evaluate all the basis at a given point
232// --------------------------------------------------
233template<class T0, class T, class Container>
234static
235void
236eval_sherwin_basis (
237 const point_basic<T0>& hat_x,
238 reference_element hat_K,
239 size_t degree,
240 const T& alpha,
241 const T& beta,
242 Container& jacobi_x, // working array
243 Container& jacobi_y, // working array
244 Container& jacobi_z, // working array
245 Container& value,
246 bool do_tilde = true)
247{
249 typedef typename Container::value_type value_type;
250 typedef typename float_traits<value_type>::type float_type;
251 value.resize (reference_element::n_node(hat_K.variant(), degree));
252 switch (hat_K.variant()) {
254 value [0] = 1;
255 break;
256 }
258 // HesWar-2008 refers to [-1,1] instead of hat_K=[0,1]:
259 value_type tilde_x = do_tilde ? 2*hat_x[0] - 1: hat_x[0];
260 eval_sherwin_basis_e (tilde_x, degree, alpha, beta, jacobi_x, value);
261 break;
262 }
264 // Kir-2010 refers to [-1,1]^2 instead of [0,1]^2:
265 point_basic<value_type> tilde_x
266 = do_tilde ? point_basic<value_type>(2*hat_x[0]-1, 2*hat_x[1]-1)
267 : point_basic<value_type>( hat_x[0], hat_x[1]);
268 eval_sherwin_basis_t (tilde_x, degree, alpha, beta, jacobi_x, jacobi_y, value);
269 break;
270 }
271#ifdef TODO
273 // Dubiner, by algorithm Kir-2010, page 5:7
274 // Kir-2010 refers to [-1,1]^2 instead of hat_K subset [0,1]^2
275 point_basic<T> tilde_x
276 = do_tilde ? point_basic<T>(2*hat_x[0]-1, 2*hat_x[1]-1, 2*hat_x[2]-1) : hat_x;
277 eval_dubiner_basis_T (tilde_x, degree, perm, value, power_index);
278 break;
279 }
281 Eigen::Matrix<T,Eigen::Dynamic,1> value0, value1;
282 std::vector<size_t> id_e (degree+1);
283 for (size_type i = 0; i < id_e.size(); i++) { id_e [i] = i; }
284 std::vector<point_basic<size_t> > power_index_e;
285 eval_dubiner_basis_e (hat_x[0], degree, id_e, value0, power_index_e);
286 eval_dubiner_basis_e (hat_x[1], degree, id_e, value1, power_index_e);
287 for (size_type i = 0; i <= degree; i++) {
288 for (size_type j = 0; j <= degree; j++) {
289 size_t loc_idof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
290 size_t loc_jdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(j));
291 size_t loc_idof = reference_element_q::ilat2loc_inod (degree, point_basic<size_t>(i,j));
292 power_index[perm[loc_idof]] = point_basic<size_t>(i,j);
293 value [perm[loc_idof]] = value0[loc_idof_e]*value1[loc_jdof_e];
294 }}
295 break;
296 }
298 Eigen::Matrix<T,Eigen::Dynamic,1> value0, value1, value2;
299 std::vector<size_t> id_e (degree+1);
300 for (size_type i = 0; i < id_e.size(); i++) { id_e [i] = i; }
301 std::vector<point_basic<size_t> > power_index_e;
302 eval_dubiner_basis_e (hat_x[0], degree, id_e, value0, power_index_e);
303 eval_dubiner_basis_e (hat_x[1], degree, id_e, value1, power_index_e);
304 eval_dubiner_basis_e (hat_x[2], degree, id_e, value2, power_index_e);
305 for (size_type i = 0; i <= degree; i++) {
306 for (size_type j = 0; j <= degree; j++) {
307 for (size_type k = 0; k <= degree; k++) {
308 size_t loc_idof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
309 size_t loc_jdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(j));
310 size_t loc_kdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(k));
311 size_t loc_idof = reference_element_H::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
312 power_index[perm[loc_idof]] = point_basic<size_t>(i,j,k);
313 value [perm[loc_idof]] = value0[loc_idof_e]*value1[loc_jdof_e]*value2[loc_kdof_e];
314 }}}
315 break;
316 }
318 Eigen::Matrix<T,Eigen::Dynamic,1> value01, value2;
319 std::vector<size_t> id_e (degree+1), id_t ((degree+1)*(degree+2)/2);
320 for (size_type iloc = 0; iloc < id_e.size(); iloc ++) { id_e [iloc ] = iloc ; }
321 for (size_type iloc = 0; iloc < id_t.size(); iloc ++) { id_t [iloc ] = iloc ; }
322 point_basic<T> tilde_x
323 = do_tilde ? point_basic<T>(2*hat_x[0]-1, 2*hat_x[1]-1)
324 : point_basic<T>( hat_x[0], hat_x[1]);
325 std::vector<point_basic<size_t> > power_index_e, power_index_t;
326 eval_dubiner_basis_t (tilde_x, degree, id_t, value01, power_index_t);
327 eval_dubiner_basis_e (hat_x[2], degree, id_e, value2, power_index_e);
328 for (size_type i = 0; i <= degree; i++) {
329 for (size_type j = 0; i+j <= degree; j++) {
330 for (size_type k = 0; k <= degree; k++) {
331 size_t loc_ijdof_t = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j));
332 size_t loc_kdof_e = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(k));
333 size_t loc_idof = reference_element_P::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
334 power_index[perm[loc_idof]] = point_basic<size_t>(i,j,k);
335 value [perm[loc_idof]] = value01[loc_ijdof_t]*value2[loc_kdof_e];
336 }}}
337 break;
338 }
339#endif // TODO
340 default:
341 error_macro ("unsupported element: "<<hat_K.name());
342 }
343}
344
345}// namespace rheolef
346#endif // _RHEOLEF_SHERWIN_ICC
basis_option - finite element method options
field::size_type size_type
Definition branch.cc:430
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
see the reference_element page for the full documentation
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 size_type n_node(variant_type variant, size_type order)
static const variant_type T
static const variant_type P
static const variant_type t
typename float_traits< value_type >::type float_type
result_type value_type
#define error_macro(message)
Definition dis_macros.h:49
Expr1::float_type T
Definition field_expr.h:230
This file is part of Rheolef.
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