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
moment.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): Gilles Steiner <gilles.steiner@epfl.ch>
6  Christophe Prud'homme <christophe.prudhomme@feelpp.org>
7  Date: 2005-12-05
8 
9  Copyright (C) 2005,2006 EPFL
10  Copyright (C) 2011 Universite Joseph Fourier
11 
12  This library is free software; you can redistribute it and/or
13  modify it under the terms of the GNU Lesser General Public
14  License as published by the Free Software Foundation; either
15  version 3.0 of the License, or (at your option) any later version.
16 
17  This library is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20  Lesser General Public License for more details.
21 
22  You should have received a copy of the GNU Lesser General Public
23  License along with this library; if not, write to the Free Software
24  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
25 */
31 #ifndef __Moment_H
32 #define __Moment_H 1
33 
34 #include <vector>
35 
39 #include <feel/feelalg/lu.hpp>
41 #include <feel/feelpoly/policy.hpp>
43 
44 
47 
48 namespace Feel
49 {
50 
66 template<uint16_type Dim,
67  uint16_type Degree,
68  // typename NormalizationPolicy = Normalized<false>,
69  typename TheConvex=Simplex<Dim>,
70  typename T = double,
71  template<class> class StoragePolicy = StorageUBlas>
72 class Moment
73 {
74 public:
75 
76  static const uint16_type nDim = Dim;
77  static const uint16_type nRealDim = Dim;
78  static const uint16_type nOrder = Degree;
79  static const uint16_type nConvexOrder = mpl::if_<mpl::bool_<TheConvex::is_simplex>,
80  mpl::int_<nDim+nOrder+1>,
81  mpl::int_<nOrder+2> >::type::value;
82  // static const bool is_normalized = NormalizationPolicy::is_normalized;
83  static const uint16_type convex_is_simplex = TheConvex::is_simplex;
84  static const uint16_type convex_is_hypercube = TheConvex::is_hypercube;
85  static const bool is_product = false;
86 
90 
91  typedef Moment<Dim, Degree,TheConvex,T, StoragePolicy> self_type;
92 
93  typedef self_type basis_type;
94 
95  /*
96  * numerical type
97  */
98  typedef T value_type;
99 
100  template<uint16_type order, typename V = value_type>
101  struct Convex
102  {
103  typedef typename mpl::if_<mpl::bool_<TheConvex::is_simplex>,
104  mpl::identity<Simplex<nDim, order, /*nRealDim*/nDim> >,
105  mpl::identity<Hypercube<nDim, order, /*nRealDim*/nDim> > >::type::type type;
106  typedef typename mpl::if_<mpl::bool_<TheConvex::is_simplex>,
107  mpl::identity<Reference<Simplex<nDim, order, /*nRealDim*/nDim>, nDim, order, nDim/*nRealDim*/, V > >,
108  mpl::identity<Reference<Hypercube<nDim, order, /*nRealDim*/nDim>, nDim, order, nDim/*nRealDim*/, V > > >::type::type reference_type;
109  };
110 
111  /*
112  * Geometry where the polynomials are defined and constructed
113  */
114  typedef TheConvex convex_type;
115  typedef typename Convex<nOrder>::reference_type reference_convex_type;
116  //typedef typename reference_convex_type::points_type points_type;
117 
118  /*
119  * storage policy
120  */
121  typedef StoragePolicy<value_type> storage_policy;
122  typedef typename storage_policy::matrix_type matrix_type;
123  typedef typename storage_policy::vector_matrix_type vector_matrix_type;
124  typedef typename storage_policy::matrix_node_type matrix_node_type;
125  typedef typename storage_policy::node_type node_type;
126  typedef typename storage_policy::points_type points_type;
127 
129 
133 
134  Moment()
135  :
136  M_refconvex(),
137  M_pts( M_refconvex.points() ),//M_refconvex.makePoints( nDim, 0 ) ),
138  M_coord( convex_type::polyDims( nOrder ),nDim )
139  {
140  //std::cout << "[polydims: " << convex_type::polyDims( nOrder ) << "\n";
141  if ( nDim == 1 )
142  for ( uint32_type i=0; i < M_coord.size1(); ++i )
143  M_coord( i,0 ) = i;
144 
145  else if ( nDim == 2 )
146  {
147  uint32_type G_i = 0;
148 
149  if ( convex_is_simplex )
150  {
151  for ( int32_type n = 0; n < nOrder+1; ++n )
152  for ( int32_type i = n; i >= 0; --i )
153  {
154  M_coord( G_i,0 ) = i;
155  M_coord( G_i,1 ) = n-i;
156  ++G_i;
157  }
158  }
159 
160  else if ( convex_is_hypercube )
161  {
162  for ( int32_type n = 0; n < nOrder+1; ++n )
163  for ( int32_type i = 0; i < nOrder+1; ++i )
164  {
165  M_coord( G_i,0 ) = i;
166  M_coord( G_i,1 ) = n;
167  ++G_i;
168  }
169  }
170  }
171 
172  else if ( nDim == 3 )
173  {
174  uint32_type G_i = 0;
175 
176  if ( convex_is_simplex )
177  {
178  for ( int32_type n = 0; n < nOrder+1; ++n )
179  for ( int32_type i = n; i >= 0; --i ) // x id
180  for ( int32_type j = n-i; j >= 0 ; --j )
181  {
182  M_coord( G_i,0 ) = i;
183  M_coord( G_i,1 ) = j;
184  M_coord( G_i,2 ) = n-i-j;
185  ++G_i;
186  }
187  }
188 
189  else if ( convex_is_hypercube )
190  {
191  for ( int32_type n = 0; n < nOrder+1; ++n )
192  for ( int32_type i = 0; i < nOrder+1; ++i )
193  for ( int32_type j = 0; j < nOrder+1; ++j )
194  {
195  M_coord( G_i,0 ) = j;
196  M_coord( G_i,1 ) = i;
197  M_coord( G_i,1 ) = n;
198  ++G_i;
199  }
200  }
201  }
202 
203  //std::cout << "nDim: " << nDim << "\n";
204  //std::cout << "orders: " << M_coord << "\n";
205  //std::cout << "pts: " << M_pts << "\n";
206 
207 
208  matrix_type A( ublas::trans( evaluate( M_pts ) ) );
209  //std::cout << "A=" << A << "\n";
210  matrix_type D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
211  //std::cout << "D=" << D << "\n";
212  LU<matrix_type> lu( A );
213  matrix_type C = lu.solve( D );
214  //std::cout << "C=" << C << "\n";
215  vector_matrix_type d ( derivate( M_pts ) );
216  M_D.resize( d.size() );
217 
218  for ( size_type i = 0; i < d.size(); ++i )
219  {
220  M_D[i] = ublas::prod( C, ublas::trans( d[i] ) );
221  //std::cout << "DM=" << M_D[i] << "\n";
222  }
223  }
224  Moment( Moment const & d )
225  :
226  M_refconvex(),
227  M_pts( d.M_pts ),
228  M_coord( d.M_coord ),
229  M_D( d.M_D )
230  {}
231 
232  ~Moment()
233  {}
234 
236 
241  points_type points()
242  {
243  return M_pts;
244  }
245 
246 
250 
251  self_type const& operator=( self_type const& d )
252  {
253  if ( this != &d )
254  {
255  M_pts = d.M_pts;
256  M_D = d.M_D;
257  }
258 
259  return *this;
260  }
261 
262  //ublas::matrix_column<matrix_type const> operator()( node_type const& pt ) const
263  matrix_type operator()( node_type const& pt ) const
264  {
265  points_type pts( pt.size(), 1 );
266  ublas::column( pts, 0 ) = pt;
267  return evaluate( pts );
268  }
269 
270  matrix_type operator()( points_type const& pts ) const
271  {
272  return evaluate( pts );
273  }
274 
276 
280 
281  ublas::matrix<int16_type> coord()
282  {
283  return M_coord;
284  }
285 
290  uint16_type degree() const
291  {
292  return nOrder;
293  }
294 
298  self_type const& basis() const
299  {
300  return *this;
301  }
302 
316  std::string familyName() const
317  {
318  return "moment";
319  }
320 
322 
326 
327 
329 
333 
343  matrix_type coeff() const
344  {
345 #if 0
346  std::cout << "[Moment::coeff] coeff = "
347  << ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() )
348  << "\n";
349 #endif
350  return ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() );
351  }
352 
353 
359  matrix_type evaluate( points_type const& __pts ) const
360  {
361  return evaluate( __pts, mpl::int_<nDim>() );
362  }
363 
364  template<typename AE>
365  vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts ) const
366  {
367  return derivate( __pts, mpl::int_<nDim>() );
368  }
369 
370 
377  matrix_type const& d( uint16_type i ) const
378  {
379  return M_D[i];
380  }
381 
382  template<template<uint16_type> class PolySetType = Scalar>
383  Polynomial<self_type,PolySetType> pick( int i, int c = 0 ) const
384  {
385  size_type dim_p = convex_type::polyDims( nDim );
386  int nComponents = PolySetType<Dim>::nComponents;
387  matrix_type coeff( nComponents, dim_p );
388  coeff = ublas::scalar_matrix<value_type>( coeff.size1(), coeff.size2(), 0. );
389  coeff( c, i ) = 1;
390  //std::cout << "pick.coeff(" << i << "," << c << ") = " << coeff << "\n";
391  return Polynomial<self_type, PolySetType>( *this, coeff, true );
392  //std::cout << "p1.coeff(" << i << "," << c << ") = " << p.coeff() << "\n";
393  //return p;
394  }
395 
403  value_type SP_integrate( ublas::matrix<value_type> const& __Comp )
404  {
405  return SP_integrate( __Comp, mpl::int_<nDim>() );
406  }
407 
408 
415  value_type S_integrate( ublas::matrix<value_type> const& __Comp )
416  {
417  return S_integrate( __Comp, mpl::int_<nDim>() );
418  }
419 
420 
421 
422 
424 
425 private:
426 
430  matrix_type
431  evaluate( points_type const& __pts, mpl::int_<1> ) const
432  {
433  matrix_type m ( nOrder+1 ,__pts.size2() );
434  ublas::vector<value_type> pts( ublas::row( __pts,0 ) );
435 
436  for ( uint32_type i = 0; i < m.size2(); ++i )
437  m( 0,i ) = 1;
438 
439  for ( uint32_type i = 1; i < m.size1(); ++i )
440  ublas::row( m,i ) = ublas::element_prod( ublas::row( m,i-1 ), pts );
441 
442  return m;
443  }
444 
449  matrix_type
450  evaluate( ublas::vector<value_type> const& __pts ) const
451  {
452  matrix_type m ( nOrder+1 ,__pts.size() );
453 
454  for ( uint32_type i = 0; i < m.size2(); ++i )
455  m( 0,i ) = 1;
456 
457  for ( uint32_type i = 1; i < m.size1(); ++i )
458  ublas::row( m,i ) = ublas::element_prod( ublas::row( m,i-1 ), __pts );
459 
460  return m;
461  }
462 
466  template<typename AE>
467  vector_matrix_type
468  derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<1> ) const
469  {
470  FEELPP_ASSERT( __pts().size1() == 1 )( __pts().size1() )( __pts().size2() ).error( "invalid points" );
471  vector_matrix_type D( 1 );
472 
473  D[0].resize( nOrder+1, __pts().size2() );
474  matrix_type E ( evaluate( __pts ) );
475 
476  for ( uint32_type i = 0; i < E.size2(); ++i )
477  ublas::row( D[0], 0 )( i ) = value_type( 0.0 );
478 
479  for ( uint32_type i = 1; i < E.size1(); ++i )
480  ublas::row( D[0], i ) = value_type( i )*ublas::row( E,i-1 );
481 
482  return D;
483  }
484 
488  vector_matrix_type
489  derivate( ublas::vector<value_type> const& __pts ) const
490  {
491  vector_matrix_type D( 1 );
492 
493  D[0].resize( nOrder+1, __pts.size() );
494  matrix_type E ( evaluate( __pts ) );
495 
496  for ( uint32_type i = 0; i < E.size2(); ++i )
497  ublas::row( D[0], 0 )( i ) = value_type( 0.0 );
498 
499  for ( int32_type i = 1; i < E.size1(); ++i )
500  ublas::row( D[0], i ) = value_type( i )*ublas::row( E,i-1 );
501 
502  return D;
503  }
504 
509  matrix_type evaluate( points_type const& __pts, mpl::int_<2> ) const;
510 
515  template<typename AE>
516  vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<2> ) const;
517 
522  matrix_type evaluate( points_type const& __pts, mpl::int_<3> ) const;
523 
528  template<typename AE>
529  vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<3> ) const;
530 
535  value_type SP_integrate( ublas::matrix<value_type> const& __Comp, mpl::int_<1> )
536  {
537  value_type res( 0.0 );
538 
539  if ( __Comp.size1() == 1 ) // Only 1 function to integrate : \f$ \int u \f$
540  {
541  for ( int ki = 0; ki <__Comp.size2(); ki += 2 )
542  {
543  res+= 2.0*__Comp( 0,ki ) / value_type( ki+1 );
544  }
545  }
546 
547  else if ( __Comp.size1() == 2 ) // 2 functions to integrate : \f$ \int uv \f$
548  {
549  for ( int k1 = 0; k1 <__Comp.size2(); ++k1 )
550  for ( int k2 = 0; k2 <__Comp.size2(); ++k2 )
551  {
552  if ( ( k1+k2 )%2 == 0 )
553  res+= 2.0 * __Comp( 0,k1 ) * __Comp( 1,k2 ) / value_type( k1+k2+1 );
554  }
555  }
556 
557  else if ( __Comp.size1() == 3 ) // 3 functions to integrate : \f$ \int uvw \f$
558  {
559  for ( int k1 = 0; k1 <__Comp.size2(); ++k1 )
560  for ( int k2 = 0; k2 <__Comp.size2(); ++k2 )
561  for ( int k3 = 0; k3 <__Comp.size2(); ++k3 )
562  {
563  if ( ( k1+k2+k3 )%2 == 0 )
564  res+= 2.0 * __Comp( 0,k1 ) * __Comp( 1,k2 )* __Comp( 2,k3 ) / value_type( k1+k2+k3+1 );
565  }
566  }
567 
568  return res;
569  }
570 
571 
576  value_type SP_integrate( ublas::matrix<value_type> const& __Comp, mpl::int_<2> )
577  {
578  value_type res( 0.0 );
579 
580  if ( __Comp.size1() == 1 ) // Only 1 function to integrate : \f$ \int u \f$
581  {
582  for ( int ki = 0; ki <__Comp.size2(); ++ki )
583  {
584  int a( this->coord()( ki,0 ) );
585  int b( this->coord()( ki,1 ) );
586 
587  if ( ( a%2 == 0 ) && ( b%2 == 0 ) )
588  res+= 4.0*__Comp( 0,ki ) / value_type( a+1 ) / value_type( b+1 );
589  }
590  }
591 
592  else if ( __Comp.size1() == 2 ) // 2 functions to integrate : \f$ \int uv \f$
593  {
594  for ( int k1 = 0; k1 <__Comp.size2(); ++k1 )
595  for ( int k2 = 0; k2 <__Comp.size2(); ++k2 )
596  {
597  int a( this->coord()( k1,0 ) + this->coord()( k2,0 ) );
598  int b( this->coord()( k2,1 ) + this->coord()( k2,1 ) );
599 
600  if ( ( a%2 == 0 ) && ( b%2 == 0 ) )
601  res+= 4.0 * __Comp( 0,k1 ) * __Comp( 1,k2 ) / value_type( a+1 ) / value_type( b+1 );
602  }
603  }
604 
605  else if ( __Comp.size1() == 3 ) // 3 functions to integrate : \f$ \int uvw \f$
606  {
607  for ( int k1 = 0; k1 <__Comp.size2(); ++k1 )
608  for ( int k2 = 0; k2 <__Comp.size2(); ++k2 )
609  for ( int k3 = 0; k3 <__Comp.size2(); ++k3 )
610  {
611  int a( this->coord()( k1,0 ) + this->coord()( k2,0 ) + this->coord()( k3,0 ) );
612  int b( this->coord()( k1,1 ) + this->coord()( k2,1 ) + this->coord()( k3,1 ) );
613 
614  if ( ( a%2 == 0 ) && ( b%2 == 0 ) )
615  res+= 4.0 * __Comp( 0,k1 ) * __Comp( 1,k2 ) * __Comp( 2,k3 ) / value_type( a+1 ) / value_type( b+1 );
616  }
617  }
618 
619  return res;
620  }
625  value_type SP_integrate( ublas::matrix<value_type> const& __Comp, mpl::int_<3> )
626  {
627  value_type res( 0.0 );
628 
629  if ( __Comp.size1() == 1 ) // Only 1 function to integrate : \f$ \int u \f$
630  {
631  for ( int ki = 0; ki <__Comp.size2(); ++ki )
632  {
633  int a( this->coord()( ki,0 ) );
634  int b( this->coord()( ki,1 ) );
635  int c( this->coord()( ki,2 ) );
636 
637  if ( ( a%2 == 0 ) && ( b%2 == 0 ) && ( c%2 == 0 ) )
638  res+= 8.0*__Comp( 0,ki ) / value_type( a+1 ) / value_type( b+1 ) / value_type( c+1 );
639  }
640  }
641 
642  else if ( __Comp.size1() == 2 ) // 2 functions to integrate : \f$ \int uv \f$
643  {
644  for ( int k1 = 0; k1 <__Comp.size2(); ++k1 )
645  for ( int k2 = 0; k2 <__Comp.size2(); ++k2 )
646  {
647  int a( this->coord()( k1,0 ) + this->coord()( k2,0 ) );
648  int b( this->coord()( k1,1 ) + this->coord()( k2,1 ) );
649  int c( this->coord()( k1,2 ) + this->coord()( k2,2 ) );
650 
651  if ( ( a%2 == 0 ) && ( b%2 == 0 ) && ( c%2 == 0 ) )
652  res+= 8.0*__Comp( 0,k1 )* __Comp( 1,k2 ) / value_type( a+1 ) / value_type( b+1 ) / value_type( c+1 );
653  }
654  }
655 
656  else if ( __Comp.size1() == 3 ) // 3 functions to integrate : \f$ \int uvw \f$
657  {
658  for ( int k1 = 0; k1 <__Comp.size2(); ++k1 )
659  for ( int k2 = 0; k2 <__Comp.size2(); ++k2 )
660  for ( int k3 = 0; k3 <__Comp.size2(); ++k3 )
661  {
662  int a( this->coord()( k1,0 ) + this->coord()( k2,0 ) + this->coord()( k3,0 ) );
663  int b( this->coord()( k1,1 ) + this->coord()( k2,1 ) + this->coord()( k3,1 ) );
664  int c( this->coord()( k1,2 ) + this->coord()( k2,2 ) + this->coord()( k3,2 ) );
665 
666  if ( ( a%2 == 0 ) && ( b%2 == 0 ) && ( c%2 == 0 ) )
667  res+= 8.0*__Comp( 0,k1 )* __Comp( 1,k2 )* __Comp( 2,k3 )/value_type( a+1 )/value_type( b+1 )/value_type( c+1 );
668  }
669  }
670 
671  return res;
672  }
673 
678  value_type S_integrate( ublas::matrix<value_type> const& __Comp, mpl::int_<1> )
679  {
680  return SP_integrate( __Comp, mpl::int_<nDim>() );
681  }
682 
687  value_type S_integrate( ublas::matrix<value_type> const& __Comp, mpl::int_<2> )
688  {
689  value_type res( 0.0 );
690 
691  if ( __Comp.size1() == 1 ) // Only 1 function to integrate : \f$ \int u \f$
692  {
693  for ( uint32_type ki = 0; ki <__Comp.size2(); ++ki )
694  {
695  int i( this->coord()( ki,0 ) );
696  int j( this->coord()( ki,1 ) );
697 
698  // i and j parity
699 
700  int pi = i%2;
701  int pj = j%2;
702 
703  if ( ( i%2 == 0 ) || ( ( i+j )%2 != 0 ) )
704  {
705  if ( pi==0 && pj==0 ) // (-1)/(j+1)*(-2/(i+1))
706  res+= __Comp( 0,ki )/value_type( j+1 ) * ( 2.0/value_type( i+1 ) ) ;
707 
708  else if ( pi==0 && pj==1 ) // 1/(j+1)*(2/(i+j+1)-2/(i+1))
709  res+= __Comp( 0,ki )/value_type( j+1 ) * ( 2.0/value_type( i+j+2 )-2.0/value_type( i+1 ) ) ;
710 
711  else if ( pi==1 && pj==0 ) // (-1)/(j+1)*(2/(i+j+1))
712  res-= __Comp( 0,ki )/value_type( j+1 ) * ( 2.0/value_type( i+j+2 ) ) ;
713  }
714  }
715  }
716 
717  else if ( __Comp.size1() == 2 ) // 2 function2 to integrate : \f$ \int uv \f$
718  {
719  for ( uint32_type k1 = 0; k1 <__Comp.size2(); ++k1 )
720  for ( uint32_type k2 = 0; k2 <__Comp.size2(); ++k2 )
721  {
722  int i( this->coord()( k1,0 ) + this->coord()( k2,0 ) );
723  int j( this->coord()( k1,1 ) + this->coord()( k2,1 ) );
724 
725  // i and j parity
726 
727  int pi = i%2;
728  int pj = j%2;
729 
730  if ( ( i%2 == 0 ) || ( ( i+j )%2 != 0 ) )
731  {
732  if ( pi==0 && pj==0 ) // (-1)/(j+1)*(-2/(i+1))
733  res+= __Comp( 0,k1 ) * __Comp( 1,k2 )/value_type( j+1 ) * ( 2.0/value_type( i+1 ) ) ;
734 
735  else if ( pi==0 && pj==1 ) // 1/(j+1)*(2/(i+j+1)-2/(i+1))
736  res+= __Comp( 0,k1 ) * __Comp( 1,k2 )/value_type( j+1 ) * ( 2.0/value_type( i+j+2 )-2.0/value_type( i+1 ) ) ;
737 
738  else if ( pi==1 && pj==0 ) // (-1)/(j+1)*(2/(i+j+1))
739  res-= __Comp( 0,k1 ) * __Comp( 1,k2 )/value_type( j+1 ) * ( 2.0/value_type( i+j+2 ) ) ;
740  }
741  }
742  }
743 
744  else if ( __Comp.size1() == 3 ) // 3 functions to integrate : \f$ \int uvw \f$
745  {
746  for ( uint32_type k1 = 0; k1 <__Comp.size2(); ++k1 )
747  for ( uint32_type k2 = 0; k2 <__Comp.size2(); ++k2 )
748  for ( uint32_type k3 = 0; k3 <__Comp.size2(); ++k3 )
749  {
750  int i( this->coord()( k1,0 ) + this->coord()( k2,0 )+ this->coord()( k3,0 ) );
751  int j( this->coord()( k1,1 ) + this->coord()( k2,1 )+ this->coord()( k3,1 ) );
752 
753  // i and j parity
754 
755  int pi = i%2;
756  int pj = j%2;
757 
758  if ( ( i%2 == 0 ) || ( ( i+j )%2 != 0 ) )
759  {
760  if ( pi==0 && pj==0 ) // (-1)/(j+1)*(-2/(i+1))
761  res+= __Comp( 0,k1 ) * __Comp( 1,k2 ) * __Comp( 2,k3 )/value_type( j+1 ) * ( 2.0/value_type( i+1 ) ) ;
762 
763  else if ( pi==0 && pj==1 ) // 1/(j+1)*(2/(i+j+1)-2/(i+1))
764  res+= __Comp( 0,k1 ) * __Comp( 1,k2 ) * __Comp( 2,k3 )/value_type( j+1 ) * ( 2.0/value_type( i+j+2 )-2.0/value_type( i+1 ) ) ;
765 
766  else if ( pi==1 && pj==0 ) // (-1)/(j+1)*(2/(i+j+1))
767  res-= __Comp( 0,k1 ) * __Comp( 1,k2 ) * __Comp( 2,k3 )/value_type( j+1 ) * ( 2.0/value_type( i+j+2 ) ) ;
768  }
769  }
770  }
771 
772  return res;
773  }
774 
779  value_type S_integrate( ublas::matrix<value_type> const& __Comp, mpl::int_<3> )
780  {
781  return 0.0;
782  }
783 
784 
785 
786 private:
787  reference_convex_type M_refconvex;
788  points_type M_pts;
789 
795  ublas::matrix<int16_type> M_coord;
796 
800  std::vector<matrix_type> M_D;
801 };
802 
803 template<uint16_type Dim,
804  uint16_type Degree,
805  typename TheConvex,
806  typename T,
807  template<class> class StoragePolicy>
808 typename Moment<Dim, Degree,TheConvex, T, StoragePolicy>::matrix_type
809 Moment<Dim, Degree,TheConvex, T, StoragePolicy>::evaluate( points_type const& __pts, mpl::int_<2> ) const
810 {
811  matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
812 
813  ublas::vector<value_type> pts_x = ublas::row( __pts, 0 );
814  ublas::vector<value_type> pts_y = ublas::row( __pts, 1 );
815 
816  matrix_type E_x( evaluate( pts_x ) );
817  matrix_type E_y( evaluate( pts_y ) );
818 
819  // std::cout << "\nDimensions E_x = " << E_x.size1() << " x " << E_x.size2() << std::endl;
820  // std::cout << "\nDimensions E_y = " << E_y.size1() << " x " << E_y.size2() << std::endl;
821  // std::cout << "Order = " << nOrder+1 << std::endl;
822 
823  uint32_type G_i=0; // Global indice
824 
825  if ( convex_is_hypercube )
826  {
827  for ( int32_type n = 0; n < nOrder+1; ++n )
828  for ( int32_type i = 0; i < nOrder+1; ++i )
829 
830  {
831  ublas::row( res, G_i )=ublas::element_prod( ublas::row( E_x, i ) , ublas::row( E_y, n ) );
832  ++G_i;
833  }
834  }
835  else if ( convex_is_simplex )
836  {
837  for ( int32_type n = 0; n < nOrder+1; ++n )
838  for ( int32_type i = n; i >= 0; --i )
839  {
840  ublas::row( res, G_i )=ublas::element_prod( ublas::row( E_x, i ) , ublas::row( E_y, n-i ) );
841  ++G_i;
842  }
843  }
844 
845  return res;
846 }
847 
848 template<uint16_type Dim,
849  uint16_type Degree,
850  typename TheConvex,
851  typename T,
852  template<class> class StoragePolicy>
853 template<typename AE>
854 typename Moment<Dim, Degree, TheConvex, T, StoragePolicy>::vector_matrix_type
855 Moment<Dim, Degree, TheConvex, T, StoragePolicy>::derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<2> ) const
856 {
857  vector_matrix_type res( 2 );
858  res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
859  res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
860 
861  ublas::vector<value_type> pts_x = ublas::row( __pts(), 0 );
862  ublas::vector<value_type> pts_y = ublas::row( __pts(), 1 );
863 
864  matrix_type E_x( evaluate( pts_x ) );
865  matrix_type E_y( evaluate( pts_y ) );
866 
867  vector_matrix_type D_x( derivate( pts_x ) );
868  vector_matrix_type D_y( derivate( pts_y ) );
869 
870  uint32_type G_i=0;
871 
872  if ( convex_is_simplex )
873  {
874  for ( int32_type n = 0; n < nOrder+1; ++n )
875  for ( int32_type i = n; i >= 0; --i )
876  {
877  ublas::row( res[0], G_i )=ublas::element_prod( ublas::row( D_x[0], i ) , ublas::row( E_y, n-i ) );
878  ublas::row( res[1], G_i )=ublas::element_prod( ublas::row( E_x, i ) , ublas::row( D_y[0], n-i ) );
879  ++G_i;
880  }
881  }
882  else if ( convex_is_hypercube )
883  {
884  for ( int32_type n = 0; n < nOrder+1; ++n )
885  for ( int32_type i = 0; i < nOrder+1; ++i )
886  {
887  ublas::row( res[0], G_i )=ublas::element_prod( ublas::row( D_x[0], i ) , ublas::row( E_y, n ) );
888  ublas::row( res[1], G_i )=ublas::element_prod( ublas::row( E_x, i ) , ublas::row( D_y[0], n ) );
889  ++G_i;
890  }
891  }
892 
893  return res;
894 }
895 
896 template<uint16_type Dim,
897  uint16_type Degree,
898  typename TheConvex,
899  typename T,
900  template<class> class StoragePolicy>
901 typename Moment<Dim, Degree, TheConvex, T, StoragePolicy>::matrix_type
902 Moment<Dim, Degree, TheConvex, T, StoragePolicy>::evaluate( points_type const& __pts, mpl::int_<3> ) const
903 {
904  matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
905 
906  FEELPP_ASSERT( __pts.size1() == 3 )( __pts.size1() ).error( "invalid space dimension" );
907 
908  ublas::vector<value_type> pts_x = ublas::row( __pts, 0 );
909  ublas::vector<value_type> pts_y = ublas::row( __pts, 1 );
910  ublas::vector<value_type> pts_z = ublas::row( __pts, 2 );
911 
912  matrix_type E_x( evaluate( pts_x ) );
913  matrix_type E_y( evaluate( pts_y ) );
914  matrix_type E_z( evaluate( pts_z ) );
915 
916  uint32_type G_i=0; // Global indice
917 
918  for ( int32_type n = 0; n < nOrder+1; ++n )
919  for ( int32_type i = n; i >= 0; --i ) // x id
920  for ( int32_type j = n-i; j >= 0 ; --j )
921  {
922  ublas::vector<value_type> tmp( ublas::element_prod( ublas::row( E_x, i ) , ublas::row( E_y, j ) ) );
923  ublas::row( res, G_i )=ublas::element_prod( tmp , ublas::row( E_z, n-i-j ) );
924  ++G_i;
925  }
926 
927  return res;
928 }
929 
930 template<uint16_type Dim,
931  uint16_type Degree,
932  typename TheConvex,
933  typename T,
934  template<class> class StoragePolicy>
935 template<typename AE>
936 typename Moment<Dim, Degree, TheConvex, T, StoragePolicy>::vector_matrix_type
937 Moment<Dim, Degree, TheConvex, T, StoragePolicy>::derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<3> ) const
938 {
939  vector_matrix_type res( 3 );
940  res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
941  res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
942  res[2].resize( convex_type::polyDims( nOrder ), __pts().size2() );
943 
944  FEELPP_ASSERT( __pts().size1() == 3 )( __pts().size1() ).error( "invalid space dimension" );
945 
946  ublas::vector<value_type> pts_x = ublas::row( __pts(), 0 );
947  ublas::vector<value_type> pts_y = ublas::row( __pts(), 1 );
948  ublas::vector<value_type> pts_z = ublas::row( __pts(), 2 );
949 
950  matrix_type E_x( evaluate( pts_x ) );
951  matrix_type E_y( evaluate( pts_y ) );
952  matrix_type E_z( evaluate( pts_z ) );
953 
954  vector_matrix_type D_x( derivate( pts_x ) );
955  vector_matrix_type D_y( derivate( pts_y ) );
956  vector_matrix_type D_z( derivate( pts_z ) );
957 
958  uint32_type G_i=0;
959 
960  for ( int32_type n = 0; n < nOrder+1; ++n )
961  for ( int32_type i = n; i >= 0; --i ) // x id
962  for ( int32_type j = n-i; j >=0 ; --j )
963  {
964  ublas::vector<value_type> tmp( ublas::element_prod( ublas::row( D_x[0], i ) , ublas::row( E_y, j ) ) );
965  ublas::row( res[0], G_i )=ublas::element_prod( tmp , ublas::row( E_z, n-i-j ) );
966 
967  tmp = ublas::element_prod( ublas::row( E_x, i ) , ublas::row( D_y[0], j ) );
968  ublas::row( res[1], G_i )=ublas::element_prod( tmp , ublas::row( E_z, n-i-j ) );
969 
970  tmp = ublas::element_prod( ublas::row( E_x, i ) , ublas::row( E_y, j ) );
971  ublas::row( res[2], G_i )=ublas::element_prod( tmp , ublas::row( D_z[0], n-i-j ) );
972 
973  ++G_i;
974  }
975 
976  return res;
977 }
978 
979 namespace fem
980 {
982 namespace detail
983 {
992 template<uint16_type Dim,
993  uint16_type Order,
994  uint16_type RealDim,
995  template<uint16_type> class PolySetType = Scalar,
996  typename T = double,
997  template<uint16_type,uint16_type,uint16_type> class Convex = Simplex>
998 class MomentPolynomialSet
999  :
1000 public PolynomialSet<Moment<Dim, Order, Convex<Dim,1,Dim>, T, StorageUBlas>, PolySetType >
1001 {
1002  typedef PolynomialSet<Moment<Dim, Order, Convex<Dim,1,Dim>, T, StorageUBlas>, PolySetType > super;
1003 public:
1004 
1005  static const uint16_type nDim = Dim;
1006  static const uint16_type nOrder = Order;
1007  static const uint16_type nRealDim = RealDim;
1008  static const bool isTransformationEquivalent = true;
1009  typedef MomentPolynomialSet<Dim, Order,RealDim, PolySetType, T, Simplex> self_type;
1010  typedef self_type component_basis_type;
1011 
1012  typedef typename super::polyset_type polyset_type;
1013  static const bool is_tensor2 = polyset_type::is_tensor2;
1014  static const bool is_vectorial = polyset_type::is_vectorial;
1015  static const bool is_scalar = polyset_type::is_scalar;
1016  static const bool is_continuous = false;
1017  static const bool is_modal = true;
1018  static const uint16_type nComponents = polyset_type::nComponents;
1019  static const bool is_product = true;
1020  static const bool isContinuous = false;
1021  typedef Discontinuous continuity_type;
1022 
1023  typedef typename super::component_type component_type;
1024 
1025  typedef T value_type;
1026  typedef Moment<Dim, Order, Convex<Dim,1,Dim>, T, StorageUBlas> basis_type;
1027  typedef Convex<Dim, Order, /*RealDim*/Dim> convex_type;
1028  template<int O>
1029  struct convex
1030  {
1031  typedef Convex<Dim, O, /*RealDim*/Dim> type;
1032  };
1033  typedef Reference<convex_type, nDim, nOrder, nDim/*nRealDim*/, value_type> reference_convex_type;
1034 
1035  typedef typename super::polynomial_type polynomial_type;
1036 
1038  static const uint16_type nDofPerVertex = convex_type::nbPtsPerVertex;
1040  static const uint16_type nDofPerEdge = convex_type::nbPtsPerEdge;
1042  static const uint16_type nDofPerFace = convex_type::nbPtsPerFace;
1043 
1045  static const uint16_type nDofPerVolume = convex_type::nbPtsPerVolume;
1046 
1047  static const uint16_type nLocalDof = convex_type::numPoints;
1048 
1049  static const uint16_type nDof = nLocalDof;
1050  static const uint16_type nNodes = nDof;
1051  static const uint16_type nDofGrad = super::nDim*nDof;
1052  static const uint16_type nDofHess = super::nDim*super::nDim*nDof;
1053  //typedef typename matrix_node<value_type>::type points_type;
1054  typedef StorageUBlas<value_type> storage_policy;
1055  typedef typename storage_policy::matrix_type matrix_type;
1056  typedef typename storage_policy::vector_matrix_type vector_matrix_type;
1057  typedef typename storage_policy::matrix_node_type matrix_node_type;
1058  typedef typename storage_policy::points_type points_type;
1059  typedef typename storage_policy::node_type node_type;
1060  MomentPolynomialSet()
1061  :
1062  super( basis_type() )
1063 
1064  {
1065  ublas::matrix<value_type> m( ublas::identity_matrix<value_type>( nComponents*convex_type::polyDims( nOrder ) ) );
1066 
1067  if ( is_tensor2 )
1068  std::cout << "[orthonormalpolynomialset] m = " << m << "\n";
1069 
1070  if ( !( ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) -
1071  m ) < 1e-10 ) )
1072  std::cout << "m1=" << m << "\n"
1073  << "m2=" << polyset_type::toMatrix( polyset_type::toType( m ) ) << "\n"
1074  << ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) - m ) << "\n";
1075 
1076  FEELPP_ASSERT( ublas::norm_frobenius( polyset_type::toMatrix( polyset_type::toType( m ) ) -
1077  m ) < 1e-10 )( m ).warn ( "invalid transformation" );
1078  this->setCoefficient( polyset_type::toType( m ), true );
1079  }
1080 
1081  MomentPolynomialSet<Dim, Order, RealDim, Scalar,T, Simplex > toScalar() const
1082  {
1083  return MomentPolynomialSet<Dim, Order, RealDim, Scalar,T, Simplex >();
1084  }
1085 
1089  std::string familyName() const
1090  {
1091  return "dubiner";
1092  }
1093 
1094 
1095  points_type points() const
1096  {
1097  return points_type();
1098  }
1099  points_type points( int f ) const
1100  {
1101  return points_type();
1102  }
1103 };
1104 
1105 template<uint16_type Dim,
1106  uint16_type Order,
1107  uint16_type RealDim,
1108  template<uint16_type> class PolySetType,
1109  typename T,
1110  template<uint16_type,uint16_type,uint16_type> class Convex>
1111 const uint16_type MomentPolynomialSet<Dim, Order, RealDim, PolySetType,T, Convex>::nLocalDof;
1112 
1113 } // detail
1115 
1116 template<uint16_type Order,
1117  template<uint16_type Dim> class PolySetType = Scalar>
1118 class MomentPolynomialSet
1119 {
1120 public:
1121  template<uint16_type N,
1122  typename T = double,
1123  typename Convex = Simplex<N> >
1124  struct apply
1125  {
1126  typedef typename mpl::if_<mpl::bool_<Convex::is_simplex>,
1127  mpl::identity<detail::MomentPolynomialSet<N,Order,N,PolySetType,T,Simplex> >,
1128  mpl::identity<detail::MomentPolynomialSet<N,Order,N,PolySetType,T,Hypercube> > >::type::type result_type;
1129  typedef result_type type;
1130  };
1131 
1132  typedef MomentPolynomialSet<Order,Scalar> component_basis_type;
1133 };
1134 }
1135 
1136 
1137 }
1138 #endif /* __Moment_H */

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