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
matvec.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: 2007-07-20
7 
8  Copyright (C) 2007 Universite Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __VFVec_H
30 #define __VFVec_H 1
31 
32 #include <boost/fusion/container/vector.hpp>
33 #include <boost/fusion/algorithm/iteration.hpp>
34 #include <boost/mpl/bitor.hpp>
35 #include <boost/mpl/bitwise.hpp>
36 #include <boost/mpl/transform.hpp>
37 
38 #include <boost/mpl/plus.hpp>
39 #include <boost/mpl/arithmetic.hpp>
40 
41 
42 namespace Feel
43 {
44 namespace fusion = boost::fusion;
45 namespace mpl = boost::mpl;
46 namespace vf
47 {
49 namespace detail
50 {
51 struct GetContext
52 {
53  template<typename Sig>
54  struct result;
55 
56  template<typename Lhs, typename Rhs>
57  struct result<GetContext( Lhs,Rhs )>
58 #if BOOST_VERSION < 104200
59 : boost::remove_reference<Lhs>
60 #else
61 : boost::remove_reference<Rhs>
62 #endif
63  {
64 
65  typedef typename boost::remove_reference<Lhs>::type lhs_noref_type;
66  typedef typename boost::remove_reference<Rhs>::type rhs_noref_type;
67 #if BOOST_VERSION < 104200
68  typedef typename mpl::bitor_<mpl::size_t<lhs_noref_type::context>, rhs_noref_type >::type type;
69 #else
70  typedef typename mpl::bitor_<mpl::size_t<rhs_noref_type::context>, lhs_noref_type >::type type;
71 #endif
72  };
73 #if 0
74  template<typename Expr1, typename Expr2>
75  struct result
76  {
77  typedef typename mpl::bitor_<mpl::size_t<Expr1::context>, Expr2 >::type type;
78  };
79 #endif
80  template<typename Lhs, typename Rhs>
81 #if BOOST_VERSION < 104200
82  Lhs
83 #else
84  Rhs
85 #endif
86  operator()( const Lhs& lhs, const Rhs& rhs ) const
87  {
88  return lhs | rhs;
89  }
90 };
91 struct GetImOrder
92 {
93  template<typename Sig>
94  struct result;
95 
96  template<typename Lhs, typename Rhs>
97  struct result<GetImOrder( Lhs,Rhs )>
98 #if BOOST_VERSION < 104200
99 : boost::remove_reference<Lhs>
100 #else
101 : boost::remove_reference<Rhs>
102 #endif
103  {
104 
105  typedef typename boost::remove_reference<Lhs>::type lhs_noref_type;
106  typedef typename boost::remove_reference<Rhs>::type rhs_noref_type;
107 #if BOOST_VERSION < 104200
108  typedef typename boost::mpl::max< boost::mpl::size_t<lhs_noref_type::imorder>, rhs_noref_type >::type type;
109 #else
110  typedef typename boost::mpl::max< boost::mpl::size_t<rhs_noref_type::imorder>, lhs_noref_type >::type type;
111 #endif
112  };
113 
114  template<typename Lhs, typename Rhs>
115 #if BOOST_VERSION < 104200
116  Lhs
117 #else
118  Rhs
119 #endif
120  operator()( const Lhs& lhs, const Rhs& rhs ) const
121  {
122  }
123 };
124 struct GetImIsPoly
125 {
126  template<typename Sig>
127  struct result;
128 
129  template<typename Lhs, typename Rhs>
130  struct result<GetImIsPoly( Lhs,Rhs )>
131 #if BOOST_VERSION < 104200
132 : boost::remove_reference<Lhs>
133 #else
134 : boost::remove_reference<Rhs>
135 #endif
136  {
137 
138  typedef typename boost::remove_reference<Lhs>::type lhs_noref_type;
139  typedef typename boost::remove_reference<Rhs>::type rhs_noref_type;
140 #if BOOST_VERSION < 104200
141  typedef typename boost::mpl::and_< boost::mpl::bool_<lhs_noref_type::imIsPoly>, rhs_noref_type >::type type;
142 #else
143  typedef typename boost::mpl::and_< boost::mpl::bool_<rhs_noref_type::imIsPoly>, lhs_noref_type >::type type;
144 #endif
145  };
146 
147  template<typename Lhs, typename Rhs>
148 #if BOOST_VERSION < 104200
149  Lhs
150 #else
151  Rhs
152 #endif
153  operator()( const Lhs& lhs, const Rhs& rhs ) const
154  {
155  }
156 };
157 template<typename Func>
158 struct ExprHasTestFunction
159 {
160  template<typename Sig>
161  struct result;
162 
163  template<typename Lhs, typename Rhs>
164  struct result<ExprHasTestFunction( Lhs,Rhs )>
165 #if BOOST_VERSION < 104200
166 : boost::remove_reference<Lhs>
167 #else
168 : boost::remove_reference<Rhs>
169 #endif
170  {
171 
172  typedef typename boost::remove_reference<Lhs>::type lhs_noref_type;
173  typedef typename boost::remove_reference<Rhs>::type rhs_noref_type;
174 #if BOOST_VERSION < 104200
175  typedef typename mpl::or_<mpl::bool_<lhs_noref_type::template HasTestFunction<Func>::result>, rhs_noref_type >::type type;
176 #else
177  typedef typename mpl::or_<mpl::bool_<rhs_noref_type::template HasTestFunction<Func>::result>, lhs_noref_type >::type type;
178 #endif
179  };
180  template<typename Lhs, typename Rhs>
181 #if BOOST_VERSION < 104200
182  Lhs
183 #else
184  Rhs
185 #endif
186  operator()( const Lhs& lhs, const Rhs& rhs ) const
187  {
188 #if 0
189  typedef typename boost::remove_reference<Lhs>::type lhs_noref_type;
190  typedef typename boost::remove_reference<Rhs>::type rhs_noref_type;
191  return ( lhs_noref_type::template ExprHasTestFunction<Func>::result ||
192  lhs_noref_type::template ExprHasTestFunction<Func>::result );
193 #endif
194  }
195 };
196 template<typename Func>
197 struct ExprHasTrialFunction
198 {
199  template<typename Sig>
200  struct result;
201 
202  template<typename Lhs, typename Rhs>
203  struct result<ExprHasTrialFunction( Lhs,Rhs )>
204 #if BOOST_VERSION < 104200
205 : boost::remove_reference<Lhs>
206 #else
207 : boost::remove_reference<Rhs>
208 #endif
209  {
210 
211  typedef typename boost::remove_reference<Lhs>::type lhs_noref_type;
212  typedef typename boost::remove_reference<Rhs>::type rhs_noref_type;
213 #if BOOST_VERSION < 104200
214  typedef typename mpl::or_<mpl::bool_<lhs_noref_type::template HasTrialFunction<Func>::result>, rhs_noref_type >::type type;
215 #else
216  typedef typename mpl::or_<mpl::bool_<rhs_noref_type::template HasTrialFunction<Func>::result>, lhs_noref_type >::type type;
217 #endif
218  };
219  template<typename Lhs, typename Rhs>
220 #if BOOST_VERSION < 104200
221  Lhs
222 #else
223  Rhs
224 #endif
225  operator()( const Lhs& lhs, const Rhs& rhs ) const
226  {
227 #if 0
228  return ( Lhs::template HasTrialFunction<Func>::result ||
229  Rhs::template HasTrialFunction<Func>::result );
230 #endif
231  }
232 };
233 template<typename Geo_t, typename Basis_i_t, typename Basis_j_t>
234 struct initialize_expression_gij
235 {
236  template<typename Sig>
237  struct result;
238 
239  template<typename ExprT>
240  struct result<initialize_expression_gij( ExprT )>
241  {
242  typedef typename boost::remove_reference<ExprT>::type::template tensor<Geo_t, Basis_i_t, Basis_j_t> type;
243  };
244 
245  initialize_expression_gij( Geo_t const& geom , Basis_i_t const& fev, Basis_j_t const& feu )
246  :
247  M_geom( geom ),
248  M_fev( fev ),
249  M_feu( feu )
250  {}
251 
252  template <typename ExprT>
253  typename ExprT::template tensor<Geo_t, Basis_i_t, Basis_j_t>
254  operator()( ExprT& expr ) const
255  {
256  return typename ExprT::template tensor<Geo_t, Basis_i_t, Basis_j_t>( expr, M_geom, M_fev, M_feu );
257  }
258  const Geo_t& M_geom;
259  const Basis_i_t& M_fev;
260  const Basis_j_t& M_feu;
261 };
262 template<typename Geo_t, typename Basis_i_t>
263 struct initialize_expression_gi
264 {
265  template<typename Sig>
266  struct result;
267 
268  template<typename ExprT>
269  struct result<initialize_expression_gi( ExprT )>
270  {
271  typedef typename boost::remove_reference<ExprT>::type::template tensor<Geo_t, Basis_i_t> type;
272  };
273  initialize_expression_gi( Geo_t const& geom , Basis_i_t const& fev )
274  :
275  M_geom( geom ),
276  M_fev( fev )
277  {}
278 
279  template <typename ExprT>
280  typename ExprT::template tensor<Geo_t, Basis_i_t>
281  operator()( ExprT& expr ) const
282  {
283  return typename ExprT::template tensor<Geo_t, Basis_i_t>( expr, M_geom, M_fev );
284  }
285 
286  const Geo_t& M_geom;
287  const Basis_i_t& M_fev;
288 };
289 
290 template<typename Geo_t>
291 struct initialize_expression_g
292 {
293  template<typename Sig>
294  struct result;
295 
296  template<typename ExprT>
297  struct result<initialize_expression_g( ExprT )>
298  {
299  typedef typename boost::remove_reference<ExprT>::type::template tensor<Geo_t> type;
300  };
301  initialize_expression_g( Geo_t const& geom )
302  :
303  M_geom( geom )
304  {}
305 
306  template <typename ExprT>
307  typename ExprT::template tensor<Geo_t>
308  operator()( ExprT& expr ) const
309  {
310  return typename ExprT::template tensor<Geo_t>( expr, M_geom );
311  //return typename ExprT::template tensor<Geo_t>( expr, M_geom );
312  }
313 
314  const Geo_t& M_geom;
315 };
316 
317 
318 template<typename Geo_t, typename Basis_i_t, typename Basis_j_t>
319 struct update_expression_gij
320 {
321  update_expression_gij( Geo_t const& geom , Basis_i_t const& fev, Basis_j_t const& feu )
322  :
323  M_geom( geom ),
324  M_fev( fev ),
325  M_feu( feu )
326  {}
327 
328  template <typename ExprT>
329  void operator()( ExprT& expr ) const
330  {
331  expr.update( M_geom, M_fev, M_feu );
332  }
333  const Geo_t& M_geom;
334  const Basis_i_t& M_fev;
335  const Basis_j_t& M_feu;
336 };
337 template<typename Geo_t, typename Basis_i_t>
338 struct update_expression_gi
339 {
340  update_expression_gi( Geo_t const& geom , Basis_i_t const& fev )
341  :
342  M_geom( geom ),
343  M_fev( fev )
344  {}
345 
346  template <typename ExprT>
347  void operator()( ExprT& expr ) const
348  {
349  expr.update( M_geom, M_fev );
350  }
351  const Geo_t& M_geom;
352  const Basis_i_t& M_fev;
353 };
354 
355 template<typename Geo_t>
356 struct update_expression_g
357 {
358  update_expression_g( Geo_t const& geom )
359  :
360  M_geom( geom )
361  {}
362 
363  template <typename ExprT>
364  void operator()( ExprT& expr ) const
365  {
366  expr.update( M_geom );
367  }
368  const Geo_t& M_geom;
369 };
370 
371 
372 template<typename Geo_t>
373 struct update_expression_face_g
374 {
375  update_expression_face_g( Geo_t const& geom, uint16_type face )
376  :
377  M_geom( geom ),
378  M_face( face )
379  {}
380 
381  template <typename ExprT>
382  void operator()( ExprT& expr ) const
383  {
384  expr.update( M_geom, M_face );
385  }
386  const Geo_t& M_geom;
387  const uint16_type M_face;
388 };
389 
390 template<typename IM_t>
391 struct init_expression
392 {
393  init_expression( IM_t const& im )
394  :
395  M_im( im )
396  {}
397 
398  template <typename ExprT>
399  void operator()( ExprT& expr ) const
400  {
401  expr.init( M_im );
402  }
403  const IM_t& M_im;
404 };
405 
406 template<typename T>
407 struct evaluate_expression_gij
408 {
409  template<typename Sig>
410  struct result;
411 
412  template<typename ExprT, typename V>
413 #if BOOST_VERSION < 104200
414  struct result<evaluate_expression_gij<T>( ExprT,V )>
415 #else
416  struct result<evaluate_expression_gij<T>( V,ExprT )>
417 #endif
418 :
419  boost::remove_reference<V>
420  {};
421  evaluate_expression_gij( uint16_type indexi , uint16_type indexj, uint16_type c1, uint16_type q )
422  :
423  M_indexi( indexi ),
424  M_indexj( indexj ),
425  M_c1( c1 ),
426  M_c2( 0 ),
427  M_q( q ),
428  M_index( M_c1 ),
429  M_current( 0 )
430  {}
431 
432  evaluate_expression_gij( int n, uint16_type indexi, uint16_type indexj, uint16_type c1, uint16_type c2, uint16_type q )
433  :
434  M_indexi( indexi ),
435  M_indexj( indexj ),
436  M_c1( c1 ),
437  M_c2( c2 ),
438  M_q( q ),
439  M_index( M_c1*n+M_c2 ),
440  M_current( 0 )
441 
442  {}
443 
444  template <typename ExprT>
445  T
446 #if BOOST_VERSION < 104200
447  operator()( ExprT const& expr, T const& res ) const
448 #else
449  operator()( T const& res, ExprT const& expr ) const
450 #endif
451  {
452  T ret = res;
453 
454  if ( M_current == M_index )
455  ret = expr.evalijq( M_indexi, M_indexj, 0, 0, M_q );
456 
457  ++M_current;
458  return ret;
459  }
460  const uint16_type& M_indexi;
461  const uint16_type& M_indexj;
462  const uint16_type M_c1;
463  const uint16_type M_c2;
464  const uint16_type M_q;
465  const uint16_type M_index;
466  mutable int M_current;
467 
468 };
469 template<typename T=double>
470 struct evaluate_expression_gi
471 {
472  template<typename Sig>
473  struct result;
474 
475  template<typename ExprT,typename V>
476 #if BOOST_VERSION < 104200
477  struct result<evaluate_expression_gi<T>( ExprT,V )>
478 #else
479  struct result<evaluate_expression_gi<T>( V,ExprT )>
480 #endif
481 :
482  boost::remove_reference<V>
483  {};
484 
485  evaluate_expression_gi( uint16_type indexi, uint16_type c1, uint16_type q )
486  :
487  M_indexi( indexi ),
488  M_c1( c1 ),
489  M_c2( 0 ),
490  M_q( q ),
491  M_index( M_c1 ),
492  M_current( 0 )
493  {}
494 
495  evaluate_expression_gi( int n, uint16_type indexi, uint16_type c1, uint16_type c2, uint16_type q )
496  :
497  M_indexi( indexi ),
498  M_c1( c1 ),
499  M_c2( c2 ),
500  M_q( q ),
501  M_index( M_c1*n+M_c2 ),
502  M_current( 0 )
503 
504  {}
505  template <typename ExprT>
506  T
507 #if BOOST_VERSION < 104200
508  operator()( ExprT const& expr, T const& res ) const
509 #else
510  operator()( T const& res, ExprT const& expr ) const
511 #endif
512  {
513  T ret = res;
514 
515  if ( M_current == M_index )
516  ret = expr.evaliq( M_indexi, 0, 0, M_q );
517 
518  ++M_current;
519  return ret;
520 
521  }
522  const uint16_type& M_indexi;
523  const uint16_type M_c1;
524  const uint16_type M_c2;
525  const uint16_type M_q;
526  const uint16_type M_index;
527  mutable int M_current;
528 };
529 
530 template<typename T=double>
531 struct evaluate_expression_g
532 {
533  template<typename Sig>
534  struct result;
535 
536  template<typename ExprT,typename V>
537 #if BOOST_VERSION < 104200
538  struct result<evaluate_expression_g<T>( ExprT,V )>
539 #else
540  struct result<evaluate_expression_g<T>( V,ExprT )>
541 #endif
542 :
543  boost::remove_reference<V>
544  {};
545 
546  evaluate_expression_g( uint16_type c1, uint16_type q )
547  :
548  M_c1( c1 ),
549  M_c2( 0 ),
550  M_q( q ),
551  M_index( M_c1 ),
552  M_current( 0 )
553  {}
554 
555  evaluate_expression_g( int n, uint16_type c1, uint16_type c2, uint16_type q )
556  :
557  M_c1( c1 ),
558  M_c2( c2 ),
559  M_q( q ),
560  M_index( M_c1*n+M_c2 ),
561  M_current( 0 )
562 
563  {}
564 
565  template <typename ExprT>
566  T
567 #if BOOST_VERSION < 104200
568  operator()( ExprT const& expr, T const& res ) const
569 #else
570  operator()( T const& res, ExprT const& expr ) const
571 #endif
572  {
573  T ret = res;
574 
575  if ( M_current == M_index )
576  ret = expr.evalq( 0, 0, M_q );
577 
578  ++M_current;
579  return ret;
580  }
581  const uint16_type M_c1;
582  const uint16_type M_c2;
583  const uint16_type M_q;
584  const uint16_type M_index;
585  mutable int M_current;
586 };
594 template<typename VectorExpr>
595 class Vec
596 {
597 public:
598 
599 
603  static const size_type context = fusion::result_of::accumulate<VectorExpr,mpl::size_t<0>,GetContext>::type::value;
604  //static const size_type context = fusion::accumulate( VectorExpr(),size_type(0),GetContext() );
605  static const bool is_terminal = false;
606  static const uint16_type imorder = fusion::result_of::accumulate<VectorExpr,mpl::size_t<0>,GetImOrder>::type::value;
607  static const bool imIsPoly = fusion::result_of::accumulate<VectorExpr,mpl::bool_<true>,GetImIsPoly>::type::value;
608 
609  template<typename Func>
610  struct HasTestFunction
611  {
612  static const bool result = fusion::result_of::accumulate<VectorExpr,mpl::bool_<false>,ExprHasTestFunction<Func> >::type::value;
613  };
614  template<typename Func>
615  struct HasTrialFunction
616  {
617  static const bool result = fusion::result_of::accumulate<VectorExpr,mpl::bool_<false>,ExprHasTrialFunction<Func> >::type::value;
618  };
619 
620  typedef VectorExpr expression_vector_type;
621  typedef Vec<expression_vector_type> this_type;
622 
623  typedef double value_type;
624 
625  static const uint16_type vector_size = fusion::result_of::size<expression_vector_type>::type::value;
626 
628 
632 
633  Vec( VectorExpr const& expr )
634  :
635  M_expr( expr )
636  {}
637  Vec( Vec const & expr )
638  :
639  M_expr( expr.M_expr )
640  {}
641  ~Vec()
642  {}
643 
645 
649 
650 
652 
656 
657  expression_vector_type const& expression() const
658  {
659  return M_expr;
660  }
661  expression_vector_type & expression()
662  {
663  return M_expr;
664  }
665 
667 
671 
672 
674 
678 
679 
681 
682  template<typename Geo_t, typename Basis_i_t = boost::none_t, typename Basis_j_t = Basis_i_t>
683  struct tensor
684  {
685  typedef this_type expression_type;
686  typedef typename expression_type::expression_vector_type expression_vector_type;
687  typedef Shape<expression_type::vector_size, Vectorial, false, false> shape;
688  static const bool theshape = ( shape::M == expression_type::vector_size && shape::N == 1 );
689  BOOST_MPL_ASSERT_MSG( theshape,
690  INVALID_TENSOR_SHAPE_SHOULD_BE_RANK_1,
691  ( mpl::int_<shape::M>, mpl::int_<shape::N> ) );
692 
693  typedef typename boost::remove_reference<typename fusion::result_of::at<expression_vector_type, boost::mpl::int_<0> >::type>::type::value_type value_type;
694  //typedef typename expression_type::value_type value_type;
695 
696  template<typename ExprT>
697  struct ExpressionToTensor
698  {
699  typedef typename boost::remove_reference<ExprT>::type expr_type;
700  typedef typename expr_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> type;
701  };
702  typedef typename mpl::transform<expression_vector_type, ExpressionToTensor<mpl::_1>, mpl::back_inserter<fusion::vector<> > >::type tensor_vector_type;
703 
704 
705  struct is_zero
706  {
707  static const bool value = false;
708  };
709 
710  tensor( expression_type const& expr,
711  Geo_t const& geom,
712  Basis_i_t const& fev,
713  Basis_j_t const& feu )
714  :
715  M_expr( fusion::transform( expr.expression(), initialize_expression_gij<Geo_t,Basis_i_t,Basis_j_t>( geom, fev, feu ) ) )
716  {
717 
718  update( geom, fev, feu );
719  }
720  tensor( expression_type const& expr,
721  Geo_t const& geom,
722  Basis_i_t const& fev )
723  :
724  M_expr( fusion::transform( expr.expression(), initialize_expression_gi<Geo_t,Basis_i_t>( geom, fev ) ) )
725  {
726  update( geom, fev );
727  }
728  tensor( expression_type const& expr,
729  Geo_t const& geom )
730  :
731  M_expr( fusion::transform( expr.expression(), initialize_expression_g<Geo_t>( geom ) ) )
732  {
733  update( geom );
734  }
735  template<typename IM>
736  void init( IM const& im )
737  {
738  fusion::for_each( M_expr,vf::detail::init_expression<IM>( im ) );
739  }
740  void update( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
741  {
742  fusion::for_each( M_expr,vf::detail::update_expression_gij<Geo_t, Basis_i_t, Basis_j_t>( geom, fev, feu ) );
743 
744  }
745  void update( Geo_t const& geom, Basis_i_t const& fev )
746  {
747  fusion::for_each( M_expr,vf::detail::update_expression_gi<Geo_t, Basis_i_t>( geom, fev ) );
748  }
749  void update( Geo_t const& geom )
750  {
751  fusion::for_each( M_expr,vf::detail::update_expression_g<Geo_t>( geom ) );
752  }
753  void update( Geo_t const& geom, uint16_type face )
754  {
755  fusion::for_each( M_expr,vf::detail::update_expression_face_g<Geo_t>( geom, face ) );
756  }
757 
758  value_type
759  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type /*c2*/, uint16_type q ) const
760  {
761  return fusion::accumulate( M_expr, value_type( 0 ),vf::detail::evaluate_expression_gij<value_type>( i, j, c1, q ) );
762  }
763  template<int PatternContext>
764  value_type
765  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type /*c2*/, uint16_type q,
766  mpl::int_<PatternContext> ) const
767  {
768  return fusion::accumulate( M_expr, value_type( 0 ),
769  vf::detail::evaluate_expression_gij<value_type>( i, j, c1, q ) );
770 
771  }
772 
773  value_type
774  evaliq( uint16_type i, uint16_type c1, uint16_type /*c2*/, uint16_type q ) const
775  {
776  return fusion::accumulate( M_expr, value_type( 0 ),vf::detail::evaluate_expression_gi<value_type>( i, c1, q ) );
777  }
778  value_type
779  evalq( uint16_type c1, uint16_type /*c2*/, uint16_type q ) const
780  {
781  return fusion::accumulate( M_expr, value_type( 0 ),vf::detail::evaluate_expression_g<value_type>( c1, q ) );
782  }
783  tensor_vector_type M_expr;
784  };
785 
786 
787 
788 protected:
789 
790 private:
791  Vec();
792 
793  expression_vector_type M_expr;
794 };
795 } // detail
797 
801 template<typename Expr1>
802 inline
803 Expr<vf::detail::Vec<fusion::vector<Expr1> > >
804 vec( Expr1 expr1 )
805 {
806  typedef vf::detail::Vec<fusion::vector<Expr1> > expr_t;
807  return Expr<expr_t>( expr_t( fusion::vector<Expr1>( expr1 ) ) );
808 }
812 template<typename Expr1, typename Expr2>
813 inline
814 Expr<vf::detail::Vec<fusion::vector<Expr1, Expr2> > >
815 vec( Expr1 expr1, Expr2 expr2 )
816 {
817  typedef vf::detail::Vec<fusion::vector<Expr1, Expr2> > expr_t;
818  return Expr<expr_t>( expr_t( fusion::vector<Expr1, Expr2>( expr1, expr2 ) ) );
819 }
823 template<typename Expr1, typename Expr2, typename Expr3>
824 inline
825 Expr<vf::detail::Vec<fusion::vector<Expr1, Expr2, Expr3> > >
826 vec( Expr1 expr1, Expr2 expr2, Expr3 expr3 )
827 {
828  typedef vf::detail::Vec<fusion::vector<Expr1, Expr2, Expr3> > expr_t;
829  return Expr<expr_t>( expr_t( fusion::vector<Expr1, Expr2, Expr3>( expr1, expr2, expr3 ) ) );
830 }
831 
833 namespace detail
834 {
835 
843 template<int M, int N, typename MatrixExpr>
844 class Mat
845 {
846 public:
847 
848 
852 
853  static const size_type context = fusion::result_of::accumulate<MatrixExpr,mpl::size_t<0>,GetContext>::type::value;
854  //static const size_type context = fusion::accumulate( MatrixExpr(),size_type(0),GetContext() );
855  static const bool is_terminal = false;
856  static const uint16_type imorder = fusion::result_of::accumulate<MatrixExpr,mpl::size_t<0>,GetImOrder>::type::value;
857  static const bool imIsPoly = fusion::result_of::accumulate<MatrixExpr,mpl::bool_<true>,GetImIsPoly>::type::value;
858 
859  template<typename Func>
860  struct HasTestFunction
861  {
862  static const bool result = fusion::result_of::accumulate<MatrixExpr,mpl::bool_<false>,ExprHasTestFunction<Func> >::type::value;
863  };
864  template<typename Func>
865  struct HasTrialFunction
866  {
867  static const bool result = fusion::result_of::accumulate<MatrixExpr,mpl::bool_<false>,ExprHasTrialFunction<Func> >::type::value;
868  };
869 
870  typedef MatrixExpr expression_matrix_type;
871  typedef Mat<M, N, expression_matrix_type> this_type;
872 
873  typedef double value_type;
874 
875  static const uint16_type matrix_size1 = M;
876  static const uint16_type matrix_size2 = N;
877  static const uint16_type matrix_size = M*N;
878 
879 #if 0
880  BOOST_MPL_ASSERT_MSG( ( M*N == fusion::result_of::size<expression_matrix_type>::type::value ),
881  ( INVALID_MATRIX_SIZE ),
882  ( mpl::int_<M>, mpl::int_<N>, mpl::int_<M*N>,
883  mpl::int_<fusion::result_of::size<expression_matrix_type>::type::value> ) );
884 #endif
885 
886 
890 
891  Mat( MatrixExpr const& expr )
892  :
893  M_expr( expr )
894  {
895  }
896  Mat( Mat const & expr )
897  :
898  M_expr( expr.M_expr )
899  {
900  }
901  ~Mat()
902  {}
903 
905 
909 
910 
912 
916 
917  expression_matrix_type const& expression() const
918  {
919  return M_expr;
920  }
921  expression_matrix_type & expression()
922  {
923  return M_expr;
924  }
925 
927 
931 
932 
934 
938 
939 
941 
942  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t = Basis_i_t>
943  struct tensor
944  {
945  typedef this_type expression_type;
946  typedef typename expression_type::expression_matrix_type expression_matrix_type;
947  typedef Shape<expression_type::matrix_size1, Tensor2, false, false> shape;
948 
949 #if 0
950  static const bool theshape = ( shape::M == expression_type::matrix_size && shape::N == 1 );
951  BOOST_MPL_ASSERT_MSG( theshape,
952  INVALID_TENSOR_SHAPE_SHOULD_BE_RANK_1,
953  ( mpl::int_<shape::M>, mpl::int_<shape::N> ) );
954 #endif
955 
956  typedef typename boost::remove_reference<typename fusion::result_of::at<expression_matrix_type, boost::mpl::int_<0> >::type>::type::value_type value_type;
957  typedef typename boost::remove_reference<typename fusion::result_of::at<expression_matrix_type, boost::mpl::int_<0> >::type>::type::value_type first_value_type;
958  //typedef fusion::result_of::accumulate<VectorExpr,first_value_type,GetValueType>::type value_type;
959  //typedef typename expression_type::value_type value_type;
960 
961  template<typename ExprT>
962  struct ExpressionToTensor
963  {
964  typedef typename boost::remove_reference<ExprT>::type expr_type;
965  typedef typename expr_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> type;
966  };
967  typedef typename mpl::transform<expression_matrix_type, ExpressionToTensor<mpl::_1>, mpl::back_inserter<fusion::vector<> > >::type tensor_matrix_type;
968 
969 
970  struct is_zero
971  {
972  static const bool value = false;
973  };
974 
975  tensor( expression_type const& expr,
976  Geo_t const& geom,
977  Basis_i_t const& fev,
978  Basis_j_t const& feu )
979  :
980  M_expr( fusion::transform( expr.expression(), initialize_expression_gij<Geo_t,Basis_i_t,Basis_j_t>( geom, fev, feu ) ) )
981  {
982 
983  update( geom, fev, feu );
984  }
985  tensor( expression_type const& expr,
986  Geo_t const& geom,
987  Basis_i_t const& fev )
988  :
989  M_expr( fusion::transform( expr.expression(), initialize_expression_gi<Geo_t,Basis_i_t>( geom, fev ) ) )
990  {
991  update( geom, fev );
992  }
993  tensor( expression_type const& expr,
994  Geo_t const& geom )
995  :
996  M_expr( fusion::transform( expr.expression(), initialize_expression_g<Geo_t>( geom ) ) )
997  {
998  update( geom );
999  }
1000  template<typename IM>
1001  void init( IM const& im )
1002  {
1003  fusion::for_each( M_expr,vf::detail::init_expression<IM>( im ) );
1004  }
1005  void update( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
1006  {
1007  fusion::for_each( M_expr,vf::detail::update_expression_gij<Geo_t, Basis_i_t, Basis_j_t>( geom, fev, feu ) );
1008 
1009  }
1010  void update( Geo_t const& geom, Basis_i_t const& fev )
1011  {
1012  fusion::for_each( M_expr,vf::detail::update_expression_gi<Geo_t, Basis_i_t>( geom, fev ) );
1013  }
1014  void update( Geo_t const& geom )
1015  {
1016  fusion::for_each( M_expr,vf::detail::update_expression_g<Geo_t>( geom ) );
1017  }
1018  void update( Geo_t const& geom, uint16_type face )
1019  {
1020  fusion::for_each( M_expr,vf::detail::update_expression_face_g<Geo_t>( geom, face ) );
1021  }
1022 
1023  value_type
1024  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q ) const
1025  {
1026  return fusion::accumulate( M_expr, value_type( 0 ),
1027  vf::detail::evaluate_expression_gij<value_type>( expression_type::matrix_size2,
1028  i, j, c1, c2, q ) );
1029  }
1030  template<int PatternContext>
1031  value_type
1032  evalijq( uint16_type i, uint16_type j, uint16_type c1, uint16_type c2, uint16_type q,
1033  mpl::int_<PatternContext> ) const
1034  {
1035  return fusion::accumulate( M_expr, value_type( 0 ),
1036  vf::detail::evaluate_expression_gij<value_type>( expression_type::matrix_size2,
1037  i, j, c1, c2, q ) );
1038  }
1039 
1040  value_type
1041  evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q ) const
1042  {
1043  return fusion::accumulate( M_expr, value_type( 0 ),
1044  vf::detail::evaluate_expression_gi<value_type>( expression_type::matrix_size2, i, c1, c2, q ) );
1045 
1046  }
1047  value_type
1048  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const
1049  {
1050  return fusion::accumulate( M_expr, value_type( 0 ),vf::detail::evaluate_expression_g<value_type>( expression_type::matrix_size2,
1051  c1, c2, q ) );
1052  }
1053  tensor_matrix_type M_expr;
1054  };
1055 
1056 
1057 
1058 protected:
1059 
1060 private:
1061  Mat();
1062 
1063  expression_matrix_type M_expr;
1064 };
1065 } // detail
1067 
1071 template<int M, int N, typename Expr1>
1072 inline
1073 Expr<vf::detail::Mat<M, N, fusion::vector<Expr1> > >
1074 mat( Expr1 expr1 )
1075 {
1076  BOOST_MPL_ASSERT_MSG( ( M == 1 && N == 1 ), INVALID_MATRIX_SIZE_SHOULD_BE_1_1, ( mpl::int_<M>, mpl::int_<N> ) );
1077  typedef vf::detail::Mat<M, N, fusion::vector<Expr1> > expr_t;
1078  return Expr<expr_t>( expr_t( fusion::vector<Expr1>( expr1 ) ) );
1079 }
1080 
1084 template<int M, int N, typename Expr1, typename Expr2>
1085 inline
1086 Expr<vf::detail::Mat<M, N, fusion::vector<Expr1, Expr2> > >
1087 mat( Expr1 expr1, Expr2 expr2 )
1088 {
1089  BOOST_MPL_ASSERT_MSG( ( M == 2 && N == 1 || M == 1 && N == 2 ), INVALID_MATRIX_SIZE, ( mpl::int_<M>, mpl::int_<N> ) );
1090  typedef vf::detail::Mat<M,N,fusion::vector<Expr1, Expr2> > expr_t;
1091  return Expr<expr_t>( expr_t( fusion::vector<Expr1, Expr2>( expr1, expr2 ) ) );
1092 }
1093 
1097 template<int M, int N, typename Expr1, typename Expr2, typename Expr3>
1098 inline
1099 Expr<vf::detail::Mat<M, N, fusion::vector<Expr1, Expr2, Expr3> > >
1100 mat( Expr1 expr1, Expr2 expr2, Expr3 expr3 )
1101 {
1102 #if 0
1103  BOOST_MPL_ASSERT_MSG( ( M == 3 && N == 1 ||
1104  M == 1 && N == 3 ),
1105  ( INVALID_MATRIX_SIZE ),
1106  ( mpl::int_<M>, mpl::int_<N> ) );
1107 #endif
1108  typedef vf::detail::Mat<M, N, fusion::vector<Expr1, Expr2, Expr3> > expr_t;
1109  return Expr<expr_t>( expr_t( fusion::vector<Expr1, Expr2, Expr3>( expr1, expr2, expr3 ) ) );
1110 }
1111 
1115 template<int M, int N, typename Expr1, typename Expr2, typename Expr3, typename Expr4>
1116 inline
1117 Expr<vf::detail::Mat<M, N, fusion::vector<Expr1, Expr2, Expr3, Expr4> > >
1118 mat( Expr1 expr1, Expr2 expr2, Expr3 expr3, Expr4 expr4 )
1119 {
1120 #if 0
1121  BOOST_MPL_ASSERT_MSG( ( M == 4 && N == 1 ||
1122  M == 1 && N == 4 ||
1123  M == 2 && N == 2 ),
1124  ( INVALID_MATRIX_SIZE ),
1125  ( mpl::int_<M>, mpl::int_<N> ) );
1126 #endif
1127  typedef vf::detail::Mat<M, N, fusion::vector<Expr1, Expr2, Expr3, Expr4> > expr_t;
1128  return Expr<expr_t>( expr_t( fusion::vector<Expr1, Expr2, Expr3, Expr4>( expr1, expr2, expr3, expr4 ) ) );
1129 }
1130 
1134 template<int M, int N, typename Expr1, typename Expr2, typename Expr3, typename Expr4, typename Expr5>
1135 inline
1136 Expr<vf::detail::Mat<M, N, fusion::vector<Expr1, Expr2, Expr3, Expr4, Expr5> > >
1137 mat( Expr1 expr1, Expr2 expr2, Expr3 expr3, Expr4 expr4, Expr5 expr5 )
1138 {
1139 #if 0
1140  BOOST_MPL_ASSERT_MSG( ( M == 4 && N == 1 ||
1141  M == 1 && N == 4 ||
1142  M == 2 && N == 2 ),
1143  ( INVALID_MATRIX_SIZE ),
1144  ( mpl::int_<M>, mpl::int_<N> ) );
1145 #endif
1146  typedef vf::detail::Mat<M, N, fusion::vector<Expr1, Expr2, Expr3, Expr4, Expr5> > expr_t;
1147  return Expr<expr_t>( expr_t( fusion::vector<Expr1, Expr2, Expr3, Expr4, Expr5>( expr1, expr2, expr3, expr4, expr5 ) ) );
1148 }
1149 
1153 template<int M, int N, typename Expr1, typename Expr2, typename Expr3, typename Expr4, typename Expr5, typename Expr6>
1154 inline
1155 Expr<vf::detail::Mat<M, N, fusion::vector<Expr1, Expr2, Expr3, Expr4, Expr5, Expr6> > >
1156 mat( Expr1 expr1, Expr2 expr2, Expr3 expr3, Expr4 expr4, Expr5 expr5, Expr6 expr6 )
1157 {
1158 #if 0
1159  BOOST_MPL_ASSERT_MSG( ( M == 4 && N == 1 ||
1160  M == 1 && N == 4 ||
1161  M == 2 && N == 2 ),
1162  ( INVALID_MATRIX_SIZE ),
1163  ( mpl::int_<M>, mpl::int_<N> ) );
1164 #endif
1165  typedef vf::detail::Mat<M, N, fusion::vector<Expr1, Expr2, Expr3, Expr4, Expr5, Expr6> > expr_t;
1166  return Expr<expr_t>( expr_t( fusion::vector<Expr1, Expr2, Expr3, Expr4, Expr5, Expr6>( expr1, expr2, expr3, expr4, expr5, expr6 ) ) );
1167 }
1168 
1172 template<int M, int N,
1173  typename Expr1, typename Expr2, typename Expr3, typename Expr4,
1174  typename Expr5, typename Expr6, typename Expr7, typename Expr8,
1175  typename Expr9>
1176 inline
1177 Expr<vf::detail::Mat<M, N, fusion::vector<Expr1, Expr2, Expr3, Expr4, Expr5, Expr6, Expr7, Expr8, Expr9 > > >
1178 mat( Expr1 expr1, Expr2 expr2, Expr3 expr3, Expr4 expr4,
1179  Expr5 expr5, Expr6 expr6, Expr7 expr7, Expr8 expr8, Expr9 expr9 )
1180 {
1181 #if 0
1182  BOOST_MPL_ASSERT_MSG( ( M == 4 && N == 1 ||
1183  M == 1 && N == 4 ||
1184  M == 2 && N == 2 ),
1185  ( INVALID_MATRIX_SIZE ),
1186  ( mpl::int_<M>, mpl::int_<N> ) );
1187 #endif
1188  typedef vf::detail::Mat<M, N, fusion::vector<Expr1, Expr2, Expr3, Expr4,Expr5, Expr6, Expr7, Expr8, Expr9> > expr_t;
1189  return Expr<expr_t>( expr_t( fusion::vector<Expr1, Expr2, Expr3, Expr4, Expr5, Expr6, Expr7, Expr8, Expr9>( expr1, expr2, expr3, expr4, expr5, expr6, expr7, expr8, expr9 ) ) );
1190 }
1191 
1192 } // vf
1193 } // Feel
1194 #endif /* __VFVec_H */

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