30 #ifndef __Polynomial_H
31 #define __Polynomial_H 1
33 #include <boost/numeric/ublas/vector.hpp>
34 #include <boost/numeric/ublas/vector_proxy.hpp>
37 #include <feel/feelcore/feel.hpp>
39 #include <feel/feelpoly/policy.hpp>
43 namespace ublas = boost::numeric::ublas;
45 template<
typename,
template<u
int16_type>
class PolySetType >
class PolynomialSet;
64 template<
typename Poly,
65 template<u
int16_type>
class PolySetType = Scalar,
66 typename Container =
typename Poly::basis_type::matrix_type>
75 static const uint16_type nDim = Poly::nDim;
76 static const uint16_type nOrder = Poly::nOrder;
86 typedef typename Poly::value_type value_type;
87 typedef typename Poly::basis_type basis_type;
91 typedef PolySetType<nDim> polyset_type;
92 static const bool is_tensor2 = polyset_type::is_tensor2;
93 static const bool is_vectorial = polyset_type::is_vectorial;
94 static const bool is_scalar = polyset_type::is_scalar;
95 static const uint16_type nComponents = polyset_type::nComponents;
96 static const uint16_type nComponents1 = polyset_type::nComponents1;
97 static const uint16_type nComponents2 = polyset_type::nComponents2;
99 typedef typename Component<polyset_type>::type component_type;
102 typedef typename basis_type::points_type points_type;
103 typedef typename basis_type::matrix_type matrix_type;
104 typedef Container container_type;
108 BOOST_STATIC_ASSERT( ( boost::is_same<typename matrix_type::value_type, value_type>::value ) );
109 BOOST_STATIC_ASSERT( ( boost::is_same<typename matrix_type::value_type, typename points_type::value_type>::value ) );
123 M_coeff( M_basis.
coeff() )
134 M_basis( __poly.
basis() ),
135 M_coeff( M_basis.
coeff() )
146 Polynomial( Poly
const& __poly, container_type
const& __coeff,
bool __as_is =
false )
148 M_basis( __poly.
basis() ),
149 M_coeff( M_basis.
coeff() )
160 Polynomial( container_type
const& __coeff,
bool __as_is =
false )
163 M_coeff( M_basis.
coeff() )
176 Polynomial( Poly
const& __poly, ublas::matrix_expression<AE>
const& __coeff,
bool __as_is =
false )
178 M_basis( __poly.
basis() ),
179 M_coeff( M_basis.
coeff() )
186 M_basis( p.M_basis ),
199 self_type
const& operator=( self_type
const& __p )
203 M_basis = __p.M_basis;
204 M_coeff = __p.M_coeff;
210 self_type
const& operator-=( self_type
const& __p )
212 M_coeff -= __p.M_coeff;
216 self_type
const& operator()( self_type
const& __p )
const
220 M_basis = __p.M_basis;
221 M_coeff = __p.M_coeff;
235 const int ncols = M_coeff.size2();
237 ublas::slice( nComponents*i+i, 1, 1 ),
238 ublas::slice( 0, 1, ncols ) ) );
250 return ublas::prod( M_coeff, M_basis( __x ) );
261 return ublas::prod( M_coeff, M_basis( __pts ) );
276 return ublas::norm_2( M_coeff ) < Feel::type_traits<value_type>::epsilon();
318 M_coeff = ublas::prod( polyset_type::toMatrix( __c ), polyset_type::toMatrix( M_coeff ) );
319 M_coeff = polyset_type::toType( M_coeff );
342 return ublas::prod( M_coeff, M_basis( __x ) );
351 matrix_type
evaluate( points_type
const& __pts )
const
353 return ublas::prod( M_coeff, M_basis( __pts ) );
356 template<
typename AE>
357 matrix_type derivate( uint16_type i, ublas::matrix_expression<AE>
const& pts )
const
359 ublas::vector<matrix_type> der( M_basis.derivate( pts ) );
360 matrix_type res( M_coeff.size1(), pts().size2() );
361 ublas::axpy_prod( M_coeff, der[i], res );
365 template<
typename AE>
366 matrix_type derivate( uint16_type i, uint16_type j, ublas::matrix_expression<AE>
const& pts )
const
369 matrix_type eval( M_basis.evaluate( pts ) );
372 matrix_type p1 = ublas::prod( M_coeff, M_basis.d( i ) );
373 matrix_type p2 = ublas::prod( p1, M_basis.d( j ) );
374 return ublas::prod( p2, eval );
384 matrix_type
const&
d( uint16_type i )
const
386 return M_basis.d( i );
397 return self_type( Poly(), ublas::prod( M_coeff, M_basis.d( l ) ),
true );
412 matrix_type c = M_coeff-p.M_coeff;
413 std::cout <<
"c=" << c <<
"\n";
425 container_type M_coeff;
427 template<
typename Poly,
template<u
int16_type>
class PolySetType>
428 Polynomial<Poly, PolySetType> operator-( Polynomial<Poly, PolySetType>
const& p1,Polynomial<Poly, PolySetType>
const& p2 )
430 auto c = p1.coeff()-p2.coeff();
432 return Polynomial<Poly, PolySetType>( Poly(), c );
435 template<
typename Poly,
template<u
int16_type>
class PolySetType,
typename Container>
const bool Polynomial<Poly,PolySetType,Container>::is_scalar;
436 template<
typename Poly,
template<u
int16_type>
class PolySetType,
typename Container>
const bool Polynomial<Poly,PolySetType, Container>::is_vectorial;
437 template<
typename Poly,
template<u
int16_type>
class PolySetType,
typename Container>
const bool Polynomial<Poly,PolySetType, Container>::is_tensor2;