Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
monomial.icc
Go to the documentation of this file.
1#ifndef _RHEOLEF_MONOMIAL_ICC
2#define _RHEOLEF_MONOMIAL_ICC
23//
24// monomial basis, as initial local FEM basis
25// see also bernstein and dubiner that are better conditionned
26//
27// author: Pierre.Saramito@imag.fr
28//
29// date: 2 september 2017
30//
31// note: initially copied from rheolef/fem/lib/Pk.cc
32//
33#include "rheolef/reference_element.h"
34#include "basis_ordering.icc"
35
36namespace rheolef {
37
38// mx[k][i] = x[i]^k, k=0..degree, i=0..d-1
39template<class T>
40static
41void
42precompute_power_monomial (
43 reference_element hat_K,
44 size_t d,
45 const point_basic<T>& x,
46 size_t degree,
47 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& xm)
48{
49 xm.resize (degree+1);
50 for (size_t mu = 0; mu < d; mu++) {
51 xm[0][mu] = 1;
52 }
53 if (degree == 0) return;
54 xm[1] = x;
55 for (size_t i = 2; i <= degree; i++) {
56 for (size_t mu = 0; mu < d; mu++) {
57 xm[i][mu] = xm[i-1][mu]*xm[1][mu];
58 }
59 }
60}
61// ----------------------------------------------
62// monomial basis evaluation
63// result = prod_{i=0}^{d-1} x[i]^m[i]
64// use predomputed x[i]^m for efficiency
65// ----------------------------------------------
66template<class T>
67static
68T
69eval_monomial_internal (
70 reference_element hat_K,
71 size_t d,
72 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& xm,
73 const point_basic<size_t>& m) // power_index
74{
75 T result;
76 if (d == 2) {
77 result = xm[m[0]][0]
78 *xm[m[1]][1];
79 } else if (d == 3) {
80 result = xm[m[0]][0]
81 *xm[m[1]][1]
82 *xm[m[2]][2];
83 } else if (d == 1) {
84 result = xm[m[0]][0];
85 } else { // d == 0
86 result = 1;
87 }
88 return result;
89}
90
91}// namespace rheolef
92#endif // _RHEOLEF_MONOMIAL_ICC
Expr1::float_type T
Definition field_expr.h:230
This file is part of Rheolef.