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
products.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: 2011-08-04
7 
8  Copyright (C) 2011 Université 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 __Products_H
30 #define __Products_H 1
31 
32 #include <Eigen/Core>
33 
34 namespace Feel
35 {
36 namespace vf
37 {
39 
46 template<typename ExprL, typename ExprR, int Type = 1>
47 class Product
48 {
49 public:
50 
51  static const size_type context = ExprL::context | ExprR::context;
52  static const bool is_terminal = false;
53  static const int product_type = Type;
54 
55  static const uint16_type imorder = ExprL::imorder+ExprR::imorder;
56  static const bool imIsPoly = ExprL::imIsPoly && ExprR::imIsPoly;
57 
58  template<typename Func>
59  struct HasTestFunction
60  {
61  static const bool result =
62  ExprL::template HasTestFunction<Func>::result ||
63  ExprR::template HasTestFunction<Func>::result ;
64  };
65 
66  template<typename Func>
67  struct HasTrialFunction
68  {
69  static const bool result =
70  ExprL::template HasTrialFunction<Func>::result||
71  ExprR::template HasTrialFunction<Func>::result ;
72  };
73 
74 
78 
79  typedef ExprL left_expression_type;
80  typedef ExprR right_expression_type;
81  typedef typename left_expression_type::value_type value_type;
82  typedef Product<ExprL,ExprR,Type> this_type;
83 
84 
86 
90 
91  explicit Product( left_expression_type const & left_expr,
92  right_expression_type const & right_expr )
93  :
94  M_left_expr( left_expr ),
95  M_right_expr( right_expr )
96  {}
97  Product( Product const & te )
98  :
99  M_left_expr( te.M_left_expr ),
100  M_right_expr( te.M_right_expr )
101  {}
102  ~Product()
103  {}
104 
106 
110 
111 
113 
117 
118 
120 
124 
125 
127 
131 
132  left_expression_type const& left() const
133  {
134  return M_left_expr;
135  }
136  right_expression_type const& right() const
137  {
138  return M_right_expr;
139  }
140 
142 
143  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t>
144  struct tensor
145  {
146  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
147  typedef typename left_expression_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> l_tensor_expr_type;
148  typedef typename right_expression_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> r_tensor_expr_type;
149  typedef typename l_tensor_expr_type::value_type value_type;
150 
151  typedef typename l_tensor_expr_type::shape left_shape;
152  typedef typename r_tensor_expr_type::shape right_shape;
153  typedef Shape<left_shape::nDim,Scalar,false,false> shape;
154  static const bool l_is_terminal = left_expression_type::is_terminal;
155  static const bool r_is_terminal = right_expression_type::is_terminal;
156  template <class Args> struct sig
157  {
158  typedef value_type type;
159  };
160 
161  BOOST_MPL_ASSERT_MSG( (left_shape::M == right_shape::M) && (left_shape::N == right_shape::N) ,
162  INVALID_RANK_LEFT_AND_RIGHT_SHOULD_BE_THE_SAME,
163  (mpl::int_<left_shape::M>,mpl::int_<right_shape::M>,
164  mpl::int_<left_shape::N>,mpl::int_<right_shape::N>));
165 
166  struct is_zero
167  {
168  static const bool value = l_tensor_expr_type::is_zero::value || r_tensor_expr_type::is_zero::value;
169  };
170 
171  tensor( this_type const& expr,
172  Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
173  :
174  M_l_tensor_expr( expr.left(), geom, fev, feu ),
175  M_r_tensor_expr( expr.right(), geom, fev, feu )
176  {
177  }
178 
179  tensor( this_type const& expr,
180  Geo_t const& geom, Basis_i_t const& fev )
181  :
182  M_l_tensor_expr( expr.left(), geom, fev ),
183  M_r_tensor_expr( expr.right(), geom, fev )
184  {
185  }
186 
187  tensor( this_type const& expr, Geo_t const& geom )
188  :
189  M_l_tensor_expr( expr.left(), geom ),
190  M_r_tensor_expr( expr.right(), geom )
191  {
192  }
193 
194  template<typename IM>
195  void init( IM const& im )
196  {
197  M_l_tensor_expr.init( im );
198  M_r_tensor_expr.init( im );
199  }
200  void update( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
201  {
202  M_l_tensor_expr.update( geom, fev, feu );
203  M_r_tensor_expr.update( geom, fev, feu );
204  }
205  void update( Geo_t const& geom, Basis_i_t const& fev )
206  {
207  M_l_tensor_expr.update( geom, fev );
208  M_r_tensor_expr.update( geom, fev );
209  }
210  void update( Geo_t const& geom )
211  {
212  M_l_tensor_expr.update( geom );
213  M_r_tensor_expr.update( geom );
214  }
215  void update( Geo_t const& geom, uint16_type face )
216  {
217  M_l_tensor_expr.update( geom, face );
218  M_r_tensor_expr.update( geom, face );
219  }
220 
221  value_type
222  evalijq( uint16_type i, uint16_type j, uint16_type cc1, uint16_type cc2, uint16_type q ) const
223  {
224  return evalijq( i, j, cc1, cc2, q, typename mpl::and_<mpl::bool_<l_is_terminal>,mpl::bool_<l_is_terminal> >::type() );
225  }
226  value_type
227  evalijq( uint16_type i, uint16_type j, uint16_type cc1, uint16_type cc2, uint16_type q, mpl::bool_<0> ) const
228  {
229  double res = 0;
230 
231  if ( Type == 1 )
232  {
233  for ( int c2 = 0; c2 < left_shape::N; ++ c2 )
234  for ( int c1 = 0; c1 < left_shape::M; ++ c1 )
235  {
236  res += M_l_tensor_expr.evalijq( i, j, c1, c2, q )*M_r_tensor_expr.evalijq( i, j, c1, c2, q );
237  }
238  }
239 
240  return res;
241  }
242  value_type
243  evalijq( uint16_type i, uint16_type j, uint16_type cc1, uint16_type cc2, uint16_type q, mpl::bool_<1> ) const
244  {
245  //if ( Type == 1 )
246  {
247  return ( M_l_tensor_expr.evalijq( i,j,q ).adjoint()*M_r_tensor_expr.evalijq( i,j,q ) ).trace();
248  }
249  }
250  value_type
251  evaliq( uint16_type i, uint16_type cc1, uint16_type cc2, uint16_type q ) const
252  {
253  double res = 0;
254 
255  if ( Type == 1 )
256  {
257  for ( int c2 = 0; c2 < left_shape::N; ++ c2 )
258  for ( int c1 = 0; c1 < left_shape::M; ++ c1 )
259  {
260  res += M_l_tensor_expr.evaliq( i, c1, c2, q )*M_r_tensor_expr.evaliq( i, c1, c2, q );
261  }
262  }
263 
264  return res;
265  }
266  value_type
267  evalq( uint16_type cc1, uint16_type cc2, uint16_type q ) const
268  {
269  double res = 0;
270 
271  if ( Type == 1 )
272  {
273  for ( int c2 = 0; c2 < left_shape::N; ++ c2 )
274  for ( int c1 = 0; c1 < left_shape::M; ++ c1 )
275  {
276  res += M_l_tensor_expr.evalq( c1, c2, q )*M_r_tensor_expr.evalq( c1, c2, q );
277  }
278  }
279 
280  return res;
281  }
282 
283  private:
284  l_tensor_expr_type M_l_tensor_expr;
285  r_tensor_expr_type M_r_tensor_expr;
286  mutable Eigen::Matrix<value_type,left_shape::M,left_shape::N> M1,M2;
287  };
288 
289 private:
290  mutable left_expression_type M_left_expr;
291  mutable right_expression_type M_right_expr;
292 };
294 
298 template<typename ExprL, typename ExprR>
299 inline
300 Expr< Product<ExprL, ExprR,1> >
301 inner( ExprL l, ExprR r )
302 {
303  typedef Product<ExprL, ExprR,1> product_t;
304  return Expr< product_t >( product_t( l, r ) );
305 }
306 
308 
315 template<typename ExprL, typename ExprR>
316 class CrossProduct
317 {
318 public:
319 
320  static const size_type context = ExprL::context | ExprR::context;
321  static const bool is_terminal = false;
322 
323  static const uint16_type imorder = ExprL::imorder+ExprR::imorder;
324  static const bool imIsPoly = ExprL::imIsPoly && ExprR::imIsPoly;
325 
326  template<typename Func>
327  struct HasTestFunction
328  {
329  static const bool result =
330  ExprL::template HasTestFunction<Func>::result ||
331  ExprR::template HasTestFunction<Func>::result ;
332  };
333 
334  template<typename Func>
335  struct HasTrialFunction
336  {
337  static const bool result =
338  ExprL::template HasTrialFunction<Func>::result||
339  ExprR::template HasTrialFunction<Func>::result ;
340  };
341 
342 
346 
347  typedef ExprL left_expression_type;
348  typedef ExprR right_expression_type;
349  typedef typename left_expression_type::value_type value_type;
350  typedef CrossProduct<ExprL,ExprR> this_type;
351 
352 
354 
358 
359  explicit CrossProduct( left_expression_type const & left_expr,
360  right_expression_type const & right_expr )
361  :
362  M_left_expr( left_expr ),
363  M_right_expr( right_expr )
364  {}
365  CrossProduct( CrossProduct const & te )
366  :
367  M_left_expr( te.M_left_expr ),
368  M_right_expr( te.M_right_expr )
369  {}
370  ~CrossProduct()
371  {}
372 
374 
378 
379 
381 
385 
386 
388 
392 
393 
395 
399 
400  left_expression_type const& left() const
401  {
402  return M_left_expr;
403  }
404  right_expression_type const& right() const
405  {
406  return M_right_expr;
407  }
408 
410 
411  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t>
412  struct tensor
413  {
414  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
415  typedef typename left_expression_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> l_tensor_expr_type;
416  typedef typename right_expression_type::template tensor<Geo_t, Basis_i_t, Basis_j_t> r_tensor_expr_type;
417  typedef typename l_tensor_expr_type::value_type value_type;
418 
419  typedef typename l_tensor_expr_type::shape left_shape;
420  typedef typename r_tensor_expr_type::shape right_shape;
421  typedef typename mpl::if_<mpl::equal_to<mpl::int_<left_shape::nDim>,mpl::int_<2> >,
422  mpl::identity<Shape<2,Scalar,false,false> >,
423  mpl::identity<Shape<3,Vectorial,false,false> > >::type::type shape;
424  static const bool l_is_terminal = left_expression_type::is_terminal;
425  static const bool r_is_terminal = right_expression_type::is_terminal;
426 
427  BOOST_MPL_ASSERT_MSG( ( left_shape::nDim > 1 ),
428  INVALID_DIMENSION,
429  (mpl::int_<left_shape::nDim>,mpl::int_<right_shape::nDim>));
430  BOOST_MPL_ASSERT_MSG( left_shape::nDim == right_shape::nDim,
431  INVALID_DIMENSION_LEFT_AND_RIGHT_SHOULD_BE_THE_SAME,
432  (mpl::int_<left_shape::nDim>,mpl::int_<right_shape::nDim>));
433  BOOST_MPL_ASSERT_MSG( (left_shape::M == right_shape::M) && (left_shape::N == right_shape::N) ,
434  INVALID_RANK_LEFT_AND_RIGHT_SHOULD_BE_THE_SAME,
435  (mpl::int_<left_shape::M>,mpl::int_<right_shape::M>,
436  mpl::int_<left_shape::N>,mpl::int_<right_shape::N>));
437 
438 
439  template <class Args> struct sig
440  {
441  typedef value_type type;
442  };
443 
444  struct is_zero
445  {
446  static const bool value = l_tensor_expr_type::is_zero::value || r_tensor_expr_type::is_zero::value;
447  };
448 
449  tensor( this_type const& expr,
450  Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
451  :
452  M_l_tensor_expr( expr.left(), geom, fev, feu ),
453  M_r_tensor_expr( expr.right(), geom, fev, feu )
454  {
455  }
456 
457  tensor( this_type const& expr,
458  Geo_t const& geom, Basis_i_t const& fev )
459  :
460  M_l_tensor_expr( expr.left(), geom, fev ),
461  M_r_tensor_expr( expr.right(), geom, fev )
462  {
463  }
464 
465  tensor( this_type const& expr, Geo_t const& geom )
466  :
467  M_l_tensor_expr( expr.left(), geom ),
468  M_r_tensor_expr( expr.right(), geom )
469  {
470  }
471 
472  template<typename IM>
473  void init( IM const& im )
474  {
475  M_l_tensor_expr.init( im );
476  M_r_tensor_expr.init( im );
477  }
478  void update( Geo_t const& geom, Basis_i_t const& fev, Basis_j_t const& feu )
479  {
480  M_l_tensor_expr.update( geom, fev, feu );
481  M_r_tensor_expr.update( geom, fev, feu );
482  }
483  void update( Geo_t const& geom, Basis_i_t const& fev )
484  {
485  M_l_tensor_expr.update( geom, fev );
486  M_r_tensor_expr.update( geom, fev );
487  }
488  void update( Geo_t const& geom )
489  {
490  M_l_tensor_expr.update( geom );
491  M_r_tensor_expr.update( geom );
492  }
493  void update( Geo_t const& geom, uint16_type face )
494  {
495  M_l_tensor_expr.update( geom, face );
496  M_r_tensor_expr.update( geom, face );
497  }
498 
499  value_type
500  evalijq( uint16_type i, uint16_type j, uint16_type cc1, uint16_type cc2, uint16_type q ) const
501  {
502  return evalijq( i, j, cc1, cc2, q, mpl::int_<left_shape::nDim>() );
503  }
504  value_type
505  evalijq( uint16_type i, uint16_type j, uint16_type cc1, uint16_type cc2, uint16_type q, mpl::int_<2> ) const
506  {
507  double res = M_l_tensor_expr.evalijq( i, j, 0, 0, q )*M_r_tensor_expr.evalijq( i, j, 1, 0, q );
508  res -= M_l_tensor_expr.evalijq( i, j, 1, 0, q )*M_r_tensor_expr.evalijq( i, j, 0, 0, q );
509  return res;
510  }
511  value_type
512  evalijq( uint16_type i, uint16_type j, uint16_type cc1, uint16_type cc2, uint16_type q, mpl::int_<3> ) const
513  {
514  double res = M_l_tensor_expr.evalijq( i, j, (cc1+1)%3, 0, q )*M_r_tensor_expr.evalijq( i, j, (cc1+2)%3, 0, q );
515  res -= M_l_tensor_expr.evalijq( i, j, (cc1+2)%3, 0, q )*M_r_tensor_expr.evalijq( i, j, (cc1+1)%3, 0, q );
516  return res;
517  }
518  value_type
519  evaliq( uint16_type i, uint16_type cc1, uint16_type cc2, uint16_type q ) const
520  {
521  return evaliq( i, cc1, cc2, q, mpl::int_<left_shape::nDim>() );
522  }
523  value_type
524  evaliq( uint16_type i, uint16_type cc1, uint16_type cc2, uint16_type q, mpl::int_<2> ) const
525  {
526  double res = M_l_tensor_expr.evaliq( i, 0, 0, q )*M_r_tensor_expr.evaliq( i, 1, 0, q );
527  res -= M_l_tensor_expr.evaliq( i, 1, 0, q )*M_r_tensor_expr.evaliq( i, 0, 0, q );
528  return res;
529  }
530  value_type
531  evaliq( uint16_type i, uint16_type cc1, uint16_type cc2, uint16_type q, mpl::int_<3> ) const
532  {
533  double res = M_l_tensor_expr.evaliq( i, (cc1+1)%3, 0, q )*M_r_tensor_expr.evaliq( i, (cc1+2)%3, 0, q );
534  res -= M_l_tensor_expr.evaliq( i, (cc1+2)%3, 0, q )*M_r_tensor_expr.evaliq( i, (cc1+1)%3, 0, q );
535  return res;
536  }
537  value_type
538  evalq( uint16_type cc1, uint16_type cc2, uint16_type q ) const
539  {
540  return evalq( cc1, cc2, q, mpl::int_<left_shape::nDim>() );
541  }
542  value_type
543  evalq( uint16_type cc1, uint16_type cc2, uint16_type q, mpl::int_<2> ) const
544  {
545  double res = M_l_tensor_expr.evalq( 0, 0, q )*M_r_tensor_expr.evalq( 1, 0, q );
546  res -= M_l_tensor_expr.evalq( 1, 0, q )*M_r_tensor_expr.evalq( 0, 0, q );
547  return res;
548  }
549  value_type
550  evalq( uint16_type cc1, uint16_type cc2, uint16_type q, mpl::int_<3> ) const
551  {
552  double res = M_l_tensor_expr.evalq( (cc1+1)%3, 0, q )*M_r_tensor_expr.evalq( (cc1+2)%3, 0, q );
553  res -= M_l_tensor_expr.evalq( (cc1+2)%3, 0, q )*M_r_tensor_expr.evalq( (cc1+1)%3, 0, q );
554  return res;
555  }
556 
557  private:
558  l_tensor_expr_type M_l_tensor_expr;
559  r_tensor_expr_type M_r_tensor_expr;
560  };
561 
562 private:
563  mutable left_expression_type M_left_expr;
564  mutable right_expression_type M_right_expr;
565 };
567 
571 template<typename ExprL, typename ExprR>
572 inline
573 Expr< CrossProduct<ExprL, ExprR> >
574 cross( ExprL l, ExprR r )
575 {
576  typedef CrossProduct<ExprL, ExprR> product_t;
577  return Expr< product_t >( product_t( l, r ) );
578 }
579 
580 } // vf
581 
582 
583 } // Feel
584 #endif /* __Products_H */
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319

Generated on Sun Dec 22 2013 13:11:11 for Feel++ by doxygen 1.8.5