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
legendre.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: 2006-02-20
7 
8  Copyright (C) 2006 EPFL
9  Copyright (C) 2007 Université Joseph Fourier Grenoble 1
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 __Legendre_H
31 #define __Legendre_H 1
32 
33 #include <vector>
34 
35 #include <boost/lambda/if.hpp>
36 
38 #include <feel/feelalg/glas.hpp>
39 #include <feel/feelalg/lu.hpp>
41 #include <feel/feelpoly/policy.hpp>
45 
46 namespace Feel
47 {
48 template< class Convex,
49  uint16_type Order,
50  typename T >
51 class PointSetGaussLobatto;
52 
53 template<uint16_type Dim,
54  uint16_type RealDim,
55  uint16_type Degree,
56  typename NormalizationPolicy,
57  typename T,
58  template<class> class StoragePolicy>
59 class Legendre;
60 
61 
62 template<uint16_type Dim,
63  uint16_type RealDim,
64  uint16_type Degree,
65  typename NormalizationPolicy = Normalized<true>,
66  typename T = double,
67  template<class> class StoragePolicy = StorageUBlas>
68 struct LegendreTraits
69 {
70  static const uint16_type nDim = Dim;
71  static const uint16_type nRealDim = RealDim;
72  static const uint16_type nOrder = Degree;
73  static const uint16_type nConvexOrderDiff = nOrder+2;
74  static const bool is_normalized = NormalizationPolicy::is_normalized;
75 
79 
80  /*
81  * numerical type
82  */
83  typedef T value_type;
84 
85  template<uint16_type order, typename V = value_type>
86  struct Convex
87  {
88  typedef Hypercube<nDim, order, /*nRealDim*/nDim> type;
89  typedef Reference<Hypercube<nDim, order, /*nRealDim*/nDim>, nDim, order, nDim/*nRealDim*/, V> reference_type;
90  };
91 
92  template<typename NewT>
93  struct ChangeValueType
94  {
96  typedef LegendreTraits<Dim, RealDim, Degree, NormalizationPolicy, NewT, StoragePolicy> traits_type;
97  };
98 
99  template<uint16_type NewOrder>
100  struct ChangeOrder
101  {
102  typedef Legendre<Dim, RealDim, NewOrder, NormalizationPolicy, T, StoragePolicy> type;
103  typedef LegendreTraits<Dim, RealDim, NewOrder, NormalizationPolicy, T, StoragePolicy> traits_type;
104  };
105 
106  /*
107  * Geometry where the polynomials are defined and constructed
108  */
109  typedef typename Convex<nOrder>::type convex_type;
110  typedef typename Convex<nOrder>::reference_type reference_convex_type;
111 
112  typedef typename Convex<nConvexOrderDiff>::type diff_convex_type;
113  typedef typename Convex<nConvexOrderDiff>::reference_type diff_reference_convex_type;
114 
115  typedef PointSetGaussLobatto<diff_convex_type, nConvexOrderDiff, value_type> diff_pointset_type;
116 
117  /*
118  * storage policy
119  */
120  typedef StoragePolicy<value_type> storage_policy;
121  typedef typename storage_policy::matrix_type matrix_type;
122  typedef typename storage_policy::vector_matrix_type vector_matrix_type;
123  typedef typename storage_policy::matrix_node_type matrix_node_type;
124  typedef typename storage_policy::points_type points_type;
125  typedef typename storage_policy::node_type node_type;
126 };
127 
128 template<int D, int O>
129 struct LegendreTag
130 {
131  static const int Dim = D;
132  static const int Order = O;
133 };
155 template<uint16_type Dim,
156  uint16_type RealDim,
157  uint16_type Degree,
158  typename NormalizationPolicy = Normalized<true>,
159  typename T = double,
160  template<class> class StoragePolicy = StorageUBlas>
161 class Legendre
162 {
163 public:
164 
165  typedef LegendreTraits<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy> traits_type;
166  static const uint16_type nDim = Dim;
167  static const uint16_type nRealDim = RealDim;
168  static const uint16_type nOrder = Degree;
169  static const uint16_type nConvexOrder = nOrder+2;
170  static const bool is_normalized = NormalizationPolicy::is_normalized;
171  static const bool isTransformationEquivalent = true;
172  static const bool isContinuous = false;
173  static const bool is_product = true;
174  typedef Discontinuous continuity_type;
175 
179 
180  typedef Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy> self_type;
181 
182  /*
183  * for now can be used only as a basis but we might be interested
184  * to have then expressed in other basis like in the moment basis
185  */
186  typedef self_type basis_type;
187 
188  typedef typename traits_type::value_type value_type;
189 
190  /*
191  * Geometry where the polynomials are defined and constructed
192  */
193  typedef typename traits_type::convex_type convex_type;
194  typedef typename traits_type::reference_convex_type reference_convex_type;
195 
196  typedef typename traits_type::diff_pointset_type diff_pointset_type;
197 
198  /*
199  * storage policy
200  */
201  typedef typename traits_type::storage_policy storage_policy;
202  typedef typename traits_type::matrix_type matrix_type;
203  typedef typename traits_type::vector_matrix_type vector_matrix_type;
204  typedef typename traits_type::matrix_node_type matrix_node_type;
205  typedef typename traits_type::points_type points_type;
206  typedef typename traits_type::node_type node_type;
207 
208 
210 
214 
215  Legendre()
216  :
217  M_refconvex(),
218  M_pts( M_refconvex.makePoints( nDim, 0 ) )
219  {
220  this->initDerivation();
221  }
222  Legendre( Legendre const & d )
223  :
224  M_refconvex(),
225  M_pts( d.M_pts )
226  {
227 
228  }
229 
230  ~Legendre()
231  {}
232 
234 
238 
239  self_type const& operator=( self_type const& d )
240  {
241  if ( this != &d )
242  {
243  M_pts = d.M_pts;
244  _S_D = d._S_D;
245  }
246 
247  return *this;
248  }
249 
250  //ublas::matrix_column<matrix_type const> operator()( node_type const& pt ) const
251  matrix_type operator()( node_type const& pt ) const
252  {
253  points_type pts( pt.size(), 1 );
254  ublas::column( pts, 0 ) = pt;
255  return evaluate( pts );
256  }
257 
258  matrix_type operator()( points_type const& pts ) const
259  {
260  return evaluate( pts );
261  }
262 
264 
268 
272  size_type size() const
273  {
274  return convex_type::polyDims( nOrder );
275  }
276 
281  uint16_type degree() const
282  {
283  return nOrder;
284  }
285 
289  self_type const& basis() const
290  {
291  return *this;
292  }
293 
298  bool isNormalized() const
299  {
300  return is_normalized;
301  }
302 
306  std::string familyName() const
307  {
308  return "legendre";
309  }
310 
312 
316 
317 
319 
323 
333  matrix_type coeff() const
334  {
335 #if 0
336  std::cout << "[Legendre::coeff] coeff = "
337  << ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() )
338  << "\n";
339 #endif
340  return ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() );
341  }
342 
343 
349  static matrix_type evaluate( points_type const& __pts )
350  {
351  return evaluate( __pts, mpl::int_<nDim>() );
352  }
353 
354  template<typename AE>
355  static vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts )
356  {
357  return derivate( __pts, mpl::int_<nDim>() );
358  }
359 
366  static matrix_type const& d( uint16_type i )
367  {
368  return _S_D[i];
369  }
370 
377  static matrix_type const& derivate( uint16_type i )
378  {
379  return _S_D[i];
380  }
381 
382 
384 
385 private:
386 private:
387 
388  static value_type normalization( int i )
389  {
390  return ( is_normalized?math::sqrt( value_type( i ) + 0.5 ) : value_type( 1 ) );
391  }
392  static value_type normalization( int i, int j )
393  {
394  return ( is_normalized?math::sqrt( ( value_type( i ) + 0.5 ) *
395  ( value_type( j ) + 0.5 ) ) : value_type( 1 ) );
396  }
397  static value_type normalization( int i, int j, int k )
398  {
399  return ( is_normalized?math::sqrt( ( value_type( i ) + 0.5 ) *
400  ( value_type( j ) + 0.5 ) *
401  ( value_type( k ) + 0.5 ) ) : value_type( 1 ) );
402  }
407  static matrix_type
408  evaluate( points_type const& __pts, mpl::int_<1> )
409  {
410  matrix_type m ( JacobiBatchEvaluation<nOrder,value_type>( 0.0, 0.0, ublas::row( __pts, 0 ) ) );
411 
412  if ( is_normalized )
413  {
414  for ( uint16_type i = 0; i < m.size1(); ++i )
415  ublas::row( m, i ) *= normalization( i );
416  }
417 
418  return m;
419  }
420 
425  template<typename AE>
426  static vector_matrix_type
427  derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<1> )
428  {
429  FEELPP_ASSERT( __pts().size1() == 1 )( __pts().size1() )( __pts().size2() ).error( "invalid points" );
430  // VLOG(1) << "Expansion::derivate<1>] number of points " << __pts().size2() << "\n";
431 
432  vector_matrix_type D( 1 );
433  D[0].resize( nOrder+1, __pts().size2() );
434  D[0] = JacobiBatchDerivation<nOrder,value_type>( 0.0, 0.0, ublas::row( __pts(),0 ) );
435 
436  if ( is_normalized )
437  for ( uint16_type i = 0; i < nOrder+1; ++i )
438  ublas::row( D[0], i ) *= normalization( i );
439 
440  return D;
441  }
442 
447  static matrix_type evaluate( points_type const& __pts, mpl::int_<2> );
448 
453  template<typename AE>
454  static vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<2> );
455 
460  static matrix_type evaluate( points_type const& __pts, mpl::int_<3> );
461 
466  template<typename AE>
467  static vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<3> );
468 
469  static void initDerivation();
470 private:
471  reference_convex_type M_refconvex;
472  points_type M_pts;
477  static bool _S_has_derivation;
478 
483  static std::vector<matrix_type> _S_D;
484 
485 }; // class Dubiner
486 
487 template<uint16_type Dim,
488  uint16_type RealDim,
489  uint16_type Degree,
490  typename NormalizationPolicy,
491  typename T,
492  template<class> class StoragePolicy>
493 bool Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::_S_has_derivation = false;
494 
495 template<uint16_type Dim,
496  uint16_type RealDim,
497  uint16_type Degree,
498  typename NormalizationPolicy,
499  typename T,
500  template<class> class StoragePolicy>
501 std::vector<typename Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::matrix_type>
502 Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::_S_D;
503 
504 template<uint16_type Dim,
505  uint16_type RealDim,
506  uint16_type Degree,
507  typename NormalizationPolicy,
508  typename T,
509  template<class> class StoragePolicy>
510 void
511 Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::initDerivation()
512 {
513 #if 0
514  typedef typename traits_type::convex_type convex_type;
515  typedef typename traits_type::reference_convex_type reference_convex_type;
516 
517  typedef typename traits_type::diff_pointset_type diff_pointset_type;
518 
519  typedef typename traits_type::storage_policy storage_policy;
520  typedef typename traits_type::matrix_type matrix_type;
521  typedef typename traits_type::vector_matrix_type vector_matrix_type;
522  typedef typename traits_type::matrix_node_type matrix_node_type;
523  typedef typename traits_type::points_type points_type;
524  typedef typename traits_type::node_type node_type;
525 #endif // 0
526 
527  if ( _S_has_derivation == false )
528  {
529  _S_has_derivation = true;
530 
531  reference_convex_type refconvex;
532  // constructor pointset for differentiation only in
533  // the interior(1)
534  diff_pointset_type diff_pts( 1 );
535  matrix_type A( evaluate( diff_pts.points() ) );
536 
537 #if 1
538  matrix_type D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
539  LU<matrix_type> lu( A );
540  matrix_type C = lu.solve( D );
541 
542  vector_matrix_type d ( derivate( diff_pts.points() ) );
543  _S_D.resize( d.size() );
544 
545  for ( size_type i = 0; i < d.size(); ++i )
546  {
547  _S_D[i] = ublas::prod( d[i], C );
548  glas::clean( _S_D[i] );
549  }
550 
551 #endif
552  }
553 }
554 
555 template<uint16_type Dim,
556  uint16_type RealDim,
557  uint16_type Degree,
558  typename NormalizationPolicy,
559  typename T,
560  template<class> class StoragePolicy>
561 typename Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::matrix_type
563 {
564  matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
565 
566  ublas::vector<value_type> eta1s = ublas::row( __pts, 0 );
567  ublas::vector<value_type> eta2s = ublas::row( __pts, 1 );
568 
569  matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
570  matrix_type bs( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta2s ) );
571 
572  for ( uint16_type cur = 0, i = 0; i < nOrder+1; ++i )
573  {
574  for ( uint16_type j = 0; j < nOrder+1; ++j,++cur )
575  {
576  ublas::row( res, cur ) = normalization( i, j ) * ublas::element_prod( ublas::row( as, i ),
577  ublas::row( bs, j ) );
578  }
579  }
580 
581  return res;
582 }
583 
584 template<uint16_type Dim,
585  uint16_type RealDim,
586  uint16_type Degree,
587  typename NormalizationPolicy,
588  typename T,
589  template<class> class StoragePolicy>
590 template<typename AE>
591 typename Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::vector_matrix_type
592 Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<2> )
593 {
594  vector_matrix_type res( 2 );
595  res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
596  res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
597 
598  // evaluate Legendre polynomials components
599  matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 0 ) ) );
600  matrix_type das( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 0 ) ) );
601  matrix_type bs( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 1 ) ) );
602  matrix_type dbs( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 1 ) ) );
603 
604  for ( uint16_type cur = 0, i = 0; i < nOrder+1; ++i )
605  {
606  for ( uint16_type j = 0; j < nOrder+1; ++j,++cur )
607  {
608  ublas::row( res[0], cur ) = normalization( i, j ) * ublas::element_prod( ublas::row( das, i ),
609  ublas::row( bs, j ) );
610  ublas::row( res[1], cur ) = normalization( i, j ) * ublas::element_prod( ublas::row( as, i ),
611  ublas::row( dbs, j ) );
612  }
613  }
614 
615  return res;
616 }
617 
618 template<uint16_type Dim,
619  uint16_type RealDim,
620  uint16_type Degree,
621  typename NormalizationPolicy,
622  typename T,
623  template<class> class StoragePolicy>
624 typename Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::matrix_type
626 {
627  matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
628 
629  ublas::vector<value_type> eta1s = ublas::row( __pts, 0 );
630  ublas::vector<value_type> eta2s = ublas::row( __pts, 1 );
631  ublas::vector<value_type> eta3s = ublas::row( __pts, 2 );
632 
633  matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
634  matrix_type bs( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta2s ) );
635  matrix_type cs( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta3s ) );
636 
637  for ( uint16_type cur = 0, i = 0; i < nOrder+1; ++i )
638  {
639  for ( uint16_type j = 0; j < nOrder+1; ++j )
640  {
641  for ( uint16_type k = 0; k < nOrder+1; ++k,++cur )
642  {
643  ublas::row( res, cur ) =
644  normalization( i, j, k )*
645  ublas::element_prod( ublas::element_prod( ublas::row( as, i ),
646  ublas::row( bs, j ) ),
647  ublas::row( cs, k ) );
648  }
649  }
650  }
651 
652  return res;
653 }
654 
655 template<uint16_type Dim,
656  uint16_type RealDim,
657  uint16_type Degree,
658  typename NormalizationPolicy,
659  typename T,
660  template<class> class StoragePolicy>
661 template<typename AE>
662 typename Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::vector_matrix_type
663 Legendre<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<3> )
664 {
665  vector_matrix_type res( 3 );
666  res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
667  res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
668  res[2].resize( convex_type::polyDims( nOrder ), __pts().size2() );
669 
670  // evaluate Legendre polynomials components
671  matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 0 ) ) );
672  matrix_type das( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 0 ) ) );
673  matrix_type bs( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 1 ) ) );
674  matrix_type dbs( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 1 ) ) );
675  matrix_type cs( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 2 ) ) );
676  matrix_type dcs( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, ublas::row( __pts(), 2 ) ) );
677 
678  for ( uint16_type cur = 0, i = 0; i < nOrder+1; ++i )
679  {
680  for ( uint16_type j = 0; j < nOrder+1; ++j )
681  {
682  for ( uint16_type k = 0; k < nOrder+1; ++k,++cur )
683  {
684  ublas::row( res[0], cur ) = ( normalization( i, j, k ) *
685  ublas::element_prod( ublas::element_prod( ublas::row( das, i ),
686  ublas::row( bs, j ) ),
687  ublas::row( cs, k ) ) );
688 
689  ublas::row( res[1], cur ) = ( normalization( i, j, k ) *
690  ublas::element_prod( ublas::element_prod( ublas::row( as, i ),
691  ublas::row( dbs, j ) ),
692  ublas::row( cs, k ) ) );
693 
694  ublas::row( res[2], cur ) = ( normalization( i, j, k ) *
695  ublas::element_prod( ublas::element_prod( ublas::row( as, i ),
696  ublas::row( bs, j ) ),
697  ublas::row( dcs, k ) ) );
698 
699 
700  }
701  }
702  }
703 
704  return res;
705 }
706 
707 
708 }
709 #endif /* __Legendre_H */

Generated on Fri Oct 25 2013 14:24:17 for Feel++ by doxygen 1.8.4