Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
polynomial.hpp
Go to the documentation of this file.
1 /* -*- 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=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-10-06
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2011 Universite de Grenoble
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 #ifndef __Polynomial_H
31 #define __Polynomial_H 1
32 
33 #include <boost/numeric/ublas/vector.hpp>
34 #include <boost/numeric/ublas/vector_proxy.hpp>
35 
36 
37 #include <feel/feelcore/feel.hpp>
38 #include <feel/feelalg/glas.hpp>
39 #include <feel/feelpoly/policy.hpp>
40 
41 namespace Feel
42 {
43 namespace ublas = boost::numeric::ublas;
44 
45 template<typename, template<uint16_type> class PolySetType > class PolynomialSet;
46 
64 template<typename Poly,
65  template<uint16_type> class PolySetType = Scalar,
66  typename Container = typename Poly::basis_type::matrix_type>
68 {
69 public:
70 
74 
75  static const uint16_type nDim = Poly::nDim;
76  static const uint16_type nOrder = Poly::nOrder;
77 
79 
80 
84 
86  typedef typename Poly::value_type value_type;
87  typedef typename Poly::basis_type basis_type;
88 
89 
90 
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;
98 
99  typedef typename Component<polyset_type>::type component_type;
101 
102  typedef typename basis_type::points_type points_type;
103  typedef typename basis_type::matrix_type matrix_type;
104  typedef Container container_type;
105 
106  typedef typename node<value_type>::type node_type;
107 
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 ) );
110 
112 
115  //@
116 
121  :
122  M_basis(),
123  M_coeff( M_basis.coeff() )
124  {
125  }
126 
127 
132  Polynomial( Poly const& __poly )
133  :
134  M_basis( __poly.basis() ),
135  M_coeff( M_basis.coeff() )
136  {
137  }
138 
146  Polynomial( Poly const& __poly, container_type const& __coeff, bool __as_is = false )
147  :
148  M_basis( __poly.basis() ),
149  M_coeff( M_basis.coeff() )
150  {
151  setCoefficient( __coeff, __as_is );
152  }
153 
160  Polynomial( container_type const& __coeff, bool __as_is = false )
161  :
162  M_basis(),
163  M_coeff( M_basis.coeff() )
164  {
165  setCoefficient( __coeff, __as_is );
166  }
167 
175  template<class AE>
176  Polynomial( Poly const& __poly, ublas::matrix_expression<AE> const& __coeff, bool __as_is = false )
177  :
178  M_basis( __poly.basis() ),
179  M_coeff( M_basis.coeff() )
180  {
181  setCoefficient( __coeff, __as_is );
182  }
183 
184  Polynomial( Polynomial const & p )
185  :
186  M_basis( p.M_basis ),
187  M_coeff( p.M_coeff )
188  {}
189 
190  ~Polynomial()
191  {}
192 
194 
198 
199  self_type const& operator=( self_type const& __p )
200  {
201  if ( this != &__p )
202  {
203  M_basis = __p.M_basis;
204  M_coeff = __p.M_coeff;
205  }
206 
207  return *this;
208  }
209 
210  self_type const& operator-=( self_type const& __p )
211  {
212  M_coeff -= __p.M_coeff;
213  return *this;
214  }
215 
216  self_type const& operator()( self_type const& __p ) const
217  {
218  if ( this != &__p )
219  {
220  M_basis = __p.M_basis;
221  M_coeff = __p.M_coeff;
222  }
223 
224  return *this;
225  }
226 
227 
233  component_type operator[]( int i ) const
234  {
235  const int ncols = M_coeff.size2();
236  return component_type( Poly(), ublas::project( M_coeff,
237  ublas::slice( nComponents*i+i, 1, 1 ),
238  ublas::slice( 0, 1, ncols ) ) );
239  }
240 
247  matrix_type
248  operator()( node_type const& __x ) const
249  {
250  return ublas::prod( M_coeff, M_basis( __x ) );
251  }
252 
259  matrix_type operator()( points_type const& __pts ) const
260  {
261  return ublas::prod( M_coeff, M_basis( __pts ) );
262  }
263 
265 
269 
270 
274  bool isZero() const
275  {
276  return ublas::norm_2( M_coeff ) < Feel::type_traits<value_type>::epsilon();
277 
278  }
282  matrix_type const& coeff() const
283  {
284  return M_coeff;
285  }
286 
290  matrix_type const& coefficients() const
291  {
292  return M_coeff;
293  }
294 
298  basis_type const& basis() const
299  {
300  return M_basis;
301  }
302 
304 
308 
312  void setCoefficient( matrix_type const& __c, bool __as_is = false )
313  {
314  //FEELPP_ASSERT( __c.size1() == nComponents*nComponents && __c.size2() == M_coeff.size2() )
315  // ( is_scalar )( is_vectorial )( is_tensor2 )( __c )( M_coeff ).error( "invalid polynomial coefficients" );
316  if ( !__as_is )
317  {
318  M_coeff = ublas::prod( polyset_type::toMatrix( __c ), polyset_type::toMatrix( M_coeff ) );
319  M_coeff = polyset_type::toType( M_coeff );
320  }
321 
322  else
323  M_coeff = __c;
324  }
325 
326 
328 
332 
339  matrix_type
340  evaluate( node_type const& __x ) const
341  {
342  return ublas::prod( M_coeff, M_basis( __x ) );
343  }
344 
351  matrix_type evaluate( points_type const& __pts ) const
352  {
353  return ublas::prod( M_coeff, M_basis( __pts ) );
354  }
355 
356  template<typename AE>
357  matrix_type derivate( uint16_type i, ublas::matrix_expression<AE> const& pts ) const
358  {
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 );
362  return res;
363  }
364 
365  template<typename AE>
366  matrix_type derivate( uint16_type i, uint16_type j, ublas::matrix_expression<AE> const& pts ) const
367  {
368  //std::cout << "[derivate2] M_coeff = " << M_coeff << "\n";
369  matrix_type eval( M_basis.evaluate( pts ) );
370  //matrix_type res( M_coeff.size1(), pts().size2() );
371  //ublas::axpy_prod( M_coeff, der[i], res );
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 );
375  }
376 
377 
384  matrix_type const& d( uint16_type i ) const
385  {
386  return M_basis.d( i );
387  }
388 
395  self_type derivative( uint16_type l ) const
396  {
397  return self_type( Poly(), ublas::prod( M_coeff, M_basis.d( l ) ), true );
398  }
399 
405  PolynomialSet<Poly,PolySetType> toSet( bool asis = false ) const
406  {
407  return PolynomialSet<Poly,PolySetType>( Poly(), M_coeff, asis );
408  }
409 #if 0
411  {
412  matrix_type c = M_coeff-p.M_coeff;
413  std::cout << "c=" << c << "\n";
414  return Polynomial<Poly, PolySetType>( Poly(), c );
415  }
416 #endif
417 
418 
419 
420 protected:
421 
422 private:
423 
424  basis_type M_basis;
425  container_type M_coeff;
426 };
427 template<typename Poly, template<uint16_type> class PolySetType>
428 Polynomial<Poly, PolySetType> operator-( Polynomial<Poly, PolySetType> const& p1,Polynomial<Poly, PolySetType> const& p2 )
429 {
430  auto c = p1.coeff()-p2.coeff();
431  //std::cout << "operator- c=" << c << "\n";
432  return Polynomial<Poly, PolySetType>( Poly(), c );
433 }
434 
435 template<typename Poly,template<uint16_type> class PolySetType, typename Container> const bool Polynomial<Poly,PolySetType,Container>::is_scalar;
436 template<typename Poly,template<uint16_type> class PolySetType, typename Container> const bool Polynomial<Poly,PolySetType, Container>::is_vectorial;
437 template<typename Poly,template<uint16_type> class PolySetType, typename Container> const bool Polynomial<Poly,PolySetType, Container>::is_tensor2;
438 
439 }
440 #endif /* __Polynomial_H */

Generated on Sun Oct 20 2013 08:25:03 for Feel++ by doxygen 1.8.4