29 #ifndef __operations_H
30 #define __operations_H 1
48 template<u
int16_type>
class PsetType
50 typename Poly::value_type
53 return p.
coeff()( 0, 0 );
61 template<
typename Poly,
template<u
int16_type>
class Type>
62 Polynomial<Poly, Type>
72 template<
typename Poly,
template<u
int16_type>
class Type>
73 PolynomialSet<Poly, Type>
83 template<
typename Poly,
template<u
int16_type>
class Type>
84 Polynomial<Poly, Type>
95 template<
typename Poly,
template<u
int16_type>
class Type>
96 PolynomialSet<Poly, Type>
107 template<
typename Poly,
template<u
int16_type>
class Type>
108 Polynomial<Poly, Type>
119 template<
typename Poly,
template<u
int16_type>
class Type>
120 PolynomialSet<Poly, Type>
131 template<
typename Poly>
132 Polynomial<Poly, Vectorial>
137 typedef typename Poly::value_type value_type;
138 const int ndof = p.
coeff().size2();
139 ublas::matrix<value_type> c ( nDim*nDim, ndof );
142 for (
int i = 0; i <nDim; ++i )
143 ublas::row( c, nDim*i+i ) = ublas::row( ublas::prod( p.
coeff(), p.
basis().d( i ) ), 0 );
152 template<
typename Poly>
153 Polynomial<Poly, Tensor2>
159 typedef typename Poly::value_type value_type;
160 const int ndof = p.
coeff().size2();
161 ublas::matrix<value_type> c ( nComponents, ndof );
164 for (
int i = 0; i < nDim; ++i )
165 for (
int j = 0; j < nDim; ++j )
166 ublas::row( c, nComponents*i+nDim*j ) = ublas::row( ublas::prod( p[j].coeff(), p.
basis().d( i ) ), 0 );
175 template<
typename Poly>
176 PolynomialSet<Poly, Vectorial>
181 typedef typename Poly::value_type value_type;
182 const int n1 = p.
coeff().size1();
183 const int n2 = p.
coeff().size2();
184 ublas::matrix<value_type> c ( nDim*nDim*n1, n2 );
187 for (
int i = 0; i <nDim; ++i )
190 ublas::slice( nDim*n1*i+i, nDim, n1 ),
191 ublas::slice( 0, 1, n2 ) ) = ublas::prod( p.
coeff(), p.
basis().d( i ) );
201 template<
typename Poly>
202 PolynomialSet<Poly, Tensor2>
208 typedef typename Poly::value_type value_type;
209 const int n1 = p.
coeff().size1();
210 const int n2 = p.
coeff().size2();
211 ublas::matrix<value_type> c ( nComponents*n1, n2 );
214 for (
int i = 0; i <nDim; ++i )
215 for (
int j = 0; j < nDim; ++j )
218 ublas::slice( nDim*n1*i+i, nDim, n1 ),
219 ublas::slice( 0, 1, n2 ) ) = ublas::prod( p.
coeff(), p.
basis().d( i ) );
229 template<
typename Poly>
230 PolynomialSet<Poly, Tensor2>
236 typedef typename Poly::value_type value_type;
237 const int n1 = p.
coeff().size1();
238 const int n2 = p.
coeff().size2();
239 ublas::matrix<value_type> c ( nComponents*nComponents*n1, n2 );
242 for (
int i = 0; i <nDim; ++i )
243 for (
int j = 0; j < nDim; ++j )
246 ublas::slice( nComponents*n1*( nDim*j+i ), nComponents, n1 ),
247 ublas::slice( 0, 1, n2 ) ) = ublas::prod( p.
coeff(), p.
d( i, j ) );
258 template<
typename Poly>
259 Polynomial<Poly, Scalar>
264 typedef typename Poly::value_type value_type;
266 ublas::matrix<value_type> c ( ublas::zero_matrix<value_type>( 1, p.
coeff().size2() ) );
268 for (
int i = 0; i <nDim; ++i )
270 ublas::noalias( c ) += ublas::prod( p[i].coeff(), p.
basis().d( i ) );
280 template<
typename Poly>
281 PolynomialSet<Poly, Scalar>
286 typedef typename Poly::value_type value_type;
288 ublas::matrix<value_type> c ( ublas::zero_matrix<value_type>( p.
coeff().size1()/nComponents,
289 p.
coeff().size2() ) );
291 for (
int i = 0; i <nComponents; ++i )
293 ublas::noalias( c ) += ublas::prod( p[i].coeff(), p.
basis().d( i ) );
303 template<
typename Poly>
304 Polynomial<Poly, Vectorial>
311 typedef typename Poly::value_type value_type;
313 ublas::matrix<value_type> c ( p.
coeff().size1(), p.
coeff().size2() );
315 const uint16_type dim_p = p.
coeff().size1()/nComponents;
316 const uint16_type ndof = p.
coeff().size2();
319 ublas::slice( 0, nComponents, dim_p/nComponents ),
320 ublas::slice( 0, 1, ndof ) ) = p[2].derivate( 1 ).coeff() - p[1].derivate( 2 ).coeff();
322 ublas::slice( dim_p+1, nComponents, dim_p/nComponents ),
323 ublas::slice( 0, 1, ndof ) ) = p[0].derivate( 2 ).coeff() - p[2].derivate( 0 ).coeff();
325 ublas::slice( 2*dim_p+2, nComponents, dim_p/nComponents ),
326 ublas::slice( 0, 1, ndof ) ) = p[1].derivate( 0 ).coeff() - p[0].derivate( 1 ).coeff();
335 template<
typename Poly>
336 PolynomialSet<Poly, Vectorial>
343 typedef typename Poly::value_type value_type;
345 ublas::matrix<value_type> c ( p.
coeff().size1(), p.
coeff().size2() );
347 const uint16_type dim_p = p.
coeff().size1()/nComponents;
348 const uint16_type ndof = p.
coeff().size2();
351 ublas::slice( 0, nComponents, dim_p/nComponents ),
352 ublas::slice( 0, 1, ndof ) ) = p[2].
derivate( 1 ).coeff() - p[1].
derivate( 2 ).coeff();
354 ublas::slice( dim_p+1, nComponents, dim_p/nComponents ),
355 ublas::slice( 0, 1, ndof ) ) = p[0].
derivate( 2 ).coeff() - p[2].
derivate( 0 ).coeff();
357 ublas::slice( 2*dim_p+2, nComponents, dim_p/nComponents ),
358 ublas::slice( 0, 1, ndof ) ) = p[1].
derivate( 0 ).coeff() - p[0].
derivate( 1 ).coeff();
367 template<
typename Poly>
368 Polynomial<Poly, Vectorial>
375 typedef typename Poly::value_type value_type;
377 ublas::matrix<value_type> c ( p.
coeff().size1(), p.
coeff().size2() );
378 const uint16_type dim_p = p.
coeff().size1()/nComponents;
379 const uint16_type ndof = p.
coeff().size2();
382 ublas::slice( 0, nComponents, dim_p/nComponents ),
383 ublas::slice( 0, 1, ndof ) ) = - p[2].derivate( 1 ).coeff() + p[1].derivate( 2 ).coeff();
385 ublas::slice( dim_p+1, nComponents, dim_p/nComponents ),
386 ublas::slice( 0, 1, ndof ) ) = - p[0].derivate( 2 ).coeff() + p[2].derivate( 0 ).coeff();
388 ublas::slice( 2*dim_p+2, nComponents, dim_p/nComponents ),
389 ublas::slice( 0, 1, ndof ) ) = - p[1].derivate( 0 ).coeff() + p[0].derivate( 1 ).coeff();
398 template<
typename Poly>
399 PolynomialSet<Poly, Vectorial>
406 typedef typename Poly::value_type value_type;
408 ublas::matrix<value_type> c ( p.
coeff().size1(), p.
coeff().size2() );
409 const uint16_type dim_p = p.
coeff().size1()/nComponents;
410 const uint16_type ndof = p.
coeff().size2();
413 ublas::slice( 0, nComponents, dim_p/nComponents ),
414 ublas::slice( 0, 1, ndof ) ) = - p[2].
derivate( 1 ).coeff() + p[1].
derivate( 2 ).coeff();
416 ublas::slice( dim_p+1, nComponents, dim_p/nComponents ),
417 ublas::slice( 0, 1, ndof ) ) = - p[0].
derivate( 2 ).coeff() + p[2].
derivate( 0 ).coeff();
419 ublas::slice( 2*dim_p+2, nComponents, dim_p/nComponents ),
420 ublas::slice( 0, 1, ndof ) ) = - p[1].
derivate( 0 ).coeff() + p[0].
derivate( 1 ).coeff();
432 template<
typename Pset,
435 Polynomial<Pset, Scalar>
436 project( Pset
const& pset, Func
const& f, IM
const& im )
438 typedef typename Func::value_type value_type;
441 ublas::vector<value_type> fq( f( im.points() ) );
444 ublas::matrix<value_type> psetq( pset.evaluate( im.points() ) );
446 ublas::matrix<value_type> c( 1, psetq.size1() );
447 ublas::row( c, 0 ) = ublas::prod( psetq, ublas::element_prod( fq, im.weights() ) );
457 template<
typename P,
template<u
int16_type>
class Type,
458 template<
class,
template<u
int16_type>
class>
class Poly1,
459 template<class, template<uint16_type> class> class Poly2 >
460 PolynomialSet<P, Type>
462 Poly2<P, Type>
const& pset2 )
465 typedef typename res_type::value_type value_type;
467 FEELPP_ASSERT( pset1.coeff().size2() == pset2.coeff().size2() )
468 ( pset1.coeff().size2() )( pset2.coeff().size2() ).error(
"incompatible size" );
470 ublas::matrix<value_type> M( pset1.coeff().size1()+pset2.coeff().size1(), pset1.coeff().size2() );
472 ublas::range( 0, pset1.coeff().size1() ),
473 ublas::range( 0,pset1.coeff().size2() ) ) = pset1.coeff();
475 ublas::range( pset1.coeff().size1(), pset1.coeff().size1()+pset2.coeff().size1() ),
476 ublas::range( 0,pset1.coeff().size2() ) ) = pset2.coeff();
478 M = res_type::polyset_type::toMatrix( M );
480 std::cout <<
"before SVD M = " << M <<
"\n";
481 Eigen::MatrixXd A ( Eigen::MatrixXd::Zero( M.size1(), M.size2() ) );
483 for (
int i = 0; i < M.size1(); ++i )
484 for (
int j = 0; j < M.size2(); ++j )
489 Eigen::SVD<Eigen::MatrixXd> svdOfA( A );
490 std::cout <<
"SVDofA.S() = " << svdOfA.singularValues() <<
"\n";
491 std::cout <<
"SVDofA.V() = " << svdOfA.matrixV().transpose() <<
"\n";
496 ublas::matrix<value_type> m ( svdOfA.singularValues().size(), M.size2() );
498 for (
int i = 0; i < m.size1(); ++i )
499 for (
int j = 0; j < m.size2(); ++j )
501 m( i,j ) = svdOfA.matrixV()( j, i );
505 ublas::matrix<value_type> m ( ublas::subrange( svd.
V(), 0, svd.
S().
size(), 0, M.size2() ) );
507 return res_type( P(), res_type::polyset_type::toType( m ),
true );
self_type derivate(uint16_type l) const
Derivate with respect to the l-th direction.
Definition: polynomialset.hpp:554
Poly::value_type integrate(Polynomial< Poly, PsetType > const &p)
integrate a polynomial over a convex
Definition: operations.hpp:51
Polynomial< Pset, Scalar > project(Pset const &pset, Func const &f, IM const &im)
Definition: operations.hpp:436
basis_type const & basis() const
Definition: polynomialset.hpp:253
Polynomial< Poly, Type > dx(Polynomial< Poly, Type > const &p)
compute of a polynomial p
Definition: operations.hpp:63
PolynomialSet< Poly, Tensor2 > hessian(PolynomialSet< Poly, Scalar > const &p)
compute the gradient of a vectorial polynomial set
Definition: operations.hpp:231
polynomial class
Definition: polynomial.hpp:67
Polynomial< Poly, Vectorial > gradient(Polynomial< Poly, Scalar > const &p)
compute the gradient of a scalar polynomial p
Definition: operations.hpp:133
matrix_type const & coeff() const
Definition: polynomialset.hpp:245
matrix_type const & d(uint16_type i) const
derivatives of Dubiner polynomials the derivatives are computed at the nodes of the lattice ...
Definition: polynomialset.hpp:539
Polynomial< Poly, Type > dz(Polynomial< Poly, Type > const &p)
compute of a polynomial p
Definition: operations.hpp:109
PolynomialSet< P, Type > unite(Poly1< P, Type > const &pset1, Poly2< P, Type > const &pset2)
union of two polynomial sets P1 and P2 :
Definition: operations.hpp:461
matrix_type const & coeff() const
Definition: polynomial.hpp:282
Polynomial< Poly, Vectorial > curl(Polynomial< Poly, Vectorial > const &p)
compute the curl of a vectorial polynomial p
Definition: operations.hpp:305
vector_type const & S() const
Definition: svd.hpp:121
Polynomial< Poly, Vectorial > curlt(Polynomial< Poly, Vectorial > const &p)
compute the transpose of the curl of a vectorial polynomial p
Definition: operations.hpp:369
a Set of polynomials
Definition: fe.hpp:37
Singular Value Decomposition of a rectangular matrix.
Definition: svd.hpp:73
matrix_type const & V() const
Definition: svd.hpp:115
Polynomial< Poly, Scalar > divergence(Polynomial< Poly, Vectorial > const &p)
compute the divergence of a vectorial polynomial p
Definition: operations.hpp:260
virtual size_type size() const
Definition: vector.hpp:290
basis_type const & basis() const
Definition: polynomial.hpp:298
Polynomial< Poly, Type > dy(Polynomial< Poly, Type > const &p)
compute of a polynomial p
Definition: operations.hpp:85