29 #ifndef __expansions_H
30 #define __expansions_H 1
38 enum Shape { LINE = 1, TRIANGLE = 2, TETRAHEDRON = 3 };
42 static const uint16_type value = mpl::if_<mpl::equal_to<mpl::int_<sh>, mpl::int_<LINE> >,
44 typename mpl::if_<mpl::equal_to<mpl::int_<sh>, mpl::int_<TRIANGLE> >,
46 mpl::int_<3> >::type>::type::value;
59 template<Shape sh,
typename T =
double>
65 struct xi<TRIANGLE, T>
68 typedef typename node<value_type>::type node_type;
75 xi( node_type
const& eta )
79 M_xi[0] = 0.5*( 1.0 + eta[0] )*( 1.0 - eta[1] ) - 1.0;
83 node_type
const& operator()( node_type
const& eta )
85 M_xi[0] = 0.5*( 1.0 + eta[0] )*( 1.0 - eta[1] ) - 1.0;
89 node_type
const& operator()()
const
97 struct xi<TETRAHEDRON, T>
100 typedef typename node<value_type>::type node_type;
107 xi( node_type
const& eta )
111 M_xi[0] = 0.25*( 1.0 + eta[0] )*( 1.0 - eta[1] )*( 1.0 - eta[2] ) - 1.0;
112 M_xi[1] = 0.5*( 1.0 + eta[1] )*( 1.0 - eta[2] ) - 1.0;
116 node_type
const& operator()( node_type
const& eta )
118 M_xi[0] = 0.25*( 1.0 + eta[0] )*( 1.0 - eta[1] )*( 1.0 - eta[2] ) - 1.0;
119 M_xi[1] = 0.5*( 1.0 + eta[1] )*( 1.0 - eta[2] ) - 1.0;
123 node_type
const& operator()()
const
134 template<Shape sh,
typename T =
double>
138 template<Shape sh,
typename T =
double>
144 struct eta<TRIANGLE, T>
146 typedef T value_type;
147 typedef typename node<value_type>::type node_type;
154 eta( node_type
const& xi )
162 M_eta[0] = 2.0 * ( 1.0 + xi[0] ) / ( 1.0 - xi[1] ) - 1.0;
167 node_type
const& operator()( node_type
const& xi )
173 M_eta[0] = 2.0 * ( 1.0 + xi[0] ) / ( 1.0 - xi[1] ) - 1.0;
178 node_type
const& operator()()
const
185 struct etas<TRIANGLE, T>
187 typedef T value_type;
188 typedef typename matrix_node<value_type>::type matrix_node_type;
195 etas( matrix_node_type
const& xi )
197 M_eta( xi.size1(), xi.size2() )
199 for (
size_type i = 0; i < xi.size2(); ++i )
201 if ( xi( 1, i ) == 1.0 )
202 M_eta( 0, i ) = -1.0;
205 M_eta( 0, i ) = 2.0 * ( 1.0 + xi( 0, i ) ) / ( 1.0 - xi( 1, i ) ) - 1.0;
207 M_eta( 1, i ) = xi( 1, i );
212 matrix_node_type
const& operator()( matrix_node_type
const& xi )
214 M_eta.resize( xi.size1(), xi.size2() );
216 for (
size_type i = 0; i < xi.size2(); ++i )
218 if ( xi( 1, i ) == 1.0 )
219 M_eta( 0, i ) = -1.0;
222 M_eta( 0, i ) = 2.0 * ( 1.0 + xi( 0, i ) ) / ( 1.0 - xi( 1, i ) ) - 1.0;
224 M_eta( 1, i ) = xi( 1, i );
229 matrix_node_type
const& operator()()
const
233 matrix_node_type M_eta;
237 struct etas<TETRAHEDRON, T>
239 typedef T value_type;
240 typedef typename matrix_node<value_type>::type matrix_node_type;
247 etas( matrix_node_type
const& xi )
249 M_eta( xi.size1(), xi.size2() )
251 for (
size_type i = 0; i < xi.size2(); ++i )
253 if ( xi( 1, i ) + xi( 2, i ) == 0. )
257 M_eta( 0, i ) = -2. * ( 1. + xi( 0, i ) ) / ( xi( 1, i ) + xi( 2, i ) ) - 1.;
259 if ( xi( 2, i ) == 1. )
263 M_eta( 1, i ) = 2. * ( 1. + xi( 1, i ) ) / ( 1. - xi( 2, i ) ) - 1.;
265 M_eta( 2, i ) = xi( 2, i );
269 matrix_node_type
const& operator()( matrix_node_type
const& xi )
271 M_eta.resize( xi.size1(), xi.size2() );
273 for (
size_type i = 0; i < xi.size2(); ++i )
275 if ( xi( 1, i ) + xi( 2, i ) == 0. )
279 M_eta( 0, i ) = -2. * ( 1. + xi( 0, i ) ) / ( xi( 1, i ) + xi( 2, i ) ) - 1.;
281 if ( xi( 2, i ) == 1. )
285 M_eta( 1, i ) = 2. * ( 1. + xi( 1, i ) ) / ( 1. - xi( 2, i ) ) - 1.;
287 M_eta( 2, i ) = xi( 2, i );
292 matrix_node_type
const& operator()()
const
296 matrix_node_type M_eta;
315 typedef T value_type;
328 psitilde(
int i,
int j )
334 psitilde(
int i,
int j,
int k )
336 p( k, 2*( i+j+1 ), 0.0 ),
345 value_type a( value_type
const& x )
const
354 value_type b( value_type
const& x )
const
356 return math::pow( 0.5 *( 1.0-x ), M_b )*p( x );
363 value_type c( value_type
const& x )
const
365 return math::pow( 0.5 *( 1.0-x ), M_b )*p( x );
377 template<u
int16_type N,
typename T =
double>
380 typedef T value_type;
381 typedef ublas::matrix<value_type> matrix_type;
382 typedef typename matrix_node<value_type>::type matrix_node_type;
389 scalings( ublas::vector<value_type>
const& pts )
391 M_s( N+1, pts.size() )
394 ublas::vector<value_type> one( ublas::scalar_vector<value_type>( pts.size(), 1.0 ) );
395 ublas::row( M_s, 0 ) = one;
397 ublas::row( M_s, 0 ) = ublas::scalar_vector<value_type>( pts.size(), 1.0 );
402 ublas::row( M_s, 1 ) = 0.5 * ( ublas::row( M_s, 0 ) - pts );
404 for ( uint16_type k = 2; k < N+1; ++k )
406 ublas::row( M_s, k ) = ublas::element_prod( ublas::row( M_s, k-1 ),
407 ublas::row( M_s, 1 ) );
411 matrix_type
const& operator()()
const
1D Jacobi polynomial
Definition: jacobi.hpp:270
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319