30 #define __principal_H 1
37 #include <feel/feelpoly/policy.hpp>
57 template< uint16_type Degree,
typename T = double,
58 template<
class>
class StoragePolicy = StorageUBlas>
63 static const uint16_type nOrder = Degree;
79 typedef StoragePolicy<value_type> storage_policy;
80 typedef typename storage_policy::matrix_type matrix_type;
81 typedef typename storage_policy::vector_matrix_type vector_matrix_type;
82 typedef typename storage_policy::matrix_node_type matrix_node_type;
83 typedef typename storage_policy::node_type node_type;
84 typedef typename storage_policy::vector_vector_matrix_type vector_vector_matrix_type;
85 typedef typename storage_policy::vector_type vector_type;
94 : M_a( value_type( 1.0 ) ), M_b( value_type( 1.0 ) )
123 uint16_type degree()
const
151 matrix_type
evaluate_1( vector_type
const& __pts )
const;
155 vector_matrix_type
evaluate_2( vector_type
const& __pts )
const;
160 vector_vector_matrix_type
evaluate_3( vector_type
const& __pts )
const;
164 matrix_type derivate_1( vector_type
const& __pts )
const;
166 vector_matrix_type derivate_2( vector_type
const& __pts )
const;
168 vector_vector_matrix_type derivate_3( vector_type
const& __pts )
const;
178 template<uint16_type Degree,
180 template<
class>
class StoragePolicy>
181 typename Principal<Degree, T, StoragePolicy>::matrix_type
185 matrix_type J ( JacobiBatchEvaluation<nOrder-1,value_type>( M_a, M_b, __pts ) );
186 matrix_type D( nOrder+1,__pts.size() );
188 vector_type ones( ublas::scalar_vector<value_type>( __pts.size(), value_type( 1.0 ) ) );
190 ublas::row( D,0 ) = value_type( 0.5 )*( ones - __pts );
191 ublas::row( D,nOrder ) = value_type( 0.5 )*( ones + __pts );
193 vector_type tmp( ublas::element_prod( ublas::row( D,0 ),ublas::row( D,D.size1()-1 ) ) );
195 for ( uint16_type i = 1; i < D.size1()-1; ++i )
196 ublas::row( D, i ) = ublas::element_prod( tmp, ublas::row( J,i-1 ) );
202 template<uint16_type Degree,
204 template<
class>
class StoragePolicy>
205 typename Principal<Degree, T, StoragePolicy>::vector_matrix_type
209 matrix_type psi_1( evaluate_1( __pts ) );
211 vector_matrix_type m( nOrder+1 );
212 vector_matrix_type J( nOrder-1 );
214 m[0].resize( nOrder+1 ,__pts.size() );
218 vector_type ones( ublas::scalar_vector<value_type>( __pts.size(), value_type( 1.0 ) ) );
220 vector_type tmp1 = ( ones - __pts )/value_type( 2.0 );
221 vector_type tmp2 = ( ones + __pts )/value_type( 2.0 );
222 vector_type tmp_i = tmp1;
224 for ( uint16_type i= 1; i <= J.size() ; ++i )
226 m[i].resize( nOrder ,__pts.size() );
227 J[i-1] = JacobiBatchEvaluation<nOrder-1,value_type>( value_type( 2.0*i+1.0 ), M_b, __pts );
229 tmp_i=ublas::element_prod( tmp_i,tmp1 );
231 ublas::row( m[i],0 ) = tmp_i;
232 vector_type tmp3( ublas::element_prod( tmp2,tmp_i ) );
234 for ( int16_type j=1; j < nOrder; ++j )
236 ublas::row( m[i],j ) = ublas::element_prod( tmp3, ublas::row( J[i-1],j-1 ) );
244 template<uint16_type Degree,
246 template<
class>
class StoragePolicy>
247 typename Principal<Degree, T, StoragePolicy>::vector_vector_matrix_type
251 vector_vector_matrix_type m( nOrder+1 );
253 vector_matrix_type psi_2( evaluate_2( __pts ) );
256 m[0].resize( nOrder+1 );
260 m[nOrder].resize( nOrder+1 );
263 vector_type ones( ublas::scalar_vector<value_type>( __pts.size(), value_type( 1.0 ) ) );
264 vector_type tmp1 = ( ones - __pts )/value_type( 2.0 );
265 vector_type tmp2 = ( ones + __pts )/value_type( 2.0 );
266 vector_type tmp_i = tmp1;
268 vector_vector_matrix_type J( nOrder+1 );
270 for ( int16_type i=1; i < nOrder; ++i )
272 m[i].resize( nOrder+1 );
273 J[i].resize( nOrder+1 );
276 m[i][0].resize( psi_2[i].size1(), psi_2[i].size2() );
280 m[i][nOrder].resize( psi_2[i].size1(),psi_2[i].size2() );
281 m[i][nOrder] = psi_2[i];
283 tmp_i=ublas::element_prod( tmp_i,tmp1 );
285 vector_type tmp_i_j( tmp_i );
287 for ( int16_type j=1; j < nOrder; ++j )
290 m[i][j].resize( nOrder,__pts.size() );
291 J[i][j] = JacobiBatchEvaluation<nOrder-1,value_type>( ( 2.0*i+2.0*j+1.0 ), 1.0, __pts );
292 tmp_i_j = ublas::element_prod( tmp_i_j,tmp1 );
294 ublas::row( m[i][j],0 ) = tmp_i_j;
296 for ( int16_type k=1; k < nOrder; ++k )
298 vector_type tmp3( ublas::element_prod( tmp2,tmp_i_j ) );
299 ublas::row( m[i][j],k ) = ublas::element_prod( tmp3, ublas::row( J[i][j],k-1 ) );
308 template<uint16_type Degree,
310 template<
class>
class StoragePolicy>
311 typename Principal<Degree, T, StoragePolicy>::matrix_type
315 matrix_type D( nOrder+1, __pts.
size() );
316 matrix_type J ( JacobiBatchDerivation<nOrder-1,value_type>( M_a, M_b, __pts ) );
317 vector_type demi( ublas::scalar_vector<value_type>( __pts.
size(), value_type( 0.5 ) ) );
319 ublas::row( D,0 ) = -demi;
320 ublas::row( D,nOrder ) = demi;
322 vector_type ones( ublas::scalar_vector<value_type>( __pts.
size(), value_type( 1.0 ) ) );
325 matrix_type E ( JacobiBatchEvaluation<nOrder-1,value_type>( M_a, M_b, __pts ) );
327 vector_type tmp( ublas::element_prod( 0.5*( ones - __pts ),0.5*( ones + __pts ) ) );
329 for ( uint16_type i = 1; i < D.size1()-1; ++i )
330 ublas::row( D, i ) = ublas::element_prod( tmp, ublas::row( J,i-1 ) ) - 0.5*ublas::element_prod( __pts, ublas::row( E,i-1 ) );
335 template<uint16_type Degree,
337 template<
class>
class StoragePolicy>
338 typename Principal<Degree, T, StoragePolicy>::vector_matrix_type
339 Principal<Degree, T, StoragePolicy>::derivate_2( vector_type
const& __pts )
const
342 matrix_type dpsi_1( derivate_1( __pts ) );
344 vector_matrix_type m( nOrder+1 );
345 vector_matrix_type J( nOrder-1 );
346 vector_matrix_type dJ( nOrder-1 );
348 m[0].resize( nOrder+1 ,__pts.size() );
352 vector_type ones( ublas::scalar_vector<value_type>( __pts.size(), value_type( 1.0 ) ) );
354 vector_type tmp1 = ( ones - __pts )/value_type( 2.0 );
355 vector_type tmp2 = ( ones + __pts )/value_type( 2.0 );
357 vector_type tmp_prod = ublas::element_prod( tmp1,tmp2 );
359 vector_type tmp_i = tmp1;
361 for ( uint16_type i= 1; i <= J.size() ; ++i )
363 m[i].resize( nOrder ,__pts.size() );
364 J[i-1] = JacobiBatchEvaluation<nOrder-1,value_type>( ( 2.0*i+1.0 ), M_b, __pts );
365 dJ[i-1] = JacobiBatchDerivation<nOrder-1,value_type>( ( 2.0*i+1.0 ), M_b, __pts );
367 ublas::row( m[i],0 ) = ( - value_type( i ) - 1.0 ) * tmp_i / 2.0;
369 for ( int16_type j=1; j < nOrder; ++j )
373 ublas::element_prod( tmp1, ublas::row( J[i-1],j-1 ) )
374 + 2.0*element_prod( tmp_prod, ublas::row( dJ[i-1],j-1 ) )
375 - ( value_type( i ) + 1.0 ) * ublas::element_prod( tmp2,ublas::row( J[i-1],j-1 ) )
377 ublas::row( m[i],j ) = ( ublas::element_prod( tmp_i, tmp3 ) ) / value_type( 2.0 ) ;
379 vector_type A( ublas::element_prod( ublas::row( m[i],0 ), tmp2 ) );
380 A = ublas::element_prod( A, ublas::row( J[i-1],j-1 ) );
382 vector_type B( ublas::element_prod( 0.5*tmp_i , ublas::row( J[i-1],j-1 ) ) );
384 vector_type C( ublas::element_prod( tmp_i , tmp_prod ) );
385 C = ublas::element_prod( C ,ublas::row( dJ[i-1],j-1 ) );
387 ublas::row( m[i],j ) = A + B + C;
391 tmp_i=ublas::element_prod( tmp_i,tmp1 );
398 template<uint16_type Degree,
400 template<
class>
class StoragePolicy>
401 typename Principal<Degree, T, StoragePolicy>::vector_vector_matrix_type
402 Principal<Degree, T, StoragePolicy>::derivate_3( vector_type
const& __pts )
const
405 vector_vector_matrix_type m( nOrder+1 );
407 vector_matrix_type dpsi_2( derivate_2( __pts ) );
410 m[0].resize( nOrder+1 );
414 m[nOrder].resize( nOrder+1 );
417 vector_type ones( ublas::scalar_vector<value_type>( __pts.size(), value_type( 1.0 ) ) );
418 vector_type tmp1 = ( ones - __pts )/value_type( 2.0 );
419 vector_type tmp2 = ( ones + __pts )/value_type( 2.0 );
420 vector_type tmp_i = ones;
421 vector_type tmp_prod = ublas::element_prod( tmp1,tmp2 );
423 vector_vector_matrix_type J( nOrder+1 );
424 vector_vector_matrix_type dJ( nOrder+1 );
426 for ( int16_type i=1; i < nOrder; ++i )
428 m[i].resize( nOrder+1 );
429 J[i].resize( nOrder+1 );
430 dJ[i].resize( nOrder+1 );
433 m[i][0].resize( dpsi_2[i].size1(), dpsi_2[i].size2() );
437 m[i][nOrder].resize( dpsi_2[i].size1(),dpsi_2[i].size2() );
438 m[i][nOrder] = dpsi_2[i];
440 tmp_i=ublas::element_prod( tmp_i,tmp1 );
442 vector_type tmp_i_j( tmp_i );
444 for ( int16_type j=1; j < nOrder; ++j )
447 m[i][j].resize( nOrder,__pts.size() );
448 J[i][j] = JacobiBatchEvaluation<nOrder-1,value_type>( ( 2.0*i+2.0*j+1.0 ), 1.0, __pts );
449 dJ[i][j] = JacobiBatchDerivation<nOrder-1,value_type>( ( 2.0*i+2.0*j+1.0 ), 1.0, __pts );
451 tmp_i_j=ublas::element_prod( tmp_i_j,tmp1 );
453 ublas::row( m[i][j],0 ) = -value_type( i+j+1.0 )*tmp_i_j / value_type( 2.0 );
455 for ( int16_type k=1; k < nOrder; ++k )
457 vector_type tmp3( ublas::element_prod( tmp1, ublas::row( J[i][j],k-1 ) )
458 + 2.0*element_prod( tmp_prod, ublas::row( dJ[i][j],k-1 ) )
459 - ( value_type( i+j+1.0 ) ) * ublas::element_prod( tmp2,ublas::row( J[i][j],k-1 ) ) );
461 ublas::row( m[i][j],k ) = ( ublas::element_prod( tmp_i_j, tmp3 ) ) / 2.0 ;