https://github.com/feelpp/feelpp
Tip revision: 77ff318a3d2f026f0e706896aecf570465004838 authored by Christophe Prud'homme on 23 February 2022, 16:15:49 UTC
electric does not support P3 @vincentchabannes
electric does not support P3 @vincentchabannes
Tip revision: 77ff318
polynomial_set.hpp
/**
\file Polynomial_set.hpp
\author Abdoulaye Samake <samakeablo@gmail.com>
\date 2010-05-28
*/
/* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=cpp:et:sw=4:ts=4:sts=4
This file is part of the Feel library
Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
Date: 2010-07-02
Copyright (C) 2010 Université Joseph Fourier Grenoble 1
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 3.0 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
/**
\file Polynomial_set.hpp
\author Abdoulaye Samake <samakeablo@gmail.com>
\date 2010-07-02
*/
# ifndef DEF_Polynomial_set_H
# define DEF_Polynomial_set_H
#include <vector>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
namespace Feel
{
/**
* \class Polynomial_set
* \brief polynomial_set class
*
* The polynomial is expressed in the basis from \p Poly. The
* coefficients of the polynomial in this basis are represented by a
* matrix whose lines are the polymomial components coefficients (1 if
* \code is_scalar == true \endcode, \p nDim if \code is_vectorial ==
* true\endcode and columns are the basis
*
* Evaluating the polynomial at a set of points(or just one point) is
* then simply a matrix-matrix product.
*
* \ingroup Polynomial_set
* @author Abdoulaye Samake < samakeablo@gmail.com>
* @see
*/
template< uint nDim ,uint nOrder , uint dim ,uint number >
class Polynomialset
{
public:
/** @name Typedefs
*/
//@{
typedef boost::numeric::ublas::matrix<double> matrix_type ;
typedef boost::numeric::ublas::vector<double> vector_type ;
typedef boost::numeric::ublas::identity_matrix<double> matrix_id ;
typedef boost::numeric::ublas::zero_vector<double> vector_zero ;
typedef boost::numeric::ublas::vector<double> Points_type ;
typedef boost::numeric::ublas::vector<Points_type> vectors_type ;
//@}
//@}
/** @name Constructors, destructor
*/
//@
/**
* default constructor
*/
Polynomialset() ;
Polynomialset ( matrix_type const& m , vector_type const& b ) ;
Polynomialset ( matrix_type const& m ) ;
Polynomialset ( const Polynomialset<nDim,nOrder ,dim ,number> & P ) ;
~Polynomialset() ;
matrix_type evaluates_Points( vectors_type const& P );
private :
matrix_type M ;
vector_type B ;
}; // Polynomial_set
// Implementation
template< uint nDim ,uint nOrder , uint dim ,uint number >
Polynomialset <nDim , nOrder ,dim ,number>::Polynomialset() :M( matrix_id( nOrder,dim ) ),B( vector_zero( nDim ) )
{
};
template< uint nDim ,uint nOrder ,uint dim ,uint number >
Polynomialset <nDim , nOrder , dim ,number>::Polynomialset ( matrix_type const& m , vector_type const& b )
:
M( m ),
B( b )
{
}
template< uint nDim , uint nOrder ,uint dim ,uint number >
Polynomialset <nDim , nOrder , dim ,number>::Polynomialset ( matrix_type const& m )
:
M( m ),
B( ublas::scalar_vector<double>( nDim, 1. ) )
{
}
template< uint nDim , uint nOrder ,uint dim ,uint number >
Polynomialset<nDim ,nOrder ,dim ,number >::Polynomialset ( const Polynomialset< nDim, nOrder ,dim ,number> & P )
:
M( P.M ),
B( P.B )
{
}
template< uint nDim , uint nOrder ,uint dim ,uint number >
Polynomialset<nDim ,nOrder ,dim ,number>::~Polynomialset()
{
}
template< uint nDim , uint nOrder ,uint dim ,uint number >
typename Polynomialset<nDim ,nOrder ,dim ,number>::matrix_type
Polynomialset<nDim ,nOrder ,dim ,number>::evaluates_Points( vectors_type const& P )
{
matrix_type matrice( nDim , number ) ;
matrix_type matrice_valeur( nDim , number ) ;
if ( P( 0 ).size() == 2 )
{
for ( uint j = 0 ; j< matrice.size2() ; ++j )
{
matrice( 0 ,j ) = 1. ;
matrice( 1 ,j ) = B( 1 )*P( j )( 0 ) ;
matrice( 2 ,j ) = B( 2 )*P( j )( 1 ) ;
matrice( 3 ,j ) = B( 3 )*P( j )( 0 )*P( j )( 1 ) ;
matrice( 4 ,j ) = B( 4 )*pow( P( j )( 0 ),2 ) ;
matrice( 5 ,j ) = B( 5 )*pow( P( j )( 1 ),2 ) ;
matrice( 6 ,j ) = B( 6 )*pow( P( j )( 0 ),2 )*P( j )( 1 ) ;
matrice( 7 ,j ) = B( 7 )*pow( P( j )( 1 ),2 )*P( j )( 0 ) ;
matrice( 8 ,j ) = B( 8 )*pow( P( j )( 0 ),3 ) ;
matrice( 9 ,j ) = B( 9 )*pow( P( j )( 1 ),3 ) ;
}
}
else
{
for ( uint j = 0 ; j< matrice.size2() ; ++j )
{
matrice( 0 ,j ) = 1. ;
matrice( 1 ,j ) = B( 1 )*P( j )( 0 ) ;
matrice( 2 ,j ) = B( 2 )*P( j )( 1 ) ;
matrice( 3 ,j ) = B( 3 )*P( j )( 2 ) ;
matrice( 4 ,j ) = B( 4 )*P( j )( 0 )*P( j )( 1 ) ;
matrice( 5 ,j ) = B( 5 )*P( j )( 0 )*P( j )( 2 ) ;
matrice( 6 ,j ) = B( 6 )*P( j )( 1 )*P( j )( 2 ) ;
matrice( 7 ,j ) = B( 7 )*pow( P( j )( 0 ),2 ) ;
matrice( 8 ,j ) = B( 8 )*pow( P( j )( 1 ),2 ) ;
matrice( 9 ,j ) = B( 9 )*pow( P( j )( 2 ),2 ) ;
matrice( 10 ,j ) = B( 10 )*pow( P( j )( 0 ),2 )*P( j )( 1 ) ;
matrice( 11 ,j ) = B( 11 )*pow( P( j )( 1 ),2 )*P( j )( 0 ) ;
matrice( 12 ,j ) = B( 12 )*pow( P( j )( 0 ),2 )*P( j )( 2 ) ;
matrice( 13 ,j ) = B( 13 )*pow( P( j )( 2 ),2 )*P( j )( 0 ) ;
matrice( 14 ,j ) = B( 14 )*pow( P( j )( 1 ),2 )*P( j )( 2 ) ;
matrice( 15 ,j ) = B( 15 )*pow( P( j )( 2 ),2 )*P( j )( 1 ) ;
matrice( 16 ,j ) = B( 16 )*P( j )( 0 )*P( j )( 1 )*P( j )( 2 ) ;
matrice( 17 ,j ) = B( 17 )*pow( P( j )( 0 ),3 ) ;
matrice( 18 ,j ) = B( 18 )*pow( P( j )( 1 ),3 ) ;
matrice( 19 ,j ) = B( 19 )*pow( P( j )( 2 ),3 ) ;
}
}
#if 0
for ( uint i = 0 ; i < matrice_valeur.size1() ; ++i )
{
for ( uint j = 0 ; j< matrice_valeur.size2() ; ++j )
{
matrice_valeur ( i , j ) = matrice ( i , j ) ;
}
}
#endif
return boost::numeric::ublas::prod( M , matrice ) ;
}
} //namespace Feel
#endif