34 #include <boost/version.hpp>
35 #if (BOOST_VERSION >= 103400)
36 #include <boost/none.hpp>
38 #include <boost/none_t.hpp>
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>
54 #include <feel/feelpoly/policy.hpp>
55 #include <feel/feelpoly/context.hpp>
64 typedef node<double>::type node_type;
72 template<
typename ExprT>
77 static const size_type context = ExprT::context;
78 static const bool is_terminal =
false;
80 static const uint16_type imorder = ExprT::imorder;
82 static const bool imIsPoly = ExprT::imIsPoly;
85 template<
typename Func>
86 struct HasTestFunction
88 static const bool result = ExprT::template HasTestFunction<Func>::result;
91 template<
typename Func>
92 struct HasTrialFunction
94 static const bool result = ExprT::template HasTrialFunction<Func>::result;
104 typedef ExprT expression_type;
105 typedef typename expression_type::value_type value_type;
106 typedef ComponentsExpr<ExprT> this_type;
119 explicit ComponentsExpr( expression_type
const & __expr,
int c1,
int c2 )
130 expression_type
const& expression()
const
139 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t = Basis_i_t>
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;
151 template <
class Args>
struct sig
153 typedef value_type type;
158 static const bool value = tensor_expr_type::is_zero::value;
161 tensor( this_type
const& expr,
162 Geo_t
const& geom, Basis_i_t
const& fev, Basis_j_t
const& feu )
164 M_tensor_expr( expr.expression(), geom, fev, feu ),
169 tensor( this_type
const& expr,
170 Geo_t
const& geom, Basis_i_t
const& fev )
172 M_tensor_expr( expr.expression(), geom, fev ),
177 tensor( this_type
const& expr, Geo_t
const& geom )
179 M_tensor_expr( expr.expression(), geom ),
185 template<
typename IM>
186 void init( IM
const& im )
188 M_tensor_expr.init( im );
190 void update( Geo_t
const& geom, Basis_i_t
const& fev, Basis_j_t
const& feu )
192 M_tensor_expr.update( geom, fev, feu );
194 void update( Geo_t
const& geom, Basis_i_t
const& fev )
196 M_tensor_expr.update( geom, fev );
198 void update( Geo_t
const& geom )
200 M_tensor_expr.update( geom );
202 void update( Geo_t
const& geom, uint16_type face )
204 M_tensor_expr.update( geom, face );
209 evalij( uint16_type i, uint16_type j )
const
211 return M_tensor_expr.evalij( i, j );
216 evalijq( uint16_type i, uint16_type j, uint16_type , uint16_type , uint16_type q )
const
218 return M_tensor_expr.evalijq( i, j, M_c1, M_c2, q );
221 template<
int PatternContext>
223 evalijq( uint16_type i, uint16_type j, uint16_type , uint16_type , uint16_type q,
224 mpl::int_<PatternContext> )
const
226 return M_tensor_expr.evalijq( i, j, M_c1, M_c2, q, mpl::int_<PatternContext>() );
231 evaliq( uint16_type i, uint16_type , uint16_type , uint16_type q )
const
233 return M_tensor_expr.evaliq( i, M_c1, M_c2, q );
237 evalq( uint16_type , uint16_type , uint16_type q )
const
239 return M_tensor_expr.evalq( M_c1, M_c2, q );
242 tensor_expr_type M_tensor_expr;
243 const int M_c1, M_c2;
245 expression_type M_expr;
249 class IntegratorBase {};
250 class LambdaExprBase {};
251 class LambdaExpr1 :
public LambdaExprBase
256 static const bool is_terminal =
false;
258 static const uint16_type imorder = 0;
259 static const bool imIsPoly =
false;
261 template<
typename Func>
262 struct HasTestFunction
264 static const bool result =
false;
266 template<
typename Func>
267 struct HasTrialFunction
269 static const bool result =
false;
272 typedef double value_type;
274 template<
typename TheExpr>
277 typedef typename TheExpr::expression_type type;
280 template<
typename ExprT>
281 typename Lambda<ExprT>::type
282 operator()( ExprT
const& e ) {
283 return e.expression();
286 template<
typename ExprT>
287 typename Lambda<ExprT>::type
288 operator()( ExprT
const& e )
const {
return e.expression(); }
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>
293 typedef LambdaExpr1 expression_type;
294 typedef typename LambdaExpr1::value_type value_type;
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;
302 template<
typename Indq,
typename Indi,
typename Indj>
305 typedef value_type type;
310 static const bool value =
false;
313 tensor( expression_type
const& expr,
314 Geo_t
const& , Basis_i_t
const& , Basis_j_t
const& )
317 tensor( expression_type
const& expr,
318 Geo_t
const& , Basis_i_t
const& )
321 tensor( expression_type
const& expr, Geo_t
const& )
324 template<
typename IM>
325 void init( IM
const& )
328 void update( Geo_t
const&, Basis_i_t
const& , Basis_j_t
const& )
331 void update( Geo_t
const& , Basis_i_t
const& )
334 void update( Geo_t
const& )
337 void update( Geo_t
const&, uint16_type )
342 evalij( uint16_type , uint16_type )
const
349 evalijq( uint16_type , uint16_type , uint16_type , uint16_type , uint16_type )
const
353 template<
int PatternContext>
355 evalijq( uint16_type , uint16_type , uint16_type , uint16_type , uint16_type ,
356 mpl::int_<PatternContext> )
const
362 evaliq( uint16_type , uint16_type , uint16_type , uint16_type )
const
367 evalq( uint16_type , uint16_type , uint16_type )
const
384 template<
typename ExprT>
389 static const size_type context = ExprT::context;
390 static const bool is_terminal = ExprT::is_terminal;
393 static const uint16_type imorder = ExprT::imorder;
395 static const bool imIsPoly = ExprT::imIsPoly;
398 template<
typename Func>
399 struct HasTestFunction
401 static const bool result = ExprT::template HasTestFunction<Func>::result;
404 template<
typename Func>
405 struct HasTrialFunction
407 static const bool result = ExprT::template HasTrialFunction<Func>::result;
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;
432 explicit Expr( expression_type
const & __expr )
444 Expr<ComponentsExpr<Expr<ExprT> > >
445 operator()(
int c1 = 0,
int c2 = 0 )
const
447 auto ex = ComponentsExpr<Expr<ExprT> >( Expr<ExprT>( M_expr ), c1, c2 );
448 return Expr<ComponentsExpr<Expr<ExprT> > >( ex );
451 template<
typename TheExpr>
454 typedef typename ExprT::template Lambda<TheExpr>::type expr_type;
455 typedef Expr<expr_type> type;
461 template<
typename TheExpr>
462 typename Lambda<TheExpr>::type
463 operator()( TheExpr
const& e )
468 return expr( M_expr( e ) );
471 template<
typename TheExpr>
472 typename Lambda<TheExpr>::type
473 operator()( TheExpr
const& e )
const {
return expr(M_expr(e)); }
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>
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;
483 typedef typename tensor_expr_type::shape shape;
485 template <
class Args>
struct sig
487 typedef value_type type;
492 static const bool value = tensor_expr_type::is_zero::value;
495 tensor( this_type
const& expr,
496 Geo_t
const& geom, Basis_i_t
const& fev, Basis_j_t
const& feu )
498 M_tensor_expr( expr.expression(), geom, fev, feu )
501 tensor( this_type
const& expr,
502 Geo_t
const& geom, Basis_i_t
const& fev )
504 M_tensor_expr( expr.expression(), geom, fev )
507 tensor( this_type
const& expr, Geo_t
const& geom )
509 M_tensor_expr( expr.expression(), geom )
513 template<
typename IM>
514 void init( IM
const& im )
516 M_tensor_expr.init( im );
518 void update( Geo_t
const& geom, Basis_i_t
const& fev, Basis_j_t
const& feu )
520 M_tensor_expr.update( geom, fev, feu );
522 void update( Geo_t
const& geom, Basis_i_t
const& fev )
524 M_tensor_expr.update( geom, fev );
526 void update( Geo_t
const& geom )
528 M_tensor_expr.update( geom );
530 void update( Geo_t
const& geom, uint16_type face )
532 M_tensor_expr.update( geom, face );
537 evalij( uint16_type i, uint16_type j )
const
539 return M_tensor_expr.evalij( i, j );
542 Eigen::Matrix<value_type, shape::M, shape::N>
const&
543 evalijq( uint16_type i, uint16_type j, uint16_type q )
const
545 return M_tensor_expr.evalijq( i, j, q );
549 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q )
const
551 return M_tensor_expr.evalijq( i, j, c1, c2, q );
554 template<
int PatternContext>
556 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
557 mpl::int_<PatternContext> )
const
559 return M_tensor_expr.evalijq( i, j, c1, c2, q, mpl::int_<PatternContext>() );
564 evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q )
const
566 return M_tensor_expr.evaliq( i, c1, c2, q );
570 evalq( uint16_type c1, uint16_type c2, uint16_type q )
const
572 return M_tensor_expr.evalq( c1, c2, q );
575 tensor_expr_type M_tensor_expr;
582 value_type value()
const
584 return __expression.value();
587 value_type grad(
int __ith )
const
589 return __expression.grad( __ith );
592 value_type
hessian(
int __i,
int __j )
const
594 return __expression.hessian( __i, __j );
605 bool isSymetric()
const
607 return M_expr.isSymetric();
612 expression_type
const& expression()
const
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
634 DVLOG(2) <<
"calling assemble(u,v)\n";
635 M_expr.assemble( __u, __v, __f );
636 DVLOG(2) <<
"calling assemble(u,v) done\n";
639 template<
typename Elem1,
typename FormType>
640 void assemble( boost::shared_ptr<Elem1>
const& __v,
641 FormType& __f )
const
643 DVLOG(2) <<
"calling assemble(v)\n";
644 M_expr.assemble( __v, __f );
645 DVLOG(2) <<
"calling assemble(v) done\n";
648 template<
typename P0hType>
649 typename P0hType::element_type
650 broken( boost::shared_ptr<P0hType>& P0h )
const
652 return M_expr.broken( P0h );
657 typename expression_type::value_type
658 evaluate(
bool parallel =
true, WorldComm
const& worldcomm = Environment::worldComm() )
const
660 return M_expr.evaluate( parallel,worldcomm );
663 typename expression_type::value_type
664 evaluateAndSum()
const
666 return M_expr.evaluateAndSum();
668 std::string expressionStr()
const
670 return std::string();
681 mutable expression_type M_expr;
684 template <
typename ExprT>
686 expr( ExprT
const& exprt )
688 return Expr<ExprT>( exprt );
691 template <
typename ExprT>
692 boost::shared_ptr<Expr<ExprT> >
693 exprPtr( ExprT
const& exprt )
695 return boost::shared_ptr<Expr<ExprT> >(
new Expr<ExprT>( exprt ) );
698 extern Expr<LambdaExpr1> _e1;
711 template<
typename IntElts,
typename ExprT>
712 struct ExpressionOrder
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;
720 static const uint16_type nOrderGeo = the_element_type::nOrder;
722 static const bool is_polynomial = ExprT::imIsPoly;
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;
732 static const int value = ( ExprT::imorder )?( ExprT::imorder*nOrderGeo ):( nOrderGeo );
733 static const int value_1 = ExprT::imorder;
740 template<
typename Pr
intExprT>
745 static const size_type context = PrintExprT::context;
746 static const bool is_terminal =
false;
748 static const uint16_type imorder = PrintExprT::imorder;
749 static const bool imIsPoly = PrintExprT::imIsPoly;
751 template<
typename Func>
752 struct HasTestFunction
754 static const bool result = PrintExprT::template HasTestFunction<Func>::result;
756 template<
typename Func>
757 struct HasTrialFunction
759 static const bool result = PrintExprT::template HasTrialFunction<Func>::result;
766 typedef PrintExprT expression_type;
767 typedef typename expression_type::value_type value_type;
768 typedef PrintExpr<PrintExprT> this_type;
776 explicit PrintExpr( expression_type
const & __expr,
777 std::string
const & __tag )
791 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t>
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;
797 typedef typename tensor_expr_type::shape shape;
799 template <
class Args>
struct sig
801 typedef value_type type;
805 static const bool value = tensor_expr_type::is_zero::value;
808 tensor( this_type
const& expr,
809 Geo_t
const& geom, Basis_i_t
const& fev, Basis_j_t
const& feu )
811 M_tensor_expr( expr.expression(), geom, fev, feu ),
815 tensor( this_type
const& expr,
816 Geo_t
const& geom, Basis_i_t
const& fev )
818 M_tensor_expr( expr.expression(), geom, fev ),
822 tensor( this_type
const& expr, Geo_t
const& geom )
824 M_tensor_expr( expr.expression(), geom ),
828 template<
typename IM>
829 void init( IM
const& im )
831 M_tensor_expr.init( im );
833 void update( Geo_t
const& geom, Basis_i_t
const& fev, Basis_j_t
const& feu )
835 M_tensor_expr.update( geom, fev, feu );
837 void update( Geo_t
const& geom, Basis_i_t
const& fev )
839 M_tensor_expr.update( geom, fev );
841 void update( Geo_t
const& geom )
843 M_tensor_expr.update( geom );
845 void update( Geo_t
const& geom, uint16_type face )
847 M_tensor_expr.update( geom, face );
852 evalij( uint16_type i, uint16_type j )
const
854 return M_tensor_expr.evalij( i, j );
859 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q )
const
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";
865 template<
int PatternContext>
867 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
868 mpl::int_<PatternContext> )
const
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";
878 evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q )
const
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";
886 evalq( uint16_type c1, uint16_type c2, uint16_type q )
const
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";
893 tensor_expr_type M_tensor_expr;
901 value_type value()
const
903 return __expression.value();
906 value_type grad(
int __ith )
const
908 return __expression.grad( __ith );
911 value_type
hessian(
int __i,
int __j )
const
913 return __expression.hessian( __i, __j );
924 bool isSymetric()
const
926 return M_expr.isSymetric();
929 expression_type
const& expression()
const
934 const std::string& tag()
const
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
956 DVLOG(2) <<
"calling assemble(u,v)\n";
957 M_expr.assemble( __u, __v, __f );
958 DVLOG(2) <<
"calling assemble(u,v) done\n";
961 template<
typename Elem1,
typename FormType>
962 void assemble( boost::shared_ptr<Elem1>
const& __v,
963 FormType& __f )
const
965 DVLOG(2) <<
"calling assemble(v)\n";
966 M_expr.assemble( __v, __f );
967 DVLOG(2) <<
"calling assemble(v) done\n";
971 ublas::matrix<typename expression_type::value_type>
974 return M_expr.evaluate();
985 mutable expression_type M_expr;
986 const std::string M_tag;
988 template<
typename ExprT>
990 Expr< PrintExpr<ExprT> >
991 print( ExprT v, std::string t=
"" )
993 typedef PrintExpr<ExprT> print_t;
994 return Expr< print_t >( print_t( v, t ) );
1004 template<
typename ExprT>
1009 static const size_type context = ExprT::context;
1010 static const bool is_terminal =
false;
1012 static const uint16_type imorder = ExprT::imorder;
1013 static const bool imIsPoly = ExprT::imIsPoly;
1015 template<
typename Func>
1016 struct HasTestFunction
1018 static const bool result = ExprT::template HasTestFunction<Func>::result;
1021 template<
typename Func>
1022 struct HasTrialFunction
1024 static const bool result = ExprT::template HasTrialFunction<Func>::result;
1031 typedef ExprT expression_type;
1032 typedef typename expression_type::value_type value_type;
1033 typedef Trans<ExprT> this_type;
1041 explicit Trans( expression_type
const & __expr )
1053 template<
typename TheExpr>
1056 typedef Trans<typename expression_type::template Lambda<TheExpr>::type> type;
1058 template<
typename TheExpr>
1059 typename Lambda<TheExpr>::type
1060 operator()( TheExpr
const& e ) {
return trans(M_expr(e)); }
1062 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t>
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;
1068 typedef typename Transpose<typename tensor_expr_type::shape>::type shape;
1070 template <
class Args>
struct sig
1072 typedef value_type type;
1077 static const bool value = tensor_expr_type::is_zero::value;
1080 tensor( this_type
const& expr,
1081 Geo_t
const& geom, Basis_i_t
const& fev, Basis_j_t
const& feu )
1083 M_tensor_expr( expr.expression(), geom, fev, feu )
1086 tensor( this_type
const& expr,
1087 Geo_t
const& geom, Basis_i_t
const& fev )
1089 M_tensor_expr( expr.expression(), geom, fev )
1092 tensor( this_type
const& expr, Geo_t
const& geom )
1094 M_tensor_expr( expr.expression(), geom )
1097 template<
typename IM>
1098 void init( IM
const& im )
1100 M_tensor_expr.init( im );
1103 void update( Geo_t
const& geom, Basis_i_t
const& fev, Basis_j_t
const& feu )
1105 M_tensor_expr.update( geom, fev, feu );
1107 void update( Geo_t
const& geom, Basis_i_t
const& fev )
1109 M_tensor_expr.update( geom, fev );
1111 void update( Geo_t
const& geom )
1113 M_tensor_expr.update( geom );
1115 void update( Geo_t
const& geom, uint16_type face )
1117 M_tensor_expr.update( geom, face );
1122 evalij( uint16_type i, uint16_type j )
const
1124 return M_tensor_expr.evalij( i, j );
1129 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q )
const
1131 return M_tensor_expr.evalijq( i, j, c2, c1, q );
1133 template<
int PatternContext>
1135 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
1136 mpl::int_<PatternContext> )
const
1138 return M_tensor_expr.evalijq( i, j, c2, c1, q, mpl::int_<PatternContext>() );
1143 evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q )
const
1145 return M_tensor_expr.evaliq( i, c2, c1, q );
1149 evalq( uint16_type c1, uint16_type c2, uint16_type q )
const
1151 return M_tensor_expr.evalq( c2, c1, q );
1154 tensor_expr_type M_tensor_expr;
1163 bool isSymetric()
const
1165 return M_expr.isSymetric();
1168 expression_type
const& expression()
const
1191 mutable expression_type M_expr;
1194 template<
typename ExprT>
1196 Expr< Trans<ExprT> >
1199 typedef Trans<ExprT> trans_t;
1200 return Expr< trans_t >( trans_t( v ) );
1205 class Cst :
public CstBase
1211 static const size_type context = vm::JACOBIAN;
1212 static const bool is_terminal =
false;
1214 static const uint16_type imorder = 0;
1215 static const bool imIsPoly =
true;
1217 template<
typename Func>
1218 struct HasTestFunction
1220 static const bool result =
false;
1222 template<
typename Func>
1223 struct HasTrialFunction
1225 static const bool result =
false;
1228 typedef typename mpl::if_<boost::is_reference_wrapper<T>,
1230 mpl::identity<mpl::identity<T> > >::type::type::type value_type;
1232 typedef Cst<T> expression_type;
1234 constexpr
explicit Cst(
const T& value )
1240 Cst( Cst
const& __cst )
1242 M_constant( __cst.M_constant )
1250 M_constant = c.M_constant;
1254 constexpr value_type value()
const
1259 constexpr value_type evaluate()
const
1263 constexpr value_type evaluate(
bool )
const
1268 template<
typename TheExpr>
1271 typedef expression_type type;
1273 template<
typename TheExpr>
1274 typename Lambda<TheExpr>::type
1275 operator()( TheExpr
const& e ) {
return typename Lambda<TheExpr>::type(M_constant); }
1277 template<
typename Geo_t,
typename Basis_i_t=mpl::
void_,
typename Basis_j_t = Basis_i_t>
1280 typedef typename Cst<T>::expression_type expression_type;
1281 typedef typename Cst<T>::value_type value_type;
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;
1289 template<
typename Indq,
typename Indi,
typename Indj>
1292 typedef value_type type;
1297 static const bool value =
false;
1300 tensor( expression_type
const& expr,
1301 Geo_t
const& , Basis_i_t
const& , Basis_j_t
const& )
1303 M_constant( expr.value() )
1306 tensor( expression_type
const& expr,
1307 Geo_t
const& , Basis_i_t
const& )
1309 M_constant( expr.value() )
1312 tensor( expression_type
const& expr, Geo_t
const& )
1314 M_constant( expr.value() )
1317 template<
typename IM>
1318 void init( IM
const& )
1321 void update( Geo_t
const&, Basis_i_t
const& , Basis_j_t
const& )
1324 void update( Geo_t
const& , Basis_i_t
const& )
1327 void update( Geo_t
const& )
1330 void update( Geo_t
const&, uint16_type )
1334 constexpr value_type
1335 evalij( uint16_type , uint16_type )
const
1341 constexpr value_type
1342 evalijq( uint16_type , uint16_type , uint16_type , uint16_type , uint16_type )
const
1346 template<
int PatternContext>
1347 constexpr value_type
1348 evalijq( uint16_type , uint16_type , uint16_type , uint16_type , uint16_type ,
1349 mpl::int_<PatternContext> )
const
1354 constexpr value_type
1355 evaliq( uint16_type , uint16_type , uint16_type , uint16_type )
const
1359 constexpr value_type
1360 evalq( uint16_type , uint16_type , uint16_type )
const
1364 const value_type M_constant;
1368 Cst() : M_constant( 0 )
1376 template<
typename T>
1381 typedef Cst<T> cst_t;
1382 return Expr< cst_t >( cst_t( v ) );
1385 template<
typename T>
1390 typedef Cst<T> cst_t;
1391 return Expr< cst_t >( cst_t( v ) );
1393 template<
typename T>
1395 Expr< Cst<boost::reference_wrapper<T> > >
1398 typedef Cst<boost::reference_wrapper<T> > cst_t;
1399 return Expr< cst_t >( cst_t( boost::ref( v ) ) );
1401 template<
typename T>
1403 Expr< Cst<boost::reference_wrapper<T> > >
1404 constant_ref( T& v )
1406 return cst_ref( v );
1416 static const bool is_terminal =
false;
1418 static const uint16_type imorder = 0;
1419 static const bool imIsPoly =
true;
1422 template<
typename Func>
1423 struct HasTestFunction
1425 static const bool result =
false;
1427 template<
typename Func>
1428 struct HasTrialFunction
1430 static const bool result =
false;
1434 typedef One<CType> this_type;
1436 typedef double value_type;
1439 One( One
const& ) {}
1441 template<
typename TheExpr>
1444 typedef this_type type;
1446 template<
typename TheExpr>
1447 typename Lambda<TheExpr>::type
1448 operator()( TheExpr
const& e ) {
return this_type(); }
1451 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t = Basis_i_t>
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> ) );
1466 typedef typename expression_type::value_type value_type;
1468 static const uint16_type nComponents = gmc_type::nDim;
1469 static const int16_type vector_comp = ( CType==-1 )?1:CType;
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;
1477 static const bool value =
false;
1480 tensor( expression_type
const& ,
1485 M_one( nComponents, vector_comp )
1489 tensor( expression_type
const& ,
1493 M_one( nComponents, vector_comp )
1496 tensor( expression_type
const& ,
1499 M_one( nComponents, vector_comp )
1505 template<
typename IM>
1506 void init( IM
const& )
1509 void update( Geo_t
const& , Basis_i_t
const& , Basis_j_t
const& )
1512 void update( Geo_t
const& , Basis_i_t
const& )
1515 void update( Geo_t
const& )
1518 void update( Geo_t
const& , uint16_type )
1522 FEELPP_STRONG_INLINE value_type
1523 evalijq( uint16_type , uint16_type , uint16_type c1, uint16_type , uint16_type )
const
1525 return ( gmc_type::nDim>=c1 )&&( ( c1==(uint16_type)CType ) || ( CType==-1 ) );
1528 template<
int PatternContext>
1529 FEELPP_STRONG_INLINE value_type
1530 evalijq( uint16_type , uint16_type , uint16_type c1, uint16_type , uint16_type ,
1531 mpl::int_<PatternContext> )
const
1533 return ( gmc_type::nDim>=c1 )&&( ( c1==(uint16_type)CType ) || ( CType==-1 ) );
1537 FEELPP_STRONG_INLINE value_type
1538 evaliq( uint16_type , uint16_type c1, uint16_type , uint16_type )
const
1540 return ( gmc_type::nDim>=c1 )&&( ( c1==(uint16_type)CType ) || ( CType==-1 ) );
1543 FEELPP_STRONG_INLINE value_type
1544 evalq( uint16_type c1, uint16_type , uint16_type )
const
1546 return ( gmc_type::nDim>=c1 )&&( ( c1==(uint16_type)CType ) || ( CType==-1 ) );
1558 return Expr< One<-1> >( One<-1>() );
1565 return Expr< One<0> >( One<0>() );
1572 return Expr< One<1> >( One<1>() );
1579 return Expr< One<2> >( One<2>() );
1585 return Expr< One<0> >( One<0>() );
1592 return Expr< One<1> >( One<1>() );
1599 return Expr< One<2> >( One<2>() );
1608 template <
class ExprT >
1613 static const size_type context = ExprT::context;
1614 static const bool is_terminal =
false;
1616 static const uint16_type imorder = ExprT::imorder;
1617 static const bool imIsPoly = ExprT::imIsPoly;
1619 template<
typename Func>
1620 struct HasTestFunction
1622 static const bool result = ExprT::template HasTestFunction<Func>::result;
1625 template<
typename Func>
1626 struct HasTrialFunction
1628 static const bool result = ExprT::template HasTrialFunction<Func>::result;
1631 typedef ExprT expression_type;
1632 typedef typename ExprT::value_type value_type;
1633 typedef UnaryPlus<ExprT> this_type;
1635 UnaryPlus(
const ExprT& expr )
1640 UnaryPlus( UnaryPlus
const& e )
1646 expression_type
const& expression()
const
1651 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t = Basis_i_t>
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;
1660 static const bool value = tensor_expr_type::is_zero::value;
1663 tensor( this_type
const& expr,
1665 Basis_i_t
const& fev,
1666 Basis_j_t
const& feu )
1668 M_t_expr( expr.expression(), geom, fev, feu )
1670 tensor( this_type
const& expr,
1672 Basis_i_t
const& fev )
1674 M_t_expr( expr.expression(), geom, fev )
1676 tensor( this_type
const& expr,
1679 M_t_expr( expr.expression(), geom )
1681 template<
typename IM>
1682 void init( IM
const& im )
1684 M_t_expr.init( im );
1686 void update( Geo_t
const& geom, Basis_i_t
const& fev , Basis_j_t
const& feu )
1688 M_t_expr.update( geom, fev, feu );
1690 void update( Geo_t
const& geom, Basis_i_t
const& fev )
1692 M_t_expr.update( geom, fev );
1694 void update( Geo_t
const& geom )
1696 M_t_expr.update( geom );
1698 void update( Geo_t
const& geom, uint16_type face )
1700 M_t_expr.update( geom, face );
1704 evalij( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2 )
const
1706 return M_t_expr.evalij( i, j, c1, c2 );
1710 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q )
const
1712 return M_t_expr.evalijq( i, j, c1, c2, q );
1714 template<
int PatternContext>
1716 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
1717 mpl::int_<PatternContext> )
const
1719 return M_t_expr.evalijq( i, j, c1, c2, q, mpl::int_<PatternContext>() );
1723 evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q )
const
1725 return M_t_expr.evaliq( i, c1, c2, q );
1728 evalq( uint16_type c1, uint16_type c2, uint16_type q )
const
1730 return M_t_expr.evalq( c1, c2, q );
1732 tensor_expr_type M_t_expr;
1737 expression_type M_expr;
1739 template <
class T>
inline
1740 Expr< UnaryPlus< Expr<T> > >
1741 operator + (
const Expr<T>& expr )
1743 typedef UnaryPlus< Expr<T> > expr_t;
1745 return Expr< expr_t >( expr_t( expr ) );
1752 template <
class ExprT >
1757 static const size_type context = ExprT::context;
1758 static const bool is_terminal =
false;
1760 static const uint16_type imorder = ExprT::imorder;
1761 static const bool imIsPoly = ExprT::imIsPoly;
1763 template<
typename Func>
1764 struct HasTestFunction
1766 static const bool result = ExprT::template HasTestFunction<Func>::result;
1769 template<
typename Func>
1770 struct HasTrialFunction
1772 static const bool result = ExprT::template HasTrialFunction<Func>::result;
1775 typedef ExprT expression_type;
1776 typedef typename ExprT::value_type value_type;
1777 typedef UnaryMinus<ExprT> this_type;
1779 UnaryMinus(
const ExprT& expr )
1786 UnaryMinus( UnaryMinus
const& u )
1793 UnaryMinus( UnaryMinus&& u )
1807 expression_type
const& expression()
const
1812 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t = Basis_i_t>
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;
1821 static const bool value = tensor_expr_type::is_zero::value;
1824 tensor( this_type
const& expr,
1826 Basis_i_t
const& fev,
1827 Basis_j_t
const& feu )
1829 M_t_expr( expr.expression(), geom, fev, feu )
1831 tensor( this_type
const& expr,
1833 Basis_i_t
const& fev )
1835 M_t_expr( expr.expression(), geom, fev )
1837 tensor( this_type
const& expr,
1840 M_t_expr( expr.expression(), geom )
1842 template<
typename IM>
1843 void init( IM
const& im )
1845 M_t_expr.init( im );
1847 void update( Geo_t
const& geom, Basis_i_t
const& fev , Basis_j_t
const& feu )
1849 M_t_expr.update( geom, fev, feu );
1851 void update( Geo_t
const& geom, Basis_i_t
const& fev )
1853 M_t_expr.update( geom, fev );
1855 void update( Geo_t
const& geom )
1857 M_t_expr.update( geom );
1859 void update( Geo_t
const& geom, uint16_type face )
1861 M_t_expr.update( geom, face );
1865 evalij( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2 )
const
1867 return -M_t_expr.evalij( i, j, c1, c2 );
1871 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q )
const
1873 return -M_t_expr.evalijq( i, j, c1, c2, q );
1875 template<
int PatternContext>
1877 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
1878 mpl::int_<PatternContext> )
const
1880 return -M_t_expr.evalijq( i, j, c1, c2, q, mpl::int_<PatternContext>() );
1884 evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q )
const
1886 return -M_t_expr.evaliq( i, c1, c2, q );
1889 evalq( uint16_type c1, uint16_type c2, uint16_type q )
const
1891 return -M_t_expr.evalq( c1, c2, q );
1893 tensor_expr_type M_t_expr;
1898 expression_type M_expr;
1900 template <
class T>
inline
1901 Expr< UnaryMinus< Expr<T> > >
1902 operator - (
const Expr<T>& expr )
1904 typedef UnaryMinus< Expr<T> > expr_t;
1906 return Expr< expr_t >( expr_t( expr ) );
1909 template <
typename ExprT1,
typename ExprT2 >
1913 static const size_type context = ExprT1::context | ExprT2::context;
1914 static const bool is_terminal =
false;
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;
1919 template<
typename Func>
1920 struct HasTestFunction
1922 static const bool result =
false;
1924 template<
typename Func>
1925 struct HasTrialFunction
1927 static const bool result =
false;
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 )
1937 M_expr_1( __expr1 ),
1940 DVLOG(2) <<
"OpMax::OpMax default constructor\n";
1943 OpMax( OpMax
const& __vfp )
1945 M_expr_1( __vfp.M_expr_1 ),
1946 M_expr_2( __vfp.M_expr_2 )
1948 DVLOG(2) <<
"OpMax::OpMax copy constructor\n";
1951 expression_1_type
const& left()
const
1955 expression_2_type
const& right()
const
1960 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t = Basis_i_t>
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;
1967 typedef typename strongest_numeric_type<
typename l_type::value_type,
1968 typename r_type::value_type>::type value_type;
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;
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;
1983 static const bool value =
false;
1986 tensor( expression_type
const& expr, Geo_t
const& geom,
1987 Basis_i_t
const& fev, Basis_j_t
const& feu )
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 )
1994 tensor( expression_type
const& expr, Geo_t
const& geom,
1995 Basis_i_t
const& fev )
1997 M_gmc( fusion::at_key<key_type>( geom ).get() ),
1998 M_left( expr.left(), geom, fev ),
1999 M_right( expr.right(), geom, fev )
2002 tensor( expression_type
const& expr, Geo_t
const& geom )
2004 M_gmc( fusion::at_key<key_type>( geom ).get() ),
2005 M_left( expr.left(), geom ),
2006 M_right( expr.right(), geom )
2009 template<
typename IM>
2010 void init( IM
const& im )
2015 void update( Geo_t
const& geom, Basis_i_t
const& , Basis_j_t
const& )
2019 void update( Geo_t
const& geom, Basis_i_t
const& )
2023 void update( Geo_t
const& geom )
2025 M_gmc = fusion::at_key<key_type>( geom ).
get();
2026 M_left.update( geom );
2027 M_right.update( geom );
2030 void update( Geo_t
const& geom, uint16_type face )
2032 M_gmc = fusion::at_key<key_type>( geom ).
get();
2033 M_left.update( geom, face );
2034 M_right.update( geom, face );
2039 evalijq( uint16_type , uint16_type , uint16_type c1, uint16_type c2, uint16_type q )
const
2041 return evalq( c1, c2, q );
2043 template<
int PatternContext>
2045 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
2046 mpl::int_<PatternContext> )
const
2048 Feel::detail::ignore_unused_variable_warning( i );
2049 Feel::detail::ignore_unused_variable_warning( j );
2050 return evalq( c1, c2, q );
2054 evaliq( uint16_type , uint16_type c1, uint16_type c2, uint16_type q )
const
2056 return evalq( c1, c2, q );
2061 evalq( uint16_type c1, uint16_type c2, uint16_type q )
const
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 );
2077 expression_1_type M_expr_1;
2078 expression_2_type M_expr_2;
2081 template<
typename ExprT1,
typename ExprT2>
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
2090 max( ExprT1
const& __e1, ExprT2
const& __e2 )
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 ) ) );
2102 template <
typename ExprT1,
typename ExprT2 >
2106 static const size_type context = ExprT1::context | ExprT2::context;
2107 static const bool is_terminal =
false;
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;
2112 template<
typename Func>
2113 struct HasTestFunction
2115 static const bool result =
false;
2118 template<
typename Func>
2119 struct HasTrialFunction
2121 static const bool result =
false;
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;
2130 explicit OpMin( expression_1_type
const& __expr1, expression_2_type
const& __expr2 )
2132 M_expr_1( __expr1 ),
2135 DVLOG(2) <<
"OpMin::OpMin default constructor\n";
2138 OpMin( OpMin
const& __vfp )
2140 M_expr_1( __vfp.M_expr_1 ),
2141 M_expr_2( __vfp.M_expr_2 )
2143 DVLOG(2) <<
"OpMin::OpMin copy constructor\n";
2146 expression_1_type
const& left()
const
2150 expression_2_type
const& right()
const
2155 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t = Basis_i_t>
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;
2162 typedef typename strongest_numeric_type<
typename l_type::value_type,
2163 typename r_type::value_type>::type value_type;
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;
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;
2178 static const bool value =
false;
2181 tensor( expression_type
const& expr, Geo_t
const& geom,
2182 Basis_i_t
const& fev, Basis_j_t
const& feu )
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 )
2189 tensor( expression_type
const& expr, Geo_t
const& geom,
2190 Basis_i_t
const& fev )
2192 M_gmc( fusion::at_key<key_type>( geom ).get() ),
2193 M_left( expr.left(), geom, fev ),
2194 M_right( expr.right(), geom, fev )
2197 tensor( expression_type
const& expr, Geo_t
const& geom )
2199 M_gmc( fusion::at_key<key_type>( geom ).get() ),
2200 M_left( expr.left(), geom ),
2201 M_right( expr.right(), geom )
2204 template<
typename IM>
2205 void init( IM
const& im )
2210 void update( Geo_t
const& geom, Basis_i_t
const& , Basis_j_t
const& )
2214 void update( Geo_t
const& geom, Basis_i_t
const& )
2218 void update( Geo_t
const& geom )
2220 M_gmc = fusion::at_key<key_type>( geom ).
get();
2221 M_left.update( geom );
2222 M_right.update( geom );
2225 void update( Geo_t
const& geom, uint16_type face )
2227 M_gmc = fusion::at_key<key_type>( geom ).
get();
2228 M_left.update( geom, face );
2229 M_right.update( geom, face );
2234 evalijq( uint16_type , uint16_type , uint16_type c1, uint16_type c2, uint16_type q )
const
2236 return evalq( c1, c2, q );
2238 template<
int PatternContext>
2240 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
2241 mpl::int_<PatternContext> )
const
2243 return evalq( c1, c2, q );
2247 evaliq( uint16_type , uint16_type c1, uint16_type c2, uint16_type q )
const
2249 return evalq( c1, c2, q );
2254 evalq( uint16_type c1, uint16_type c2, uint16_type q )
const
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 );
2270 expression_1_type M_expr_1;
2271 expression_2_type M_expr_2;
2274 template<
typename ExprT1,
typename ExprT2>
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 )
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 ) ) );
2294 template <
typename ExprT1,
typename ExprT2 >
2299 static const size_type context = ExprT1::context|ExprT2::context;
2300 static const bool is_terminal =
false;
2306 static const uint16_type imorder = ExprT1::imorder;
2307 static const bool imIsPoly = ExprT1::imIsPoly && ExprT2::imIsPoly;
2309 template<
typename Func>
2310 struct HasTestFunction
2312 static const bool result =
false;
2315 template<
typename Func>
2316 struct HasTrialFunction
2318 static const bool result =
false;
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;
2329 BOOST_STATIC_ASSERT( ::boost::is_arithmetic<value_1_type>::value &&
2330 ::boost::is_arithmetic<value_2_type>::value );
2332 explicit Pow( expression_1_type
const& __expr1, expression_2_type
const& __expr2 )
2334 M_expr_1( __expr1 ),
2337 DVLOG(2) <<
"Pow::Pow default constructor\n";
2340 Pow( Pow
const& __vfp )
2342 M_expr_1( __vfp.M_expr_1 ),
2343 M_expr_2( __vfp.M_expr_2 )
2345 DVLOG(2) <<
"Pow::Pow copy constructor\n";
2348 bool isSymetric()
const
2353 expression_1_type
const& left()
const
2357 expression_2_type
const& right()
const
2362 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t = Basis_i_t>
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;
2369 typedef typename strongest_numeric_type<
typename l_type::value_type,
2370 typename r_type::value_type>::type value_type;
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;
2380 typedef typename Eigen::Matrix<value_type,shape::M,shape::N> loc_type;
2384 static const bool value = l_type::is_zero::value;
2387 tensor( expression_type
const& expr, Geo_t
const& geom,
2388 Basis_i_t
const& fev, Basis_j_t
const& feu )
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()] )
2397 tensor( expression_type
const& expr, Geo_t
const& geom,
2398 Basis_i_t
const& fev )
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()] )
2407 tensor( expression_type
const& expr, Geo_t
const& geom )
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()] )
2416 template<
typename IM>
2417 void init( IM
const& im )
2422 void update( Geo_t
const& geom, Basis_i_t
const& , Basis_j_t
const& )
2426 void update( Geo_t
const& geom, Basis_i_t
const& )
2430 void update( Geo_t
const& geom )
2432 const int npts = fusion::at_key<key_type>( geom ).
get()->nPoints();
2433 M_left.update( geom );
2434 M_right.update( geom );
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 )
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 );
2445 void update( Geo_t
const& geom, uint16_type face )
2447 M_gmc = fusion::at_key<key_type>( geom ).
get();
2448 M_left.update( geom, face );
2449 M_right.update( geom, face );
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 )
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 );
2462 evalijq( uint16_type , uint16_type , uint16_type c1, uint16_type c2, uint16_type q )
const
2464 return evalq( c1, c2, q );
2466 template<
int PatternContext>
2468 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
2469 mpl::int_<PatternContext> )
const
2471 Feel::detail::ignore_unused_variable_warning( i );
2472 Feel::detail::ignore_unused_variable_warning( j );
2473 return evalq( c1, c2, q );
2477 evaliq( uint16_type , uint16_type c1, uint16_type c2, uint16_type q )
const
2479 return evalq( c1, c2, q );
2483 evalq( uint16_type c1, uint16_type c2, uint16_type q )
const
2485 return evalq( c1, c2, q, mpl::int_<shape::rank>() );
2488 evalq( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<0> )
const
2490 Feel::detail::ignore_unused_variable_warning( c1 );
2491 Feel::detail::ignore_unused_variable_warning( c2 );
2492 return M_loc[q]( 0,0 );
2495 evalq( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<1> )
const
2497 if ( shape::M > shape::N )
2498 return M_loc[q]( c1,0 );
2500 return M_loc[q]( 0,c2 );
2503 evalq( uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<2> )
const
2505 return M_loc[q]( c1,c2 );
2513 boost::multi_array<loc_type,1> M_loc;
2519 expression_1_type M_expr_1;
2520 expression_2_type M_expr_2;
2523 template<
typename ExprT1,
typename ExprT2>
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 )
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 ) ) );
2544 template<
typename ExprT1,
typename ExprT2>
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 )
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 ) ) );
2576 template<
int GeoId,
typename ExprT>
2581 static const size_type context = ExprT::context;
2582 static const bool is_terminal =
false;
2584 static const uint16_type imorder = ExprT::imorder;
2585 static const bool imIsPoly = ExprT::imIsPoly;
2590 typedef ExprT expression_type;
2591 typedef typename expression_type::value_type value_type;
2592 typedef EvalFace<GeoId,ExprT> this_type;
2600 explicit EvalFace( expression_type
const & __expr )
2613 template<
typename VecGeo_t,
typename Basis_i_t = boost::none_t,
typename Basis_j_t = Basis_i_t>
2616 typedef typename fusion::result_of::at_c<VecGeo_t,GeoId>::type Geo_t;
2618 typedef typename expression_type::template tensor<Geo_t,
2620 Basis_j_t> tensor_expr_type;
2621 typedef typename tensor_expr_type::value_type value_type;
2623 template <
class Args>
struct sig
2625 typedef value_type type;
2628 tensor( this_type
const& expr,
2629 VecGeo_t
const& geom, Basis_i_t
const& fev, Basis_j_t
const& feu )
2631 M_tensor_expr( expr.expression(), fusion::at_c<GeoId>( geom ), fev, feu )
2634 tensor( this_type
const& expr,
2635 VecGeo_t
const& geom, Basis_i_t
const& fev )
2637 M_tensor_expr( expr.expression(), fusion::at_c<GeoId>( geom ), fev )
2640 tensor( this_type
const& expr,
2641 VecGeo_t
const& geom )
2643 M_tensor_expr( expr.expression(), fusion::at_c<GeoId>( geom ) )
2645 template<
typename IM>
2646 void init( IM
const& im )
2648 M_tensor_expr.init( im );
2650 void update( VecGeo_t
const& geom, Basis_i_t
const& fev, Basis_j_t
const& feu )
2652 M_tensor_expr.update( fusion::at_c<GeoId>( geom ), fev, feu );
2654 void update( VecGeo_t
const& geom, Basis_i_t
const& fev )
2656 M_tensor_expr.update( fusion::at_c<GeoId>( geom ), fev );
2658 void update( VecGeo_t
const& geom )
2660 M_tensor_expr.update( fusion::at_c<GeoId>( geom ) );
2665 evalij( uint16_type i, uint16_type j )
const
2667 return M_tensor_expr.evalij( i, j );
2672 evalijq( uint16_type i, uint16_type j,
int q )
const
2674 return M_tensor_expr.evalijq( i, j, q );
2679 evaliq( uint16_type i,
int q )
const
2681 return M_tensor_expr.evaliq( i, q );
2685 evalq(
int q,
int c )
const
2687 return M_tensor_expr.evalq( q, c );
2690 tensor_expr_type M_tensor_expr;
2699 expression_type
const& expression()
const
2721 mutable expression_type M_expr;
2724 template<
int GeoId,
typename ExprT>
2726 Expr< EvalFace<GeoId, ExprT> >
2727 evalface( ExprT
const& v )
2729 typedef EvalFace<GeoId, ExprT> eval_t;
2730 return Expr< eval_t >( eval_t( v ) );
2733 template <
class Element,
int Type>
2738 static const size_type context = vm::JACOBIAN |vm::POINT;
2739 static const bool is_terminal =
false;
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;
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;
2757 static const uint16_type imorder = element_type::functionspace_type::basis_type::nOrder;
2758 static const bool imIsPoly =
true;
2760 template<
typename Func>
2761 struct HasTestFunction
2763 static const bool result = ( Type==0 );
2766 template<
typename Func>
2767 struct HasTrialFunction
2769 static const bool result = ( Type==1 );
2773 GElem ( std::map<
size_type,std::vector<element_ptrtype> >
const& v )
2777 typename basis_type::iterator it = M_basis.begin();
2778 typename basis_type::iterator en = M_basis.end();
2780 for ( ; it != en; ++it )
2781 for ( uint16_type i = 0; i < it->second.size(); ++i )
2782 it->second[i]->updateGlobalValues();
2784 GElem( GElem
const& op )
2786 M_basis ( op.M_basis )
2791 basis_type
const& basis()
const
2797 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t = Basis_i_t>
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;
2807 typedef typename mpl::if_<mpl::equal_to<mpl::int_<rank>,
2809 mpl::identity<Shape<gmc_type::NDim, Scalar, false> >,
2810 typename mpl::if_<mpl::equal_to<mpl::int_<rank>,
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;
2819 typedef typename expression_type::value_type value_type;
2823 static const bool value =
false;
2826 tensor( expression_type
const& 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() )
2837 for ( uint16_type i = 0; i < M_loc.size(); ++i )
2838 M_loc[i].resize( boost::extents[M_pc.nPoints()][nComponents1][nComponents2] );
2840 tensor( expression_type
const& 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() )
2850 for ( uint16_type i = 0; i < M_loc.size(); ++i )
2851 M_loc[i].resize( boost::extents[M_pc.nPoints()] );
2853 tensor( expression_type
const& 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() )
2862 for ( uint16_type i = 0; i < M_loc.size(); ++i )
2863 M_loc[i].resize( boost::extents[M_pc.nPoints()] );
2865 template<
typename IM>
2866 void init( IM
const& )
2870 void update( Geo_t
const& geom, Basis_i_t
const& , Basis_j_t
const& )
2874 void update( Geo_t
const& geom, Basis_i_t
const& )
2878 void update( Geo_t
const& geom )
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();
2884 FEELPP_ASSERT( it != en )( fusion::at_key<key_type>( geom )->
id() ).error (
"invalid basis function to integrate" );
2886 for ( uint16_type i = 0; i < M_loc.size(); ++i )
2893 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q )
const
2896 return M_loc[i]( c1,c2,q );
2898 return M_loc[j]( c1,c2,q );
2900 template<
int PatternContext>
2902 evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q, mpl::int_<PatternContext> )
const
2905 return M_loc[i]( c1,c2,q );
2907 return M_loc[j]( c1,c2,q );
2911 evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q )
const
2913 return M_loc[i]( c1,c2,q );
2916 evalq( uint16_type , uint16_type , uint16_type )
const
2921 this_type
const& M_expr;
2924 std::vector<typename element_type::id_type> M_loc;
2931 template<
typename Elem>
2933 Expr< GElem<Elem,1> >
2934 basist( std::map<
size_type,std::vector<boost::shared_ptr<Elem> > >
const& v )
2936 typedef GElem<Elem,1> expr_t;
2937 return Expr< expr_t >( expr_t( v ) );
2939 template<
typename Elem>
2941 Expr< GElem<Elem,0> >
2942 basis( std::map<
size_type,std::vector<boost::shared_ptr<Elem> > >
const& v )
2944 typedef GElem<Elem,0> expr_t;
2945 return Expr< expr_t >( expr_t( v ) );