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
polynomialset.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  Gilles Steiner <gilles.steiner@epfl.ch>
7  Date: 2005-10-11
8 
9  Copyright (C) 2005,2006 EPFL
10  Copyright (C) 2006-2011 Université Joseph Fourier Grenoble 1
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 __PolynomialSet_H
32 #define __PolynomialSet_H 1
33 
34 #include <vector>
35 #include <boost/multi_array.hpp>
36 #include <boost/multi_array/extent_gen.hpp>
37 #include <boost/foreach.hpp>
38 #include <boost/optional.hpp>
39 #include <boost/mpl/min_max.hpp>
40 #include <Eigen/Core>
41 
42 #include <feel/feelpoly/context.hpp>
43 
44 #include <feel/feelalg/svd.hpp>
48 #include <feel/feelpoly/tensorisedboundadapted.hpp>
52 
53 
54 namespace Feel
55 {
56 
74 template<typename Poly, template<uint16_type> class PolySetType = Scalar >
75 class PolynomialSet
76 {
77 public:
78 
82 
83  static const uint16_type nDim = Poly::nDim;
84  static const uint16_type nRealDim = Poly::nRealDim;
85  static const uint16_type nOrder = Poly::nOrder;
86 
88 
92  typedef PolynomialSet<Poly, PolySetType> self_type;
93  typedef boost::shared_ptr<self_type> self_ptrtype;
94  typedef typename Poly::value_type value_type;
95  typedef typename Poly::basis_type basis_type;
96  static const bool is_product = Poly::is_product;
97 
98 
99  typedef PolySetType<nRealDim> polyset_type;
100  static const bool is_tensor2 = polyset_type::is_tensor2;
101  static const bool is_vectorial = polyset_type::is_vectorial;
102  static const bool is_scalar = polyset_type::is_scalar;
103  static const uint16_type nComponents = polyset_type::nComponents;
104  static const uint16_type nComponents1 = polyset_type::nComponents1;
105  static const uint16_type nComponents2 = polyset_type::nComponents2;
106  static const uint16_type rank = polyset_type::rank;
107 
108  typedef PolynomialSet<Poly, Scalar> component_type;
109  typedef Polynomial<Poly, PolySetType> polynomial_type;
110  //typedef Polynomial<Poly, PolySetType, ublas::vector_range<ublas::vector<value_type> > > polynomial_type;
111  typedef polynomial_type polynomial_view_type;
112 
113  typedef typename Poly::convex_type convex_type;
114  typedef typename basis_type::matrix_type matrix_type;
115  typedef typename basis_type::points_type points_type;
116 
117 
118  typedef PolynomialSet<Poly, Vectorial> gradient_polynomialset_type;
119 
120  BOOST_STATIC_ASSERT( ( boost::is_same<typename matrix_type::value_type, value_type>::value ) );
121  BOOST_STATIC_ASSERT( ( boost::is_same<typename matrix_type::value_type, typename points_type::value_type>::value ) );
122 
124 
128 
129  PolynomialSet()
130  :
131  M_basis(),
132  M_coeff(),
133  M_fname( "pset" )
134  {
135  //std::cout << "[PolynomialSet::default] dim = " << nDim << " order = " << nOrder << "\n";
136  }
137  PolynomialSet( Poly const& p )
138  :
139  M_basis( p.basis() ),
140  M_coeff( p.coeff() ),
141  M_fname( p.familyName() )
142  {}
145  //template<typename AE>
146  //PolynomialSet( Poly const& p, ublas::matrix_expression<AE> const& c )
147  PolynomialSet( Poly const& p, matrix_type const& c, bool __as_is = false )
148  :
149  M_basis( p.basis() ),
150  M_coeff( p.coeff() ),
151  M_fname( p.familyName() )
152  {
153  setCoefficient( c, __as_is );
154  //FEELPP_ASSERT( c.size2() == p.coeff().size1() )( c.size2() )( p.coeff().size1() ).error( "invalid dimension\n" );
155  //std::cout << "[PolynomialSet] dim = " << nDim << " order = " << nOrder << "\n";
156  //std::cout << "[PolynomialSet] c = " << c << "\n";
157  //std::cout << "[PolynomialSet] p.coeff = " << p.coeff() << "\n";
158  //std::cout << "[PolynomialSet] coeff = " << M_coeff << "\n";
159  }
160 
163  //template<typename AE>
164  //PolynomialSet( Poly const& p, ublas::matrix_expression<AE> const& c )
165  PolynomialSet( matrix_type const& c, bool __as_is = false )
166  :
167  M_basis(),
168  M_coeff( c ),
169  M_fname( "pset" )
170  {
171  setCoefficient( c, __as_is );
172  //FEELPP_ASSERT( c.size2() == p.coeff().size1() )( c.size2() )( p.coeff().size1() ).error( "invalid dimension\n" );
173  //std::cout << "[PolynomialSet] dim = " << nDim << " order = " << nOrder << "\n";
174  //std::cout << "[PolynomialSet] c = " << c << "\n";
175  //std::cout << "[PolynomialSet] p.coeff = " << p.coeff() << "\n";
176  //std::cout << "[PolynomialSet] coeff = " << M_coeff << "\n";
177  }
178 
179 
180  PolynomialSet( PolynomialSet const & p )
181  :
182  M_basis( p.M_basis ),
183  M_coeff( p.M_coeff ),
184  M_fname( p.M_fname )
185  {
186  //std::cout << "[PolynomialSet::copy] dim = " << nDim << " order = " << nOrder << "\n";
187  //std::cout << "[PolynomialSet::copy] p.coeff = " << p.coeff() << "\n";
188  //std::cout << "[PolynomialSet::copy] coeff = " << M_coeff << "\n";
189  }
190 
191  virtual ~PolynomialSet()
192  {}
193 
195 
199 
200  self_type& operator=( self_type const& pset )
201  {
202  if ( this != &pset )
203  {
204  M_basis = pset.M_basis;
205  M_coeff = pset.M_coeff;
206  }
207 
208  return *this;
209  }
210 
216  component_type operator[]( uint16_type i ) const
217  {
218  BOOST_STATIC_ASSERT( is_vectorial );
219  FEELPP_ASSERT( i < nComponents )( i )( nComponents ).error ( "invalid component index" );
220  const int nrows = M_coeff.size1()/nComponents;
221  const int ncols = M_coeff.size2();
222  return component_type( Poly(), ublas::project( M_coeff,
223  ublas::slice( nrows*i+i, nComponents, nrows/nComponents ),
224  ublas::slice( 0, 1, ncols ) ), true );
225  }
226 
228 
232 
236  uint16_type degree() const
237  {
238  return nOrder;
239  }
240 
245  matrix_type const& coeff() const
246  {
247  return M_coeff;
248  }
249 
253  basis_type const& basis() const
254  {
255  return M_basis;
256  }
257 
261  static bool isScalar()
262  {
263  return is_scalar;
264  }
265 
269  static bool isVectorial()
270  {
271  return is_vectorial;
272  }
273 
277  static uint16_type numberOfComponents()
278  {
279  return nComponents;
280  }
281 
286  {
287  return M_coeff.size1()/nComponents;
288  }
289 
294  {
295  return M_coeff.size2();
296  }
297 
301  bool isZero() const
302  {
303  return ublas::norm_frobenius( M_coeff ) < Feel::type_traits<value_type>::epsilon();
304  }
305 
310  virtual std::string familyName() const
311  {
312  return M_fname;
313  }
314 
326  std::string name( std::string sep = "." ) const
327  {
328  std::ostringstream os;
329  os << this->familyName() << sep << nDim << sep << nOrder;
330  return os.str();
331  }
332 
334 
338 
342  void setCoefficient( matrix_type const& __c, bool __as_is = false )
343  {
344  if ( is_scalar )
345  {
346  if ( !__as_is )
347  M_coeff = ublas::prod( __c, M_coeff );
348 
349  else
350  M_coeff = __c;
351  }
352 
353  else
354  {
355  if ( !__as_is )
356  {
357  M_coeff = ublas::prod( __c, polyset_type::toMatrix( M_coeff ) );
358  M_coeff = polyset_type::toType( M_coeff );
359  }
360 
361  else
362  M_coeff = __c;
363  }
364  }
365 
367 
371 
378  PolynomialSet<Poly, PolySetType> polynomials( std::vector<uint16_type> const& list_p ) const
379  {
380  size_type dim_p = this->polynomialDimension();
381  size_type new_dim_p = nComponents*list_p.size();
382  matrix_type coeff( nComponents*nComponents*list_p.size(), M_coeff.size2() );
383  int j = 0;
384  BOOST_FOREACH( uint16_type i, list_p )
385  {
386  for ( int c = 0; c < nComponents; ++c )
387  {
389  ublas::range( c*new_dim_p+nComponents*j, c*dim_p+nComponents*j+nComponents ),
390  ublas::range( 0, M_coeff.size2() ) ) =
391  ublas::project( M_coeff,
392  ublas::range( c*dim_p+nComponents*i, c*dim_p+nComponents*i+nComponents ),
393  ublas::range( 0, M_coeff.size2() ) );
394  }
395 
396  ++j;
397  }
398  return PolynomialSet<Poly, PolySetType>( Poly(), coeff, true );
399  }
400 
408  {
409  matrix_type coeff( nComponents*nComponents*dim_p, M_coeff.size2() );
410 
411  for ( int c = 0; c < nComponents; ++c )
412  {
413  size_type nc = c*this->polynomialDimension();
415  ublas::range( c*nComponents*dim_p, ( c+1 )*nComponents*dim_p ),
416  ublas::range( 0, M_coeff.size2() ) ) =
417  ublas::project( M_coeff,
418  ublas::range( nc, nc+nComponents*dim_p ),
419  ublas::range( 0, M_coeff.size2() ) );
420  }
421 
422  return PolynomialSet<Poly, PolySetType>( Poly(), coeff, true );
423  }
424 
433  polynomialsRange( uint16_type dim_bot, uint16_type dim_top ) const
434  {
435  uint16_type dim_p = dim_top-dim_bot;
436  matrix_type coeff( nComponents*nComponents*dim_p,
437  M_coeff.size2() );
438 
439  for ( int c = 0; c < nComponents; ++c )
440  {
441  size_type nc = c*this->polynomialDimension();
443  ublas::range( c*nComponents*dim_bot, ( c+1 )*nComponents*dim_top ),
444  ublas::range( 0, M_coeff.size2() ) ) =
445  ublas::project( M_coeff,
446  ublas::range( nc+dim_bot, nc+nComponents*dim_top ),
447  ublas::range( 0, M_coeff.size2() ) );
448  }
449 
450  return PolynomialSet<Poly, PolySetType>( Poly(), coeff, true );
451  }
452 
460  {
461  size_type dim_p = this->polynomialDimension();
462  matrix_type coeff( nComponents, M_coeff.size2() );
463  //for ( int c = 0; c < nComponents; ++c )
464  {
466  ublas::range( 0, nComponents ),
467  ublas::range( 0, M_coeff.size2() ) ) =
468  ublas::project( M_coeff,
469  //ublas::range( c*dim_p+nComponents*i, c*dim_p+nComponents*i+nComponents ),
470  ublas::range( nComponents*i, nComponents*( i+1 ) ),
471  ublas::range( 0, M_coeff.size2() ) );
472  }
473  return Polynomial<Poly, PolySetType> ( Poly(), coeff, true );
474  }
475 
476 
483  template<typename AE>
484  ublas::vector<value_type> evaluate( uint16_type i, ublas::vector_expression<AE> const& __pt ) const
485  {
486  return ublas::row( ublas::prod( M_coeff, M_basis( __pt ) ), i );
487  }
488 
493  template<typename AE>
494  ublas::matrix<value_type> evaluate( ublas::vector_expression<AE> const& __pt ) const
495  {
496  return ublas::prod( M_coeff, M_basis( __pt ) );
497  }
498 
499 
510  template<typename AE>
511  matrix_type evaluate( ublas::matrix_expression<AE> const& __pts ) const
512  {
513  matrix_type m ( M_basis.evaluate( __pts ) );
514  FEELPP_ASSERT( M_coeff.size2() == m.size1() )( M_coeff.size2() )( m.size1() ).error( "invalid size" );
515  return ublas::prod( M_coeff, m );
516  }
517 
539  matrix_type const& d( uint16_type i ) const
540  {
541  return M_basis.d( i );
542  }
543 
544  matrix_type d( uint16_type i, uint16_type j ) const
545  {
546  return ublas::prod( M_basis.d( i ), M_basis.d( j ) );
547  }
548 
554  self_type derivate( uint16_type l ) const
555  {
556  return self_type( Poly(), ublas::prod( M_coeff, M_basis.d( l ) ), true );
557  }
558 
559  template<typename AE>
560  ublas::vector<matrix_type> derivate( ublas::matrix_expression<AE> const& pts ) const
561  {
562  ublas::vector<matrix_type> der( M_basis.derivate( pts ) );
563  ublas::vector<matrix_type> res( nDim );
564 
565  for ( uint16_type i = 0; i < nDim; ++i )
566  {
567  res[i].resize( M_coeff.size1(), pts().size2() );
568  ublas::axpy_prod( M_coeff, der[i], res[i] );
569  }
570 
571  return res;
572  }
573 
574  template<typename AE>
575  matrix_type derivate( uint16_type i, ublas::matrix_expression<AE> const& pts ) const
576  {
577  ublas::vector<matrix_type> der( M_basis.derivate( pts ) );
578  matrix_type res( M_coeff.size1(), pts().size2() );
579  ublas::axpy_prod( M_coeff, der[i], res );
580  return res;
581  }
582 
583  template<typename AE>
584  matrix_type derivate( uint16_type i, uint16_type j, ublas::matrix_expression<AE> const& pts ) const
585  {
586  //std::cout << "[derivate2] M_coeff = " << M_coeff << "\n";
587  matrix_type eval( M_basis.evaluate( pts ) );
588  //matrix_type res( M_coeff.size1(), pts().size2() );
589  //ublas::axpy_prod( M_coeff, der[i], res );
590  matrix_type p1 = ublas::prod( M_coeff, M_basis.d( i ) );
591  matrix_type p2 = ublas::prod( p1, M_basis.d( j ) );
592  return ublas::prod( p2, eval );
593  }
599  gradient_polynomialset_type
600  gradient() const
601  {
602  return gradient( mpl::int_<polyset_type::rank>() );
603  }
604 
605  gradient_polynomialset_type
606  gradient( mpl::int_<0> ) const
607  {
608  const int n1 = M_coeff.size1();
609  const int n2 = M_coeff.size2();
610  ublas::matrix<value_type> c ( nDim*nDim*n1, n2 );
611  c.clear();
612 
613  for ( int i = 0; i <nDim; ++i )
614  {
615  ublas::project( c,
616  ublas::slice( nDim*n1*i+i, nDim, n1 ),
617  ublas::slice( 0, 1, n2 ) ) = ublas::prod( M_coeff, M_basis.d( i ) );
618  }
619 
620  return gradient_polynomialset_type( c, true );
621  }
622 
623  gradient_polynomialset_type
624  gradient( mpl::int_<1> ) const
625  {
626  // we deal with a VECTORIAL polynomial set: nComponents = nDim*nDim
627  // number of basis functions times the numb er of components(was vertorial function)
628  const int n1 = M_coeff.size1();
629  const int n2 = M_coeff.size2();
630  ublas::matrix<value_type> c ( nComponents*nComponents*n1, n2 );
631  c.clear();
632 #if 0
633  std::cout << "M_coeff = " << M_coeff << "\n"
634  << "c = " << c << "\n";
635 #endif
636 
637  for ( int i = 0; i <nDim; ++i )
638  {
639 #if 0
640  std::cout << "i=" << i << "\n"
641  << " prod = "
642  << ublas::prod( M_coeff, M_basis.d( i ) ) << "\n"
643  << " slice = "
644  << ublas::project( c,
645  ublas::slice( nDim*n1*i+i, nDim, n1 ),
646  ublas::slice( 0, 1, n2 ) ) << "\n";
647  // c[i][c1][c2][q]
648 #endif
649  ublas::project( c,
650  ublas::slice( nDim*n1*i+i, nDim, n1 ),
651  ublas::slice( 0, 1, n2 ) ) = ublas::prod( M_coeff, M_basis.d( i ) );
652  }
653 
654  //std::cout << "coeff = " << c << "\n";
655  return gradient_polynomialset_type( c, true );
656  }
657 
658 
659 
660 
664  uint16_type nbDof() const
665  {
666  return M_coeff.size1();
667  }
668 
669 
673  void insert( PolynomialSet<Poly,PolySetType> const& p, bool erase = false )
674  {
675  FEELPP_ASSERT( p.coeff().size2() == coeff().size2() )( p.coeff().size2() )( coeff().size2() ).warn( "invalid polynomial set" );
676 
677 #if 0
678  std::cout << "coeff = " << M_coeff << "\n"
679  << "p = " << p.coeff() << "\n";
680 #endif
681 
682  if ( erase )
683  {
684  M_coeff = p.M_coeff;
685  //std::cout << "insert after erase =" << M_coeff << "\n";
686  return;
687  }
688 
689  matrix_type oldcoeff = M_coeff;
690  M_coeff.resize( M_coeff.size1() + p.coeff().size1(), p.coeff().size2(), false );
691 
692 
693  if ( oldcoeff.size1() > 0 )
694  {
695  ublas::project( M_coeff,
696  ublas::range( 0, oldcoeff.size1() ),
697  ublas::range( 0, oldcoeff.size2() ) ) = oldcoeff;
698  }
699 
700  ublas::project( M_coeff,
701  ublas::range( oldcoeff.size1(), oldcoeff.size1()+p.coeff().size1() ),
702  ublas::range( 0, oldcoeff.size2() ) ) = p.coeff();
703  //std::cout << "insert =" << M_coeff << "\n";
704 
705  }
706 
707 
709 
710 
711  /******************************************************/
715  class PreCompute
716  {
717  public:
718  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
719  typedef PolynomialSet<Poly, PolySetType > reference_element_type;
720  typedef boost::shared_ptr<reference_element_type> reference_element_ptrtype;
721 
722  typedef typename reference_element_type::value_type value_type;
723 
724  static const uint16_type nDim = reference_element_type::nDim;
725  static const uint16_type nComponents1 = reference_element_type::nComponents1;
726  static const uint16_type nComponents2 = reference_element_type::nComponents2;
727  //static const uint16_type nComponents3 = reference_element_type::nComponents3;
728  static const uint16_type nComponents = reference_element_type::nComponents;
729 
730  typedef typename reference_element_type::points_type matrix_node_t_type;
731  typedef ublas::matrix<value_type> matrix_type;
732  typedef Eigen::Matrix<value_type,nComponents1,1> id_type;
733  typedef Eigen::Matrix<value_type,nComponents1,nDim> g_type;
734  typedef Eigen::Matrix<value_type,nDim,nDim> h_type;
735  typedef boost::multi_array<id_type,2> functionvalue_type;
736  typedef boost::multi_array<g_type,2> grad_type;
737  typedef boost::multi_array<h_type,2> hessian_type;
738  PreCompute() {}
739 
740 #if 0
741 
745  PreCompute( matrix_node_t_type const& __pts )
746  :
747  M_ref_ele( new reference_element_type() ),
748  M_nodes( __pts ),
749  M_phi(),
750  M_grad(),
751  M_hessian()
752  {
753  init( M_ref_ele, __pts, mpl::int_<rank>() );
754  }
755 #endif
756 
761  PreCompute( reference_element_ptrtype const& __ref_ele,
762  matrix_node_t_type const& __pts )
763  :
764  M_ref_ele( __ref_ele ),
765  M_nodes( __pts ),
766  M_phi(),
767  M_grad(),
768  M_hessian()
769  {
770  init( M_ref_ele, __pts, mpl::int_<rank>() );
771  }
772 
774  PreCompute( PreCompute const& __pc )
775  :
776  M_ref_ele( __pc.M_ref_ele ),
777  M_nodes( __pc.M_nodes ),
778  M_phi( __pc.M_phi ),
779  M_grad( __pc.M_grad ),
780  M_hessian( __pc.M_hessian )
781  {}
782 
784  ~PreCompute()
785  {}
786 
788  PreCompute& operator=( PreCompute const& __pc )
789  {
790  if ( this != &__pc )
791  {
792  M_ref_ele = __pc.M_ref_ele;
793  M_nodes = __pc.M_nodes;
794  M_phi = __pc.M_phi;
795  M_grad = __pc.M_grad;
796  M_hessian = __pc.M_hessian;
797  }
798 
799  return *this;
800  }
801 
802  void update( matrix_node_t_type const& __pts )
803  {
804  M_nodes = __pts;
805  init( M_ref_ele, __pts, mpl::int_<rank>() );
806  }
807 
811  uint16_type dim() const
812  {
813  return reference_element_type::nDim;
814  }
815 
820  uint16_type nComputedNodes() const
821  {
822  return M_nodes.size2();
823  }
824 
829  uint16_type nPoints() const
830  {
831  return M_nodes.size2();
832  }
833 
838  matrix_node_t_type const& nodes() const
839  {
840  return M_nodes;
841  }
842 
847  ublas::matrix_column<matrix_node_t_type const> node( uint16_type __i ) const
848  {
849  return ublas::column( M_nodes, __i );
850  }
851 
858  functionvalue_type const& phi() const
859  {
860  return M_phi;
861  }
862 
867  value_type phi( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q ) const
868  {
869  FEELPP_ASSERT( q < nComputedNodes() )( q )( nComputedNodes() ).error( "invalid node index" );
870  FEELPP_ASSERT( i < M_ref_ele->nbDof() )( i )( M_ref_ele->nbDof() ).error( "invalid dof index" );
871 
872  return M_phi[i][q]( c1,c2 );
873  }
874 
875  grad_type const& grad() const
876  {
877  return M_grad;
878  }
879 
880  value_type grad( size_type i, uint16_type c1, uint16_type c2, uint16_type q ) const
881  {
882  return M_grad[i][q]( c1,c2 );
883  }
884 
885  hessian_type const& hessian() const
886  {
887  return M_hessian;
888  }
889 
890  value_type hessian( size_type i, uint16_type c1, uint16_type c2, uint16_type q ) const
891  {
892  return M_hessian[i][q]( c1,c2 );
893  }
894 
895  private:
896 
897  void
898  init( reference_element_ptrtype const& M_ref_ele,
899  matrix_node_t_type const& __pts,
900  mpl::int_<0> )
901  {
902  M_phi.resize( boost::extents[M_ref_ele->nbDof()][__pts.size2()] );
903  M_grad.resize( boost::extents[M_ref_ele->nbDof()][__pts.size2()] );
904  M_hessian.resize( boost::extents[M_ref_ele->nbDof()][__pts.size2()] );
905 
906  matrix_type phiv = M_ref_ele->evaluate( __pts );
907  ublas::vector<matrix_type> __grad( M_ref_ele->derivate( __pts ) );
908  matrix_type __hessian( M_ref_ele->gradient().gradient().evaluate( __pts ) );
909 
910  typedef typename grad_type::index index;
911  const index I = M_ref_ele->nbDof();
912  const index Q = __pts.size2();
913 
914  for ( index i = 0; i < I; ++i )
915  {
916  for ( index q = 0; q < Q; ++q )
917  M_phi[i][q]( 0,0 ) = phiv( i, q );
918 
919  for ( index q = 0; q < Q; ++q )
920  for ( index j = 0; j < nDim; ++j )
921  M_grad[i][q]( 0,j ) = __grad[j]( i, q );
922 
923  for ( index q = 0; q < Q; ++q )
924  for ( index j = 0; j < nDim; ++j )
925  for ( index k = j; k < nDim; ++k )
926  {
927  value_type t = __hessian( nDim*nDim*I*( nDim*k+j )+nDim*nDim*i+nDim*j+k, q );
928  M_hessian[i][q]( j,k ) = t;
929  M_hessian[i][q]( k,j ) = t;
930  }
931  }
932 
933  }
934  void
935  init( reference_element_ptrtype const& M_ref_ele,
936  matrix_node_t_type const& __pts,
937  mpl::int_<1> )
938  {
939  typedef typename grad_type::index index;
940 #if 0
941  std::cout << "family = " << M_ref_ele->familyName() << std::endl;
942  std::cout << "is_product = " << reference_element_type::is_product << std::endl;
943  std::cout << "nbDof = " << M_ref_ele->nbDof() << std::endl;
944 #endif
945  const index I = ( reference_element_type::is_product?M_ref_ele->nbDof()/nRealDim/nRealDim:M_ref_ele->nbDof()/nRealDim );
946  //std::cout << "I = " << I << std::endl;
947  //std::cout << "nbDof = " << M_ref_ele->nbDof() << std::endl;
948 
949  const index Q = __pts.size2();
950  //std::cout << "Q = " << I << std::endl;
951 
952  const int ncdof= ( reference_element_type::is_product?nRealDim:1 );
953  int nldof= ( reference_element_type::is_product?I*nRealDim:I );
954  //std::cout << "ncdof = " << ncdof << ", nldof = " << nldof << "\n";
955  M_phi.resize( boost::extents[nldof][__pts.size2()] );
956  M_grad.resize( boost::extents[nldof][__pts.size2()] );
957 
958  matrix_type phiv = M_ref_ele->evaluate( __pts );
959  ublas::vector<matrix_type> __grad( M_ref_ele->derivate( __pts ) );
960 
961  for ( index i = 0; i < I; ++i )
962  {
963 
964  for ( index c1 = 0; c1 < ncdof; ++c1 )
965  {
966  for ( index q = 0; q < Q; ++q )
967  for ( index j = 0; j < nRealDim; ++j )
968  M_phi[I*c1+i][q]( j,0 ) = phiv( nldof*c1+nRealDim*i+j, q );
969 
970  //M_phi[I*c1+i][q](j,0) = phiv( nDim*I*c1+nDim*i+j, q );
971 
972  for ( index q = 0; q < Q; ++q )
973  for ( index j = 0; j < nRealDim; ++j )
974  for ( index l = 0; l < nDim; ++l )
975  {
976  //M_grad[I*c1+i][j](l,q) = __grad[l]( nDim*I*c1+nDim*i+j, q );
977  M_grad[I*c1+i][q]( j,l ) = __grad[l]( nldof*c1+nRealDim*i+j, q );
978  //M_grad[I*c1+i][j][nRealDim-1][q] = __grad[l]( nldof*c1+nRealDim*i+j, q );
979  //std::cout << "grad(" << i << "," << c1 << "," << j << "," << l << "," << q << ")=" << M_grad[I*c1+i][j](l,q) << "\n";
980  }
981  }
982  }
983  }
984  private:
985 
986  reference_element_ptrtype M_ref_ele;
987  matrix_node_t_type M_nodes;
988  functionvalue_type M_phi;
989  grad_type M_grad;
990  hessian_type M_hessian;
991  };
993  typedef PreCompute precompute_type;
994  typedef boost::shared_ptr<precompute_type> precompute_ptrtype;
995 
996  precompute_ptrtype
997  preCompute( self_ptrtype p, points_type const& P )
998  {
999  return precompute_ptrtype( new PreCompute( p, P ) );
1000  }
1001 
1002  typedef std::vector<std::map<typename convex_type::permutation_type, precompute_ptrtype> > faces_precompute_type;
1003 
1004  std::vector<std::map<typename convex_type::permutation_type, precompute_ptrtype> >
1005  preComputeOnFaces( self_ptrtype p, points_type const& P )
1006  {
1007 #if 0
1008  QuadMapped<pointset_type> qm;
1009  typedef typename QuadMapped<pointset_type>::permutation_type permutation_type;
1010  typename QuadMapped<pointset_type>::permutation_points_type ppts( qm( P ) );
1011 #endif
1012  typedef typename convex_type::permutation_type permutation_type;
1013  std::vector<std::map<permutation_type, precompute_ptrtype> > geopc( convex_type::numTopologicalFaces );
1014 
1015  for ( uint16_type __f = 0; __f < convex_type::numTopologicalFaces; ++__f )
1016  {
1017  for ( permutation_type __p( permutation_type::IDENTITY );
1018  __p < permutation_type( permutation_type::N_PERMUTATIONS ); ++__p )
1019  {
1020  //FEELPP_ASSERT( ppts[__f].find(__p)->second.size2() != 0 ).warn( "invalid quadrature type" );
1021  geopc[__f][__p] = precompute_ptrtype( new precompute_type( p, P ) );
1022  }
1023  }
1024 
1025  return geopc;
1026  }
1027 
1028 
1029  template<size_type context_v, typename Basis_t, typename Geo_t, typename ElementType, size_type context_g = context_v>
1030  class Context
1031  {
1032  public:
1033  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
1034  static const int context = context_v;
1035  // reference space dimension
1036  static const uint16_type PDim = ElementType::nDim;
1037  // real space dimension
1038  static const uint16_type NDim = ElementType::nRealDim;
1039  static const uint16_type nDof = Basis_t::nLocalDof;
1040  static const bool is_product = Basis_t::is_product;
1041  static const bool is_scalar = Basis_t::is_scalar;
1042  static const bool is_vectorial = Basis_t::is_vectorial;
1043  static const bool is_tensor2 = Basis_t::is_tensor2;
1044  static const uint16_type nComponents = Basis_t::nComponents;
1045  static const uint16_type nComponents1 = Basis_t::nComponents1;
1046  static const uint16_type nComponents2 = Basis_t::nComponents2;
1047 
1048  typedef typename Basis_t::polyset_type polyset_type;
1049  static const uint16_type rank = polyset_type::rank;
1050 
1051  typedef Basis_t reference_element_type;
1052  typedef boost::shared_ptr<Basis_t> reference_element_ptrtype;
1053 
1054  typedef typename reference_element_type::value_type value_type;
1055 
1056  typedef ElementType geometric_element_type;
1057  typedef typename Geo_t::template Context<context_g, ElementType> geometric_mapping_context_type;
1058  typedef boost::shared_ptr<geometric_mapping_context_type> geometric_mapping_context_ptrtype;
1059 
1060  typedef typename node<value_type>::type node_type;
1061 
1062  typedef ublas::matrix<value_type> matrix_type;
1063  typedef points_type matrix_node_t_type;
1064 #if 0
1065  typedef value_type dn_type;
1066  typedef value_type grad_type;
1067  typedef value_type div_type;
1068 
1069 
1070  typedef typename matrix_type::value_type phi_type;
1071  typedef typename matrix_type::value_type dphi_type;
1072  typedef typename matrix_type::value_type id_type;
1073 #else
1074  typedef Eigen::Matrix<value_type,nComponents1,1> id_type;
1075  typedef Eigen::Matrix<value_type,nComponents1,nDim> ref_grad_type;
1076  typedef Eigen::Matrix<value_type,nComponents1,NDim> grad_type;
1077  typedef Eigen::Matrix<value_type,NDim,NDim> hess_type;
1078  typedef Eigen::Matrix<value_type,1,1> div_type;
1079  typedef Eigen::Matrix<value_type,nComponents1,1> dn_type;
1080  typedef Eigen::Matrix<value_type,3,1> curl_type;
1081  typedef Eigen::Matrix<value_type,nComponents1,1> dx_type;
1082  typedef Eigen::Matrix<value_type,nComponents1,1> dy_type;
1083  typedef Eigen::Matrix<value_type,nComponents1,1> dz_type;
1084 
1085 #if 0
1086  typedef boost::multi_array<id_type,2> functionvalue_type;
1087  typedef boost::multi_array<g_type,2> grad_type;
1088  typedef boost::multi_array<h_type,2> hessian_type;
1089  typedef boost::multi_array<n_type,2> dn_type;
1090  typedef boost::multi_array<d_type,2> div_type;
1091 #endif
1092 #endif
1093 
1094  typedef geometric_mapping_context_type gmc_type;
1095  typedef Eigen::Matrix<value_type,Eigen::Dynamic, Eigen::Dynamic> matrix_eigen_type;
1096  typedef Eigen::Matrix<value_type,gmc_type::NDim,gmc_type::PDim> matrix_eigen_NP_type;
1097  typedef Eigen::Matrix<value_type,gmc_type::PDim,gmc_type::NDim> matrix_eigen_PN_type;
1098  typedef Eigen::Matrix<value_type,gmc_type::NDim,gmc_type::NDim> matrix_eigen_NN_type;
1099  typedef Eigen::Matrix<value_type,nComponents1,NDim> matrix_eigen_grad_type;
1100  typedef typename Eigen::Map<const Eigen::Matrix<value_type,gmc_type::NDim,gmc_type::PDim,Eigen::RowMajor> > matrix_eigen_ublas_NP_type;
1101  typedef typename Eigen::Map<const Eigen::Matrix<value_type,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> > matrix_eigen_ublas_type;
1102 
1103 
1104  template<uint16_type TheRank = polyset_type::rank+2>
1105  struct Index
1106  {
1107  static const uint16_type rank = TheRank;
1108  typedef boost::array<size_type,rank> index_type;
1109  typedef boost::detail::multi_array::extent_gen<rank> extents_type;
1110  template<uint16_type> friend class Rank;
1111  Index()
1112  {
1113  //std::cout << "Index::rank = " << rank << "\n";
1114  init( mpl::int_<rank>() );
1115  }
1116  Index( Index const& i )
1117  :
1118  M_index( i.M_index ),
1119  M_extents( i.M_extents ),
1120  M_comp( 0 )
1121  {
1122  //std::cout << "Index::rank = " << rank << "\n";
1123  }
1124 
1125  // rank + 1
1126  Index( Index<rank-1> const& __index )
1127  :
1128  M_extents( __index.extents()[RankUp<polyset_type>::type::nComponentsLast] )
1129  {
1130  std::copy( __index.beginIndex(), __index.endIndex(), M_index.begin() );
1131  //std::cout << "[rankup] index<" << rank << ">: ";
1132  //std::for_each( M_index.begin(), M_index.end(), std::cout << lambda::_1 << " " );
1133  //std::cout << "\n";
1134  }
1135 #if 0
1136  // rank reduction wrt the last extent
1137  Index( Index<rank+1> const& index )
1138  :
1139  M_extents()
1140  {
1141 #if 0
1142  // copy only up to rank
1143  std::copy( index.M_extents.begin(),
1144  boost::prior( index.M_extents.end() ),
1145  M_extents.begin() );
1146 #endif
1147  }
1148 #endif
1149  ~Index() {}
1150 
1151  Index const& operator=( Index const& i )
1152  {
1153  if ( this != &i )
1154  {
1155  M_index = i.M_index;
1156  M_extents = i.M_extents;
1157  }
1158 
1159  return *this;
1160  }
1161 
1162  typename index_type::iterator beginIndex()
1163  {
1164  return M_index.begin();
1165  }
1166  typename index_type::const_iterator beginIndex() const
1167  {
1168  return M_index.begin();
1169  }
1170  typename index_type::iterator endIndex()
1171  {
1172  return M_index.end();
1173  }
1174  typename index_type::const_iterator endIndex() const
1175  {
1176  return M_index.end();
1177  }
1178 
1179  template<typename Tuple>
1180  void setIndex( Tuple const& tu )
1181  {
1182  setIndex( tu, mpl::int_<rank>() );
1183  }
1184 
1185  void setIndex( uint16_type c, size_type i )
1186  {
1187  M_index[c] = i;
1188  }
1189 
1190 
1191  size_type index() const
1192  {
1193  return index( mpl::int_<rank>() );
1194  }
1195 
1196  size_type div() const
1197  {
1198  return nComponents*nDof*M_index[1] + nComponents*M_index[0] + M_index[1];
1199  }
1200 
1201  uint16_type component() const
1202  {
1203  return M_comp;
1204  }
1205 
1206  Index<rank+1> rankUp() const
1207  {
1208  return Index<rank+1>( M_extents[RankUp<polyset_type>::type::nComponentsLast] );
1209  }
1210 
1211  operator size_type() const
1212  {
1213  return index( mpl::int_<rank>() );
1214  }
1215 
1216  extents_type extents() const
1217  {
1218  return M_extents;
1219  }
1220 
1221  private:
1222  void init( mpl::int_<2> )
1223  {
1224  M_extents = boost::extents[nDof][nComponents];
1225  M_index[0] = size_type( -1 );
1226  }
1227  void init( mpl::int_<3> )
1228  {
1229  M_extents = boost::extents[nDof][nComponents][nComponents];
1230  M_index[0] = size_type( -1 );
1231  M_index[1] = size_type( -1 );
1232  M_index[2] = size_type( -1 );
1233  }
1234  void init( mpl::int_<4> )
1235  {
1236  M_extents = boost::extents[nDof][nComponents][nComponents][nComponents];
1237  M_index[0] = size_type( -1 );
1238  M_index[1] = size_type( -1 );
1239  M_index[2] = size_type( -1 );
1240  }
1241  size_type index( mpl::int_<2> ) const
1242  {
1243  return M_index[0];
1244  }
1245  size_type index( mpl::int_<3> ) const
1246  {
1247  return nDof*M_index[1] + M_index[0];// + M_index[2];
1248  }
1249  size_type index( mpl::int_<4> ) const
1250  {
1251  const size_type d = M_extents.ranges_[1].size();
1252  const size_type d2 = d*d;
1253  const size_type n = M_extents.ranges_[0].size();
1254 #if 0
1255  return ( d2*n*( d*M_index[1]+M_index[2] ) +
1256  d2*M_index[0] +
1257  d*M_index[1] + M_index[2] );
1258 #else
1259  return ( d2*n*M_index[1]*M_index[2] +
1260  d2*M_index[0]*M_index[2] +
1261  d*M_index[1]*M_index[2] +
1262  d2*n*M_index[1]+
1263  d2*M_index[0] +
1264  d*M_index[1] +
1265  M_index[2] );
1266 #endif
1267  }
1268  template<typename Tuple>
1269  void setIndex( Tuple const& tu, mpl::int_<2> )
1270  {
1271  M_index[ 0 ] = boost::get<0>( tu );
1272  M_comp = 0;
1273  }
1274  template<typename Tuple>
1275  void setIndex( Tuple const& tu, mpl::int_<3> )
1276  {
1277  M_index[ 0 ] = boost::get<0>( tu );
1278  M_index[ 1 ] = boost::get<1>( tu );
1279  M_index[ 2 ] = boost::get<2>( tu );
1280  M_comp = M_index[ 1 ];
1281  }
1282  template<typename Tuple>
1283  void setIndex( Tuple const& tu, mpl::int_<4> )
1284  {
1285  M_index[ 0 ] = boost::get<0>( tu );
1286  M_index[ 1 ] = boost::get<1>( tu );
1287  M_index[ 2 ] = boost::get<2>( tu );
1288  }
1289 
1290  private:
1291  index_type M_index;
1292 
1293  extents_type M_extents;
1294 
1295  uint16_type M_comp;
1296  };
1297 
1298  Context( reference_element_ptrtype const& __RefEle,
1299  geometric_mapping_context_ptrtype const& __gmc,
1300  precompute_ptrtype const& __pc )
1301  :
1302  M_pc( __pc ),
1303  M_npoints( __pc->nPoints() ),
1304 
1305  M_ipt( 0 ),
1306  M_ref_ele( __RefEle ),
1307 
1308  M_gmc( __gmc ),
1309  M_phi( __pc->phi() ),
1310  M_gradphi( __pc->grad() ),
1311  M_dn(),
1312  M_grad(),
1313  M_dx(),
1314  M_dy(),
1315  M_dz()
1316  {
1317  if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1318  {
1319  const int ntdof = nDof*nComponents1;
1320  M_grad.resize( boost::extents[ntdof][M_npoints] );
1321  M_dx.resize( boost::extents[ntdof][M_npoints] );
1322  M_dy.resize( boost::extents[ntdof][M_npoints] );
1323  M_dz.resize( boost::extents[ntdof][M_npoints] );
1324 
1325  if ( vm::has_first_derivative_normal<context>::value )
1326  {
1327  M_dn.resize( boost::extents[ntdof][M_npoints] );
1328  }
1329 
1330  if ( vm::has_div<context>::value )
1331  {
1332  M_div.resize( boost::extents[ntdof][M_npoints] );
1333  }
1334 
1335  if ( vm::has_curl<context>::value )
1336  {
1337  M_curl.resize( boost::extents[ntdof][M_npoints] );
1338  }
1339 
1340  if ( vm::has_hessian<context>::value || vm::has_second_derivative<context>::value )
1341  {
1342  M_hessian.resize( boost::extents[ntdof][M_npoints] );
1343  }
1344  }
1345 
1346  update( __gmc );
1347  }
1348 
1353  void transformationEquivalence( geometric_mapping_context_ptrtype const& __gmc,
1354  mpl::bool_<true> )
1355  {
1356  // M_phi = phi;
1357  // if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1358  // M_gradphi = M_pc->grad();
1359  // if ( vm::has_hessian<context>::value || vm::has_second_derivative<context>::value )
1360  // M_hessphi = M_pc->hessian();
1361  }
1362 
1369  void transformationEquivalence( geometric_mapping_context_ptrtype const& __gmc,
1370  mpl::bool_<false> )
1371  {
1372  //M_ref_ele->transform( __gmc, phi, M_phi, M_gradphi, M_hessphi );
1373  M_ref_ele->transform( __gmc, M_pc.get(), M_phi,
1374  M_gradphi, ( vm::has_div<context>::value || vm::has_curl<context>::value || vm::has_grad<context>::value || vm::has_first_derivative<context>::value ),
1375  M_hessphi, ( vm::has_hessian<context>::value || vm::has_second_derivative<context>::value )
1376  );
1377  }
1378 
1379  void update( geometric_mapping_context_ptrtype const& __gmc,
1380  precompute_ptrtype const& __pc )
1381  {
1382  M_pc = __pc;
1383  M_gmc = __gmc ;
1384 
1385  if ( ( M_npoints != M_gmc->nPoints() ) ||
1386  ( M_npoints != __pc->nPoints() ) ||
1387  ( M_npoints != M_phi.shape()[1] ) )
1388  {
1389  //std::cout << "M_npoints = " << M_npoints << "\n";
1390  //std::cout << "pc->npoints = " << __pc->nPoints() << "\n";
1391  //std::cout << "phi->npoints = " << M_phi.shape()[1] << "\n";
1392  //std::cout << "gmc->npoints = " << M_gmc->nPoints() << "\n";
1393  M_npoints = __pc->nPoints();
1394 
1395  const int ntdof = nDof*nComponents1;
1396  M_phi.resize( boost::extents[ntdof][M_npoints] );
1397  M_gradphi.resize( boost::extents[ntdof][M_npoints] );
1398 
1399  if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1400  {
1401  M_grad.resize( boost::extents[ntdof][M_npoints] );
1402  M_dx.resize( boost::extents[ntdof][M_npoints] );
1403  M_dy.resize( boost::extents[ntdof][M_npoints] );
1404  M_dz.resize( boost::extents[ntdof][M_npoints] );
1405 
1406  if ( vm::has_div<context>::value )
1407  {
1408  M_div.resize( boost::extents[ntdof][M_npoints] );
1409  }
1410 
1411  if ( vm::has_curl<context>::value )
1412  {
1413  M_curl.resize( boost::extents[ntdof][M_npoints] );
1414  }
1415 
1416  if ( vm::has_first_derivative_normal<context>::value )
1417  {
1418  M_dn.resize( boost::extents[ntdof][M_npoints] );
1419  }
1420 
1421  if ( vm::has_hessian<context>::value || vm::has_second_derivative<context>::value )
1422  {
1423  M_hessian.resize( boost::extents[ntdof][M_npoints] );
1424  }
1425  }
1426  }
1427 
1428  M_phi = M_pc.get()->phi();
1429  M_gradphi = M_pc.get()->grad();
1430 
1431  update( __gmc );
1432  }
1433  void update( geometric_mapping_context_ptrtype const& __gmc )
1434  {
1435  //M_phi = M_pc->get()->phi();
1436  //M_gradphi = M_pc->get()->grad();
1437  static const bool do_opt= ( nOrder<=1 ) && ( Geo_t::nOrder==1 ) && ( convex_type::is_simplex );
1438  transformationEquivalence( __gmc, mpl::bool_<Basis_t::isTransformationEquivalent>() );
1439 #if 0
1440 
1441  for ( int i = 0; i < M_gradphi.num_elements(); ++i )
1442  std::cout << "M_gradphi[" << i << "]=" << *( M_gradphi.data()+i ) << "\n";
1443 
1444 #endif
1445  update( __gmc, mpl::int_<rank>(), mpl::bool_<do_opt>() );
1446  //update( __gmc, mpl::int_<rank>(), mpl::bool_<false>() );
1447  }
1448  void update( geometric_mapping_context_ptrtype const& __gmc, mpl::int_<0>, mpl::bool_<true> )
1449  {
1450 
1451 
1452  //#pragma omp parallel
1453  //precompute_type* __pc = M_pc.get().get();
1454  geometric_mapping_context_type* thegmc = __gmc.get();
1455 
1456  if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1457  {
1458  const uint16_type Q = M_npoints;//__gmc->nPoints();//M_grad.size2();
1459  const uint16_type I = nDof; //M_ref_ele->nbDof();
1460 
1461 
1462  matrix_eigen_ublas_type Bt ( thegmc->B( 0 ).data().begin(), gmc_type::NDim, gmc_type::PDim );
1463  matrix_eigen_grad_type grad_real = matrix_eigen_grad_type::Zero();
1464 
1465  for ( uint16_type i = 0; i < I; ++i )
1466  {
1467  grad_real.noalias() = M_gradphi[i][0]*Bt.transpose();
1468  for ( uint16_type q = 0; q < Q; ++q )
1469  {
1470  M_grad[i][q] = grad_real;
1471  M_dx[i][q] = M_grad[i][q].col( 0 );
1472 
1473  if ( NDim == 2 )
1474  M_dy[i][q] = M_grad[i][q].col( 1 );
1475 
1476  if ( NDim == 3 )
1477  M_dz[i][q] = M_grad[i][q].col( 2 );
1478 
1479 
1480  }
1481  }
1482 
1483  // we need the normal derivative
1484  if ( vm::has_first_derivative_normal<context>::value )
1485  {
1486  std::fill( M_dn.data(), M_dn.data()+M_dn.num_elements(), dn_type::Zero() );
1487  const uint16_type I = M_ref_ele->nbDof()*nComponents;
1488  const uint16_type Q = nPoints();
1489 
1490  for ( int i = 0; i < I; ++i )
1491  {
1492  for ( uint16_type q = 0; q < Q; ++q )
1493  {
1494  for ( uint16_type l = 0; l < NDim; ++l )
1495  {
1496  value_type n = thegmc->unitNormal( l, 0 );
1497  value_type gn = M_grad[i][0]( 0,l ) * n;
1498  M_dn[i][q]( 0,0 ) += gn;
1499  }
1500  }
1501 
1502  }
1503  }
1504 
1505  } // grad
1506 
1507  if ( vm::has_hessian<context>::value || vm::has_second_derivative<context>::value )
1508  {
1509  precompute_type* __pc = M_pc.get().get();
1510  const uint16_type Q = M_npoints;//__gmc->nPoints();//M_grad.size2();
1511  const uint16_type I = nDof; //M_ref_ele->nbDof();
1512 
1513  // hessian only for P1 geometric mappings
1514  boost::multi_array<value_type,4> const& B3 = thegmc->B3();
1515 
1516  std::fill( M_hessian.data(), M_hessian.data()+M_hessian.num_elements(), hess_type::Zero() );
1517 
1518  //#pragma omp for
1519  for ( uint16_type i = 0; i < I; ++i )
1520  {
1521  for ( uint16_type q = 0; q < Q; ++q )
1522  {
1523 
1524  for ( uint16_type l = 0; l < NDim; ++l )
1525  {
1526  for ( uint16_type j = 0; j < NDim; ++j )
1527  {
1528  for ( uint16_type p = 0; p < PDim; ++p )
1529  {
1530  for ( uint16_type r = 0; r < PDim; ++r )
1531  {
1532  // we have twice the same contibution thanks to the symmetry
1533  value_type h = B3[l][j][p][r] * __pc->hessian( i, p, r, q );
1534  M_hessian[i][q]( j,l ) += h;
1535  } // q
1536  } // r
1537  } // p
1538  } // j
1539  } // l
1540  } // i
1541  }
1542  }
1543 
1544  void update( geometric_mapping_context_ptrtype const& __gmc, mpl::int_<0>, mpl::bool_<false> )
1545  {
1546  //#pragma omp parallel
1547  //precompute_type* __pc = M_pc.get().get();
1548  geometric_mapping_context_type* thegmc = __gmc.get();
1549 
1550  if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1551  {
1552  const uint16_type Q = M_npoints;//__gmc->nPoints();//M_grad.size2();
1553  const uint16_type I = nDof; //M_ref_ele->nbDof();
1554 
1555  for ( uint16_type i = 0; i < I; ++i )
1556  {
1557  for ( uint16_type q = 0; q < Q; ++q )
1558  {
1559  matrix_eigen_ublas_type Bt ( thegmc->B( 0 ).data().begin(), gmc_type::NDim, gmc_type::PDim );
1560  M_grad[i][q] = M_gradphi[i][q]*Bt.transpose();
1561  M_dx[i][q] = M_grad[i][q].col( 0 );
1562 
1563  if ( NDim == 2 )
1564  M_dy[i][q] = M_grad[i][q].col( 1 );
1565 
1566  if ( NDim == 3 )
1567  M_dz[i][q] = M_grad[i][q].col( 2 );
1568  }
1569  }
1570 
1571  // we need the normal derivative
1572  if ( vm::has_first_derivative_normal<context>::value )
1573  {
1574  std::fill( M_dn.data(), M_dn.data()+M_dn.num_elements(), dn_type::Zero() );
1575  const uint16_type I = M_ref_ele->nbDof()*nComponents;
1576  const uint16_type Q = nPoints();
1577 
1578  for ( int i = 0; i < I; ++i )
1579  {
1580  for ( uint16_type q = 0; q < Q; ++q )
1581  {
1582  for ( uint16_type l = 0; l < NDim; ++l )
1583  {
1584  M_dn[i][q]( 0,0 ) += M_grad[i][q]( 0,l ) * thegmc->unitNormal( l, q );
1585  }
1586  }
1587  }
1588  }
1589 
1590  } // grad
1591 
1592  if ( vm::has_hessian<context>::value || vm::has_second_derivative<context>::value )
1593  {
1594  precompute_type* __pc = M_pc.get().get();
1595  const uint16_type Q = M_npoints;//__gmc->nPoints();//M_grad.size2();
1596  const uint16_type I = nDof; //M_ref_ele->nbDof();
1597 
1598  // hessian only for P1 geometric mappings
1599  boost::multi_array<value_type,4> const& B3 = thegmc->B3();
1600 
1601  std::fill( M_hessian.data(), M_hessian.data()+M_hessian.num_elements(), hess_type::Zero() );
1602 
1603  //#pragma omp for
1604  for ( uint16_type i = 0; i < I; ++i )
1605  {
1606  for ( uint16_type q = 0; q < Q; ++q )
1607  {
1608  for ( uint16_type l = 0; l < NDim; ++l )
1609  {
1610  //for ( uint16_type j = l; j < NDim; ++j )
1611  for ( uint16_type j = 0; j < NDim; ++j )
1612  {
1613  for ( uint16_type p = 0; p < PDim; ++p )
1614  {
1615  // deal with _extra_ diagonal contributions
1616  //for ( uint16_type r = p+1; r < PDim; ++r )
1617  for ( uint16_type r = 0; r < PDim; ++r )
1618  {
1619  // we have twice the same contibution thanks to the symmetry
1620  value_type h = B3[l][j][p][r] * __pc->hessian( i, p, r, q );
1621  M_hessian[i][q]( l,j ) += h;
1622  } // q
1623  } // r
1624  } // p
1625  } // j
1626  } // l
1627  } // i
1628  }
1629  }
1630 
1631  void update( geometric_mapping_context_ptrtype const& __gmc, mpl::int_<1>, mpl::bool_<false> )
1632  {
1633  //precompute_type* __pc = M_pc.get().get();
1634  geometric_mapping_context_type* thegmc = __gmc.get();
1635 
1636  if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1637  {
1638  const uint16_type Q = M_npoints;//__gmc->nPoints();//M_grad.size2();
1639  const uint16_type I = nDof; //M_ref_ele->nbDof();
1640 
1641  typedef typename boost::multi_array<value_type,4>::index_range range;
1642 
1643  for ( uint16_type ii = 0; ii < I; ++ii )
1644  {
1645  int ncomp= ( reference_element_type::is_product?nComponents1:1 );
1646 
1647  for ( uint16_type c = 0; c < ncomp; ++c )
1648  {
1649  uint16_type i = I*c + ii;
1650 
1651  for ( uint16_type q = 0; q < Q; ++q )
1652  {
1653  //uint16_type c1 = c;
1654  matrix_eigen_ublas_type Bt ( thegmc->B( q ).data().begin(), gmc_type::NDim, gmc_type::PDim );
1655  //matrix_type const& Bq = thegmc->B( q );
1656  M_grad[i][q] = M_gradphi[i][q]*Bt.transpose();
1657  M_dx[i][q] = M_grad[i][q].col( 0 );
1658 
1659  if ( NDim == 2 )
1660  M_dy[i][q] = M_grad[i][q].col( 1 );
1661 
1662  if ( NDim == 3 )
1663  M_dz[i][q] = M_grad[i][q].col( 2 );
1664 
1665  // update divergence if needed
1666  if ( vm::has_div<context>::value )
1667  {
1668  M_div[i][q]( 0,0 ) = M_grad[i][q].trace();
1669  }
1670 
1671  // update curl if needed
1672  if ( vm::has_curl<context>::value )
1673  {
1674  if ( NDim == 2 )
1675  {
1676  M_curl[i][q]( 0 ) = M_grad[i][q]( 1,0 ) - M_grad[i][q]( 0,1 );
1677  M_curl[i][q]( 1 ) = M_curl[i][q]( 0 );
1678  M_curl[i][q]( 2 ) = M_curl[i][q]( 0 );
1679  }
1680 
1681  else if ( NDim == 3 )
1682  {
1683  M_curl[i][q]( 0 ) = M_grad[i][q]( 2,1 ) - M_grad[i][q]( 1,2 );
1684  M_curl[i][q]( 1 ) = M_grad[i][q]( 0,2 ) - M_grad[i][q]( 2,0 );
1685  M_curl[i][q]( 2 ) = M_grad[i][q]( 1,0 ) - M_grad[i][q]( 0,1 );
1686  }
1687  }
1688  }
1689  }
1690  }
1691 
1692  // we need the normal derivative
1693  if ( vm::has_first_derivative_normal<context>::value )
1694  {
1695  std::fill( M_dn.data(), M_dn.data()+M_dn.num_elements(), dn_type::Zero() );
1696  const uint16_type I = nDof*nComponents1;
1697  const uint16_type Q = nPoints();
1698 
1699  for ( int i = 0; i < I; ++i )
1700  for ( uint16_type q = 0; q < Q; ++q )
1701  {
1702  for ( uint16_type c1 = 0; c1 < NDim; ++c1 )
1703  {
1704  for ( uint16_type l = 0; l < NDim; ++l )
1705  {
1706  M_dn[i][q]( c1,0 ) += M_grad[i][q]( c1,l ) * thegmc->unitNormal( l, q );
1707  }
1708  }
1709  }
1710  }
1711  }
1712  }
1713  void update( geometric_mapping_context_ptrtype const& __gmc, mpl::int_<1>, mpl::bool_<true> )
1714  {
1715  //precompute_type* __pc = M_pc.get().get();
1716  geometric_mapping_context_type* thegmc = __gmc.get();
1717 
1718  if ( vm::has_grad<context>::value || vm::has_first_derivative<context>::value )
1719  {
1720  const uint16_type Q = M_npoints;//__gmc->nPoints();//M_grad.size2();
1721  const uint16_type I = nDof; //M_ref_ele->nbDof();
1722 
1723  typedef typename boost::multi_array<value_type,4>::index_range range;
1724 
1725  matrix_eigen_ublas_type Bt ( thegmc->B( 0 ).data().begin(), gmc_type::NDim, gmc_type::PDim );
1726  matrix_eigen_grad_type grad_real = matrix_eigen_grad_type::Zero();
1727  //matrix_eigen_PN_type B=Bt.transpose();
1728 
1729  for ( uint16_type ii = 0; ii < I; ++ii )
1730  {
1731  int ncomp= ( reference_element_type::is_product?nComponents1:1 );
1732 
1733  for ( uint16_type c = 0; c < ncomp; ++c )
1734  {
1735  //std::cout << "component " << c << "\n";
1736  uint16_type i = I*c + ii;
1737  //uint16_type c1 = c;
1738  grad_real.noalias() = M_gradphi[i][0]*Bt.transpose();
1739 
1740  for ( uint16_type q = 0; q < Q; ++q )
1741  {
1742  M_grad[i][q] = grad_real;
1743  M_dx[i][q] = M_grad[i][q].col( 0 );
1744 
1745  if ( NDim == 2 )
1746  M_dy[i][q] = M_grad[i][q].col( 1 );
1747 
1748  if ( NDim == 3 )
1749  M_dz[i][q] = M_grad[i][q].col( 2 );
1750 
1751  // update divergence if needed
1752  if ( vm::has_div<context>::value )
1753  {
1754  M_div[i][q]( 0,0 ) = M_grad[i][q].trace();
1755  }
1756 
1757  // update curl if needed
1758  if ( vm::has_curl<context>::value )
1759  {
1760  if ( NDim == 2 )
1761  {
1762 #if 0
1763  M_curl[i][q]( 0 ) = 0;
1764  M_curl[i][q]( 1 ) = 0;
1765 #else
1766  M_curl[i][q]( 0 ) = M_grad[i][q]( 1,0 ) - M_grad[i][q]( 0,1 );
1767  M_curl[i][q]( 1 ) = M_grad[i][q]( 1,0 ) - M_grad[i][q]( 0,1 );
1768 #endif
1769  M_curl[i][q]( 2 ) = M_grad[i][q]( 1,0 ) - M_grad[i][q]( 0,1 );
1770  }
1771 
1772  else if ( NDim == 3 )
1773  {
1774  M_curl[i][q]( 0 ) = M_grad[i][q]( 2,1 ) - M_grad[i][q]( 1,2 );
1775  M_curl[i][q]( 1 ) = M_grad[i][q]( 0,2 ) - M_grad[i][q]( 2,0 );
1776  M_curl[i][q]( 2 ) = M_grad[i][q]( 1,0 ) - M_grad[i][q]( 0,1 );
1777  }
1778  }
1779  } // q
1780  } // c
1781  } // ii
1782 
1783  // we need the normal derivative
1784  if ( vm::has_first_derivative_normal<context>::value )
1785  {
1786  std::fill( M_dn.data(), M_dn.data()+M_dn.num_elements(), dn_type::Zero() );
1787  const uint16_type I = nDof*nComponents1;
1788  const uint16_type Q = nPoints();
1789 
1790  for ( int i = 0; i < I; ++i )
1791  for ( uint16_type q = 0; q < Q; ++q )
1792  {
1793  for ( uint16_type c1 = 0; c1 < NDim; ++c1 )
1794  {
1795  for ( uint16_type l = 0; l < NDim; ++l )
1796  {
1797  M_dn[i][q]( c1,0 ) += M_grad[i][q]( c1,l ) * thegmc->unitNormal( l, q );
1798  }
1799  }
1800  }
1801  }
1802 
1803  } // grad
1804  }
1805 
1810  uint16_type nPoints() const
1811  {
1812  return M_npoints;
1813  }
1814 
1818  geometric_mapping_context_ptrtype const& gmContext() const
1819  {
1820  return M_gmc;
1821  }
1822 
1824  size_type eId() const
1825  {
1826  return M_gmc->id();
1827  }
1828 
1830  matrix_node_t_type const& xRefs() const
1831  {
1832  return M_gmc->xRefs();
1833  }
1834 
1836  precompute_ptrtype const& pc() const
1837  {
1838  return M_pc.get();
1839  }
1840 
1841  boost::multi_array<value_type,4> const& id() const
1842  {
1843  return M_phi;
1844  }
1845 
1846  value_type const& id( uint32_type i,
1847  uint16_type c1,
1848  uint16_type c2,
1849  uint32_type q ) const
1850  {
1851  return id( i, c1, c2, q, mpl::int_<rank>() );
1852  }
1853 
1854  value_type const& id( uint32_type i,
1855  uint16_type /*c1*/,
1856  uint16_type /*c2*/,
1857  uint32_type q,
1858  mpl::int_<0> ) const
1859  {
1860  return M_phi[i][q]( 0,0 );
1861  }
1862 
1863  value_type const& id( uint32_type i,
1864  uint16_type c1,
1865  uint16_type c2,
1866  uint32_type q,
1867  mpl::int_<1> ) const
1868  {
1869  detail::ignore_unused_variable_warning( c2 );
1870  return M_phi[i][q]( c1,0 );
1871  }
1872 
1873  id_type const& id( uint32_type i, uint32_type q ) const
1874  {
1875  return M_phi[i][q];
1876  }
1877 
1878  value_type const& id( uint32_type i,
1879  uint16_type c1,
1880  uint16_type c2,
1881  uint32_type q,
1882  mpl::int_<2> ) const
1883  {
1884  return M_phi[i][q]( c1,c2 );
1885  }
1886 
1887  value_type d( uint32_type i, uint16_type c1, uint16_type c2, uint32_type q ) const
1888  {
1889  return M_grad[i][q]( c1,c2 );
1890  }
1891 
1892  value_type dx( uint32_type i, uint16_type c1, uint16_type c2, uint32_type q ) const
1893  {
1894  return dx( i, c1, c2, q, mpl::int_<rank>() );
1895  }
1896  value_type dx( uint32_type i, uint16_type /*c1*/, uint16_type /*c2*/, uint32_type q, mpl::int_<0> ) const
1897  {
1898  BOOST_MPL_ASSERT_MSG( nDim >= 1, INVALID_DIM, ( mpl::int_<nDim>, mpl::int_<0> ) );
1899  return M_grad[i][q]( 0,0 );
1900  }
1901  value_type dx( uint32_type i, uint16_type c1, uint16_type /*c2*/, uint32_type q, mpl::int_<1> ) const
1902  {
1903  BOOST_MPL_ASSERT_MSG( nDim >= 1, INVALID_DIM, ( mpl::int_<nDim>, mpl::int_<0> ) );
1904  return M_grad[i][q]( c1,0 );
1905  }
1906  value_type dy( uint32_type i, uint16_type c1, uint16_type c2, uint32_type q ) const
1907  {
1908  return dy( i, c1, c2, q, mpl::int_<rank>() );
1909  }
1910  value_type dy( uint32_type i, uint16_type /*c1*/, uint16_type /*c2*/, uint32_type q, mpl::int_<0> ) const
1911  {
1912  BOOST_MPL_ASSERT_MSG( nDim >= 2, INVALID_DIM, ( mpl::int_<nDim>, mpl::int_<1> ) );
1913  return M_grad[i][q]( 0,1 );
1914  }
1915  value_type dy( uint32_type i, uint16_type c1, uint16_type /*c2*/, uint32_type q, mpl::int_<1> ) const
1916  {
1917  BOOST_MPL_ASSERT_MSG( nDim >= 2, INVALID_DIM, ( mpl::int_<nDim>, mpl::int_<1> ) );
1918  return M_grad[i][q]( c1,1 );
1919  }
1920  value_type dz( uint32_type i, uint16_type c1, uint16_type c2, uint32_type q ) const
1921  {
1922  return dz( i, c1, c2, q, mpl::int_<rank>() );
1923  }
1924  value_type dz( uint32_type i, uint16_type /*c1*/, uint16_type /*c2*/, uint32_type q, mpl::int_<0> ) const
1925  {
1926  BOOST_MPL_ASSERT_MSG( nDim >= 3, INVALID_DIM, ( mpl::int_<nDim>, mpl::int_<2> ) );
1927  return M_grad[i][q]( 0,2 );
1928  }
1929  value_type dz( uint32_type i, uint16_type c1, uint16_type /*c2*/, uint32_type q, mpl::int_<1> ) const
1930  {
1931  BOOST_MPL_ASSERT_MSG( nDim >= 3, INVALID_DIM, ( mpl::int_<nDim>, mpl::int_<2> ) );
1932  return M_grad[i][q]( c1,2 );
1933  }
1934 
1935 
1940  //matrix_type const& dn() const { return M_dn; }
1941 
1946  value_type const& dn( uint32_type i,
1947  uint16_type c1,
1948  uint16_type c2,
1949  uint32_type q ) const
1950  {
1951  return M_dn[i][q]( c1,c2 );
1952  }
1953 
1954  dn_type const& dn( uint16_type i, uint32_type q ) const
1955  {
1956  return M_dn[i][q];
1957  }
1958 
1959  value_type grad( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q ) const
1960  {
1961  return grad( i, c1, c2, q, mpl::int_<rank>() );
1962  }
1963  value_type grad( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> ) const
1964  {
1965  detail::ignore_unused_variable_warning( c1 );
1966  detail::ignore_unused_variable_warning( c2 );
1967  return M_grad[i][q]( 0,c2 );
1968  }
1969  value_type grad( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<1> ) const
1970  {
1971  return M_grad[i][q]( c1,c2 );
1972  }
1973 
1974  grad_type const& grad( uint16_type i, uint32_type q ) const
1975  {
1976  return M_grad[i][q];
1977  }
1978  dx_type const& dx( uint16_type i, uint32_type q ) const
1979  {
1980  return M_dx[i][q];
1981  }
1982  dy_type const& dy( uint16_type i, uint32_type q ) const
1983  {
1984  return M_dy[i][q];
1985  }
1986  dz_type const& dz( uint16_type i, uint32_type q ) const
1987  {
1988  return M_dz[i][q];
1989  }
1990 
1999  value_type div( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q ) const
2000  {
2001  detail::ignore_unused_variable_warning( c1 );
2002  detail::ignore_unused_variable_warning( c2 );
2003  return div( i, c1, c2, q, mpl::int_<rank>() );
2004  }
2005 
2006  value_type div( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> ) const
2007  {
2008  detail::ignore_unused_variable_warning( i );
2009  detail::ignore_unused_variable_warning( q );
2010  detail::ignore_unused_variable_warning( c1 );
2011  detail::ignore_unused_variable_warning( c2 );
2015  FEELPP_ASSERT( 0 ).error( "divergence of a scalar function is undefined." );
2016  return 0;
2017  }
2018 
2019  value_type div( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<1> ) const
2020  {
2021  detail::ignore_unused_variable_warning( c1 );
2022  detail::ignore_unused_variable_warning( c2 );
2023  return M_div[i][q]( 0,0 );
2024 #if 0
2025  value_type res = value_type( 0 );
2026 
2027  for ( int ii = 0; ii < nDim; ++ii )
2028  res += M_grad[i][ii]( ii,q );
2029 
2030  return res;
2031 #endif
2032  }
2033 
2042  value_type curl( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q ) const
2043  {
2044  detail::ignore_unused_variable_warning( c1 );
2045  detail::ignore_unused_variable_warning( c2 );
2046  return curl( i, c1, c2, q, mpl::int_<rank>() );
2047  }
2048 
2049  value_type curl( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> ) const
2050  {
2051  detail::ignore_unused_variable_warning( i );
2052  detail::ignore_unused_variable_warning( q );
2053  detail::ignore_unused_variable_warning( c1 );
2054  detail::ignore_unused_variable_warning( c2 );
2055  throw std::logic_error( "invalid use of curl operator, field must be vectorial" );
2056  return 0;
2057  }
2058 
2059  value_type curl( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<1> ) const
2060  {
2061  detail::ignore_unused_variable_warning( c2 );
2062  return M_curl[i][q]( c1 );
2063  }
2064  value_type curlx( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q ) const
2065  {
2066  detail::ignore_unused_variable_warning( c1 );
2067  detail::ignore_unused_variable_warning( c2 );
2068  return curlx( i, c1, c2, q, mpl::int_<rank>() );
2069  }
2070  value_type curlx( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> ) const
2071  {
2072  detail::ignore_unused_variable_warning(i);
2073  detail::ignore_unused_variable_warning(q);
2074  detail::ignore_unused_variable_warning(c1);
2075  detail::ignore_unused_variable_warning(c2);
2076  throw std::logic_error("invalid use of curl operator, field must be vectorial");
2077  return 0;
2078  }
2079  value_type curlx( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<1> ) const
2080  {
2081  detail::ignore_unused_variable_warning( c1 );
2082  detail::ignore_unused_variable_warning( c2 );
2083  return M_curl[i][q]( 0 );
2084  }
2085  value_type curly( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q ) const
2086  {
2087  detail::ignore_unused_variable_warning( c1 );
2088  detail::ignore_unused_variable_warning( c2 );
2089  return curly( i, c1, c2, q, mpl::int_<rank>() );
2090  }
2091  value_type curly( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> ) const
2092  {
2093  detail::ignore_unused_variable_warning(i);
2094  detail::ignore_unused_variable_warning(q);
2095  detail::ignore_unused_variable_warning(c1);
2096  detail::ignore_unused_variable_warning(c2);
2097  throw std::logic_error("invalid use of curl operator, field must be vectorial");
2098  return 0;
2099  }
2100  value_type curly( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<1> ) const
2101  {
2102  detail::ignore_unused_variable_warning( c1 );
2103  detail::ignore_unused_variable_warning( c2 );
2104  return M_curl[i][q]( 1 );
2105  }
2106 
2107  value_type curlz( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q ) const
2108  {
2109  detail::ignore_unused_variable_warning( c1 );
2110  detail::ignore_unused_variable_warning( c2 );
2111  return curlz( i, c1, c2, q, mpl::int_<rank>() );
2112  }
2113  value_type curlz( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> ) const
2114  {
2115  detail::ignore_unused_variable_warning(i);
2116  detail::ignore_unused_variable_warning(q);
2117  detail::ignore_unused_variable_warning(c1);
2118  detail::ignore_unused_variable_warning(c2);
2119  throw std::logic_error("invalid use of curl operator, field must be vectorial");
2120  return 0;
2121  }
2122  value_type curlz( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<1> ) const
2123  {
2124  detail::ignore_unused_variable_warning( c1 );
2125  detail::ignore_unused_variable_warning( c2 );
2126  return M_curl[i][q]( 2 );
2127  }
2128 
2129  value_type hess( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q ) const
2130  {
2131  return hess( i, c1, c2, q, mpl::int_<rank>() );
2132  }
2133  value_type hess( uint16_type i, uint16_type c1, uint16_type c2, uint32_type q, mpl::int_<0> ) const
2134  {
2135  return M_hessian[i][q]( c1,c2 );
2136  }
2137 
2138  // private:
2139  Context() {}
2140  Context( Context const& ) {}
2141  private:
2142 
2143  boost::optional<precompute_ptrtype> M_pc;
2144  uint16_type M_npoints;
2145 
2146  uint16_type M_ipt;
2147 
2148  reference_element_ptrtype M_ref_ele;
2149 
2150  geometric_mapping_context_ptrtype M_gmc;
2151 
2152  boost::multi_array<id_type,2> M_phi;
2153  boost::multi_array<ref_grad_type,2> M_gradphi;
2154  boost::multi_array<hess_type,2> M_hessphi;
2155  boost::multi_array<dn_type,2> M_dn;
2156  boost::multi_array<grad_type,2> M_grad;
2157  boost::multi_array<dx_type,2> M_dx;
2158  boost::multi_array<dy_type,2> M_dy;
2159  boost::multi_array<dz_type,2> M_dz;
2160  boost::multi_array<div_type,2> M_div;
2161  boost::multi_array<curl_type,2> M_curl;
2162  boost::multi_array<hess_type,2> M_hessian;
2163  };
2164 
2165  template<size_type context_v, size_type context_g, typename BasisType, typename GeoType, typename ElementType>
2166  boost::shared_ptr<Context<context_v,BasisType, GeoType, ElementType> >
2167  context( boost::shared_ptr<BasisType> b, boost::shared_ptr<GeoType> gm, precompute_ptrtype& pc )
2168  {
2169  return boost::shared_ptr<Context<context_v,BasisType, GeoType, ElementType> >(
2170  new Context<context_v, BasisType, GeoType, ElementType>( context_v,
2171  b,
2172  gm,
2173  pc ) );
2174  }
2175 
2176  template<int contextv, int contextg, typename BasisType, typename GeoType,typename ElementType>
2177  boost::shared_ptr<Context<contextv,BasisType, GeoType, ElementType> >
2178  ctx( boost::shared_ptr<BasisType> const& b,
2179  boost::shared_ptr<typename GeoType::template Context<contextg, ElementType> > const& gm,
2180  precompute_ptrtype pc, ElementType& e )
2181  {
2182  typedef Context<contextv,BasisType, GeoType, ElementType> ctx_type;
2183  return boost::shared_ptr<ctx_type>( new ctx_type( b, gm, pc ) );
2184 
2185  }
2186  template<int contextv, typename BasisType, typename GeoType,typename ElementType>
2187  boost::shared_ptr<Context<contextv,BasisType, GeoType, ElementType> >
2188  ctx( boost::shared_ptr<BasisType> const& b,
2189  boost::shared_ptr<typename GeoType::template Context<contextv, ElementType> > const& gm,
2190  precompute_ptrtype pc, ElementType& e )
2191  {
2192  typedef Context<contextv,BasisType, GeoType, ElementType> ctx_type;
2193  return boost::shared_ptr<ctx_type>( new ctx_type( b, gm, pc ) );
2194 
2195  }
2196 
2197 protected:
2198 
2199 private:
2200 
2201 private:
2202 
2203  basis_type M_basis;
2204  matrix_type M_coeff;
2205  std::string M_fname;
2206 };
2207 
2208 template<typename Poly,template<uint16_type> class PolySetType> const bool PolynomialSet<Poly,PolySetType>::is_scalar;
2209 template<typename Poly,template<uint16_type> class PolySetType> const bool PolynomialSet<Poly,PolySetType>::is_vectorial;
2210 template<typename Poly,template<uint16_type> class PolySetType> const bool PolynomialSet<Poly,PolySetType>::is_tensor2;
2211 template<typename Poly,template<uint16_type> class PolySetType> const uint16_type PolynomialSet<Poly,PolySetType>::nComponents;
2212 template<typename Poly,template<uint16_type> class PolySetType> const uint16_type PolynomialSet<Poly,PolySetType>::nComponents1;
2213 template<typename Poly,template<uint16_type> class PolySetType> const uint16_type PolynomialSet<Poly,PolySetType>::nComponents2;
2214 
2215 } // Feel
2216 
2218 
2219 #endif /* __PolynomialSet_H */

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