31 #define __Legendre_H 1
35 #include <boost/lambda/if.hpp>
41 #include <feel/feelpoly/policy.hpp>
48 template<
class Convex,
51 class PointSetGaussLobatto;
53 template<uint16_type Dim,
56 typename NormalizationPolicy,
58 template<
class>
class StoragePolicy>
62 template<uint16_type Dim,
70 static const uint16_type nDim = Dim;
71 static const uint16_type nRealDim = RealDim;
72 static const uint16_type nOrder = Degree;
73 static const uint16_type nConvexOrderDiff = nOrder+2;
74 static const bool is_normalized = NormalizationPolicy::is_normalized;
85 template<u
int16_type order,
typename V = value_type>
88 typedef Hypercube<nDim, order, nDim> type;
89 typedef Reference<Hypercube<nDim, order, nDim>, nDim, order, nDim, V> reference_type;
92 template<
typename NewT>
93 struct ChangeValueType
96 typedef LegendreTraits<Dim, RealDim, Degree, NormalizationPolicy, NewT, StoragePolicy> traits_type;
99 template<u
int16_type NewOrder>
102 typedef Legendre<Dim, RealDim, NewOrder, NormalizationPolicy, T, StoragePolicy> type;
103 typedef LegendreTraits<Dim, RealDim, NewOrder, NormalizationPolicy, T, StoragePolicy> traits_type;
109 typedef typename Convex<nOrder>::type convex_type;
110 typedef typename Convex<nOrder>::reference_type reference_convex_type;
112 typedef typename Convex<nConvexOrderDiff>::type diff_convex_type;
113 typedef typename Convex<nConvexOrderDiff>::reference_type diff_reference_convex_type;
115 typedef PointSetGaussLobatto<diff_convex_type, nConvexOrderDiff, value_type> diff_pointset_type;
120 typedef StoragePolicy<value_type> storage_policy;
121 typedef typename storage_policy::matrix_type matrix_type;
122 typedef typename storage_policy::vector_matrix_type vector_matrix_type;
123 typedef typename storage_policy::matrix_node_type matrix_node_type;
124 typedef typename storage_policy::points_type points_type;
125 typedef typename storage_policy::node_type node_type;
128 template<
int D,
int O>
131 static const int Dim = D;
132 static const int Order = O;
155 template<uint16_type Dim,
158 typename NormalizationPolicy = Normalized<true>,
160 template<
class>
class StoragePolicy = StorageUBlas>
165 typedef LegendreTraits<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy> traits_type;
166 static const uint16_type nDim = Dim;
167 static const uint16_type nRealDim = RealDim;
168 static const uint16_type nOrder = Degree;
169 static const uint16_type nConvexOrder = nOrder+2;
170 static const bool is_normalized = NormalizationPolicy::is_normalized;
171 static const bool isTransformationEquivalent =
true;
172 static const bool isContinuous =
false;
173 static const bool is_product =
true;
174 typedef Discontinuous continuity_type;
180 typedef Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy> self_type;
186 typedef self_type basis_type;
188 typedef typename traits_type::value_type value_type;
193 typedef typename traits_type::convex_type convex_type;
194 typedef typename traits_type::reference_convex_type reference_convex_type;
196 typedef typename traits_type::diff_pointset_type diff_pointset_type;
201 typedef typename traits_type::storage_policy storage_policy;
202 typedef typename traits_type::matrix_type matrix_type;
203 typedef typename traits_type::vector_matrix_type vector_matrix_type;
204 typedef typename traits_type::matrix_node_type matrix_node_type;
205 typedef typename traits_type::points_type points_type;
206 typedef typename traits_type::node_type node_type;
218 M_pts( M_refconvex.makePoints( nDim, 0 ) )
220 this->initDerivation();
222 Legendre( Legendre
const &
d )
239 self_type
const& operator=( self_type
const&
d )
251 matrix_type operator()( node_type
const& pt )
const
253 points_type pts( pt.size(), 1 );
254 ublas::column( pts, 0 ) = pt;
258 matrix_type operator()( points_type
const& pts )
const
274 return convex_type::polyDims( nOrder );
300 return is_normalized;
336 std::cout <<
"[Legendre::coeff] coeff = "
337 << ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() )
340 return ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() );
349 static matrix_type
evaluate( points_type
const& __pts )
351 return evaluate( __pts, mpl::int_<nDim>() );
354 template<
typename AE>
355 static vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts )
357 return derivate( __pts, mpl::int_<nDim>() );
366 static matrix_type
const&
d( uint16_type i )
377 static matrix_type
const&
derivate( uint16_type i )
388 static value_type normalization(
int i )
390 return ( is_normalized?math::sqrt( value_type( i ) + 0.5 ) : value_type( 1 ) );
392 static value_type normalization(
int i,
int j )
394 return ( is_normalized?math::sqrt( ( value_type( i ) + 0.5 ) *
395 ( value_type( j ) + 0.5 ) ) : value_type( 1 ) );
397 static value_type normalization(
int i,
int j,
int k )
399 return ( is_normalized?math::sqrt( ( value_type( i ) + 0.5 ) *
400 ( value_type( j ) + 0.5 ) *
401 ( value_type( k ) + 0.5 ) ) : value_type( 1 ) );
408 evaluate( points_type
const& __pts, mpl::int_<1> )
410 matrix_type m ( JacobiBatchEvaluation<nOrder,value_type>( 0.0, 0.0, ublas::row( __pts, 0 ) ) );
414 for ( uint16_type i = 0; i < m.size1(); ++i )
415 ublas::row( m, i ) *= normalization( i );
425 template<
typename AE>
426 static vector_matrix_type
427 derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<1> )
429 FEELPP_ASSERT( __pts().size1() == 1 )( __pts().size1() )( __pts().size2() ).error(
"invalid points" );
432 vector_matrix_type D( 1 );
433 D[0].resize( nOrder+1, __pts().size2() );
434 D[0] = JacobiBatchDerivation<nOrder,value_type>( 0.0, 0.0, ublas::row( __pts(),0 ) );
437 for ( uint16_type i = 0; i < nOrder+1; ++i )
438 ublas::row( D[0], i ) *= normalization( i );
447 static matrix_type
evaluate( points_type
const& __pts, mpl::int_<2> );
453 template<
typename AE>
454 static vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<2> );
460 static matrix_type
evaluate( points_type
const& __pts, mpl::int_<3> );
466 template<
typename AE>
467 static vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<3> );
469 static void initDerivation();
471 reference_convex_type M_refconvex;
477 static bool _S_has_derivation;
483 static std::vector<matrix_type> _S_D;
487 template<uint16_type Dim,
490 typename NormalizationPolicy,
492 template<
class>
class StoragePolicy>
493 bool Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::_S_has_derivation =
false;
495 template<uint16_type Dim,
498 typename NormalizationPolicy,
500 template<
class>
class StoragePolicy>
501 std::vector<typename Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::matrix_type>
502 Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::_S_D;
504 template<uint16_type Dim,
507 typename NormalizationPolicy,
509 template<
class>
class StoragePolicy>
511 Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::initDerivation()
514 typedef typename traits_type::convex_type convex_type;
515 typedef typename traits_type::reference_convex_type reference_convex_type;
517 typedef typename traits_type::diff_pointset_type diff_pointset_type;
519 typedef typename traits_type::storage_policy storage_policy;
520 typedef typename traits_type::matrix_type matrix_type;
521 typedef typename traits_type::vector_matrix_type vector_matrix_type;
522 typedef typename traits_type::matrix_node_type matrix_node_type;
523 typedef typename traits_type::points_type points_type;
524 typedef typename traits_type::node_type node_type;
527 if ( _S_has_derivation ==
false )
529 _S_has_derivation =
true;
531 reference_convex_type refconvex;
534 diff_pointset_type diff_pts( 1 );
535 matrix_type A( evaluate( diff_pts.points() ) );
538 matrix_type D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
539 LU<matrix_type> lu( A );
540 matrix_type C = lu.solve( D );
542 vector_matrix_type d ( derivate( diff_pts.points() ) );
543 _S_D.resize( d.size() );
545 for (
size_type i = 0; i < d.size(); ++i )
547 _S_D[i] = ublas::prod( d[i], C );
548 glas::clean( _S_D[i] );
555 template<uint16_type Dim,
558 typename NormalizationPolicy,
560 template<
class>
class StoragePolicy>
561 typename Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::matrix_type
564 matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
566 ublas::vector<value_type> eta1s = ublas::row( __pts, 0 );
567 ublas::vector<value_type> eta2s = ublas::row( __pts, 1 );
569 matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
570 matrix_type bs( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta2s ) );
572 for ( uint16_type cur = 0, i = 0; i < nOrder+1; ++i )
574 for ( uint16_type j = 0; j < nOrder+1; ++j,++cur )
576 ublas::row( res, cur ) = normalization( i, j ) * ublas::element_prod( ublas::row( as, i ),
577 ublas::row( bs, j ) );
584 template<uint16_type Dim,
587 typename NormalizationPolicy,
589 template<
class>
class StoragePolicy>
590 template<
typename AE>
591 typename Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::vector_matrix_type
592 Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<2> )
594 vector_matrix_type res( 2 );
595 res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
596 res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
599 matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 0 ) ) );
600 matrix_type das( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 0 ) ) );
601 matrix_type bs( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 1 ) ) );
602 matrix_type dbs( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 1 ) ) );
604 for ( uint16_type cur = 0, i = 0; i < nOrder+1; ++i )
606 for ( uint16_type j = 0; j < nOrder+1; ++j,++cur )
608 ublas::row( res[0], cur ) = normalization( i, j ) * ublas::element_prod( ublas::row( das, i ),
609 ublas::row( bs, j ) );
610 ublas::row( res[1], cur ) = normalization( i, j ) * ublas::element_prod( ublas::row( as, i ),
611 ublas::row( dbs, j ) );
618 template<uint16_type Dim,
621 typename NormalizationPolicy,
623 template<
class>
class StoragePolicy>
624 typename Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::matrix_type
627 matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
629 ublas::vector<value_type> eta1s = ublas::row( __pts, 0 );
630 ublas::vector<value_type> eta2s = ublas::row( __pts, 1 );
631 ublas::vector<value_type> eta3s = ublas::row( __pts, 2 );
633 matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
634 matrix_type bs( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta2s ) );
635 matrix_type cs( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta3s ) );
637 for ( uint16_type cur = 0, i = 0; i < nOrder+1; ++i )
639 for ( uint16_type j = 0; j < nOrder+1; ++j )
641 for ( uint16_type k = 0; k < nOrder+1; ++k,++cur )
643 ublas::row( res, cur ) =
644 normalization( i, j, k )*
645 ublas::element_prod( ublas::element_prod( ublas::row( as, i ),
646 ublas::row( bs, j ) ),
647 ublas::row( cs, k ) );
655 template<uint16_type Dim,
658 typename NormalizationPolicy,
660 template<
class>
class StoragePolicy>
661 template<
typename AE>
662 typename Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::vector_matrix_type
663 Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<3> )
665 vector_matrix_type res( 3 );
666 res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
667 res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
668 res[2].resize( convex_type::polyDims( nOrder ), __pts().size2() );
671 matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 0 ) ) );
672 matrix_type das( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 0 ) ) );
673 matrix_type bs( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 1 ) ) );
674 matrix_type dbs( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 1 ) ) );
675 matrix_type cs( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 2 ) ) );
676 matrix_type dcs( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 2 ) ) );
678 for ( uint16_type cur = 0, i = 0; i < nOrder+1; ++i )
680 for ( uint16_type j = 0; j < nOrder+1; ++j )
682 for ( uint16_type k = 0; k < nOrder+1; ++k,++cur )
684 ublas::row( res[0], cur ) = ( normalization( i, j, k ) *
685 ublas::element_prod( ublas::element_prod( ublas::row( das, i ),
686 ublas::row( bs, j ) ),
687 ublas::row( cs, k ) ) );
689 ublas::row( res[1], cur ) = ( normalization( i, j, k ) *
690 ublas::element_prod( ublas::element_prod( ublas::row( as, i ),
691 ublas::row( dbs, j ) ),
692 ublas::row( cs, k ) ) );
694 ublas::row( res[2], cur ) = ( normalization( i, j, k ) *
695 ublas::element_prod( ublas::element_prod( ublas::row( as, i ),
696 ublas::row( bs, j ) ),
697 ublas::row( dcs, k ) ) );