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
dubiner.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) 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 __Dubiner_H
31 #define __Dubiner_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>
46 
47 namespace Feel
48 {
49 
50 template<uint16_type Dim,
51  uint16_type RealDim,
52  uint16_type Degree,
53  typename NormalizationPolicy,
54  typename T,
55  template<class> class StoragePolicy>
56 class Dubiner;
57 
58 
59 template<uint16_type Dim,
60  uint16_type RealDim,
61  uint16_type Degree,
62  typename NormalizationPolicy = Normalized<true>,
63  typename T = double,
64  template<class> class StoragePolicy = StorageUBlas>
65 struct DubinerTraits
66 {
67  static const uint16_type nDim = Dim;
68  static const uint16_type nRealDim = RealDim;
69  static const uint16_type nOrder = Degree;
70  static const uint16_type nConvexOrderDiff = nDim+nOrder+1;
71  static const bool is_normalized = NormalizationPolicy::is_normalized;
72 
76 
77  /*
78  * numerical type
79  */
80  typedef T value_type;
81 
82  template<uint16_type order, typename V = value_type>
83  struct Convex
84  {
85  typedef Simplex<nDim, order, nDim/*nRealDim*/> type;
86  typedef Reference<Simplex<nDim, order, nDim/*nRealDim*/>, nDim, order, nDim/*nRealDim*/, V> reference_type;
87  };
88 
89  template<typename NewT>
90  struct ChangeValueType
91  {
92  typedef Dubiner<Dim, RealDim, Degree, NormalizationPolicy, NewT, StoragePolicy> type;
93  typedef DubinerTraits<Dim, RealDim, Degree, NormalizationPolicy, NewT, StoragePolicy> traits_type;
94  };
95 
96  template<uint16_type NewOrder>
97  struct ChangeOrder
98  {
99  typedef Dubiner<Dim, RealDim, NewOrder, NormalizationPolicy, T, StoragePolicy> type;
100  typedef DubinerTraits<Dim, RealDim, NewOrder, NormalizationPolicy, T, StoragePolicy> traits_type;
101  };
102 
103  /*
104  * Geometry where the polynomials are defined and constructed
105  */
106  typedef typename Convex<nOrder>::type convex_type;
107  typedef typename Convex<nOrder>::reference_type reference_convex_type;
108 
109  typedef typename Convex<nConvexOrderDiff>::type diff_convex_type;
110  typedef typename Convex<nConvexOrderDiff>::reference_type diff_reference_convex_type;
111 
112  typedef typename mpl::if_<mpl::equal_to<mpl::int_<nDim>, mpl::int_<2> >,
113  mpl::identity<PointSetWarpBlend<diff_convex_type, nConvexOrderDiff, value_type> >,
114  mpl::identity<PointSetEquiSpaced<diff_convex_type, nConvexOrderDiff, value_type> > >::type::type diff_pointset_type;
115 
116  /*
117  * storage policy
118  */
119  typedef StoragePolicy<value_type> storage_policy;
120  typedef typename storage_policy::matrix_type matrix_type;
121  typedef typename storage_policy::vector_matrix_type vector_matrix_type;
122  typedef typename storage_policy::matrix_node_type matrix_node_type;
123  typedef typename storage_policy::points_type points_type;
124  typedef typename storage_policy::node_type node_type;
125 }; // class DubinerTraits
126 
127 template<int D, int O>
128 struct DubinerTag
129 {
130  static const int Dim = D;
131  static const int Order = O;
132 };
161 template<uint16_type Dim,
162  uint16_type RealDim,
163  uint16_type Degree,
164  typename NormalizationPolicy = Normalized<true>,
165  typename T = double,
166  template<class> class StoragePolicy = StorageUBlas>
167 class Dubiner
168 
169 {
170 
171 public:
172  typedef DubinerTraits<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy> traits_type;
173 
174  static const uint16_type nDim = traits_type::nDim;
175  static const uint16_type nRealDim = traits_type::nRealDim;
176  static const uint16_type nOrder = traits_type::nOrder;
177  static const uint16_type nConvexOrderDiff = traits_type::nConvexOrderDiff;
178  static const bool is_normalized = traits_type::is_normalized;
179  static const bool isTransformationEquivalent = true;
180  static const bool isContinuous = false;
181  static const bool is_product = true;
182  typedef Discontinuous continuity_type;
183 
187 
188  typedef Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy> self_type;
189 
190  /*
191  * for now can be used only as a basis but we might be interested
192  * to have then expressed in other basis like in the moment basis
193  */
194  typedef self_type basis_type;
195 
196  typedef typename traits_type::value_type value_type;
197 
198  /*
199  * Geometry where the polynomials are defined and constructed
200  */
201  typedef typename traits_type::convex_type convex_type;
202  typedef typename traits_type::reference_convex_type reference_convex_type;
203 
204  typedef typename traits_type::diff_pointset_type diff_pointset_type;
205 
206  /*
207  * storage policy
208  */
209  typedef typename traits_type::storage_policy storage_policy;
210  typedef typename traits_type::matrix_type matrix_type;
211  typedef typename traits_type::vector_matrix_type vector_matrix_type;
212  typedef typename traits_type::matrix_node_type matrix_node_type;
213  typedef typename traits_type::points_type points_type;
214  typedef typename traits_type::node_type node_type;
215 
217 
221 
222  Dubiner()
223  :
224  M_refconvex(),
225  M_pts( M_refconvex.makePoints( nDim, 0 ) )
226  {
227  this->initDerivation();
228  }
229  Dubiner( Dubiner const & d )
230  :
231  M_refconvex(),
232  M_pts( d.M_pts )
233  {
234 
235  }
236 
237  ~Dubiner()
238  {}
239 
241 
245 
246  self_type const& operator=( self_type const& d )
247  {
248  if ( this != &d )
249  {
250  M_pts = d.M_pts;
251  _S_D = d._S_D;
252  }
253 
254  return *this;
255  }
256 
257  //ublas::matrix_column<matrix_type const> operator()( node_type const& pt ) const
258  matrix_type operator()( node_type const& pt ) const
259  {
260  points_type pts( pt.size(), 1 );
261  ublas::column( pts, 0 ) = pt;
262  return evaluate( pts );
263  }
264 
265  matrix_type operator()( points_type const& pts ) const
266  {
267  return evaluate( pts );
268  }
269 
271 
275 
279  size_type size() const
280  {
281  return convex_type::polyDims( nOrder );
282  }
283 
288  uint16_type degree() const
289  {
290  return nOrder;
291  }
292 
296  self_type const& basis() const
297  {
298  return *this;
299  }
300 
305  bool isNormalized() const
306  {
307  return is_normalized;
308  }
309 
313  std::string familyName() const
314  {
315  return "dubiner";
316  }
317 
319 
323 
324 
326 
330 
331 
332 
342  matrix_type coeff() const
343  {
344 #if 0
345  std::cout << "[Dubiner::coeff] coeff = "
346  << ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() )
347  << "\n";
348 #endif
349  return ublas::identity_matrix<value_type>( reference_convex_type::polyDims( nOrder ), M_pts.size2() );
350  }
351 
352 
358  static matrix_type evaluate( points_type const& __pts )
359  {
360  return evaluate( __pts, mpl::int_<nDim>() );
361  }
362 
363  template<typename AE>
364  static vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts )
365  {
366  return derivate( __pts, mpl::int_<nDim>() );
367  }
368 
375  static matrix_type const& d( uint16_type i )
376  {
377  return _S_D[i];
378  }
379 
386  static matrix_type const& derivate( uint16_type i )
387  {
388  return _S_D[i];
389  }
390 
392 
393 private:
394 private:
395 
400  static matrix_type
401  evaluate( points_type const& __pts, mpl::int_<1> )
402  {
403  matrix_type m ( JacobiBatchEvaluation<nOrder,value_type>( 0.0, 0.0, ublas::row( __pts, 0 ) ) );
404 
405  if ( is_normalized )
406  {
407  for ( uint16_type i = 0; i < m.size1(); ++i )
408  ublas::row( m, i ) *= math::sqrt( value_type( i )+0.5 );
409  }
410 
411  return m;
412  }
413 
418  template<typename AE>
419  static vector_matrix_type
420  derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<1> )
421  {
422  FEELPP_ASSERT( __pts().size1() == 1 )( __pts().size1() )( __pts().size2() ).warn( "invalid points" );
423  // VLOG(1) << "Expansion::derivate<1>] number of points " << __pts().size2() << "\n";
424 
425  vector_matrix_type D( 1 );
426  D[0].resize( nOrder+1, __pts().size2() );
427  D[0] = JacobiBatchDerivation<nOrder,value_type>( 0.0, 0.0, ublas::row( __pts(),0 ) );
428 
429  if ( is_normalized )
430  for ( uint16_type i = 0; i < nOrder+1; ++i )
431  ublas::row( D[0], i ) *= math::sqrt( value_type( i )+0.5 );
432 
433  return D;
434  }
435 
440  static matrix_type evaluate( points_type const& __pts, mpl::int_<2> );
441 
446  template<typename AE>
447  static vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<2> );
448 
453  static matrix_type evaluate( points_type const& __pts, mpl::int_<3> );
454 
459  template<typename AE>
460  static vector_matrix_type derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<3> );
461 
462  static void initDerivation();
463 private:
464  reference_convex_type M_refconvex;
465  points_type M_pts;
466 
471  static bool _S_has_derivation;
472 
477  static std::vector<matrix_type> _S_D;
478 
479 }; // class Dubiner
480 
481 template<uint16_type Dim,
482  uint16_type RealDim,
483  uint16_type Degree,
484  typename NormalizationPolicy,
485  typename T,
486  template<class> class StoragePolicy>
487 bool Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::_S_has_derivation = false;
488 
489 template<uint16_type Dim,
490  uint16_type RealDim,
491  uint16_type Degree,
492  typename NormalizationPolicy,
493  typename T,
494  template<class> class StoragePolicy>
495 std::vector<typename Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::matrix_type>
496 Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::_S_D;
497 
498 template<uint16_type Dim,
499  uint16_type RealDim,
500  uint16_type Degree,
501  typename NormalizationPolicy,
502  typename T,
503  template<class> class StoragePolicy>
504 void
505 Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::initDerivation()
506 {
507 #if 0
508  typedef typename traits_type::convex_type convex_type;
509  typedef typename traits_type::reference_convex_type reference_convex_type;
510 
511  typedef typename traits_type::diff_pointset_type diff_pointset_type;
512 
513  typedef typename traits_type::storage_policy storage_policy;
514  typedef typename traits_type::matrix_type matrix_type;
515  typedef typename traits_type::vector_matrix_type vector_matrix_type;
516  typedef typename traits_type::matrix_node_type matrix_node_type;
517  typedef typename traits_type::points_type points_type;
518  typedef typename traits_type::node_type node_type;
519 #endif // 0
520 
521  if ( _S_has_derivation == false )
522  {
523  _S_has_derivation = true;
524 
525  reference_convex_type refconvex;
526  // constructor pointset for differentiation only in
527  // the interior(1)
528  diff_pointset_type diff_pts( 1 );
529  matrix_type A( evaluate( diff_pts.points() ) );
530 
531 #if 1
532  matrix_type D = ublas::identity_matrix<value_type>( A.size1(), A.size2() );
533  LU<matrix_type> lu( A );
534  matrix_type C = lu.solve( D );
535 
536  vector_matrix_type d ( derivate( diff_pts.points() ) );
537  _S_D.resize( d.size() );
538 
539  for ( size_type i = 0; i < d.size(); ++i )
540  {
541  _S_D[i] = ublas::prod( d[i], C );
542  glas::clean( _S_D[i] );
543  }
544 
545 #endif
546  }
547 }
548 
549 template<uint16_type Dim,
550  uint16_type RealDim,
551  uint16_type Degree,
552  typename NormalizationPolicy,
553  typename T,
554  template<class> class StoragePolicy>
555 typename Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::matrix_type
557 {
558  matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
559 
560  details::etas<TRIANGLE, value_type> etas( __pts );
561  ublas::vector<value_type> eta1s = ublas::row( etas(), 0 );
562  ublas::vector<value_type> eta2s = ublas::row( etas(), 1 );
563 
564  //std::cout << "xis = " << __pts << "\n";
565  //std::cout << "etas = " << etas() << "\n";
566 
567  matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
568  std::vector<matrix_type> bs( nOrder+1 );
569 
570  for ( int i = 0; i < nOrder+1; ++i )
571  {
572  bs[ i ].resize( nOrder-i, eta2s.size() );
573  bs[ i ] = dyna::JacobiBatchEvaluation( nOrder-i, value_type( 2*i+1 ), value_type( 0.0 ), eta2s );
574  }
575 
576 
577  details::scalings<nOrder, T> scalings( eta2s );
578 
579 
580  for ( uint16_type cur = 0, k = 0; k < nOrder+1; ++k )
581  {
582  for ( uint16_type i = 0; i < k+1; ++i,++cur )
583  {
584  uint16_type ii = k-i;
585  uint16_type jj = i;
586 
587  if ( is_normalized )
588  {
589  value_type normalization = math::sqrt( ( value_type( ii )+0.5 )*( value_type( ii+jj )+1.0 ) );
590 
591  for ( uint16_type l = 0; l < as.size2(); ++l )
592  res( cur, l ) = normalization*as( ii,l )*scalings()( ii,l )*bs[ii]( jj,l );
593  }
594 
595  else
596  {
597  for ( uint16_type l = 0; l < as.size2(); ++l )
598  res( cur, l ) = as( ii,l )*scalings()( ii,l )*bs[ii]( jj,l );
599  }
600  }
601  }
602 
603  return res;
604 }
605 
606 template<uint16_type Dim,
607  uint16_type RealDim,
608  uint16_type Degree,
609  typename NormalizationPolicy,
610  typename T,
611  template<class> class StoragePolicy>
612 template<typename AE>
613 typename Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::vector_matrix_type
614 Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<2> )
615 {
616  vector_matrix_type res( 2 );
617  res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
618  res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
619 
620  // transform wrapped coordinates in cartesian coordinates
621  details::etas<TRIANGLE, value_type> etas( __pts );
622  ublas::vector<value_type> eta1s = ublas::row( etas(), 0 );
623  ublas::vector<value_type> eta2s = ublas::row( etas(), 1 );
624 
625  //std::cout << "xis = " << __pts << "\n";
626  //std::cout << "etas = " << etas() << "\n";
627 
628  // evaluate Dubiner polynomials components
629  matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
630  matrix_type das( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
631  //std::cout << "das= " << das << "\n";
632  std::vector<matrix_type> bs( nOrder+1 );
633  std::vector<matrix_type> dbs( nOrder+1 );
634 
635  for ( uint16_type i = 0; i < nOrder+1; ++i )
636  {
637  bs[ i ].resize( nOrder-i, eta2s.size() );
638  dbs[ i ].resize( nOrder-i, eta2s.size() );
639  bs[ i ] = dyna::JacobiBatchEvaluation( nOrder-i, value_type( 2*i+1 ), value_type( 0.0 ), eta2s );
640  dbs[ i ] = dyna::JacobiBatchDerivation( nOrder-i, value_type( 2*i+1 ), value_type( 0.0 ), eta2s );
641 
642  //std::cout << "dbs["<< i << "]= " << dbs[i] << "\n";
643  }
644 
645  details::scalings<nOrder, T> scalings( eta2s );
646  //std::cout << "scalings = " << scalings() << "\n";
647  ublas::vector<value_type> one( ublas::scalar_vector<value_type>( eta1s.size(), 1.0 ) );
648  ublas::vector<value_type> tmp( ublas::scalar_vector<value_type>( eta1s.size(), 1.0 ) );
649 
650  // assemble Dubiner polynomials components
651  for ( uint16_type k = 0, cur = 0; k < nOrder+1; ++k )
652  {
653  for ( uint16_type i = 0; i < k+1; ++i, ++cur )
654  {
655  uint16_type ii = k-i;
656  uint16_type jj = i;
657 
658 
659  // x derivation
660  ublas::row( res[0], cur ) = ublas::element_prod( ublas::row( das, ii ),
661  ublas::row( bs[ii], jj ) );
662 
663  if ( ii > 0 )
664  ublas::row( res[0], cur ) = element_prod( ublas::row( res[0], cur ),
665  ublas::row( scalings(), ii-1 ) );
666 
667  // y derivation
668  ublas::row( res[1], cur ) = ublas::element_prod( ublas::row( das, ii ),
669  ublas::row( bs[ii], jj ) );
670  ublas::row( res[1], cur ) = 0.5 * element_prod( ublas::row( res[1], cur ), ( one+eta1s ) );
671 
672  if ( ii > 0 )
673  ublas::row( res[1], cur ) = element_prod( ublas::row( res[1], cur ),
674  ublas::row( scalings(), ii-1 ) );
675 
676  // derivate (1-x)^ii
677  tmp = ublas::element_prod( ublas::row( scalings(), ii ),
678  ublas::row( dbs[ii], jj ) );
679 
680  if ( ii > 0 )
681  tmp -= 0.5 * ii * ublas::element_prod( ublas::row( scalings(), ii-1 ),
682  ublas::row( bs[ii], jj ) );
683 
684  // add contrib to y derivation
685  ublas::row( res[1], cur ) += ublas::element_prod( ublas::row( as, ii ), tmp );
686 
687  // orthonormalize if required
688  if ( is_normalized )
689  {
690  value_type normalization = math::sqrt( ( value_type( ii )+0.5 )*( value_type( ii+jj )+1.0 ) );
691  ublas::row( res[0], cur ) *= normalization;
692  ublas::row( res[1], cur ) *= normalization;
693  }
694  }
695  }
696 
697  return res;
698 }
699 
700 template<uint16_type Dim,
701  uint16_type RealDim,
702  uint16_type Degree,
703  typename NormalizationPolicy,
704  typename T,
705  template<class> class StoragePolicy>
706 typename Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::matrix_type
708 {
709  matrix_type res( convex_type::polyDims( nOrder ), __pts.size2() );
710 
711  FEELPP_ASSERT( __pts.size1() == 3 )( __pts.size1() ).error( "invalid space dimension" );
712 
713  details::etas<TETRAHEDRON, value_type> etas( __pts );
714  ublas::vector<value_type> eta1s = ublas::row( etas(), 0 );
715  ublas::vector<value_type> eta2s = ublas::row( etas(), 1 );
716  ublas::vector<value_type> eta3s = ublas::row( etas(), 2 );
717 
718  //std::cout << "xis = " << __pts << "\n";
719  //std::cout << "etas = " << etas() << "\n";
720 
721  matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
722  std::vector<matrix_type> bs( nOrder+1 );
723  ublas::matrix<matrix_type> cs( nOrder+1, nOrder+1 );
724 
725  for ( int i = 0; i < nOrder+1; ++i )
726  {
727  bs[ i ].resize( nOrder-i, eta2s.size() );
728  bs[ i ] = dyna::JacobiBatchEvaluation( nOrder-i, value_type( 2*i+1 ), value_type( 0.0 ), eta2s );
729 
730  for ( int j = 0; j < nOrder+1-i; ++j )
731  {
732  cs( i, j ).resize( nOrder-i-j, eta3s.size() );
733  cs( i, j ) = dyna::JacobiBatchEvaluation( nOrder-i-j,
734  value_type( 2*( i+j+1 ) ), value_type( 0.0 ), eta3s );
735  }
736  }
737 
738 
739  details::scalings<nOrder, T> scalings2( eta2s );
740  details::scalings<nOrder, T> scalings3( eta3s );
741 
742  for ( uint16_type cur = 0, k = 0; k < nOrder+1; ++k )
743  {
744  for ( uint16_type i = 0; i < k+1; ++i )
745  {
746  for ( uint16_type j = 0; j < k+1-i; ++j,++cur )
747  {
748  uint16_type ii = k-i-j;
749  uint16_type jj = j;
750  uint16_type kk = i;
751 
752 
753  if ( is_normalized )
754  {
755  value_type normalization = math::sqrt( ( value_type( ii )+0.5 )*
756  ( value_type( ii+jj )+1.0 )*
757  ( value_type( ii+jj+kk )+1.5 ) );
758 
759  for ( uint16_type l = 0; l < as.size2(); ++l )
760  res( cur, l ) = normalization*( as( ii,l )*
761  scalings2()( ii,l )*bs[ii]( jj,l )*
762  scalings3()( ii+jj,l )*cs( ii, jj )( kk,l ) );
763  }
764 
765  else
766  {
767  for ( uint16_type l = 0; l < as.size2(); ++l )
768  res( cur, l ) = as( ii,l )*
769  scalings2()( ii,l )*bs[ii]( jj,l )*
770  scalings3()( ii+jj,l )*cs( ii, jj )( kk,l );
771  }
772  }
773  }
774  }
775 
776  return res;
777 }
778 
779 template<uint16_type Dim,
780  uint16_type RealDim,
781  uint16_type Degree,
782  typename NormalizationPolicy,
783  typename T,
784  template<class> class StoragePolicy>
785 template<typename AE>
786 typename Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::vector_matrix_type
787 Dubiner<Dim, RealDim, Degree, NormalizationPolicy, T, StoragePolicy>::derivate( ublas::matrix_expression<AE> const& __pts, mpl::int_<3> )
788 {
789  vector_matrix_type res( 3 );
790  res[0].resize( convex_type::polyDims( nOrder ), __pts().size2() );
791  res[1].resize( convex_type::polyDims( nOrder ), __pts().size2() );
792  res[2].resize( convex_type::polyDims( nOrder ), __pts().size2() );
793 
794  FEELPP_ASSERT( __pts().size1() == 3 )( __pts().size1() ).error( "invalid space dimension" );
795 
796  details::etas<TETRAHEDRON, value_type> etas( __pts );
797  ublas::vector<value_type> eta1s = ublas::row( etas(), 0 );
798  ublas::vector<value_type> eta2s = ublas::row( etas(), 1 );
799  ublas::vector<value_type> eta3s = ublas::row( etas(), 2 );
800 
801  //std::cout << "xis = " << __pts << "\n";
802  //std::cout << "etas = " << etas() << "\n";
803 
804  matrix_type as( JacobiBatchEvaluation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
805  matrix_type das( JacobiBatchDerivation<nOrder, value_type>( 0.0, 0.0, eta1s ) );
806  std::vector<matrix_type> bs( nOrder+1 );
807  std::vector<matrix_type> dbs( nOrder+1 );
808  ublas::matrix<matrix_type> cs( nOrder+1, nOrder+1 );
809  ublas::matrix<matrix_type> dcs( nOrder+1, nOrder+1 );
810 
811  for ( int i = 0; i < nOrder+1; ++i )
812  {
813  bs[ i ].resize( nOrder-i, eta2s.size() );
814  dbs[ i ].resize( nOrder-i, eta2s.size() );
815  bs[ i ] = dyna::JacobiBatchEvaluation( nOrder-i, value_type( 2*i+1 ), value_type( 0.0 ), eta2s );
816  dbs[ i ] = dyna::JacobiBatchDerivation( nOrder-i, value_type( 2*i+1 ), value_type( 0.0 ), eta2s );
817 
818  for ( int j = 0; j < nOrder+1-i; ++j )
819  {
820  cs( i, j ).resize( nOrder-i-j, eta3s.size() );
821  dcs( i, j ).resize( nOrder-i-j, eta3s.size() );
822  cs( i, j ) = dyna::JacobiBatchEvaluation( nOrder-i-j,
823  value_type( 2*( i+j+1 ) ), value_type( 0.0 ), eta3s );
824  dcs( i, j ) = dyna::JacobiBatchDerivation( nOrder-i-j,
825  value_type( 2*( i+j+1 ) ), value_type( 0.0 ), eta3s );
826  }
827  }
828 
829 
830  details::scalings<nOrder, T> scalings2( eta2s );
831  details::scalings<nOrder, T> scalings3( eta3s );
832 
833  //std::cout << "scalings = " << scalings() << "\n";
834  ublas::vector<value_type> one( ublas::scalar_vector<value_type>( eta1s.size(), 1.0 ) );
835  ublas::vector<value_type> tmp( ublas::scalar_vector<value_type>( eta1s.size(), 1.0 ) );
836 
837 
838  for ( uint16_type cur = 0, k = 0; k < nOrder+1; ++k )
839  {
840  for ( uint16_type i = 0; i < k+1; ++i )
841  {
842  for ( uint16_type j = 0; j < k+1-i; ++j,++cur )
843  {
844  uint16_type ii = k-i-j;
845  uint16_type jj = j;
846  uint16_type kk = i;
847 
848 
849  // x derivation
850  ublas::row( res[0], cur ) = ublas::element_prod( ublas::row( das, ii ),
851  ublas::row( bs[ii], jj ) );
852  ublas::row( res[0], cur ) = element_prod( ublas::row( res[0], cur ),
853  ublas::row( cs( ii, jj ), kk ) );
854 
855  if ( ii > 0 )
856  ublas::row( res[0], cur ) = element_prod( ublas::row( res[0], cur ),
857  ublas::row( scalings2(), ii-1 ) );
858 
859  if ( ii+jj > 0 )
860  ublas::row( res[0], cur ) = element_prod( ublas::row( res[0], cur ),
861  ublas::row( scalings3(), ii+jj-1 ) );
862 
863  // y derivation
864  ublas::row( res[1], cur ) = ublas::element_prod( ublas::row( das, ii ),
865  ublas::row( bs[ii], jj ) );
866  ublas::row( res[1], cur ) = element_prod( ublas::row( res[1], cur ),
867  ublas::row( cs( ii, jj ), kk ) );
868  ublas::row( res[1], cur ) = 0.5 * element_prod( ublas::row( res[1], cur ),
869  ( one+eta1s ) );
870 
871  if ( ii > 0 )
872  ublas::row( res[1], cur ) = element_prod( ublas::row( res[1], cur ),
873  ublas::row( scalings2(), ii-1 ) );
874 
875  if ( ii+jj > 0 )
876  ublas::row( res[1], cur ) = element_prod( ublas::row( res[1], cur ),
877  ublas::row( scalings3(), ii+jj-1 ) );
878 
879  // derivate (1-x)^ii
880  tmp = ublas::element_prod( ublas::row( scalings2(), ii ),
881  ublas::row( dbs[ii], jj ) );
882 
883  if ( ii > 0 )
884  tmp -= 0.5 * ii * ublas::element_prod( ublas::row( scalings2(), ii-1 ),
885  ublas::row( bs[ii], jj ) );
886 
887  tmp = ublas::element_prod( tmp,
888  ublas::row( as, ii ) );
889  tmp = ublas::element_prod( tmp,
890  ublas::row( cs( ii, jj ), kk ) );
891 
892  if ( ii+jj > 0 )
893  tmp = ublas::element_prod( tmp,
894  ublas::row( scalings3(), ii+jj-1 ) );
895 
896  // add contrib to y derivation
897  ublas::row( res[1], cur ) += tmp;
898 
899  // z derivation
900  ublas::row( res[2], cur ) = ublas::element_prod( ublas::row( das, ii ),
901  ublas::row( bs[ii], jj ) );
902  ublas::row( res[2], cur ) = element_prod( ublas::row( res[2], cur ),
903  ublas::row( cs( ii, jj ), kk ) );
904  ublas::row( res[2], cur ) = 0.5 * element_prod( ublas::row( res[2], cur ),
905  ( one+eta1s ) );
906 
907  if ( ii > 0 )
908  ublas::row( res[2], cur ) = element_prod( ublas::row( res[2], cur ),
909  ublas::row( scalings2(), ii-1 ) );
910 
911  if ( ii+jj > 0 )
912  ublas::row( res[2], cur ) = element_prod( ublas::row( res[2], cur ),
913  ublas::row( scalings3(), ii+jj-1 ) );
914 
915  // derivate (1-x)^ii
916  tmp = ublas::element_prod( ublas::row( scalings2(), ii ),
917  ublas::row( dbs[ii], jj ) );
918 
919  if ( ii > 0 )
920  tmp -= 0.5 * ii * ublas::element_prod( ublas::row( scalings2(), ii-1 ),
921  ublas::row( bs[ii], jj ) );
922 
923  tmp = ublas::element_prod( tmp,
924  ublas::row( as, ii ) );
925  tmp = ublas::element_prod( tmp,
926  ublas::row( cs( ii, jj ), kk ) );
927  tmp = 0.5 * element_prod( tmp, ( one+eta2s ) );
928 
929  if ( ii+jj > 0 )
930  tmp = ublas::element_prod( tmp,
931  ublas::row( scalings3(), ii+jj-1 ) );
932 
933  // add contrib to z derivation
934  ublas::row( res[2], cur ) += tmp;
935 
936  // derivate (1-x)^ii
937  tmp = ublas::element_prod( ublas::row( scalings3(), ii+jj ),
938  ublas::row( dcs( ii, jj ), kk ) );
939 
940  if ( ii+jj > 0 )
941  tmp -= 0.5*( ii+jj )*ublas::element_prod( ublas::row( cs( ii, jj ), kk ),
942  ublas::row( scalings3(), ii+jj-1 ) );
943 
944  tmp = ublas::element_prod( tmp,
945  ublas::row( as, ii ) );
946  tmp = ublas::element_prod( tmp,
947  ublas::row( bs[ ii ], jj ) );
948  tmp = ublas::element_prod( tmp,
949  ublas::row( scalings2(), ii ) );
950 
951  // add contrib to z derivation
952  ublas::row( res[2], cur ) += tmp;
953 
954  if ( is_normalized )
955  {
956  value_type normalization = math::sqrt( ( value_type( ii )+0.5 )*
957  ( value_type( ii+jj )+1.0 )*
958  ( value_type( ii+jj+kk )+1.5 ) );
959 
960  ublas::row( res[0], cur ) *= normalization;
961  ublas::row( res[1], cur ) *= normalization;
962  ublas::row( res[2], cur ) *= normalization;
963  }
964  }
965  }
966  }
967 
968  return res;
969 }
970 
971 
972 
973 }
974 #endif /* __Dubiner_H */

Generated on Sun Oct 20 2013 08:24:56 for Feel++ by doxygen 1.8.4