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
expr.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-01-17
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2006,2007 Université Joseph Fourier (Grenoble I)
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 __Expr_H
31 #define __Expr_H 1
32 
33 #undef max
34 #include <boost/version.hpp>
35 #if (BOOST_VERSION >= 103400)
36 #include <boost/none.hpp>
37 #else
38 #include <boost/none_t.hpp>
39 #endif /* BOOST_VERSION >= 103400 */
40 
41 #include <algorithm>
42 #include <boost/type_traits/is_base_of.hpp>
43 #include <boost/type_traits/is_same.hpp>
44 #include <boost/type_traits/is_class.hpp>
45 #include <boost/static_assert.hpp>
46 #include <boost/foreach.hpp>
47 #include <boost/fusion/sequence.hpp>
48 #include <boost/fusion/support/pair.hpp>
49 #include <boost/multi_array.hpp>
50 
51 #include <Eigen/Core>
52 
54 #include <feel/feelpoly/policy.hpp>
55 #include <feel/feelpoly/context.hpp>
56 #include <feel/feelvf/shape.hpp>
57 
58 namespace Feel
59 {
60 namespace vf
61 {
62 
64 typedef node<double>::type node_type;
65 
66 enum
67 {
68  CONTEXT_1 = ( 1<<0 ),
69  CONTEXT_2 = ( 1<<1 )
70 };
71 
72 template<typename ExprT>
73 class ComponentsExpr
74 {
75 public:
76 
77  static const size_type context = ExprT::context;
78  static const bool is_terminal = false;
79  //integration order
80  static const uint16_type imorder = ExprT::imorder;
81  //the expression is a polynomial type?
82  static const bool imIsPoly = ExprT::imIsPoly;
83 
84 
85  template<typename Func>
86  struct HasTestFunction
87  {
88  static const bool result = ExprT::template HasTestFunction<Func>::result;
89  };
90 
91  template<typename Func>
92  struct HasTrialFunction
93  {
94  static const bool result = ExprT::template HasTrialFunction<Func>::result;
95  };
96 
97 
98 
99 
103 
104  typedef ExprT expression_type;
105  typedef typename expression_type::value_type value_type;
106  typedef ComponentsExpr<ExprT> this_type;
107 
109 
113 
114  ComponentsExpr()
115  :
116  M_expr()
117  {}
118 
119  explicit ComponentsExpr( expression_type const & __expr, int c1, int c2 )
120  :
121  M_expr( __expr ),
122  M_c1( c1 ),
123  M_c2( c2 )
124  {}
125  ~ComponentsExpr()
126  {}
127 
129 
130  expression_type const& expression() const
131  {
132  return M_expr;
133  }
134 
138 
139  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t = Basis_i_t>
140  struct tensor
141  {
142 
143  typedef typename expression_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> tensor_expr_type;
144  typedef typename tensor_expr_type::value_type value_type;
145  typedef typename mpl::if_<fusion::result_of::has_key<Geo_t, vf::detail::gmc<0> >,
146  mpl::identity<vf::detail::gmc<0> >,
147  mpl::identity<vf::detail::gmc<1> > >::type::type key_type;
148  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type gmc_type;
149  typedef Shape<gmc_type::NDim, Scalar, false> shape;
150 
151  template <class Args> struct sig
152  {
153  typedef value_type type;
154  };
155 
156  struct is_zero
157  {
158  static const bool value = tensor_expr_type::is_zero::value;
159  };
160 
161  tensor( this_type const& expr,
162  Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
163  :
164  M_tensor_expr( expr.expression(), geom, fev, feu ),
165  M_c1( expr.M_c1 ),
166  M_c2( expr.M_c2 )
167  {}
168 
169  tensor( this_type const& expr,
170  Geo_t const& geom, Basis_i_t const& fev )
171  :
172  M_tensor_expr( expr.expression(), geom, fev ),
173  M_c1( expr.M_c1 ),
174  M_c2( expr.M_c2 )
175  {}
176 
177  tensor( this_type const& expr, Geo_t const& geom )
178  :
179  M_tensor_expr( expr.expression(), geom ),
180  M_c1( expr.M_c1 ),
181  M_c2( expr.M_c2 )
182  {
183  }
184 
185  template<typename IM>
186  void init( IM const& im )
187  {
188  M_tensor_expr.init( im );
189  }
190  void update( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
191  {
192  M_tensor_expr.update( geom, fev, feu );
193  }
194  void update( Geo_t const& geom, Basis_i_t const& fev )
195  {
196  M_tensor_expr.update( geom, fev );
197  }
198  void update( Geo_t const& geom )
199  {
200  M_tensor_expr.update( geom );
201  }
202  void update( Geo_t const& geom, uint16_type face )
203  {
204  M_tensor_expr.update( geom, face );
205  }
206 
207 
208  value_type
209  evalij( uint16_type i, uint16_type j ) const
210  {
211  return M_tensor_expr.evalij( i, j );
212  }
213 
214 
215  value_type
216  evalijq( uint16_type i, uint16_type j, uint16_type /*c1*/, uint16_type /*c2*/, uint16_type q ) const
217  {
218  return M_tensor_expr.evalijq( i, j, M_c1, M_c2, q );
219  }
220 
221  template<int PatternContext>
222  value_type
223  evalijq( uint16_type i, uint16_type j, uint16_type /*c1*/, uint16_type /*c2*/, uint16_type q,
224  mpl::int_<PatternContext> ) const
225  {
226  return M_tensor_expr.evalijq( i, j, M_c1, M_c2, q, mpl::int_<PatternContext>() );
227  }
228 
229 
230  value_type
231  evaliq( uint16_type i, uint16_type /*c1*/, uint16_type /*c2*/, uint16_type q ) const
232  {
233  return M_tensor_expr.evaliq( i, M_c1, M_c2, q );
234  }
235 
236  value_type
237  evalq( uint16_type /*c1*/, uint16_type /*c2*/, uint16_type q ) const
238  {
239  return M_tensor_expr.evalq( M_c1, M_c2, q );
240  }
241 
242  tensor_expr_type M_tensor_expr;
243  const int M_c1, M_c2;
244  };
245  expression_type M_expr;
246  int M_c1, M_c2;
247 };
248 class CstBase {};
249 class IntegratorBase {};
250 class LambdaExprBase {};
251 class LambdaExpr1 : public LambdaExprBase
252 {
253 public:
254 
255  static const size_type context = 0;
256  static const bool is_terminal = false;
257 
258  static const uint16_type imorder = 0;
259  static const bool imIsPoly = false;
260 
261  template<typename Func>
262  struct HasTestFunction
263  {
264  static const bool result = false;
265  };
266  template<typename Func>
267  struct HasTrialFunction
268  {
269  static const bool result = false;
270  };
271 
272  typedef double value_type;
273 
274  template<typename TheExpr>
275  struct Lambda
276  {
277  typedef typename TheExpr::expression_type type;
278  };
279 
280  template<typename ExprT>
281  typename Lambda<ExprT>::type
282  operator()( ExprT const& e ) {
283  return e.expression();
284  }
285 
286  template<typename ExprT>
287  typename Lambda<ExprT>::type
288  operator()( ExprT const& e ) const { return e.expression(); }
289 
290  template<typename Geo_t, typename Basis_i_t = fusion::map<fusion::pair<vf::detail::gmc<0>,boost::shared_ptr<vf::detail::gmc<0> > >,fusion::pair<vf::detail::gmc<1>,boost::shared_ptr<vf::detail::gmc<1> > > >, typename Basis_j_t = Basis_i_t>
291  struct tensor
292  {
293  typedef LambdaExpr1 expression_type;
294  typedef typename LambdaExpr1::value_type value_type;
295 
296  typedef typename mpl::if_<fusion::result_of::has_key<Geo_t, vf::detail::gmc<0> >, mpl::identity<vf::detail::gmc<0> >, mpl::identity<vf::detail::gmc<1> > >::type::type key_type;
297  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type* gmc_ptrtype;
298  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type gmc_type;
299  typedef Shape<gmc_type::nDim, Scalar, false, false> shape;
300 
301 
302  template<typename Indq, typename Indi, typename Indj>
303  struct expr
304  {
305  typedef value_type type;
306  };
307 
308  struct is_zero
309  {
310  static const bool value = false;
311  };
312 
313  tensor( expression_type const& expr,
314  Geo_t const& /*geom*/, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
315  {
316  }
317  tensor( expression_type const& expr,
318  Geo_t const& /*geom*/, Basis_i_t const& /*fev*/ )
319  {
320  }
321  tensor( expression_type const& expr, Geo_t const& /*geom*/ )
322  {
323  }
324  template<typename IM>
325  void init( IM const& /*im*/ )
326  {
327  }
328  void update( Geo_t const&, Basis_i_t const& , Basis_j_t const& )
329  {
330  }
331  void update( Geo_t const& , Basis_i_t const& )
332  {
333  }
334  void update( Geo_t const& )
335  {
336  }
337  void update( Geo_t const&, uint16_type )
338  {
339  }
340 
341  value_type
342  evalij( uint16_type /*i*/, uint16_type /*j*/ ) const
343  {
344  return 0;
345  }
346 
347 
348  value_type
349  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type /*c1*/, uint16_type /*c2*/, uint16_type /*q*/ ) const
350  {
351  return 0;
352  }
353  template<int PatternContext>
354  value_type
355  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type /*c1*/, uint16_type /*c2*/, uint16_type /*q*/,
356  mpl::int_<PatternContext> ) const
357  {
358  return 0;
359  }
360 
361  value_type
362  evaliq( uint16_type /*i*/, uint16_type /*c1*/, uint16_type /*c2*/, uint16_type /*q*/ ) const
363  {
364  return 0;
365  }
366  value_type
367  evalq( uint16_type /*c1*/, uint16_type /*c2*/, uint16_type /*q*/ ) const
368  {
369  return 0;
370  }
371  };
372 
373 };
374 
375 
376 
384 template<typename ExprT>
385 class Expr//: public boost::enable_shared_from_this<Expr<ExprT> >
386 {
387 public:
388 
389  static const size_type context = ExprT::context;
390  static const bool is_terminal = ExprT::is_terminal;
391 
392  //integration order
393  static const uint16_type imorder = ExprT::imorder;
394  //the expression is a polynomial type?
395  static const bool imIsPoly = ExprT::imIsPoly;
396 
397 
398  template<typename Func>
399  struct HasTestFunction
400  {
401  static const bool result = ExprT::template HasTestFunction<Func>::result;
402  };
403 
404  template<typename Func>
405  struct HasTrialFunction
406  {
407  static const bool result = ExprT::template HasTrialFunction<Func>::result;
408  };
409 
410 
411 
412 
416 
417  typedef ExprT expression_type;
418  typedef typename expression_type::value_type value_type;
419  typedef Expr<ExprT> this_type;
420  typedef boost::shared_ptr<this_type> this_ptrtype;
422 
426 
427  Expr()
428  :
429  M_expr()
430  {}
431 
432  explicit Expr( expression_type const & __expr )
433  :
434  M_expr( __expr )
435  {}
436  ~Expr()
437  {}
438 
440 
444  Expr<ComponentsExpr<Expr<ExprT> > >
445  operator()( int c1 = 0, int c2 = 0 ) const
446  {
447  auto ex = ComponentsExpr<Expr<ExprT> >( Expr<ExprT>( M_expr ), c1, c2 );
448  return Expr<ComponentsExpr<Expr<ExprT> > >( ex );
449  }
450 
451  template<typename TheExpr>
452  struct Lambda
453  {
454  typedef typename ExprT::template Lambda<TheExpr>::type expr_type;
455  typedef Expr<expr_type> type;
456  //typedef expr_type type;
457  };
458 
459 
460 
461  template<typename TheExpr>
462  typename Lambda<TheExpr>::type
463  operator()( TheExpr const& e )
464  {
465  //typename Lambda<TheExpr>::expr_type e1( M_expr(e) );
466  //typename Lambda<TheExpr>::type r( Expr(e1 ) );
467  //return r;
468  return expr( M_expr( e ) );
469  }
470 
471  template<typename TheExpr>
472  typename Lambda<TheExpr>::type
473  operator()( TheExpr const& e ) const { return expr(M_expr(e)); }
474 
475 
476  template<typename Geo_t, typename Basis_i_t = fusion::map<fusion::pair<vf::detail::gmc<0>,boost::shared_ptr<vf::detail::gmc<0> > >,fusion::pair<vf::detail::gmc<1>,boost::shared_ptr<vf::detail::gmc<1> > > >, typename Basis_j_t = Basis_i_t>
477  struct tensor
478  {
479 
480  typedef typename expression_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> tensor_expr_type;
481  typedef typename tensor_expr_type::value_type value_type;
482 
483  typedef typename tensor_expr_type::shape shape;
484 
485  template <class Args> struct sig
486  {
487  typedef value_type type;
488  };
489 
490  struct is_zero
491  {
492  static const bool value = tensor_expr_type::is_zero::value;
493  };
494 
495  tensor( this_type const& expr,
496  Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
497  :
498  M_tensor_expr( expr.expression(), geom, fev, feu )
499  {}
500 
501  tensor( this_type const& expr,
502  Geo_t const& geom, Basis_i_t const& fev )
503  :
504  M_tensor_expr( expr.expression(), geom, fev )
505  {}
506 
507  tensor( this_type const& expr, Geo_t const& geom )
508  :
509  M_tensor_expr( expr.expression(), geom )
510  {
511  }
512 
513  template<typename IM>
514  void init( IM const& im )
515  {
516  M_tensor_expr.init( im );
517  }
518  void update( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
519  {
520  M_tensor_expr.update( geom, fev, feu );
521  }
522  void update( Geo_t const& geom, Basis_i_t const& fev )
523  {
524  M_tensor_expr.update( geom, fev );
525  }
526  void update( Geo_t const& geom )
527  {
528  M_tensor_expr.update( geom );
529  }
530  void update( Geo_t const& geom, uint16_type face )
531  {
532  M_tensor_expr.update( geom, face );
533  }
534 
535 
536  value_type
537  evalij( uint16_type i, uint16_type j ) const
538  {
539  return M_tensor_expr.evalij( i, j );
540  }
541 
542  Eigen::Matrix<value_type, shape::M, shape::N> const&
543  evalijq( uint16_type i, uint16_type j, uint16_type q ) const
544  {
545  return M_tensor_expr.evalijq( i, j, q );
546  }
547 
548  value_type
549  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q ) const
550  {
551  return M_tensor_expr.evalijq( i, j, c1, c2, q );
552  }
553 
554  template<int PatternContext>
555  value_type
556  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
557  mpl::int_<PatternContext> ) const
558  {
559  return M_tensor_expr.evalijq( i, j, c1, c2, q, mpl::int_<PatternContext>() );
560  }
561 
562 
563  value_type
564  evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q ) const
565  {
566  return M_tensor_expr.evaliq( i, c1, c2, q );
567  }
568 
569  value_type
570  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const
571  {
572  return M_tensor_expr.evalq( c1, c2, q );
573  }
574 
575  tensor_expr_type M_tensor_expr;
576  };
577 #if 0
578  class Diff
579  {
580  public:
581 
582  value_type value() const
583  {
584  return __expression.value();
585  }
586 
587  value_type grad( int __ith ) const
588  {
589  return __expression.grad( __ith );
590  }
591 
592  value_type hessian( int __i, int __j ) const
593  {
594  return __expression.hessian( __i, __j );
595  }
596 
597  };
598 #endif
599 
600 
604 
605  bool isSymetric() const
606  {
607  return M_expr.isSymetric();
608  }
609 
610  //this_ptrtype ptr() { return boost::shared_from_this(); }
611 
612  expression_type const& expression() const
613  {
614  return M_expr;
615  }
616 
618 
622 
624 
628 
629  template<typename Elem1, typename Elem2, typename FormType>
630  void assemble( boost::shared_ptr<Elem1> const& __u,
631  boost::shared_ptr<Elem2> const& __v,
632  FormType& __f ) const
633  {
634  DVLOG(2) << "calling assemble(u,v)\n";
635  M_expr.assemble( __u, __v, __f );
636  DVLOG(2) << "calling assemble(u,v) done\n";
637  }
638 
639  template<typename Elem1, typename FormType>
640  void assemble( boost::shared_ptr<Elem1> const& __v,
641  FormType& __f ) const
642  {
643  DVLOG(2) << "calling assemble(v)\n";
644  M_expr.assemble( __v, __f );
645  DVLOG(2) << "calling assemble(v) done\n";
646  }
647 
648  template<typename P0hType>
649  typename P0hType::element_type
650  broken( boost::shared_ptr<P0hType>& P0h ) const
651  {
652  return M_expr.broken( P0h );
653  }
654  //__typeof__( M_expr.evaluate() )
655  //ublas::matrix<typename expression_type::value_type>
656 
657  typename expression_type::value_type
658  evaluate( bool parallel = true, WorldComm const& worldcomm = Environment::worldComm() ) const
659  {
660  return M_expr.evaluate( parallel,worldcomm );
661  }
662 
663  typename expression_type::value_type
664  evaluateAndSum() const
665  {
666  return M_expr.evaluateAndSum();
667  }
668  std::string expressionStr() const
669  {
670  return std::string();
671  //return M_expr.expressionStr();
672  }
673 
674 
676 
677 protected:
678 
679 private:
680 
681  mutable expression_type M_expr;
682 };
683 
684 template <typename ExprT>
685 Expr<ExprT>
686 expr( ExprT const& exprt )
687 {
688  return Expr<ExprT>( exprt );
689 }
690 
691 template <typename ExprT>
692 boost::shared_ptr<Expr<ExprT> >
693 exprPtr( ExprT const& exprt )
694 {
695  return boost::shared_ptr<Expr<ExprT> >( new Expr<ExprT>( exprt ) );
696 }
697 
698 extern Expr<LambdaExpr1> _e1;
699 
711 template<typename IntElts,typename ExprT>
712 struct ExpressionOrder
713 {
714 
715  typedef typename boost::tuples::template element<1, IntElts>::type element_iterator_type;
716  typedef typename boost::remove_reference<typename element_iterator_type::reference>::type const_t;
717  typedef typename boost::remove_const<const_t>::type the_face_element_type;
718  typedef typename the_face_element_type::super2::template Element<the_face_element_type>::type the_element_type;
719 
720  static const uint16_type nOrderGeo = the_element_type::nOrder;
721 
722  static const bool is_polynomial = ExprT::imIsPoly;
723 #if 0
724  static const int value = boost::mpl::if_< boost::mpl::bool_< ExprT::imIsPoly > ,
725  typename boost::mpl::if_< boost::mpl::greater< boost::mpl::int_<ExprT::imorder>,
726  boost::mpl::int_<19> > ,
727  boost::mpl::int_<19>,
728  boost::mpl::int_<ExprT::imorder> >::type,
729  boost::mpl::int_<10> >::type::value;
730 #else
731  // this is a very rough approximation
732  static const int value = ( ExprT::imorder )?( ExprT::imorder*nOrderGeo ):( nOrderGeo );
733  static const int value_1 = ExprT::imorder;
734 #endif
735 
736 
737 };
738 
739 
740 template<typename PrintExprT>
741 class PrintExpr
742 {
743 public:
744 
745  static const size_type context = PrintExprT::context;
746  static const bool is_terminal = false;
747 
748  static const uint16_type imorder = PrintExprT::imorder;
749  static const bool imIsPoly = PrintExprT::imIsPoly;
750 
751  template<typename Func>
752  struct HasTestFunction
753  {
754  static const bool result = PrintExprT::template HasTestFunction<Func>::result;
755  };
756  template<typename Func>
757  struct HasTrialFunction
758  {
759  static const bool result = PrintExprT::template HasTrialFunction<Func>::result;
760  };
761 
765 
766  typedef PrintExprT expression_type;
767  typedef typename expression_type::value_type value_type;
768  typedef PrintExpr<PrintExprT> this_type;
769 
771 
775 
776  explicit PrintExpr( expression_type const & __expr,
777  std::string const & __tag )
778  :
779  M_expr( __expr ),
780  M_tag( __tag )
781  {}
782  ~PrintExpr()
783  {}
784 
786 
790 
791  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t>
792  struct tensor
793  {
794  typedef typename expression_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> tensor_expr_type;
795  typedef typename tensor_expr_type::value_type value_type;
796 
797  typedef typename tensor_expr_type::shape shape;
798 
799  template <class Args> struct sig
800  {
801  typedef value_type type;
802  };
803  struct is_zero
804  {
805  static const bool value = tensor_expr_type::is_zero::value;
806  };
807 
808  tensor( this_type const& expr,
809  Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
810  :
811  M_tensor_expr( expr.expression(), geom, fev, feu ),
812  M_tag( expr.tag() )
813  {}
814 
815  tensor( this_type const& expr,
816  Geo_t const& geom, Basis_i_t const& fev )
817  :
818  M_tensor_expr( expr.expression(), geom, fev ),
819  M_tag( expr.tag() )
820  {}
821 
822  tensor( this_type const& expr, Geo_t const& geom )
823  :
824  M_tensor_expr( expr.expression(), geom ),
825  M_tag( expr.tag() )
826  {
827  }
828  template<typename IM>
829  void init( IM const& im )
830  {
831  M_tensor_expr.init( im );
832  }
833  void update( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
834  {
835  M_tensor_expr.update( geom, fev, feu );
836  }
837  void update( Geo_t const& geom, Basis_i_t const& fev )
838  {
839  M_tensor_expr.update( geom, fev );
840  }
841  void update( Geo_t const& geom )
842  {
843  M_tensor_expr.update( geom );
844  }
845  void update( Geo_t const& geom, uint16_type face )
846  {
847  M_tensor_expr.update( geom, face );
848  }
849 
850 
851  value_type
852  evalij( uint16_type i, uint16_type j ) const
853  {
854  return M_tensor_expr.evalij( i, j );
855  }
856 
857 
858  value_type
859  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q ) const
860  {
861  value_type res= M_tensor_expr.evalijq( i, j, c1, c2, q );
862  std::cout << "[print] " << M_tag << " shape(" << shape::M << "," << shape::N << ") evalijq( " << i << "," << j << "," << c1 << "," << c2 << "," << q << ")=" << res << "\n";
863  return res;
864  }
865  template<int PatternContext>
866  value_type
867  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
868  mpl::int_<PatternContext> ) const
869  {
870  value_type res= M_tensor_expr.evalijq( i, j, c1, c2, q, mpl::int_<PatternContext>() );
871  std::cout << "[print] " << M_tag << " shape(" << shape::M << "," << shape::N << ") evalijq( " << i << "," << j << "," << c1 << "," << c2 << "," << q << ")=" << res << "\n";
872  return res;
873  }
874 
875 
876 
877  value_type
878  evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q ) const
879  {
880  value_type res= M_tensor_expr.evaliq( i, c1, c2, q );
881  std::cout << "[print] " << M_tag << " shape(" << shape::M << "," << shape::N << ") evaliq( " << i << "," << c1 << "," << c2 << "," << q << ")=" << res << "\n";
882  return res;
883  }
884 
885  value_type
886  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const
887  {
888  value_type res= M_tensor_expr.evalq( c1, c2, q );
889  std::cout << "[print] " << M_tag << " shape(" << shape::M << "," << shape::N << ") evalq( " << c1 << "," << c2 << "," << q << ")=" << res << "\n";
890  return res;
891  }
892 
893  tensor_expr_type M_tensor_expr;
894  std::string M_tag;
895  };
896 #if 0
897  class Diff
898  {
899  public:
900 
901  value_type value() const
902  {
903  return __expression.value();
904  }
905 
906  value_type grad( int __ith ) const
907  {
908  return __expression.grad( __ith );
909  }
910 
911  value_type hessian( int __i, int __j ) const
912  {
913  return __expression.hessian( __i, __j );
914  }
915 
916  };
917 #endif
918 
919 
923 
924  bool isSymetric() const
925  {
926  return M_expr.isSymetric();
927  }
928 
929  expression_type const& expression() const
930  {
931  return M_expr;
932  }
933 
934  const std::string& tag() const
935  {
936  return M_tag;
937  }
938 
940 
944 
946 
950 
951  template<typename Elem1, typename Elem2, typename FormType>
952  void assemble( boost::shared_ptr<Elem1> const& __u,
953  boost::shared_ptr<Elem2> const& __v,
954  FormType& __f ) const
955  {
956  DVLOG(2) << "calling assemble(u,v)\n";
957  M_expr.assemble( __u, __v, __f );
958  DVLOG(2) << "calling assemble(u,v) done\n";
959  }
960 
961  template<typename Elem1, typename FormType>
962  void assemble( boost::shared_ptr<Elem1> const& __v,
963  FormType& __f ) const
964  {
965  DVLOG(2) << "calling assemble(v)\n";
966  M_expr.assemble( __v, __f );
967  DVLOG(2) << "calling assemble(v) done\n";
968  }
969 #if 0
970  //__typeof__( M_expr.evaluate() )
971  ublas::matrix<typename expression_type::value_type>
972  evaluate() const
973  {
974  return M_expr.evaluate();
975  }
976 #endif
977 
978 
980 
981 protected:
982 
983 private:
984 
985  mutable expression_type M_expr;
986  const std::string M_tag;
987 };
988 template<typename ExprT>
989 inline
990 Expr< PrintExpr<ExprT> >
991 print( ExprT v, std::string t="" )
992 {
993  typedef PrintExpr<ExprT> print_t;
994  return Expr< print_t >( print_t( v, t ) );
995 }
996 
1004 template<typename ExprT>
1005 class Trans
1006 {
1007 public:
1008 
1009  static const size_type context = ExprT::context;
1010  static const bool is_terminal = false;
1011 
1012  static const uint16_type imorder = ExprT::imorder;
1013  static const bool imIsPoly = ExprT::imIsPoly;
1014 
1015  template<typename Func>
1016  struct HasTestFunction
1017  {
1018  static const bool result = ExprT::template HasTestFunction<Func>::result;
1019  };
1020 
1021  template<typename Func>
1022  struct HasTrialFunction
1023  {
1024  static const bool result = ExprT::template HasTrialFunction<Func>::result;
1025  };
1026 
1030 
1031  typedef ExprT expression_type;
1032  typedef typename expression_type::value_type value_type;
1033  typedef Trans<ExprT> this_type;
1034 
1036 
1040 
1041  explicit Trans( expression_type const & __expr )
1042  :
1043  M_expr( __expr )
1044  {}
1045  ~Trans()
1046  {}
1047 
1049 
1053  template<typename TheExpr>
1054  struct Lambda
1055  {
1056  typedef Trans<typename expression_type::template Lambda<TheExpr>::type> type;
1057  };
1058  template<typename TheExpr>
1059  typename Lambda<TheExpr>::type
1060  operator()( TheExpr const& e ) { return trans(M_expr(e)); }
1061 
1062  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t>
1063  struct tensor
1064  {
1065  typedef typename expression_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> tensor_expr_type;
1066  typedef typename tensor_expr_type::value_type value_type;
1067 
1068  typedef typename Transpose<typename tensor_expr_type::shape>::type shape;
1069 
1070  template <class Args> struct sig
1071  {
1072  typedef value_type type;
1073  };
1074 
1075  struct is_zero
1076  {
1077  static const bool value = tensor_expr_type::is_zero::value;
1078  };
1079 
1080  tensor( this_type const& expr,
1081  Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
1082  :
1083  M_tensor_expr( expr.expression(), geom, fev, feu )
1084  {}
1085 
1086  tensor( this_type const& expr,
1087  Geo_t const& geom, Basis_i_t const& fev )
1088  :
1089  M_tensor_expr( expr.expression(), geom, fev )
1090  {}
1091 
1092  tensor( this_type const& expr, Geo_t const& geom )
1093  :
1094  M_tensor_expr( expr.expression(), geom )
1095  {
1096  }
1097  template<typename IM>
1098  void init( IM const& im )
1099  {
1100  M_tensor_expr.init( im );
1101  }
1102 
1103  void update( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
1104  {
1105  M_tensor_expr.update( geom, fev, feu );
1106  }
1107  void update( Geo_t const& geom, Basis_i_t const& fev )
1108  {
1109  M_tensor_expr.update( geom, fev );
1110  }
1111  void update( Geo_t const& geom )
1112  {
1113  M_tensor_expr.update( geom );
1114  }
1115  void update( Geo_t const& geom, uint16_type face )
1116  {
1117  M_tensor_expr.update( geom, face );
1118  }
1119 
1120 
1121  value_type
1122  evalij( uint16_type i, uint16_type j ) const
1123  {
1124  return M_tensor_expr.evalij( i, j );
1125  }
1126 
1127 
1128  value_type
1129  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q ) const
1130  {
1131  return M_tensor_expr.evalijq( i, j, c2, c1, q );
1132  }
1133  template<int PatternContext>
1134  value_type
1135  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
1136  mpl::int_<PatternContext> ) const
1137  {
1138  return M_tensor_expr.evalijq( i, j, c2, c1, q, mpl::int_<PatternContext>() );
1139  }
1140 
1141 
1142  value_type
1143  evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q ) const
1144  {
1145  return M_tensor_expr.evaliq( i, c2, c1, q );
1146  }
1147 
1148  value_type
1149  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const
1150  {
1151  return M_tensor_expr.evalq( c2, c1, q );
1152  }
1153 
1154  tensor_expr_type M_tensor_expr;
1155  };
1156 
1158 
1162 
1163  bool isSymetric() const
1164  {
1165  return M_expr.isSymetric();
1166  }
1167 
1168  expression_type const& expression() const
1169  {
1170  return M_expr;
1171  }
1172 
1174 
1178 
1180 
1184 
1186 
1187 protected:
1188 
1189 private:
1190 
1191  mutable expression_type M_expr;
1192 };
1193 
1194 template<typename ExprT>
1195 inline
1196 Expr< Trans<ExprT> >
1197 trans( ExprT v )
1198 {
1199  typedef Trans<ExprT> trans_t;
1200  return Expr< trans_t >( trans_t( v ) );
1201 }
1202 
1203 
1204 template < class T>
1205 class Cst : public CstBase
1206 {
1207 public:
1208 
1209  //BOOST_STATIC_ASSERT( ::boost::is_arithmetic<T>::value );
1210 
1211  static const size_type context = vm::JACOBIAN;
1212  static const bool is_terminal = false;
1213 
1214  static const uint16_type imorder = 0;
1215  static const bool imIsPoly = true;
1216 
1217  template<typename Func>
1218  struct HasTestFunction
1219  {
1220  static const bool result = false;
1221  };
1222  template<typename Func>
1223  struct HasTrialFunction
1224  {
1225  static const bool result = false;
1226  };
1227 
1228  typedef typename mpl::if_<boost::is_reference_wrapper<T>,
1229  mpl::identity<T>,
1230  mpl::identity<mpl::identity<T> > >::type::type::type value_type;
1231 
1232  typedef Cst<T> expression_type;
1233 
1234  constexpr explicit Cst( const T& value )
1235  :
1236  M_constant( value )
1237  {
1238  }
1239 
1240  Cst( Cst const& __cst )
1241  :
1242  M_constant( __cst.M_constant )
1243  {
1244  }
1245 
1246  Cst&
1247  operator=( Cst const& c )
1248  {
1249  if ( this != &c )
1250  M_constant = c.M_constant;
1251  return *this;
1252  }
1253 
1254  constexpr value_type value() const
1255  {
1256  return M_constant;
1257  }
1258 
1259  constexpr value_type evaluate() const
1260  {
1261  return M_constant;
1262  }
1263  constexpr value_type evaluate( bool ) const
1264  {
1265  return M_constant;
1266  }
1267 
1268  template<typename TheExpr>
1269  struct Lambda
1270  {
1271  typedef expression_type type;
1272  };
1273  template<typename TheExpr>
1274  typename Lambda<TheExpr>::type
1275  operator()( TheExpr const& e ) { return typename Lambda<TheExpr>::type(M_constant); }
1276 
1277  template<typename Geo_t, typename Basis_i_t=mpl::void_, typename Basis_j_t = Basis_i_t>
1278  struct tensor
1279  {
1280  typedef typename Cst<T>::expression_type expression_type;
1281  typedef typename Cst<T>::value_type value_type;
1282 
1283  typedef typename mpl::if_<fusion::result_of::has_key<Geo_t, vf::detail::gmc<0> >, mpl::identity<vf::detail::gmc<0> >, mpl::identity<vf::detail::gmc<1> > >::type::type key_type;
1284  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type* gmc_ptrtype;
1285  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type gmc_type;
1286  typedef Shape<gmc_type::nDim, Scalar, false, false> shape;
1287 
1288 
1289  template<typename Indq, typename Indi, typename Indj>
1290  struct expr
1291  {
1292  typedef value_type type;
1293  };
1294 
1295  struct is_zero
1296  {
1297  static const bool value = false;
1298  };
1299 
1300  tensor( expression_type const& expr,
1301  Geo_t const& /*geom*/, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
1302  :
1303  M_constant( expr.value() )
1304  {
1305  }
1306  tensor( expression_type const& expr,
1307  Geo_t const& /*geom*/, Basis_i_t const& /*fev*/ )
1308  :
1309  M_constant( expr.value() )
1310  {
1311  }
1312  tensor( expression_type const& expr, Geo_t const& /*geom*/ )
1313  :
1314  M_constant( expr.value() )
1315  {
1316  }
1317  template<typename IM>
1318  void init( IM const& /*im*/ )
1319  {
1320  }
1321  void update( Geo_t const&, Basis_i_t const& , Basis_j_t const& )
1322  {
1323  }
1324  void update( Geo_t const& , Basis_i_t const& )
1325  {
1326  }
1327  void update( Geo_t const& )
1328  {
1329  }
1330  void update( Geo_t const&, uint16_type )
1331  {
1332  }
1333 
1334  constexpr value_type
1335  evalij( uint16_type /*i*/, uint16_type /*j*/ ) const
1336  {
1337  return M_constant;
1338  }
1339 
1340 
1341  constexpr value_type
1342  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type /*c1*/, uint16_type /*c2*/, uint16_type /*q*/ ) const
1343  {
1344  return M_constant;
1345  }
1346  template<int PatternContext>
1347  constexpr value_type
1348  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type /*c1*/, uint16_type /*c2*/, uint16_type /*q*/,
1349  mpl::int_<PatternContext> ) const
1350  {
1351  return M_constant;
1352  }
1353 
1354  constexpr value_type
1355  evaliq( uint16_type /*i*/, uint16_type /*c1*/, uint16_type /*c2*/, uint16_type /*q*/ ) const
1356  {
1357  return M_constant;
1358  }
1359  constexpr value_type
1360  evalq( uint16_type /*c1*/, uint16_type /*c2*/, uint16_type /*q*/ ) const
1361  {
1362  return M_constant;
1363  }
1364  const value_type M_constant;
1365  };
1366 
1367 protected:
1368  Cst() : M_constant( 0 )
1369  {
1370  //DVLOG(2) << "Cst::Cst( default ) : constant value: " << M_constant << "\n";
1371  }
1372 
1373  const T M_constant;
1374 };
1375 
1376 template<typename T>
1377 inline
1378 Expr< Cst<T> >
1379 constant( T v )
1380 {
1381  typedef Cst<T> cst_t;
1382  return Expr< cst_t >( cst_t( v ) );
1383 }
1384 
1385 template<typename T>
1386 inline
1387 Expr< Cst<T> >
1388 cst( T v )
1389 {
1390  typedef Cst<T> cst_t;
1391  return Expr< cst_t >( cst_t( v ) );
1392 }
1393 template<typename T>
1394 inline
1395 Expr< Cst<boost::reference_wrapper<T> > >
1396 cst_ref( T& v )
1397 {
1398  typedef Cst<boost::reference_wrapper<T> > cst_t;
1399  return Expr< cst_t >( cst_t( boost::ref( v ) ) );
1400 }
1401 template<typename T>
1402 inline
1403 Expr< Cst<boost::reference_wrapper<T> > >
1404 constant_ref( T& v )
1405 {
1406  return cst_ref( v );
1407 }
1408 
1409 
1410 
1411 template<int CType>
1412 class One
1413 {
1414 public:
1415  static const size_type context = 0;
1416  static const bool is_terminal = false;
1417 
1418  static const uint16_type imorder = 0;
1419  static const bool imIsPoly = true;
1420 
1421 
1422  template<typename Func>
1423  struct HasTestFunction
1424  {
1425  static const bool result = false;
1426  };
1427  template<typename Func>
1428  struct HasTrialFunction
1429  {
1430  static const bool result = false;
1431  };
1432 
1433 
1434  typedef One<CType> this_type;
1435 
1436  typedef double value_type;
1437 
1438  One() {}
1439  One( One const& /*__vff*/ ) {}
1440 
1441  template<typename TheExpr>
1442  struct Lambda
1443  {
1444  typedef this_type type;
1445  };
1446  template<typename TheExpr>
1447  typename Lambda<TheExpr>::type
1448  operator()( TheExpr const& e ) { return this_type(); }
1449 
1450 
1451  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t = Basis_i_t>
1452  struct tensor
1453  {
1454  typedef this_type expression_type;
1455  typedef typename mpl::if_<fusion::result_of::has_key<Geo_t, vf::detail::gmc<0> >,
1456  mpl::identity<vf::detail::gmc<0> >,
1457  mpl::identity<vf::detail::gmc<1> > >::type::type key_type;
1458  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type* gmc_ptrtype;
1459  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type gmc_type;
1460  typedef Shape<gmc_type::nDim, Vectorial, false, false> shape;
1461  static const bool theshape = ( shape::M == gmc_type::nDim && shape::N == 1 );
1462  BOOST_MPL_ASSERT_MSG( theshape,
1463  INVALID_TENSOR_SHAPE_SHOULD_BE_RANK_1,
1464  ( mpl::int_<shape::M>, mpl::int_<shape::N> ) );
1465 
1466  typedef typename expression_type::value_type value_type;
1467 
1468  static const uint16_type nComponents = gmc_type::nDim;
1469  static const int16_type vector_comp = ( CType==-1 )?1:CType;
1470 
1471  typedef typename mpl::if_<mpl::equal_to<mpl::int_<CType>,mpl::int_<-1> >,
1472  mpl::identity<ublas::scalar_vector<scalar_type> >,
1473  mpl::identity<ublas::unit_vector<scalar_type> > >::type::type vector_type;
1474 
1475  struct is_zero
1476  {
1477  static const bool value = false;
1478  };
1479 
1480  tensor( expression_type const& /*expr*/,
1481  Geo_t const& /*geom*/,
1482  Basis_i_t const& /*fev*/,
1483  Basis_j_t const& /*feu*/ )
1484  :
1485  M_one( nComponents, vector_comp )
1486  {
1487  //std::cout << "one = " << M_one << "\n";
1488  }
1489  tensor( expression_type const& /*expr*/,
1490  Geo_t const& /*geom*/,
1491  Basis_i_t const& /*fev*/ )
1492  :
1493  M_one( nComponents, vector_comp )
1494  {
1495  }
1496  tensor( expression_type const& /*expr*/,
1497  Geo_t const& /*geom*/ )
1498  :
1499  M_one( nComponents, vector_comp )
1500  {
1501  // std::cout << "one = " << M_one << "\n"
1502  // << "M=" << shape::M << "\n"
1503  // << "N=" << shape::N << "\n";
1504  }
1505  template<typename IM>
1506  void init( IM const& /*im*/ )
1507  {
1508  }
1509  void update( Geo_t const& /*geom*/, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
1510  {
1511  }
1512  void update( Geo_t const& /*geom*/, Basis_i_t const& /*fev*/ )
1513  {
1514  }
1515  void update( Geo_t const& /*geom*/ )
1516  {
1517  }
1518  void update( Geo_t const& /*geom*/, uint16_type /*face*/ )
1519  {
1520  }
1521 
1522  FEELPP_STRONG_INLINE value_type
1523  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type c1, uint16_type /*c2*/, uint16_type /*q*/ ) const
1524  {
1525  return ( gmc_type::nDim>=c1 )&&( ( c1==(uint16_type)CType ) || ( CType==-1 ) );
1526  //return M_one[c1];
1527  }
1528  template<int PatternContext>
1529  FEELPP_STRONG_INLINE value_type
1530  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type c1, uint16_type /*c2*/, uint16_type /*q*/,
1531  mpl::int_<PatternContext> ) const
1532  {
1533  return ( gmc_type::nDim>=c1 )&&( ( c1==(uint16_type)CType ) || ( CType==-1 ) );
1534  //return M_one[c1];
1535  }
1536 
1537  FEELPP_STRONG_INLINE value_type
1538  evaliq( uint16_type /*i*/, uint16_type c1, uint16_type /*c2*/, uint16_type /*q*/ ) const
1539  {
1540  return ( gmc_type::nDim>=c1 )&&( ( c1==(uint16_type)CType ) || ( CType==-1 ) );
1541  //return M_one[c1];
1542  }
1543  FEELPP_STRONG_INLINE value_type
1544  evalq( uint16_type c1, uint16_type /*c2*/, uint16_type /*q*/ ) const
1545  {
1546  return ( gmc_type::nDim>=c1 )&&( ( c1==(uint16_type)CType ) || ( CType==-1 ) );
1547  //return M_one[c1];
1548  }
1549  vector_type M_one;
1550  };
1551 
1552 };
1553 
1554 inline
1555 Expr<One<-1> >
1556 one()
1557 {
1558  return Expr< One<-1> >( One<-1>() );
1559 }
1560 
1561 inline
1562 Expr<One<0> >
1563 oneX()
1564 {
1565  return Expr< One<0> >( One<0>() );
1566 }
1567 
1568 inline
1569 Expr<One<1> >
1570 oneY()
1571 {
1572  return Expr< One<1> >( One<1>() );
1573 }
1574 
1575 inline
1576 Expr<One<2> >
1577 oneZ()
1578 {
1579  return Expr< One<2> >( One<2>() );
1580 }
1581 inline
1582 Expr<One<0> >
1583 unitX()
1584 {
1585  return Expr< One<0> >( One<0>() );
1586 }
1587 
1588 inline
1589 Expr<One<1> >
1590 unitY()
1591 {
1592  return Expr< One<1> >( One<1>() );
1593 }
1594 
1595 inline
1596 Expr<One<2> >
1597 unitZ()
1598 {
1599  return Expr< One<2> >( One<2>() );
1600 }
1601 
1602 
1603 
1608 template < class ExprT >
1609 class UnaryPlus
1610 {
1611 public:
1612 
1613  static const size_type context = ExprT::context;
1614  static const bool is_terminal = false;
1615 
1616  static const uint16_type imorder = ExprT::imorder;
1617  static const bool imIsPoly = ExprT::imIsPoly;
1618 
1619  template<typename Func>
1620  struct HasTestFunction
1621  {
1622  static const bool result = ExprT::template HasTestFunction<Func>::result;
1623  };
1624 
1625  template<typename Func>
1626  struct HasTrialFunction
1627  {
1628  static const bool result = ExprT::template HasTrialFunction<Func>::result;
1629  };
1630 
1631  typedef ExprT expression_type;
1632  typedef typename ExprT::value_type value_type;
1633  typedef UnaryPlus<ExprT> this_type;
1634 
1635  UnaryPlus( const ExprT& expr )
1636  :
1637  M_expr( expr )
1638  {
1639  }
1640  UnaryPlus( UnaryPlus const& e )
1641  :
1642  M_expr( e.M_expr )
1643  {
1644  }
1645 
1646  expression_type const& expression() const
1647  {
1648  return M_expr;
1649  }
1650 
1651  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t = Basis_i_t>
1652  struct tensor
1653  {
1654  typedef typename expression_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> tensor_expr_type;
1655  typedef typename tensor_expr_type::value_type value_type;
1656  typedef typename tensor_expr_type::shape shape;
1657 
1658  struct is_zero
1659  {
1660  static const bool value = tensor_expr_type::is_zero::value;
1661  };
1662 
1663  tensor( this_type const& expr,
1664  Geo_t const& geom,
1665  Basis_i_t const& fev,
1666  Basis_j_t const& feu )
1667  :
1668  M_t_expr( expr.expression(), geom, fev, feu )
1669  {}
1670  tensor( this_type const& expr,
1671  Geo_t const& geom,
1672  Basis_i_t const& fev )
1673  :
1674  M_t_expr( expr.expression(), geom, fev )
1675  {}
1676  tensor( this_type const& expr,
1677  Geo_t const& geom )
1678  :
1679  M_t_expr( expr.expression(), geom )
1680  {}
1681  template<typename IM>
1682  void init( IM const& im )
1683  {
1684  M_t_expr.init( im );
1685  }
1686  void update( Geo_t const& geom, Basis_i_t const& fev , Basis_j_t const& feu )
1687  {
1688  M_t_expr.update( geom, fev, feu );
1689  }
1690  void update( Geo_t const& geom, Basis_i_t const& fev )
1691  {
1692  M_t_expr.update( geom, fev );
1693  }
1694  void update( Geo_t const& geom )
1695  {
1696  M_t_expr.update( geom );
1697  }
1698  void update( Geo_t const& geom, uint16_type face )
1699  {
1700  M_t_expr.update( geom, face );
1701  }
1702 
1703  value_type
1704  evalij( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2 ) const
1705  {
1706  return M_t_expr.evalij( i, j, c1, c2 );
1707  }
1708 
1709  value_type
1710  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q ) const
1711  {
1712  return M_t_expr.evalijq( i, j, c1, c2, q );
1713  }
1714  template<int PatternContext>
1715  value_type
1716  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
1717  mpl::int_<PatternContext> ) const
1718  {
1719  return M_t_expr.evalijq( i, j, c1, c2, q, mpl::int_<PatternContext>() );
1720  }
1721 
1722  value_type
1723  evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q ) const
1724  {
1725  return M_t_expr.evaliq( i, c1, c2, q );
1726  }
1727  value_type
1728  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const
1729  {
1730  return M_t_expr.evalq( c1, c2, q );
1731  }
1732  tensor_expr_type M_t_expr;
1733  };
1734 protected:
1735  UnaryPlus() {}
1736 
1737  expression_type M_expr;
1738 };
1739 template <class T> inline
1740 Expr< UnaryPlus< Expr<T> > >
1741 operator + ( const Expr<T>& expr )
1742 {
1743  typedef UnaryPlus< Expr<T> > expr_t;
1744 
1745  return Expr< expr_t >( expr_t( expr ) );
1746 }
1747 
1752 template < class ExprT >
1753 class UnaryMinus
1754 {
1755 public:
1756 
1757  static const size_type context = ExprT::context;
1758  static const bool is_terminal = false;
1759 
1760  static const uint16_type imorder = ExprT::imorder;
1761  static const bool imIsPoly = ExprT::imIsPoly;
1762 
1763  template<typename Func>
1764  struct HasTestFunction
1765  {
1766  static const bool result = ExprT::template HasTestFunction<Func>::result;
1767  };
1768 
1769  template<typename Func>
1770  struct HasTrialFunction
1771  {
1772  static const bool result = ExprT::template HasTrialFunction<Func>::result;
1773  };
1774 
1775  typedef ExprT expression_type;
1776  typedef typename ExprT::value_type value_type;
1777  typedef UnaryMinus<ExprT> this_type;
1778 
1779  UnaryMinus( const ExprT& expr )
1780  :
1781  M_expr( expr )
1782  {
1783  ;
1784  }
1785 
1786  UnaryMinus( UnaryMinus const& u )
1787  :
1788  M_expr( u.M_expr )
1789  {
1790  ;
1791  }
1792 
1793  UnaryMinus( UnaryMinus&& u )
1794  :
1795  M_expr( u.M_expr )
1796  {
1797  }
1798 
1799  UnaryMinus&
1800  operator=( UnaryMinus const& u )
1801  {
1802  if ( this != &u )
1803  M_expr = u.M_expr;
1804  return *this;
1805  }
1806 
1807  expression_type const& expression() const
1808  {
1809  return M_expr;
1810  }
1811 
1812  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t = Basis_i_t>
1813  struct tensor
1814  {
1815  typedef typename expression_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> tensor_expr_type;
1816  typedef typename tensor_expr_type::value_type value_type;
1817  typedef typename tensor_expr_type::shape shape;
1818 
1819  struct is_zero
1820  {
1821  static const bool value = tensor_expr_type::is_zero::value;
1822  };
1823 
1824  tensor( this_type const& expr,
1825  Geo_t const& geom,
1826  Basis_i_t const& fev,
1827  Basis_j_t const& feu )
1828  :
1829  M_t_expr( expr.expression(), geom, fev, feu )
1830  {}
1831  tensor( this_type const& expr,
1832  Geo_t const& geom,
1833  Basis_i_t const& fev )
1834  :
1835  M_t_expr( expr.expression(), geom, fev )
1836  {}
1837  tensor( this_type const& expr,
1838  Geo_t const& geom )
1839  :
1840  M_t_expr( expr.expression(), geom )
1841  {}
1842  template<typename IM>
1843  void init( IM const& im )
1844  {
1845  M_t_expr.init( im );
1846  }
1847  void update( Geo_t const& geom, Basis_i_t const& fev , Basis_j_t const& feu )
1848  {
1849  M_t_expr.update( geom, fev, feu );
1850  }
1851  void update( Geo_t const& geom, Basis_i_t const& fev )
1852  {
1853  M_t_expr.update( geom, fev );
1854  }
1855  void update( Geo_t const& geom )
1856  {
1857  M_t_expr.update( geom );
1858  }
1859  void update( Geo_t const& geom, uint16_type face )
1860  {
1861  M_t_expr.update( geom, face );
1862  }
1863 
1864  value_type
1865  evalij( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2 ) const
1866  {
1867  return -M_t_expr.evalij( i, j, c1, c2 );
1868  }
1869 
1870  value_type
1871  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q ) const
1872  {
1873  return -M_t_expr.evalijq( i, j, c1, c2, q );
1874  }
1875  template<int PatternContext>
1876  value_type
1877  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
1878  mpl::int_<PatternContext> ) const
1879  {
1880  return -M_t_expr.evalijq( i, j, c1, c2, q, mpl::int_<PatternContext>() );
1881  }
1882 
1883  value_type
1884  evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q ) const
1885  {
1886  return -M_t_expr.evaliq( i, c1, c2, q );
1887  }
1888  value_type
1889  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const
1890  {
1891  return -M_t_expr.evalq( c1, c2, q );
1892  }
1893  tensor_expr_type M_t_expr;
1894  };
1895 protected:
1896  UnaryMinus() {}
1897 
1898  expression_type M_expr;
1899 };
1900 template <class T> inline
1901 Expr< UnaryMinus< Expr<T> > >
1902 operator - ( const Expr<T>& expr )
1903 {
1904  typedef UnaryMinus< Expr<T> > expr_t;
1905 
1906  return Expr< expr_t >( expr_t( expr ) );
1907 }
1908 
1909 template < typename ExprT1, typename ExprT2 >
1910 class OpMax
1911 {
1912 public:
1913  static const size_type context = ExprT1::context | ExprT2::context;
1914  static const bool is_terminal = false;
1915 
1916  static const uint16_type imorder = ( ExprT1::imorder<ExprT2::imorder )*ExprT2::imorder + ( ExprT1::imorder>=ExprT2::imorder )*ExprT1::imorder;
1917  static const bool imIsPoly = ExprT1::imIsPoly && ExprT2::imIsPoly;
1918 
1919  template<typename Func>
1920  struct HasTestFunction
1921  {
1922  static const bool result = false;
1923  };
1924  template<typename Func>
1925  struct HasTrialFunction
1926  {
1927  static const bool result = false;
1928  };
1929 
1930  typedef OpMax<ExprT1, ExprT2> this_type;
1931  typedef ExprT1 expression_1_type;
1932  typedef ExprT2 expression_2_type;
1933  typedef typename strongest_numeric_type<typename expression_1_type::value_type,
1934  typename expression_2_type::value_type>::type value_type;
1935  explicit OpMax( expression_1_type const& __expr1, expression_2_type const& __expr2 )
1936  :
1937  M_expr_1( __expr1 ),
1938  M_expr_2( __expr2 )
1939  {
1940  DVLOG(2) << "OpMax::OpMax default constructor\n";
1941  }
1942 
1943  OpMax( OpMax const& __vfp )
1944  :
1945  M_expr_1( __vfp.M_expr_1 ),
1946  M_expr_2( __vfp.M_expr_2 )
1947  {
1948  DVLOG(2) << "OpMax::OpMax copy constructor\n";
1949  }
1950 
1951  expression_1_type const& left() const
1952  {
1953  return M_expr_1;
1954  }
1955  expression_2_type const& right() const
1956  {
1957  return M_expr_2;
1958  }
1959 
1960  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t = Basis_i_t>
1961  struct tensor
1962  {
1963  typedef this_type expression_type;
1964  typedef typename expression_1_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> l_type;
1965  typedef typename expression_2_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> r_type;
1966 
1967  typedef typename strongest_numeric_type<typename l_type::value_type,
1968  typename r_type::value_type>::type value_type;
1969 
1970  typedef typename mpl::if_<fusion::result_of::has_key<Geo_t, vf::detail::gmc<0> >,
1971  mpl::identity<vf::detail::gmc<0> >,
1972  mpl::identity<vf::detail::gmc<1> > >::type::type key_type;
1973  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type* gmc_ptrtype;
1974  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type gmc_type;
1975 
1976  BOOST_MPL_ASSERT_MSG( ( boost::is_same<typename l_type::shape, typename r_type::shape>::value ),
1977  INVALID_SHAPES_FOR_MIN,
1978  ( mpl::int_<l_type::shape::M>,mpl::int_<l_type::shape::N>,mpl::int_<r_type::shape::M>,mpl::int_<r_type::shape::N> ) );
1979  typedef typename l_type::shape shape;
1980 
1981  struct is_zero
1982  {
1983  static const bool value = false;
1984  };
1985 
1986  tensor( expression_type const& expr, Geo_t const& geom,
1987  Basis_i_t const& fev, Basis_j_t const& feu )
1988  :
1989  M_gmc( fusion::at_key<key_type>( geom ).get() ),
1990  M_left( expr.left(), geom, fev, feu ),
1991  M_right( expr.right(), geom, fev, feu )
1992  {
1993  }
1994  tensor( expression_type const& expr, Geo_t const& geom,
1995  Basis_i_t const& fev )
1996  :
1997  M_gmc( fusion::at_key<key_type>( geom ).get() ),
1998  M_left( expr.left(), geom, fev ),
1999  M_right( expr.right(), geom, fev )
2000  {
2001  }
2002  tensor( expression_type const& expr, Geo_t const& geom )
2003  :
2004  M_gmc( fusion::at_key<key_type>( geom ).get() ),
2005  M_left( expr.left(), geom ),
2006  M_right( expr.right(), geom )
2007  {
2008  }
2009  template<typename IM>
2010  void init( IM const& im )
2011  {
2012  M_left.init( im );
2013  M_right.init( im );
2014  }
2015  void update( Geo_t const& geom, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
2016  {
2017  update( geom );
2018  }
2019  void update( Geo_t const& geom, Basis_i_t const& /*fev*/ )
2020  {
2021  update( geom );
2022  }
2023  void update( Geo_t const& geom )
2024  {
2025  M_gmc = fusion::at_key<key_type>( geom ).get();
2026  M_left.update( geom );
2027  M_right.update( geom );
2028 
2029  }
2030  void update( Geo_t const& geom, uint16_type face )
2031  {
2032  M_gmc = fusion::at_key<key_type>( geom ).get();
2033  M_left.update( geom, face );
2034  M_right.update( geom, face );
2035 
2036  }
2037 
2038  value_type
2039  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type c1, uint16_type c2, uint16_type q ) const
2040  {
2041  return evalq( c1, c2, q );
2042  }
2043  template<int PatternContext>
2044  value_type
2045  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
2046  mpl::int_<PatternContext> ) const
2047  {
2048  Feel::detail::ignore_unused_variable_warning( i );
2049  Feel::detail::ignore_unused_variable_warning( j );
2050  return evalq( c1, c2, q );
2051  }
2052 
2053  value_type
2054  evaliq( uint16_type /*i*/, uint16_type c1, uint16_type c2, uint16_type q ) const
2055  {
2056  return evalq( c1, c2, q );
2057 
2058  }
2059 
2060  value_type
2061  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const
2062  {
2063  value_type left = M_left.evalq( c1, c2, q );
2064  value_type right = M_right.evalq( c1, c2, q );
2065  return std::max( left, right );
2066  }
2067 
2068  private:
2069  gmc_ptrtype M_gmc;
2070  l_type M_left;
2071  r_type M_right;
2072  };
2073 
2074 protected:
2075  OpMax() {}
2076 
2077  expression_1_type M_expr_1;
2078  expression_2_type M_expr_2;
2079 };
2080 
2081 template<typename ExprT1, typename ExprT2>
2082 inline
2083 Expr< OpMax<typename mpl::if_<boost::is_arithmetic<ExprT1>,
2084  mpl::identity<Cst<ExprT1> >,
2085  mpl::identity<ExprT1> >::type::type,
2086  typename mpl::if_<boost::is_arithmetic<ExprT2>,
2087  mpl::identity<Cst<ExprT2> >,
2088  mpl::identity<ExprT2> >::type::type
2089  > >
2090  max( ExprT1 const& __e1, ExprT2 const& __e2 )
2091 {
2092  typedef typename mpl::if_<boost::is_arithmetic<ExprT1>,
2093  mpl::identity<Cst<ExprT1> >,
2094  mpl::identity<ExprT1> >::type::type t1;
2095  typedef typename mpl::if_<boost::is_arithmetic<ExprT2>,
2096  mpl::identity<Cst<ExprT2> >,
2097  mpl::identity<ExprT2> >::type::type t2;
2098  typedef OpMax<t1, t2> expr_t;
2099  return Expr< expr_t >( expr_t( t1( __e1 ), t2( __e2 ) ) );
2100 }
2101 
2102 template < typename ExprT1, typename ExprT2 >
2103 class OpMin
2104 {
2105 public:
2106  static const size_type context = ExprT1::context | ExprT2::context;
2107  static const bool is_terminal = false;
2108 
2109  static const uint16_type imorder = ( ExprT1::imorder<ExprT2::imorder )*ExprT2::imorder + ( ExprT1::imorder>=ExprT2::imorder )*ExprT1::imorder;
2110  static const bool imIsPoly = ExprT1::imIsPoly && ExprT2::imIsPoly;
2111 
2112  template<typename Func>
2113  struct HasTestFunction
2114  {
2115  static const bool result = false;
2116  };
2117 
2118  template<typename Func>
2119  struct HasTrialFunction
2120  {
2121  static const bool result = false;
2122  };
2123 
2124  typedef OpMin<ExprT1, ExprT2> this_type;
2125  typedef ExprT1 expression_1_type;
2126  typedef ExprT2 expression_2_type;
2127  typedef typename strongest_numeric_type<typename expression_1_type::value_type,
2128  typename expression_2_type::value_type>::type value_type;
2129 
2130  explicit OpMin( expression_1_type const& __expr1, expression_2_type const& __expr2 )
2131  :
2132  M_expr_1( __expr1 ),
2133  M_expr_2( __expr2 )
2134  {
2135  DVLOG(2) << "OpMin::OpMin default constructor\n";
2136  }
2137 
2138  OpMin( OpMin const& __vfp )
2139  :
2140  M_expr_1( __vfp.M_expr_1 ),
2141  M_expr_2( __vfp.M_expr_2 )
2142  {
2143  DVLOG(2) << "OpMin::OpMin copy constructor\n";
2144  }
2145 
2146  expression_1_type const& left() const
2147  {
2148  return M_expr_1;
2149  }
2150  expression_2_type const& right() const
2151  {
2152  return M_expr_2;
2153  }
2154 
2155  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t = Basis_i_t>
2156  struct tensor
2157  {
2158  typedef this_type expression_type;
2159  typedef typename expression_1_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> l_type;
2160  typedef typename expression_2_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> r_type;
2161 
2162  typedef typename strongest_numeric_type<typename l_type::value_type,
2163  typename r_type::value_type>::type value_type;
2164 
2165  typedef typename mpl::if_<fusion::result_of::has_key<Geo_t, vf::detail::gmc<0> >,
2166  mpl::identity<vf::detail::gmc<0> >,
2167  mpl::identity<vf::detail::gmc<1> > >::type::type key_type;
2168  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type* gmc_ptrtype;
2169  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type gmc_type;
2170 
2171  BOOST_MPL_ASSERT_MSG( ( boost::is_same<typename l_type::shape, typename r_type::shape>::value ),
2172  INVALID_SHAPES_FOR_MIN,
2173  ( mpl::int_<l_type::shape::M>,mpl::int_<l_type::shape::N>,mpl::int_<r_type::shape::M>,mpl::int_<r_type::shape::N> ) );
2174  typedef typename l_type::shape shape;
2175 
2176  struct is_zero
2177  {
2178  static const bool value = false;
2179  };
2180 
2181  tensor( expression_type const& expr, Geo_t const& geom,
2182  Basis_i_t const& fev, Basis_j_t const& feu )
2183  :
2184  M_gmc( fusion::at_key<key_type>( geom ).get() ),
2185  M_left( expr.left(), geom, fev, feu ),
2186  M_right( expr.right(), geom, fev, feu )
2187  {
2188  }
2189  tensor( expression_type const& expr, Geo_t const& geom,
2190  Basis_i_t const& fev )
2191  :
2192  M_gmc( fusion::at_key<key_type>( geom ).get() ),
2193  M_left( expr.left(), geom, fev ),
2194  M_right( expr.right(), geom, fev )
2195  {
2196  }
2197  tensor( expression_type const& expr, Geo_t const& geom )
2198  :
2199  M_gmc( fusion::at_key<key_type>( geom ).get() ),
2200  M_left( expr.left(), geom ),
2201  M_right( expr.right(), geom )
2202  {
2203  }
2204  template<typename IM>
2205  void init( IM const& im )
2206  {
2207  M_left.init( im );
2208  M_right.init( im );
2209  }
2210  void update( Geo_t const& geom, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
2211  {
2212  update( geom );
2213  }
2214  void update( Geo_t const& geom, Basis_i_t const& /*fev*/ )
2215  {
2216  update( geom );
2217  }
2218  void update( Geo_t const& geom )
2219  {
2220  M_gmc = fusion::at_key<key_type>( geom ).get();
2221  M_left.update( geom );
2222  M_right.update( geom );
2223 
2224  }
2225  void update( Geo_t const& geom, uint16_type face )
2226  {
2227  M_gmc = fusion::at_key<key_type>( geom ).get();
2228  M_left.update( geom, face );
2229  M_right.update( geom, face );
2230 
2231  }
2232 
2233  value_type
2234  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type c1, uint16_type c2, uint16_type q ) const
2235  {
2236  return evalq( c1, c2, q );
2237  }
2238  template<int PatternContext>
2239  value_type
2240  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
2241  mpl::int_<PatternContext> ) const
2242  {
2243  return evalq( c1, c2, q );
2244  }
2245 
2246  value_type
2247  evaliq( uint16_type /*i*/, uint16_type c1, uint16_type c2, uint16_type q ) const
2248  {
2249  return evalq( c1, c2, q );
2250 
2251  }
2252 
2253  value_type
2254  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const
2255  {
2256  value_type left = M_left.evalq( c1, c2, q );
2257  value_type right = M_right.evalq( c1, c2, q );
2258  return std::min( left, right );
2259  }
2260 
2261  private:
2262  gmc_ptrtype M_gmc;
2263  l_type M_left;
2264  r_type M_right;
2265  };
2266 
2267 protected:
2268  OpMin() {}
2269 
2270  expression_1_type M_expr_1;
2271  expression_2_type M_expr_2;
2272 };
2273 
2274 template<typename ExprT1, typename ExprT2>
2275 inline
2276 Expr< OpMin<typename mpl::if_<boost::is_arithmetic<ExprT1>,
2277  mpl::identity<Cst<ExprT1> >,
2278  mpl::identity<ExprT1> >::type::type,
2279  typename mpl::if_<boost::is_arithmetic<ExprT2>,
2280  mpl::identity<Cst<ExprT2> >,
2281  mpl::identity<ExprT2> >::type::type> >
2282  min( ExprT1 const& __e1, ExprT2 const& __e2 )
2283 {
2284  typedef typename mpl::if_<boost::is_arithmetic<ExprT1>,
2285  mpl::identity<Cst<ExprT1> >,
2286  mpl::identity<ExprT1> >::type::type t1;
2287  typedef typename mpl::if_<boost::is_arithmetic<ExprT2>,
2288  mpl::identity<Cst<ExprT2> >,
2289  mpl::identity<ExprT2> >::type::type t2;
2290  typedef OpMin<t1, t2> expr_t;
2291  return Expr< expr_t >( expr_t( t1( __e1 ), t2( __e2 ) ) );
2292 }
2293 
2294 template < typename ExprT1, typename ExprT2 >
2295 class Pow
2296 {
2297 public:
2298 
2299  static const size_type context = ExprT1::context|ExprT2::context;
2300  static const bool is_terminal = false;
2301 
2306  static const uint16_type imorder = ExprT1::imorder;
2307  static const bool imIsPoly = ExprT1::imIsPoly && ExprT2::imIsPoly;
2308 
2309  template<typename Func>
2310  struct HasTestFunction
2311  {
2312  static const bool result = false;
2313  };
2314 
2315  template<typename Func>
2316  struct HasTrialFunction
2317  {
2318  static const bool result = false;
2319  };
2320 
2321  typedef Pow<ExprT1, ExprT2> this_type;
2322  typedef ExprT1 expression_1_type;
2323  typedef ExprT2 expression_2_type;
2324  typedef typename expression_1_type::value_type value_1_type;
2325  typedef typename expression_2_type::value_type value_2_type;
2326  typedef value_1_type value_type;
2327 
2328  // verify that all returning types are integral or floating types
2329  BOOST_STATIC_ASSERT( ::boost::is_arithmetic<value_1_type>::value &&
2330  ::boost::is_arithmetic<value_2_type>::value );
2331 
2332  explicit Pow( expression_1_type const& __expr1, expression_2_type const& __expr2 )
2333  :
2334  M_expr_1( __expr1 ),
2335  M_expr_2( __expr2 )
2336  {
2337  DVLOG(2) << "Pow::Pow default constructor\n";
2338  }
2339 
2340  Pow( Pow const& __vfp )
2341  :
2342  M_expr_1( __vfp.M_expr_1 ),
2343  M_expr_2( __vfp.M_expr_2 )
2344  {
2345  DVLOG(2) << "Pow::Pow copy constructor\n";
2346  }
2347 
2348  bool isSymetric() const
2349  {
2350  return false;
2351  }
2352 
2353  expression_1_type const& left() const
2354  {
2355  return M_expr_1;
2356  }
2357  expression_2_type const& right() const
2358  {
2359  return M_expr_2;
2360  }
2361 
2362  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t = Basis_i_t>
2363  struct tensor
2364  {
2365  typedef this_type expression_type;
2366  typedef typename expression_1_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> l_type;
2367  typedef typename expression_2_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> r_type;
2368 
2369  typedef typename strongest_numeric_type<typename l_type::value_type,
2370  typename r_type::value_type>::type value_type;
2371 
2372 
2373  typedef typename mpl::if_<fusion::result_of::has_key<Geo_t, vf::detail::gmc<0> >,
2374  mpl::identity<vf::detail::gmc<0> >,
2375  mpl::identity<vf::detail::gmc<1> > >::type::type key_type;
2376  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type* gmc_ptrtype;
2377  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type gmc_type;
2378  typedef typename l_type::shape shape;
2379 
2380  typedef typename Eigen::Matrix<value_type,shape::M,shape::N> loc_type;
2381 
2382  struct is_zero
2383  {
2384  static const bool value = l_type::is_zero::value;
2385  };
2386 
2387  tensor( expression_type const& expr, Geo_t const& geom,
2388  Basis_i_t const& fev, Basis_j_t const& feu )
2389  :
2390  M_gmc( fusion::at_key<key_type>( geom ).get() ),
2391  M_left( expr.left(), geom, fev, feu ),
2392  M_right( expr.right(), geom, fev, feu ),
2393  M_loc( boost::extents[M_gmc->nPoints()] )
2394  {
2395  update( geom );
2396  }
2397  tensor( expression_type const& expr, Geo_t const& geom,
2398  Basis_i_t const& fev )
2399  :
2400  M_gmc( fusion::at_key<key_type>( geom ).get() ),
2401  M_left( expr.left(), geom, fev ),
2402  M_right( expr.right(), geom, fev ),
2403  M_loc( boost::extents[M_gmc->nPoints()] )
2404  {
2405  update( geom );
2406  }
2407  tensor( expression_type const& expr, Geo_t const& geom )
2408  :
2409  M_gmc( fusion::at_key<key_type>( geom ).get() ),
2410  M_left( expr.left(), geom ),
2411  M_right( expr.right(), geom ),
2412  M_loc( boost::extents[M_gmc->nPoints()] )
2413  {
2414  update( geom );
2415  }
2416  template<typename IM>
2417  void init( IM const& im )
2418  {
2419  M_left.init( im );
2420  M_right.init( im );
2421  }
2422  void update( Geo_t const& geom, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
2423  {
2424  update( geom );
2425  }
2426  void update( Geo_t const& geom, Basis_i_t const& /*fev*/ )
2427  {
2428  update( geom );
2429  }
2430  void update( Geo_t const& geom )
2431  {
2432  const int npts = fusion::at_key<key_type>( geom ).get()->nPoints();
2433  M_left.update( geom );
2434  M_right.update( geom );
2435 
2436  for ( int q = 0; q < npts; ++q )
2437  for ( int c1 = 0; c1 < shape::M; ++c1 )
2438  for ( int c2 = 0; c2 < shape::N; ++c2 )
2439  {
2440  value_type left = M_left.evalq( c1, c2, q );
2441  value_type right = M_right.evalq( c1, c2, q );
2442  M_loc[q]( c1,c2 ) = std::pow( left, right );
2443  }
2444  }
2445  void update( Geo_t const& geom, uint16_type face )
2446  {
2447  M_gmc = fusion::at_key<key_type>( geom ).get();
2448  M_left.update( geom, face );
2449  M_right.update( geom, face );
2450 
2451  for ( int q = 0; q < M_gmc->nPoints(); ++q )
2452  for ( int c1 = 0; c1 < shape::M; ++c1 )
2453  for ( int c2 = 0; c2 < shape::N; ++c2 )
2454  {
2455  value_type left = M_left.evalq( c1, c2, q );
2456  value_type right = M_right.evalq( c1, c2, q );
2457  M_loc[q]( c1,c2 ) = std::pow( left, right );
2458  }
2459  }
2460 
2461  value_type
2462  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type c1, uint16_type c2, uint16_type q ) const
2463  {
2464  return evalq( c1, c2, q );
2465  }
2466  template<int PatternContext>
2467  value_type
2468  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
2469  mpl::int_<PatternContext> ) const
2470  {
2471  Feel::detail::ignore_unused_variable_warning( i );
2472  Feel::detail::ignore_unused_variable_warning( j );
2473  return evalq( c1, c2, q );
2474  }
2475 
2476  value_type
2477  evaliq( uint16_type /*i*/, uint16_type c1, uint16_type c2, uint16_type q ) const
2478  {
2479  return evalq( c1, c2, q );
2480  }
2481 
2482  value_type
2483  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const
2484  {
2485  return evalq( c1, c2, q, mpl::int_<shape::rank>() );
2486  }
2487  value_type
2488  evalq( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<0> ) const
2489  {
2490  Feel::detail::ignore_unused_variable_warning( c1 );
2491  Feel::detail::ignore_unused_variable_warning( c2 );
2492  return M_loc[q]( 0,0 );
2493  }
2494  value_type
2495  evalq( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<1> ) const
2496  {
2497  if ( shape::M > shape::N )
2498  return M_loc[q]( c1,0 );
2499 
2500  return M_loc[q]( 0,c2 );
2501  }
2502  value_type
2503  evalq( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<2> ) const
2504  {
2505  return M_loc[q]( c1,c2 );
2506  }
2507 
2508  private:
2509  gmc_ptrtype M_gmc;
2510  l_type M_left;
2511  r_type M_right;
2512  //ublas::vector<double> M_loc;
2513  boost::multi_array<loc_type,1> M_loc;
2514  };
2515 
2516 protected:
2517  Pow() {}
2518 
2519  expression_1_type M_expr_1;
2520  expression_2_type M_expr_2;
2521 };
2522 
2523 template<typename ExprT1, typename ExprT2>
2524 inline
2525 Expr< Pow<typename mpl::if_<boost::is_arithmetic<ExprT1>,
2526  mpl::identity<Cst<ExprT1> >,
2527  mpl::identity<ExprT1> >::type::type,
2528  typename mpl::if_<boost::is_arithmetic<ExprT2>,
2529  mpl::identity<Cst<ExprT2> >,
2530  mpl::identity<ExprT2> >::type::type> >
2531  pow( ExprT1 const& __e1, ExprT2 const& __e2 )
2532 {
2533  typedef typename mpl::if_<boost::is_arithmetic<ExprT1>,
2534  mpl::identity<Cst<ExprT1> >,
2535  mpl::identity<ExprT1> >::type::type t1;
2536  typedef typename mpl::if_<boost::is_arithmetic<ExprT2>,
2537  mpl::identity<Cst<ExprT2> >,
2538  mpl::identity<ExprT2> >::type::type t2;
2539  typedef Pow<t1, t2> expr_t;
2540  return Expr< expr_t >( expr_t( t1( __e1 ), t2( __e2 ) ) );
2541 }
2542 
2543 
2544 template<typename ExprT1, typename ExprT2>
2545 inline
2546 Expr< Pow<typename mpl::if_<boost::is_arithmetic<ExprT1>,
2547  mpl::identity<Cst<ExprT1> >,
2548  mpl::identity<Expr<ExprT1> > >::type::type,
2549  typename mpl::if_<boost::is_arithmetic<ExprT2>,
2550  mpl::identity<Cst<ExprT2> >,
2551  mpl::identity<Expr<ExprT2> > >::type::type> >
2552 operator^( typename mpl::if_<boost::is_arithmetic<ExprT1>,
2553  mpl::identity<ExprT1>,
2554  mpl::identity<Expr<ExprT1> > >::type::type const& __e1,
2555  typename mpl::if_<boost::is_arithmetic<ExprT2>,
2556  mpl::identity<ExprT2>,
2557  mpl::identity<Expr<ExprT2> > >::type::type const& __e2 )
2558 {
2559  typedef typename mpl::if_<boost::is_arithmetic<ExprT1>,
2560  mpl::identity<Cst<ExprT1> >,
2561  mpl::identity<ExprT1> >::type::type t1;
2562  typedef typename mpl::if_<boost::is_arithmetic<ExprT2>,
2563  mpl::identity<Cst<ExprT2> >,
2564  mpl::identity<ExprT2> >::type::type t2;
2565  typedef Pow<t1, t2> expr_t;
2566  return Expr< expr_t >( expr_t( t1( __e1 ), t2( __e2 ) ) );
2567 }
2568 
2576 template<int GeoId, typename ExprT>
2577 class EvalFace
2578 {
2579 public:
2580 
2581  static const size_type context = ExprT::context;
2582  static const bool is_terminal = false;
2583 
2584  static const uint16_type imorder = ExprT::imorder;
2585  static const bool imIsPoly = ExprT::imIsPoly;
2589 
2590  typedef ExprT expression_type;
2591  typedef typename expression_type::value_type value_type;
2592  typedef EvalFace<GeoId,ExprT> this_type;
2593 
2595 
2599 
2600  explicit EvalFace( expression_type const & __expr )
2601  :
2602  M_expr( __expr )
2603  {}
2604  ~EvalFace()
2605  {}
2606 
2608 
2612 
2613  template<typename VecGeo_t, typename Basis_i_t = boost::none_t, typename Basis_j_t = Basis_i_t>
2614  struct tensor
2615  {
2616  typedef typename fusion::result_of::at_c<VecGeo_t,GeoId>::type Geo_t;
2617 
2618  typedef typename expression_type::template tensor<Geo_t,
2619  Basis_i_t,
2620  Basis_j_t> tensor_expr_type;
2621  typedef typename tensor_expr_type::value_type value_type;
2622 
2623  template <class Args> struct sig
2624  {
2625  typedef value_type type;
2626  };
2627 
2628  tensor( this_type const& expr,
2629  VecGeo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
2630  :
2631  M_tensor_expr( expr.expression(), fusion::at_c<GeoId>( geom ), fev, feu )
2632  {}
2633 
2634  tensor( this_type const& expr,
2635  VecGeo_t const& geom, Basis_i_t const& fev )
2636  :
2637  M_tensor_expr( expr.expression(), fusion::at_c<GeoId>( geom ), fev )
2638  {}
2639 
2640  tensor( this_type const& expr,
2641  VecGeo_t const& geom )
2642  :
2643  M_tensor_expr( expr.expression(), fusion::at_c<GeoId>( geom ) )
2644  {}
2645  template<typename IM>
2646  void init( IM const& im )
2647  {
2648  M_tensor_expr.init( im );
2649  }
2650  void update( VecGeo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
2651  {
2652  M_tensor_expr.update( fusion::at_c<GeoId>( geom ), fev, feu );
2653  }
2654  void update( VecGeo_t const& geom, Basis_i_t const& fev )
2655  {
2656  M_tensor_expr.update( fusion::at_c<GeoId>( geom ), fev );
2657  }
2658  void update( VecGeo_t const& geom )
2659  {
2660  M_tensor_expr.update( fusion::at_c<GeoId>( geom ) );
2661  }
2662 
2663 
2664  value_type
2665  evalij( uint16_type i, uint16_type j ) const
2666  {
2667  return M_tensor_expr.evalij( i, j );
2668  }
2669 
2670 
2671  value_type
2672  evalijq( uint16_type i, uint16_type j, int q ) const
2673  {
2674  return M_tensor_expr.evalijq( i, j, q );
2675  }
2676 
2677 
2678  value_type
2679  evaliq( uint16_type i, int q ) const
2680  {
2681  return M_tensor_expr.evaliq( i, q );
2682  }
2683 
2684  value_type
2685  evalq( int q, int c ) const
2686  {
2687  return M_tensor_expr.evalq( q, c );
2688  }
2689 
2690  tensor_expr_type M_tensor_expr;
2691  };
2692 
2694 
2698 
2699  expression_type const& expression() const
2700  {
2701  return M_expr;
2702  }
2703 
2705 
2709 
2711 
2715 
2716 
2717 protected:
2718 
2719 private:
2720 
2721  mutable expression_type M_expr;
2722 };
2723 
2724 template<int GeoId, typename ExprT>
2725 inline
2726 Expr< EvalFace<GeoId, ExprT> >
2727 evalface( ExprT const& v )
2728 {
2729  typedef EvalFace<GeoId, ExprT> eval_t;
2730  return Expr< eval_t >( eval_t( v ) );
2731 }
2732 
2733 template < class Element, int Type>
2734 class GElem
2735 {
2736 public:
2737 
2738  static const size_type context = vm::JACOBIAN |vm::POINT;
2739  static const bool is_terminal = false;
2740 
2741  typedef Element element_type;
2742  typedef boost::shared_ptr<element_type> element_ptrtype;
2743  typedef GElem<element_type, Type> this_type;
2744  typedef this_type self_type;
2745 
2746  typedef typename element_type::functionspace_type functionspace_type;
2747  typedef typename functionspace_type::reference_element_type* fe_ptrtype;
2748  typedef typename functionspace_type::reference_element_type fe_type;
2749  typedef typename functionspace_type::geoelement_type geoelement_type;
2750  typedef typename functionspace_type::gm_type gm_type;
2751  typedef typename functionspace_type::value_type value_type;
2752  static const uint16_type rank = fe_type::rank;
2753  static const uint16_type nComponents1 = fe_type::nComponents1;
2754  static const uint16_type nComponents2 = fe_type::nComponents2;
2755  typedef std::map<size_type,std::vector<element_ptrtype> > basis_type;
2756 
2757  static const uint16_type imorder = element_type::functionspace_type::basis_type::nOrder;
2758  static const bool imIsPoly = true;
2759 
2760  template<typename Func>
2761  struct HasTestFunction
2762  {
2763  static const bool result = ( Type==0 );
2764  };
2765 
2766  template<typename Func>
2767  struct HasTrialFunction
2768  {
2769  static const bool result = ( Type==1 );
2770  };
2771 
2772 
2773  GElem ( std::map<size_type,std::vector<element_ptrtype> > const& v )
2774  :
2775  M_basis ( v )
2776  {
2777  typename basis_type::iterator it = M_basis.begin();
2778  typename basis_type::iterator en = M_basis.end();
2779 
2780  for ( ; it != en; ++it )
2781  for ( uint16_type i = 0; i < it->second.size(); ++i )
2782  it->second[i]->updateGlobalValues();
2783  }
2784  GElem( GElem const& op )
2785  :
2786  M_basis ( op.M_basis )
2787  {
2788 
2789  }
2790 
2791  basis_type const& basis() const
2792  {
2793  return M_basis;
2794  }
2795 
2796 
2797  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t = Basis_i_t>
2798  struct tensor
2799  {
2800  typedef this_type expression_type;
2801  typedef typename mpl::if_<fusion::result_of::has_key<Geo_t, vf::detail::gmc<0> >,
2802  mpl::identity<vf::detail::gmc<0> >,
2803  mpl::identity<vf::detail::gmc<1> > >::type::type key_type;
2804  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type* gmc_ptrtype;
2805  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type gmc_type;
2806 
2807  typedef typename mpl::if_<mpl::equal_to<mpl::int_<rank>,
2808  mpl::int_<0> >,
2809  mpl::identity<Shape<gmc_type::NDim, Scalar, false> >,
2810  typename mpl::if_<mpl::equal_to<mpl::int_<rank>,
2811  mpl::int_<1> >,
2812  mpl::identity<Shape<gmc_type::NDim, Vectorial, false> >,
2813  mpl::identity<Shape<gmc_type::NDim, Tensor2, false> > >::type>::type::type shape;
2814  typedef typename fe_type::PreCompute pc_type;
2815  typedef boost::shared_ptr<pc_type> pc_ptrtype;
2816  typedef typename fe_type::template Context<context, fe_type, gm_type,geoelement_type,gmc_type::context> ctx_type;
2817  typedef boost::shared_ptr<ctx_type> ctx_ptrtype;
2818 
2819  typedef typename expression_type::value_type value_type;
2820 
2821  struct is_zero
2822  {
2823  static const bool value = false;
2824  };
2825 
2826  tensor( expression_type const& expr,
2827  Geo_t const& geom,
2828  Basis_i_t const& /*fev*/,
2829  Basis_j_t const& /*feu*/ )
2830  :
2831  M_expr( expr ),
2832  M_pc( expr.basis().begin()->second[0]->functionSpace()->fe(), fusion::at_key<key_type>( geom )->xRefs() ),
2833  M_ctx( new ctx_type( expr.basis().begin()->second[0]->functionSpace()->fe(),
2834  fusion::at_key<key_type>( geom ), ( pc_ptrtype const& )M_pc ) ),
2835  M_loc( expr.basis().begin()->second.size() )
2836  {
2837  for ( uint16_type i = 0; i < M_loc.size(); ++i )
2838  M_loc[i].resize( boost::extents[M_pc.nPoints()][nComponents1][nComponents2] );
2839  }
2840  tensor( expression_type const& expr,
2841  Geo_t const& geom,
2842  Basis_i_t const& /*fev*/ )
2843  :
2844  M_expr( expr ),
2845  M_pc( expr.basis().begin()->second[0]->functionSpace()->fe(), fusion::at_key<key_type>( geom )->xRefs() ),
2846  M_ctx( new ctx_type( expr.basis().begin()->second[0]->functionSpace()->fe(),
2847  fusion::at_key<key_type>( geom ), ( pc_ptrtype const& )M_pc ) ),
2848  M_loc( expr.basis().begin()->second.size() )
2849  {
2850  for ( uint16_type i = 0; i < M_loc.size(); ++i )
2851  M_loc[i].resize( boost::extents[M_pc.nPoints()] );
2852  }
2853  tensor( expression_type const& expr,
2854  Geo_t const& geom )
2855  :
2856  M_expr( expr ),
2857  M_pc( expr.basis().begin()->second[0]->functionSpace()->fe(), fusion::at_key<key_type>( geom )->xRefs() ),
2858  M_ctx( new ctx_type( expr.basis().begin()->second[0]->functionSpace()->fe(),
2859  fusion::at_key<key_type>( geom ), ( pc_ptrtype const& )M_pc ) ),
2860  M_loc( expr.basis().begin()->second.size() )
2861  {
2862  for ( uint16_type i = 0; i < M_loc.size(); ++i )
2863  M_loc[i].resize( boost::extents[M_pc.nPoints()] );
2864  }
2865  template<typename IM>
2866  void init( IM const& /*im*/ )
2867  {
2868 
2869  }
2870  void update( Geo_t const& geom, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
2871  {
2872  update( geom );
2873  }
2874  void update( Geo_t const& geom, Basis_i_t const& /*fev*/ )
2875  {
2876  update( geom );
2877  }
2878  void update( Geo_t const& geom )
2879  {
2880  //VLOG(1) << "[GElem] updating element " << fusion::at_key<key_type>( geom )->id() << "\n";
2881  typename basis_type::iterator it = const_cast<basis_type&>( M_expr.basis() ).find( fusion::at_key<key_type>( geom )->id() );
2882  typename basis_type::iterator en = const_cast<basis_type&>( M_expr.basis() ).end();
2883 
2884  FEELPP_ASSERT( it != en )( fusion::at_key<key_type>( geom )->id() ).error ( "invalid basis function to integrate" );
2885 
2886  for ( uint16_type i = 0; i < M_loc.size(); ++i )
2887  {
2888  //M_loc[i] = it->second[i]->id( *M_ctx, M_pc, M_loc[i] );
2889  }
2890  }
2891 
2892  value_type
2893  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q ) const
2894  {
2895  if ( Type == 0 )
2896  return M_loc[i]( c1,c2,q );
2897 
2898  return M_loc[j]( c1,c2,q );
2899  }
2900  template<int PatternContext>
2901  value_type
2902  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<PatternContext> ) const
2903  {
2904  if ( Type == 0 )
2905  return M_loc[i]( c1,c2,q );
2906 
2907  return M_loc[j]( c1,c2,q );
2908  }
2909 
2910  value_type
2911  evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q ) const
2912  {
2913  return M_loc[i]( c1,c2,q );
2914  }
2915  value_type
2916  evalq( uint16_type /*c1*/, uint16_type /*c2*/, uint16_type /*q*/ ) const
2917  {
2918 
2919  }
2920  private:
2921  this_type const& M_expr;
2922  pc_type M_pc;
2923  ctx_ptrtype M_ctx;
2924  std::vector<typename element_type::id_type> M_loc;
2925  };
2926 private:
2927  basis_type M_basis;
2928 
2929 };
2930 
2931 template<typename Elem>
2932 inline
2933 Expr< GElem<Elem,1> >
2934 basist( std::map<size_type,std::vector<boost::shared_ptr<Elem> > > const& v )
2935 {
2936  typedef GElem<Elem,1> expr_t;
2937  return Expr< expr_t >( expr_t( v ) );
2938 }
2939 template<typename Elem>
2940 inline
2941 Expr< GElem<Elem,0> >
2942 basis( std::map<size_type,std::vector<boost::shared_ptr<Elem> > > const& v )
2943 {
2944  typedef GElem<Elem,0> expr_t;
2945  return Expr< expr_t >( expr_t( v ) );
2946 }
2947 
2949 } // vf
2950 } // feel
2951 #endif /* __Expr_H */

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