41 #include <feel/feelpoly/policy.hpp>
66 template<uint16_type Dim,
69 typename TheConvex=Simplex<Dim>,
71 template<
class>
class StoragePolicy = StorageUBlas>
76 static const uint16_type nDim = Dim;
77 static const uint16_type nRealDim = Dim;
78 static const uint16_type nOrder = Degree;
79 static const uint16_type nConvexOrder = mpl::if_<mpl::bool_<TheConvex::is_simplex>,
80 mpl::int_<nDim+nOrder+1>,
81 mpl::int_<nOrder+2> >::type::value;
83 static const uint16_type convex_is_simplex = TheConvex::is_simplex;
84 static const uint16_type convex_is_hypercube = TheConvex::is_hypercube;
85 static const bool is_product =
false;
91 typedef Moment<Dim, Degree,TheConvex,T, StoragePolicy> self_type;
93 typedef self_type basis_type;
100 template<u
int16_type order,
typename V = value_type>
103 typedef typename mpl::if_<mpl::bool_<TheConvex::is_simplex>,
104 mpl::identity<Simplex<nDim, order, nDim> >,
105 mpl::identity<Hypercube<nDim, order, nDim> > >::type::type type;
106 typedef typename mpl::if_<mpl::bool_<TheConvex::is_simplex>,
107 mpl::identity<Reference<Simplex<nDim, order, nDim>, nDim, order, nDim, V > >,
108 mpl::identity<Reference<Hypercube<nDim, order, nDim>, nDim, order, nDim, V > > >::type::type reference_type;
114 typedef TheConvex convex_type;
115 typedef typename Convex<nOrder>::reference_type reference_convex_type;
121 typedef StoragePolicy<value_type> storage_policy;
122 typedef typename storage_policy::matrix_type matrix_type;
123 typedef typename storage_policy::vector_matrix_type vector_matrix_type;
124 typedef typename storage_policy::matrix_node_type matrix_node_type;
125 typedef typename storage_policy::node_type node_type;
126 typedef typename storage_policy::points_type points_type;
137 M_pts( M_refconvex.
points() ),
138 M_coord( convex_type::polyDims( nOrder ),nDim )
142 for ( uint32_type i=0; i < M_coord.size1(); ++i )
145 else if ( nDim == 2 )
149 if ( convex_is_simplex )
151 for ( int32_type n = 0; n < nOrder+1; ++n )
152 for ( int32_type i = n; i >= 0; --i )
154 M_coord( G_i,0 ) = i;
155 M_coord( G_i,1 ) = n-i;
160 else if ( convex_is_hypercube )
162 for ( int32_type n = 0; n < nOrder+1; ++n )
163 for ( int32_type i = 0; i < nOrder+1; ++i )
165 M_coord( G_i,0 ) = i;
166 M_coord( G_i,1 ) = n;
172 else if ( nDim == 3 )
176 if ( convex_is_simplex )
178 for ( int32_type n = 0; n < nOrder+1; ++n )
179 for ( int32_type i = n; i >= 0; --i )
180 for ( int32_type j = n-i; j >= 0 ; --j )
182 M_coord( G_i,0 ) = i;
183 M_coord( G_i,1 ) = j;
184 M_coord( G_i,2 ) = n-i-j;
189 else if ( convex_is_hypercube )
191 for ( int32_type n = 0; n < nOrder+1; ++n )
192 for ( int32_type i = 0; i < nOrder+1; ++i )
193 for ( int32_type j = 0; j < nOrder+1; ++j )
195 M_coord( G_i,0 ) = j;
196 M_coord( G_i,1 ) = i;
197 M_coord( G_i,1 ) = n;
208 matrix_type A( ublas::trans( evaluate( M_pts ) ) );
210 matrix_type D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
212 LU<matrix_type> lu( A );
213 matrix_type C = lu.solve( D );
215 vector_matrix_type d ( derivate( M_pts ) );
216 M_D.resize( d.size() );
218 for (
size_type i = 0; i < d.size(); ++i )
220 M_D[i] = ublas::prod( C, ublas::trans( d[i] ) );
224 Moment( Moment
const & d )
228 M_coord( d.M_coord ),
251 self_type
const&
operator=( self_type
const& d )
263 matrix_type operator()( node_type
const& pt )
const
265 points_type pts( pt.size(), 1 );
266 ublas::column( pts, 0 ) = pt;
267 return evaluate( pts );
270 matrix_type operator()( points_type
const& pts )
const
272 return evaluate( pts );
281 ublas::matrix<int16_type> coord()
290 uint16_type degree()
const
298 self_type
const& basis()
const
316 std::string familyName()
const
343 matrix_type coeff()
const
346 std::cout <<
"[Moment::coeff] coeff = "
347 << ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() )
350 return ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() );
359 matrix_type evaluate( points_type
const& __pts )
const
361 return evaluate( __pts, mpl::int_<nDim>() );
364 template<
typename AE>
365 vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts )
const
367 return derivate( __pts, mpl::int_<nDim>() );
377 matrix_type
const& d( uint16_type i )
const
382 template<
template<u
int16_type>
class PolySetType = Scalar>
383 Polynomial<self_type,PolySetType> pick(
int i,
int c = 0 )
const
385 size_type dim_p = convex_type::polyDims( nDim );
386 int nComponents = PolySetType<Dim>::nComponents;
387 matrix_type coeff( nComponents, dim_p );
388 coeff = ublas::scalar_matrix<value_type>( coeff.size1(), coeff.size2(), 0. );
391 return Polynomial<self_type, PolySetType>( *
this, coeff, true );
403 value_type SP_integrate( ublas::matrix<value_type>
const& __Comp )
405 return SP_integrate( __Comp, mpl::int_<nDim>() );
415 value_type S_integrate( ublas::matrix<value_type>
const& __Comp )
417 return S_integrate( __Comp, mpl::int_<nDim>() );
431 evaluate( points_type
const& __pts, mpl::int_<1> )
const
433 matrix_type m ( nOrder+1 ,__pts.size2() );
434 ublas::vector<value_type> pts( ublas::row( __pts,0 ) );
436 for ( uint32_type i = 0; i < m.size2(); ++i )
439 for ( uint32_type i = 1; i < m.size1(); ++i )
440 ublas::row( m,i ) = ublas::element_prod( ublas::row( m,i-1 ), pts );
450 evaluate( ublas::vector<value_type>
const& __pts )
const
452 matrix_type m ( nOrder+1 ,__pts.size() );
454 for ( uint32_type i = 0; i < m.size2(); ++i )
457 for ( uint32_type i = 1; i < m.size1(); ++i )
458 ublas::row( m,i ) = ublas::element_prod( ublas::row( m,i-1 ), __pts );
466 template<
typename AE>
468 derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<1> )
const
470 FEELPP_ASSERT( __pts().size1() == 1 )( __pts().size1() )( __pts().size2() ).error(
"invalid points" );
471 vector_matrix_type D( 1 );
473 D[0].resize( nOrder+1, __pts().size2() );
474 matrix_type E ( evaluate( __pts ) );
476 for ( uint32_type i = 0; i < E.size2(); ++i )
477 ublas::row( D[0], 0 )( i ) = value_type( 0.0 );
479 for ( uint32_type i = 1; i < E.size1(); ++i )
480 ublas::row( D[0], i ) = value_type( i )*ublas::row( E,i-1 );
489 derivate( ublas::vector<value_type>
const& __pts )
const
491 vector_matrix_type D( 1 );
493 D[0].resize( nOrder+1, __pts.size() );
494 matrix_type E ( evaluate( __pts ) );
496 for ( uint32_type i = 0; i < E.size2(); ++i )
497 ublas::row( D[0], 0 )( i ) = value_type( 0.0 );
499 for ( int32_type i = 1; i < E.size1(); ++i )
500 ublas::row( D[0], i ) = value_type( i )*ublas::row( E,i-1 );
509 matrix_type evaluate( points_type
const& __pts, mpl::int_<2> )
const;
515 template<
typename AE>
516 vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<2> )
const;
522 matrix_type evaluate( points_type
const& __pts, mpl::int_<3> )
const;
528 template<
typename AE>
529 vector_matrix_type derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<3> )
const;
535 value_type SP_integrate( ublas::matrix<value_type>
const& __Comp, mpl::int_<1> )
537 value_type res( 0.0 );
539 if ( __Comp.size1() == 1 )
541 for (
int ki = 0; ki <__Comp.size2(); ki += 2 )
543 res+= 2.0*__Comp( 0,ki ) / value_type( ki+1 );
547 else if ( __Comp.size1() == 2 )
549 for (
int k1 = 0; k1 <__Comp.size2(); ++k1 )
550 for (
int k2 = 0; k2 <__Comp.size2(); ++k2 )
552 if ( ( k1+k2 )%2 == 0 )
553 res+= 2.0 * __Comp( 0,k1 ) * __Comp( 1,k2 ) / value_type( k1+k2+1 );
557 else if ( __Comp.size1() == 3 )
559 for (
int k1 = 0; k1 <__Comp.size2(); ++k1 )
560 for (
int k2 = 0; k2 <__Comp.size2(); ++k2 )
561 for (
int k3 = 0; k3 <__Comp.size2(); ++k3 )
563 if ( ( k1+k2+k3 )%2 == 0 )
564 res+= 2.0 * __Comp( 0,k1 ) * __Comp( 1,k2 )* __Comp( 2,k3 ) / value_type( k1+k2+k3+1 );
576 value_type SP_integrate( ublas::matrix<value_type>
const& __Comp, mpl::int_<2> )
578 value_type res( 0.0 );
580 if ( __Comp.size1() == 1 )
582 for (
int ki = 0; ki <__Comp.size2(); ++ki )
584 int a( this->coord()( ki,0 ) );
585 int b( this->coord()( ki,1 ) );
587 if ( ( a%2 == 0 ) && ( b%2 == 0 ) )
588 res+= 4.0*__Comp( 0,ki ) / value_type( a+1 ) / value_type( b+1 );
592 else if ( __Comp.size1() == 2 )
594 for (
int k1 = 0; k1 <__Comp.size2(); ++k1 )
595 for (
int k2 = 0; k2 <__Comp.size2(); ++k2 )
597 int a( this->coord()( k1,0 ) + this->coord()( k2,0 ) );
598 int b( this->coord()( k2,1 ) + this->coord()( k2,1 ) );
600 if ( ( a%2 == 0 ) && ( b%2 == 0 ) )
601 res+= 4.0 * __Comp( 0,k1 ) * __Comp( 1,k2 ) / value_type( a+1 ) / value_type( b+1 );
605 else if ( __Comp.size1() == 3 )
607 for (
int k1 = 0; k1 <__Comp.size2(); ++k1 )
608 for (
int k2 = 0; k2 <__Comp.size2(); ++k2 )
609 for (
int k3 = 0; k3 <__Comp.size2(); ++k3 )
611 int a( this->coord()( k1,0 ) + this->coord()( k2,0 ) + this->coord()( k3,0 ) );
612 int b( this->coord()( k1,1 ) + this->coord()( k2,1 ) + this->coord()( k3,1 ) );
614 if ( ( a%2 == 0 ) && ( b%2 == 0 ) )
615 res+= 4.0 * __Comp( 0,k1 ) * __Comp( 1,k2 ) * __Comp( 2,k3 ) / value_type( a+1 ) / value_type( b+1 );
625 value_type SP_integrate( ublas::matrix<value_type>
const& __Comp, mpl::int_<3> )
627 value_type res( 0.0 );
629 if ( __Comp.size1() == 1 )
631 for (
int ki = 0; ki <__Comp.size2(); ++ki )
633 int a( this->coord()( ki,0 ) );
634 int b( this->coord()( ki,1 ) );
635 int c( this->coord()( ki,2 ) );
637 if ( ( a%2 == 0 ) && ( b%2 == 0 ) && ( c%2 == 0 ) )
638 res+= 8.0*__Comp( 0,ki ) / value_type( a+1 ) / value_type( b+1 ) / value_type( c+1 );
642 else if ( __Comp.size1() == 2 )
644 for (
int k1 = 0; k1 <__Comp.size2(); ++k1 )
645 for (
int k2 = 0; k2 <__Comp.size2(); ++k2 )
647 int a( this->coord()( k1,0 ) + this->coord()( k2,0 ) );
648 int b( this->coord()( k1,1 ) + this->coord()( k2,1 ) );
649 int c( this->coord()( k1,2 ) + this->coord()( k2,2 ) );
651 if ( ( a%2 == 0 ) && ( b%2 == 0 ) && ( c%2 == 0 ) )
652 res+= 8.0*__Comp( 0,k1 )* __Comp( 1,k2 ) / value_type( a+1 ) / value_type( b+1 ) / value_type( c+1 );
656 else if ( __Comp.size1() == 3 )
658 for (
int k1 = 0; k1 <__Comp.size2(); ++k1 )
659 for (
int k2 = 0; k2 <__Comp.size2(); ++k2 )
660 for (
int k3 = 0; k3 <__Comp.size2(); ++k3 )
662 int a( this->coord()( k1,0 ) + this->coord()( k2,0 ) + this->coord()( k3,0 ) );
663 int b( this->coord()( k1,1 ) + this->coord()( k2,1 ) + this->coord()( k3,1 ) );
664 int c( this->coord()( k1,2 ) + this->coord()( k2,2 ) + this->coord()( k3,2 ) );
666 if ( ( a%2 == 0 ) && ( b%2 == 0 ) && ( c%2 == 0 ) )
667 res+= 8.0*__Comp( 0,k1 )* __Comp( 1,k2 )* __Comp( 2,k3 )/value_type( a+1 )/value_type( b+1 )/value_type( c+1 );
678 value_type S_integrate( ublas::matrix<value_type>
const& __Comp, mpl::int_<1> )
680 return SP_integrate( __Comp, mpl::int_<nDim>() );
687 value_type S_integrate( ublas::matrix<value_type>
const& __Comp, mpl::int_<2> )
689 value_type res( 0.0 );
691 if ( __Comp.size1() == 1 )
693 for ( uint32_type ki = 0; ki <__Comp.size2(); ++ki )
695 int i( this->coord()( ki,0 ) );
696 int j( this->coord()( ki,1 ) );
703 if ( ( i%2 == 0 ) || ( ( i+j )%2 != 0 ) )
705 if ( pi==0 && pj==0 )
706 res+= __Comp( 0,ki )/value_type( j+1 ) * ( 2.0/value_type( i+1 ) ) ;
708 else if ( pi==0 && pj==1 )
709 res+= __Comp( 0,ki )/value_type( j+1 ) * ( 2.0/value_type( i+j+2 )-2.0/value_type( i+1 ) ) ;
711 else if ( pi==1 && pj==0 )
712 res-= __Comp( 0,ki )/value_type( j+1 ) * ( 2.0/value_type( i+j+2 ) ) ;
717 else if ( __Comp.size1() == 2 )
719 for ( uint32_type k1 = 0; k1 <__Comp.size2(); ++k1 )
720 for ( uint32_type k2 = 0; k2 <__Comp.size2(); ++k2 )
722 int i( this->coord()( k1,0 ) + this->coord()( k2,0 ) );
723 int j( this->coord()( k1,1 ) + this->coord()( k2,1 ) );
730 if ( ( i%2 == 0 ) || ( ( i+j )%2 != 0 ) )
732 if ( pi==0 && pj==0 )
733 res+= __Comp( 0,k1 ) * __Comp( 1,k2 )/value_type( j+1 ) * ( 2.0/value_type( i+1 ) ) ;
735 else if ( pi==0 && pj==1 )
736 res+= __Comp( 0,k1 ) * __Comp( 1,k2 )/value_type( j+1 ) * ( 2.0/value_type( i+j+2 )-2.0/value_type( i+1 ) ) ;
738 else if ( pi==1 && pj==0 )
739 res-= __Comp( 0,k1 ) * __Comp( 1,k2 )/value_type( j+1 ) * ( 2.0/value_type( i+j+2 ) ) ;
744 else if ( __Comp.size1() == 3 )
746 for ( uint32_type k1 = 0; k1 <__Comp.size2(); ++k1 )
747 for ( uint32_type k2 = 0; k2 <__Comp.size2(); ++k2 )
748 for ( uint32_type k3 = 0; k3 <__Comp.size2(); ++k3 )
750 int i( this->coord()( k1,0 ) + this->coord()( k2,0 )+ this->coord()( k3,0 ) );
751 int j( this->coord()( k1,1 ) + this->coord()( k2,1 )+ this->coord()( k3,1 ) );
758 if ( ( i%2 == 0 ) || ( ( i+j )%2 != 0 ) )
760 if ( pi==0 && pj==0 )
761 res+= __Comp( 0,k1 ) * __Comp( 1,k2 ) * __Comp( 2,k3 )/value_type( j+1 ) * ( 2.0/value_type( i+1 ) ) ;
763 else if ( pi==0 && pj==1 )
764 res+= __Comp( 0,k1 ) * __Comp( 1,k2 ) * __Comp( 2,k3 )/value_type( j+1 ) * ( 2.0/value_type( i+j+2 )-2.0/value_type( i+1 ) ) ;
766 else if ( pi==1 && pj==0 )
767 res-= __Comp( 0,k1 ) * __Comp( 1,k2 ) * __Comp( 2,k3 )/value_type( j+1 ) * ( 2.0/value_type( i+j+2 ) ) ;
779 value_type S_integrate( ublas::matrix<value_type>
const& __Comp, mpl::int_<3> )
787 reference_convex_type M_refconvex;
795 ublas::matrix<int16_type> M_coord;
800 std::vector<matrix_type> M_D;
803 template<uint16_type Dim,
807 template<
class>
class StoragePolicy>
808 typename Moment<Dim, Degree,TheConvex, T, StoragePolicy>::matrix_type
809 Moment<Dim, Degree,TheConvex, T, StoragePolicy>::evaluate( points_type
const& __pts, mpl::int_<2> )
const
811 matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
813 ublas::vector<value_type> pts_x = ublas::row( __pts, 0 );
814 ublas::vector<value_type> pts_y = ublas::row( __pts, 1 );
816 matrix_type E_x( evaluate( pts_x ) );
817 matrix_type E_y( evaluate( pts_y ) );
825 if ( convex_is_hypercube )
827 for ( int32_type n = 0; n < nOrder+1; ++n )
828 for ( int32_type i = 0; i < nOrder+1; ++i )
831 ublas::row( res, G_i )=ublas::element_prod( ublas::row( E_x, i ) , ublas::row( E_y, n ) );
835 else if ( convex_is_simplex )
837 for ( int32_type n = 0; n < nOrder+1; ++n )
838 for ( int32_type i = n; i >= 0; --i )
840 ublas::row( res, G_i )=ublas::element_prod( ublas::row( E_x, i ) , ublas::row( E_y, n-i ) );
848 template<uint16_type Dim,
852 template<
class>
class StoragePolicy>
853 template<
typename AE>
854 typename Moment<Dim, Degree, TheConvex, T, StoragePolicy>::vector_matrix_type
855 Moment<Dim, Degree, TheConvex, T, StoragePolicy>::derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<2> )
const
857 vector_matrix_type res( 2 );
858 res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
859 res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
861 ublas::vector<value_type> pts_x = ublas::row( __pts(), 0 );
862 ublas::vector<value_type> pts_y = ublas::row( __pts(), 1 );
864 matrix_type E_x( evaluate( pts_x ) );
865 matrix_type E_y( evaluate( pts_y ) );
867 vector_matrix_type D_x( derivate( pts_x ) );
868 vector_matrix_type D_y( derivate( pts_y ) );
872 if ( convex_is_simplex )
874 for ( int32_type n = 0; n < nOrder+1; ++n )
875 for ( int32_type i = n; i >= 0; --i )
877 ublas::row( res[0], G_i )=ublas::element_prod( ublas::row( D_x[0], i ) , ublas::row( E_y, n-i ) );
878 ublas::row( res[1], G_i )=ublas::element_prod( ublas::row( E_x, i ) , ublas::row( D_y[0], n-i ) );
882 else if ( convex_is_hypercube )
884 for ( int32_type n = 0; n < nOrder+1; ++n )
885 for ( int32_type i = 0; i < nOrder+1; ++i )
887 ublas::row( res[0], G_i )=ublas::element_prod( ublas::row( D_x[0], i ) , ublas::row( E_y, n ) );
888 ublas::row( res[1], G_i )=ublas::element_prod( ublas::row( E_x, i ) , ublas::row( D_y[0], n ) );
896 template<uint16_type Dim,
900 template<
class>
class StoragePolicy>
901 typename Moment<Dim, Degree, TheConvex, T, StoragePolicy>::matrix_type
902 Moment<Dim, Degree, TheConvex, T, StoragePolicy>::evaluate( points_type
const& __pts, mpl::int_<3> )
const
904 matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
906 FEELPP_ASSERT( __pts.size1() == 3 )( __pts.size1() ).error(
"invalid space dimension" );
908 ublas::vector<value_type> pts_x = ublas::row( __pts, 0 );
909 ublas::vector<value_type> pts_y = ublas::row( __pts, 1 );
910 ublas::vector<value_type> pts_z = ublas::row( __pts, 2 );
912 matrix_type E_x( evaluate( pts_x ) );
913 matrix_type E_y( evaluate( pts_y ) );
914 matrix_type E_z( evaluate( pts_z ) );
918 for ( int32_type n = 0; n < nOrder+1; ++n )
919 for ( int32_type i = n; i >= 0; --i )
920 for ( int32_type j = n-i; j >= 0 ; --j )
922 ublas::vector<value_type> tmp( ublas::element_prod( ublas::row( E_x, i ) , ublas::row( E_y, j ) ) );
923 ublas::row( res, G_i )=ublas::element_prod( tmp , ublas::row( E_z, n-i-j ) );
930 template<uint16_type Dim,
934 template<
class>
class StoragePolicy>
935 template<
typename AE>
936 typename Moment<Dim, Degree, TheConvex, T, StoragePolicy>::vector_matrix_type
937 Moment<Dim, Degree, TheConvex, T, StoragePolicy>::derivate( ublas::matrix_expression<AE>
const& __pts, mpl::int_<3> )
const
939 vector_matrix_type res( 3 );
940 res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
941 res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
942 res[2].resize( convex_type::polyDims( nOrder ), __pts().size2() );
944 FEELPP_ASSERT( __pts().size1() == 3 )( __pts().size1() ).error(
"invalid space dimension" );
946 ublas::vector<value_type> pts_x = ublas::row( __pts(), 0 );
947 ublas::vector<value_type> pts_y = ublas::row( __pts(), 1 );
948 ublas::vector<value_type> pts_z = ublas::row( __pts(), 2 );
950 matrix_type E_x( evaluate( pts_x ) );
951 matrix_type E_y( evaluate( pts_y ) );
952 matrix_type E_z( evaluate( pts_z ) );
954 vector_matrix_type D_x( derivate( pts_x ) );
955 vector_matrix_type D_y( derivate( pts_y ) );
956 vector_matrix_type D_z( derivate( pts_z ) );
960 for ( int32_type n = 0; n < nOrder+1; ++n )
961 for ( int32_type i = n; i >= 0; --i )
962 for ( int32_type j = n-i; j >=0 ; --j )
964 ublas::vector<value_type> tmp( ublas::element_prod( ublas::row( D_x[0], i ) , ublas::row( E_y, j ) ) );
965 ublas::row( res[0], G_i )=ublas::element_prod( tmp , ublas::row( E_z, n-i-j ) );
967 tmp = ublas::element_prod( ublas::row( E_x, i ) , ublas::row( D_y[0], j ) );
968 ublas::row( res[1], G_i )=ublas::element_prod( tmp , ublas::row( E_z, n-i-j ) );
970 tmp = ublas::element_prod( ublas::row( E_x, i ) , ublas::row( E_y, j ) );
971 ublas::row( res[2], G_i )=ublas::element_prod( tmp , ublas::row( D_z[0], n-i-j ) );
992 template<uint16_type Dim,
995 template<u
int16_type>
class PolySetType = Scalar,
997 template<u
int16_type,u
int16_type,u
int16_type>
class Convex = Simplex>
998 class MomentPolynomialSet
1000 public PolynomialSet<Moment<Dim, Order, Convex<Dim,1,Dim>, T, StorageUBlas>, PolySetType >
1002 typedef PolynomialSet<Moment<Dim, Order, Convex<Dim,1,Dim>, T, StorageUBlas>, PolySetType > super;
1005 static const uint16_type nDim = Dim;
1006 static const uint16_type nOrder = Order;
1007 static const uint16_type nRealDim = RealDim;
1008 static const bool isTransformationEquivalent =
true;
1009 typedef MomentPolynomialSet<Dim, Order,RealDim, PolySetType, T, Simplex> self_type;
1010 typedef self_type component_basis_type;
1012 typedef typename super::polyset_type polyset_type;
1013 static const bool is_tensor2 = polyset_type::is_tensor2;
1014 static const bool is_vectorial = polyset_type::is_vectorial;
1015 static const bool is_scalar = polyset_type::is_scalar;
1016 static const bool is_continuous =
false;
1017 static const bool is_modal =
true;
1018 static const uint16_type nComponents = polyset_type::nComponents;
1019 static const bool is_product =
true;
1020 static const bool isContinuous =
false;
1021 typedef Discontinuous continuity_type;
1023 typedef typename super::component_type component_type;
1025 typedef T value_type;
1026 typedef Moment<Dim, Order, Convex<Dim,1,Dim>, T, StorageUBlas> basis_type;
1027 typedef Convex<Dim, Order, Dim> convex_type;
1031 typedef Convex<Dim, O, Dim> type;
1033 typedef Reference<convex_type, nDim, nOrder, nDim, value_type> reference_convex_type;
1035 typedef typename super::polynomial_type polynomial_type;
1038 static const uint16_type nDofPerVertex = convex_type::nbPtsPerVertex;
1040 static const uint16_type nDofPerEdge = convex_type::nbPtsPerEdge;
1042 static const uint16_type nDofPerFace = convex_type::nbPtsPerFace;
1045 static const uint16_type nDofPerVolume = convex_type::nbPtsPerVolume;
1047 static const uint16_type nLocalDof = convex_type::numPoints;
1049 static const uint16_type nDof = nLocalDof;
1050 static const uint16_type nNodes = nDof;
1051 static const uint16_type nDofGrad = super::nDim*nDof;
1052 static const uint16_type nDofHess = super::nDim*super::nDim*nDof;
1054 typedef StorageUBlas<value_type> storage_policy;
1055 typedef typename storage_policy::matrix_type matrix_type;
1056 typedef typename storage_policy::vector_matrix_type vector_matrix_type;
1057 typedef typename storage_policy::matrix_node_type matrix_node_type;
1058 typedef typename storage_policy::points_type points_type;
1059 typedef typename storage_policy::node_type node_type;
1060 MomentPolynomialSet()
1062 super( basis_type() )
1065 ublas::matrix<value_type> m( ublas::identity_matrix<value_type>( nComponents*convex_type::polyDims( nOrder ) ) );
1068 std::cout <<
"[orthonormalpolynomialset] m = " << m <<
"\n";
1070 if ( !( ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) -
1072 std::cout <<
"m1=" << m <<
"\n"
1073 <<
"m2=" << polyset_type::toMatrix( polyset_type::toType( m ) ) <<
"\n"
1074 << ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) - m ) <<
"\n";
1076 FEELPP_ASSERT( ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) -
1077 m ) < 1e-10 )( m ).warn (
"invalid transformation" );
1078 this->setCoefficient( polyset_type::toType( m ),
true );
1081 MomentPolynomialSet<Dim, Order, RealDim, Scalar,T, Simplex > toScalar()
const
1083 return MomentPolynomialSet<Dim, Order, RealDim, Scalar,T, Simplex >();
1089 std::string familyName()
const
1095 points_type
points()
const
1097 return points_type();
1099 points_type
points(
int f )
const
1101 return points_type();
1105 template<uint16_type Dim,
1107 uint16_type RealDim,
1108 template<u
int16_type>
class PolySetType,
1110 template<u
int16_type,u
int16_type,u
int16_type>
class Convex>
1111 const uint16_type MomentPolynomialSet<Dim, Order, RealDim, PolySetType,T, Convex>::nLocalDof;
1116 template<uint16_type Order,
1117 template<u
int16_type Dim>
class PolySetType = Scalar>
1118 class MomentPolynomialSet
1121 template<uint16_type N,
1122 typename T = double,
1123 typename Convex = Simplex<N> >
1126 typedef typename mpl::if_<mpl::bool_<Convex::is_simplex>,
1127 mpl::identity<detail::MomentPolynomialSet<N,Order,N,PolySetType,T,Simplex> >,
1128 mpl::identity<detail::MomentPolynomialSet<N,Order,N,PolySetType,T,Hypercube> > >::type::type result_type;
1129 typedef result_type type;
1132 typedef MomentPolynomialSet<Order,Scalar> component_basis_type;
boost::tuple< mpl::size_t< MESH_POINTS >, typename MeshTraits< MeshType >::point_const_iterator, typename MeshTraits< MeshType >::point_const_iterator > points(MeshType const &mesh)
Definition: filters.hpp:1296
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Elements & operator=(Elements const &e)
Definition: elements.hpp:335