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
ginac.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 -*-
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <prudhomme@unistra.fr>
6  Date: 2012-10-15
7 
8  Copyright (C) 2012 Université de Strasbourg
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 2.1 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 __Ginac_H
30 #define __Ginac_H 1
31 
32 #include <ginac/ginac.h>
33 #include <boost/parameter/preprocessor.hpp>
34 
35 namespace GiNaC
36 {
37 
38  matrix grad( ex const& f, std::vector<symbol> const& l );
39  ex laplacian( ex const& f, std::vector<symbol> const& l );
40  matrix grad( std::string const& s, std::vector<symbol> const& l );
41  ex laplacian( std::string const& s, std::vector<symbol> const& l );
42 
43  matrix grad( matrix const& f, std::vector<symbol> const& l );
44  matrix div( matrix const& f, std::vector<symbol> const& l );
45  matrix laplacian( matrix const& f, std::vector<symbol> const& l );
46 
47  ex diff(ex const& f, symbol const& l, const int n);
48  matrix diff(matrix const& f, symbol const& l, const int n);
49 
50  ex substitute(ex const& f, symbol const& l, const double val );
51  ex substitute(ex const& f, symbol const& l, ex const & g );
52 
53  matrix substitute(matrix const& f, symbol const& l, const double val );
54  matrix substitute(matrix const& f, symbol const& l, ex const & g );
55 
56  //ex parse( std::string const& str, std::vector<symbol> const& syms );
57  ex parse( std::string const& str, std::vector<symbol> const& syms, std::vector<symbol> const& params = std::vector<symbol>());
58 
59 } // GiNaC
60 
61 namespace Feel
62 {
63  using GiNaC::matrix;
64  using GiNaC::symbol;
65  using GiNaC::lst;
66  using GiNaC::ex;
67  using GiNaC::parser;
68 
69  template<int Dim> inline std::vector<symbol> symbols() { return {symbol("x")}; }
70  template<> inline std::vector<symbol> symbols<1>() { return {symbol("x")}; }
71  template<> inline std::vector<symbol> symbols<2>() { return {symbol("x"),symbol("y") };}
72  template<> inline std::vector<symbol> symbols<3>() { return {symbol("x"),symbol("y"),symbol("z") };}
73 
74  inline
75  std::vector<symbol>
76  symbols( std::initializer_list<std::string> l )
77  {
78  std::vector<symbol> s;
79  std::for_each( l.begin(), l.end(), [&s] ( std::string const& sym ) { s.push_back( symbol(sym) ); } );
80  return s;
81  }
82  inline
83  std::vector<symbol>
84  symbols( std::vector<std::string> l )
85  {
86  std::vector<symbol> s;
87  std::for_each( l.begin(), l.end(), [&s] ( std::string const& sym ) { s.push_back( symbol(sym) ); } );
88  return s;
89  }
90 
91  namespace vf
92  {
93 
95 
102  template<int Order = 2>
103  class GinacEx
104  {
105  public:
106 
107 
111 
112 
113  static const size_type context = vm::POINT;
114  static const bool is_terminal = false;
115  static const uint16_type imorder = Order;
116  static const bool imIsPoly = false;
117 
118  template<typename Funct>
119  struct HasTestFunction
120  {
121  static const bool result = false;
122  };
123  template<typename Funct>
124  struct HasTrialFunction
125  {
126  static const bool result = false;
127  };
128 
129  typedef GiNaC::ex expression_type;
130  typedef GinacEx<Order> this_type;
131  typedef double value_type;
132 
133  template<typename TheExpr>
134  struct Lambda
135  {
136  typedef this_type type;
137  };
138  template<typename TheExpr>
139  typename Lambda<TheExpr>::type
140  operator()( TheExpr const& e ) { return *this; }
141 
143 
147 
148  explicit GinacEx( expression_type const & fun, std::vector<GiNaC::symbol> const& syms, std::string filename="")
149  :
150  M_fun( fun ),
151  M_syms( syms),
152  M_cfun(),
153  M_filename(filename.empty()?filename:(fs::current_path()/filename).string())
154  {
155  DVLOG(2) << "Ginac constructor with expression_type \n";
156  GiNaC::lst exprs(fun);
157  GiNaC::lst syml;
158  std::for_each( M_syms.begin(),M_syms.end(), [&]( GiNaC::symbol const& s ) { syml.append(s); } );
159 
160  // If the so file already exists, no need to re-compile but only link it
161  std::string filenameWithSuffix = M_filename + ".so";
162  if( !filename.empty() && fs::exists( filenameWithSuffix ) )
163  GiNaC::link_ex(filenameWithSuffix, M_cfun);
164  else
165  GiNaC::compile_ex(exprs, syml, M_cfun, M_filename);
166  }
167 
168  GinacEx( GinacEx const & fun )
169  :
170  M_fun( fun.M_fun ),
171  M_syms( fun.M_syms),
172  M_cfun(),
173  M_filename( fun.M_filename )
174  {
175  if( !(M_fun==fun.M_fun && M_syms==fun.M_syms && M_filename==fun.M_filename) || M_filename.empty() )
176  {
177  DVLOG(2) << "Ginac copy constructor : compile object file \n";
178  GiNaC::lst exprs(M_fun);
179  GiNaC::lst syml;
180  std::for_each( M_syms.begin(),M_syms.end(), [&]( GiNaC::symbol const& s ) { syml.append(s); } );
181  GiNaC::compile_ex(exprs, syml, M_cfun, M_filename);
182  }
183  else
184  {
185  DVLOG(2) << "Ginac copy constructor : link with existing object file \n";
186  boost::mpi::communicator world;
187  // std::string pid = boost::lexical_cast<std::string>(world.rank());
188  // std::string filenameWithSuffix = M_filename + pid + ".so";
189  std::string filenameWithSuffix = M_filename + ".so";
190  GiNaC::link_ex(filenameWithSuffix, M_cfun);
191  }
192  }
193 
195 
199 
200 
202 
206 
207 
209 
213 
214 
216 
220 
221  const GiNaC::FUNCP_CUBA& fun() const
222  {
223  return M_cfun;
224  }
225 
226  std::vector<GiNaC::symbol> const& syms() const { return M_syms; }
227 
229 
230 
231  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t>
232  struct tensor
233  {
234  //typedef typename expression_type::value_type value_type;
235  typedef double value_type;
236 
237  typedef typename mpl::if_<fusion::result_of::has_key<Geo_t,vf::detail::gmc<0> >,
238  mpl::identity<vf::detail::gmc<0> >,
239  mpl::identity<vf::detail::gmc<1> > >::type::type key_type;
240  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type* gmc_ptrtype;
241  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type gmc_type;
242  // change 0 into rank
243  typedef typename mpl::if_<mpl::equal_to<mpl::int_<0>,mpl::int_<0> >,
244  mpl::identity<Shape<gmc_type::nDim, Scalar, false, false> >,
245  typename mpl::if_<mpl::equal_to<mpl::int_<0>,mpl::int_<1> >,
246  mpl::identity<Shape<gmc_type::nDim, Vectorial, false, false> >,
247  mpl::identity<Shape<gmc_type::nDim, Tensor2, false, false> > >::type >::type::type shape;
248 
249  typedef Eigen::Matrix<value_type,Eigen::Dynamic,1> vec_type;
250 
251  struct is_zero
252  {
253  static const bool value = false;
254  };
255 
256  tensor( this_type const& expr,
257  Geo_t const& geom, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
258  :
259  M_fun( expr.fun() ),
260  M_gmc( fusion::at_key<key_type>( geom ).get() ),
261  M_nsyms( expr.syms().size() ),
262  M_y( vec_type::Zero(M_gmc->nPoints()) ),
263  M_x( vec_type::Zero( M_nsyms ) )
264  {}
265 
266  tensor( this_type const& expr,
267  Geo_t const& geom, Basis_i_t const& /*fev*/ )
268  :
269  M_fun( expr.fun() ),
270  M_gmc( fusion::at_key<key_type>( geom ).get() ),
271  M_nsyms( expr.syms().size() ),
272  M_y( vec_type::Zero(M_gmc->nPoints()) ),
273  M_x( vec_type::Zero( M_nsyms ) )
274  {}
275 
276  tensor( this_type const& expr, Geo_t const& geom )
277  :
278  M_fun( expr.fun() ),
279  M_gmc( fusion::at_key<key_type>( geom ).get() ),
280  M_nsyms( expr.syms().size() ),
281  M_y( vec_type::Zero(M_gmc->nPoints()) ),
282  M_x( vec_type::Zero( M_nsyms ) )
283 
284  {
285  }
286 
287  template<typename IM>
288  void init( IM const& im )
289  {
290 
291  }
292  void update( Geo_t const& geom, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
293  {
294  update( geom );
295  }
296  void update( Geo_t const& geom, Basis_i_t const& /*fev*/ )
297  {
298  update( geom );
299  }
300  void update( Geo_t const& geom )
301  {
302  M_gmc = fusion::at_key<key_type>( geom ).get();
303 
304  int no = 1;
305  int ni = M_nsyms;
306 
307  for(int q = 0; q < M_gmc->nPoints();++q )
308  {
309  for(int k = 0;k < gmc_type::nDim;++k )
310  M_x[k]=M_gmc->xReal( q )[k];
311  for( int k = gmc_type::nDim; k < M_x.size(); ++k )
312  M_x[k] = 0;
313  M_fun(&ni,M_x.data(),&no,&M_y[q]);
314  }
315 
316  }
317 
318  void update( Geo_t const& geom, uint16_type /*face*/ )
319  {
320  M_gmc = fusion::at_key<key_type>( geom ).get();
321 
322  int no = 1;
323  int ni = M_nsyms;//gmc_type::nDim;
324  for(int q = 0; q < M_gmc->nPoints();++q )
325  {
326  for(int k = 0;k < gmc_type::nDim;++k )
327  M_x[k]=M_gmc->xReal( q )[k];
328  for( int k = gmc_type::nDim; k < M_x.size(); ++k )
329  M_x[k] = 0;
330  M_fun(&ni,M_x.data(),&no,&M_y[q]);
331  }
332  }
333 
334 
335  value_type
336  evalij( uint16_type i, uint16_type j ) const
337  {
338  return 0;
339  }
340 
341 
342  value_type
343  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type c1, uint16_type c2, uint16_type q ) const
344  {
345  return M_y[q];
346  }
347 
348 
349 
350  value_type
351  evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q ) const
352  {
353  return M_y[q];
354  }
355 
356  value_type
357  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const
358  {
359  return M_y[q];
360  }
361 
362  GiNaC::FUNCP_CUBA M_fun;
363  gmc_ptrtype M_gmc;
364 
365  int M_nsyms;
366  vec_type M_y;
367  vec_type M_x;
368 
369  };
370 
371  private:
372  mutable expression_type M_fun;
373  std::vector<GiNaC::symbol> M_syms;
374  GiNaC::FUNCP_CUBA M_cfun;
375  std::string M_filename;
376  };
377 
378  inline
379  Expr< GinacEx<2> >
380  expr( GiNaC::ex const& f, std::vector<GiNaC::symbol> const& lsym, std::string filename="" )
381  {
382  return Expr< GinacEx<2> >( GinacEx<2>( f, lsym, filename ) );
383  }
384 
385  inline
386  Expr< GinacEx<2> >
387  expr( std::string const& s, std::vector<GiNaC::symbol> const& lsym, std::string filename="" )
388  {
389  return Expr< GinacEx<2> >( GinacEx<2>( parse(s,lsym), lsym, filename) );
390  }
391 
396  template<int Order>
397  inline
398  Expr< GinacEx<Order> >
399  expr( GiNaC::ex const& f, std::vector<GiNaC::symbol> const& lsym, std::string filename="" )
400  {
401  return Expr< GinacEx<Order> >( GinacEx<Order>( f, lsym, filename ));
402  }
403 
404  template<int Order>
405  inline
406  Expr< GinacEx<Order> >
407  expr( std::string const& s, std::vector<GiNaC::symbol> const& lsym, std::string filename="" )
408  {
409  return Expr< GinacEx<Order> >( GinacEx<Order>( parse(s,lsym), lsym, filename) );
410  }
411 
412  template<int M=1, int N=1, int Order = 2>
413  class GinacMatrix
414  {
415  public:
416 
417 
421 
422 
423  static const size_type context = vm::POINT;
424  static const bool is_terminal = false;
425  static const uint16_type imorder = Order;
426  static const bool imIsPoly = false;
427 
428  template<typename Funct>
429  struct HasTestFunction
430  {
431  static const bool result = false;
432  };
433  template<typename Funct>
434  struct HasTrialFunction
435  {
436  static const bool result = false;
437  };
438 
439  typedef GiNaC::ex expression_type;
440  typedef GinacMatrix<M,N,Order> this_type;
441  typedef double value_type;
442 
443  template<typename TheExpr>
444  struct Lambda
445  {
446  typedef this_type type;
447  };
448  template<typename TheExpr>
449  typename Lambda<TheExpr>::type
450  operator()( TheExpr const& e ) { return *this; }
451 
453 
457 
458  explicit GinacMatrix( GiNaC::matrix const & fun, std::vector<GiNaC::symbol> const& syms, std::string filename="" )
459  :
460  M_fun( fun.evalm() ),
461  M_syms( syms),
462  M_cfun(),
463  M_filename(filename.empty()?filename:(fs::current_path()/filename).string())
464  {
465  GiNaC::lst exprs;
466  for( int i = 0; i < M_fun.nops(); ++i ) exprs.append( M_fun.op(i) );
467 
468  GiNaC::lst syml;
469  std::for_each( M_syms.begin(),M_syms.end(), [&]( GiNaC::symbol const& s ) { syml.append(s); } );
470  GiNaC::compile_ex(exprs, syml, M_cfun, M_filename);
471  }
472  explicit GinacMatrix( GiNaC::ex const & fun, std::vector<GiNaC::symbol> const& syms, std::string filename="" )
473  :
474  M_fun(fun.evalm()),
475  M_syms( syms),
476  M_cfun(),
477  M_filename(filename.empty()?filename:(fs::current_path()/filename).string())
478  {
479  GiNaC::lst exprs;
480  for( int i = 0; i < M_fun.nops(); ++i ) exprs.append( M_fun.op(i) );
481 
482  GiNaC::lst syml;
483  std::for_each( M_syms.begin(),M_syms.end(), [&]( GiNaC::symbol const& s ) { syml.append(s); } );
484  GiNaC::compile_ex(exprs, syml, M_cfun, M_filename);
485  }
486 
487  GinacMatrix( GinacMatrix const & fun )
488  :
489  M_fun( fun.M_fun ),
490  M_syms( fun.M_syms),
491  M_cfun(),
492  M_filename( fun.M_filename )
493  {
494  if( !(M_fun==fun.M_fun && M_syms==fun.M_syms && M_filename==fun.M_filename) || M_filename.empty() )
495  {
496  DVLOG(2) << "Ginac copy constructor : compile object file \n";
497  GiNaC::lst exprs;
498  for( int i = 0; i < fun.M_fun.nops(); ++i ) exprs.append( fun.M_fun.op(i) );
499 
500  GiNaC::lst syml;
501  std::for_each( M_syms.begin(),M_syms.end(), [&]( GiNaC::symbol const& s ) { syml.append(s); } );
502  GiNaC::compile_ex(exprs, syml, M_cfun, M_filename);
503  }
504  else
505  {
506  DVLOG(2) << "Ginac copy constructor : link with existing object file \n";
507  boost::mpi::communicator world;
508  //std::string pid = boost::lexical_cast<std::string>(world.rank());
509  std::string filenameWithSuffix = M_filename + ".so";
510  GiNaC::link_ex(filenameWithSuffix, M_cfun);
511  }
512  }
513 
514 
516 
520 
521 
523 
527 
528 
530 
534 
535 
537 
541 
542  const GiNaC::FUNCP_CUBA& fun() const
543  {
544  return M_cfun;
545  }
546 
547  std::vector<GiNaC::symbol> const& syms() const { return M_syms; }
549 
550 
551  template<typename Geo_t, typename Basis_i_t, typename Basis_j_t>
552  struct tensor
553  {
554  //typedef typename expression_type::value_type value_type;
555  typedef double value_type;
556 
557  typedef typename mpl::if_<fusion::result_of::has_key<Geo_t,vf::detail::gmc<0> >,
558  mpl::identity<vf::detail::gmc<0> >,
559  mpl::identity<vf::detail::gmc<1> > >::type::type key_type;
560  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type* gmc_ptrtype;
561  typedef typename fusion::result_of::value_at_key<Geo_t,key_type>::type::element_type gmc_type;
562 
563  typedef typename mn_to_shape<gmc_type::nDim,M,N>::type shape;
564  // be careful that the matrix passed to ginac must be Row Major,
565  // however if the number of columns is 1 then eigen3 fails with
566  // an assertion, so we have a special when N=1 and have the
567  // matrix column major which is ok in this case
568  typedef typename mpl::if_<mpl::equal_to<mpl::int_<shape::N>, mpl::int_<1>>,
569  mpl::identity<Eigen::Matrix<value_type,shape::M,1>>,
570  mpl::identity<Eigen::Matrix<value_type,shape::M,shape::N,Eigen::RowMajor>>>::type::type mat_type;
571  typedef std::vector<mat_type> loc_type;
572  typedef Eigen::Matrix<value_type,Eigen::Dynamic,1> vec_type;
573  struct is_zero
574  {
575  static const bool value = false;
576  };
577 
578  tensor( this_type const& expr,
579  Geo_t const& geom, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
580  :
581  M_fun( expr.fun() ),
582  M_gmc( fusion::at_key<key_type>( geom ).get() ),
583  M_nsyms( expr.syms().size() ),
584  M_y( M_gmc->nPoints(), mat_type::Zero() ),
585  M_x( vec_type::Zero(M_nsyms) )
586  {}
587 
588  tensor( this_type const& expr,
589  Geo_t const& geom, Basis_i_t const& /*fev*/ )
590  :
591  M_fun( expr.fun() ),
592  M_gmc( fusion::at_key<key_type>( geom ).get() ),
593  M_nsyms( expr.syms().size() ),
594  M_y( M_gmc->nPoints(), mat_type::Zero() ),
595  M_x( vec_type::Zero(M_nsyms) )
596 
597  {}
598 
599  tensor( this_type const& expr, Geo_t const& geom )
600  :
601  M_fun( expr.fun() ),
602  M_gmc( fusion::at_key<key_type>( geom ).get() ),
603  M_nsyms( expr.syms().size() ),
604  M_y( M_gmc->nPoints(), mat_type::Zero() ),
605  M_x( vec_type::Zero(M_nsyms) )
606  {
607  }
608 
609  template<typename IM>
610  void init( IM const& im )
611  {
612 
613  }
614  void update( Geo_t const& geom, Basis_i_t const& /*fev*/, Basis_j_t const& /*feu*/ )
615  {
616  update( geom );
617  }
618  void update( Geo_t const& geom, Basis_i_t const& /*fev*/ )
619  {
620  update( geom );
621  }
622  void update( Geo_t const& geom )
623  {
624  M_gmc = fusion::at_key<key_type>( geom ).get();
625 
626  int no = M*N;
627  int ni = M_nsyms;//gmc_type::nDim;
628  for(int q = 0; q < M_gmc->nPoints();++q )
629  {
630  for(int k = 0;k < gmc_type::nDim;++k )
631  M_x[k]=M_gmc->xReal( q )[k];
632  for( int k = gmc_type::nDim; k < M_x.size(); ++k )
633  M_x[k] = 0;
634  M_fun(&ni,M_x.data(),&no,M_y[q].data());
635 ;
636  }
637 
638  }
639 
640  void update( Geo_t const& geom, uint16_type /*face*/ )
641  {
642  M_gmc = fusion::at_key<key_type>( geom ).get();
643 
644  int no = M*N;
645  int ni = M_nsyms;//gmc_type::nDim;
646  for(int q = 0; q < M_gmc->nPoints();++q )
647  {
648  for(int k = 0;k < gmc_type::nDim;++k )
649  M_x[k]=M_gmc->xReal( q )[k];
650  for( int k = gmc_type::nDim; k < M_x.size(); ++k )
651  M_x[k] = 0;
652  M_fun(&ni,M_x.data(),&no,M_y[q].data());
653  }
654  }
655 
656 
657  value_type
658  evalij( uint16_type i, uint16_type j ) const
659  {
660  return 0;
661  }
662 
663 
664  value_type
665  evalijq( uint16_type /*i*/, uint16_type /*j*/, uint16_type c1, uint16_type c2, uint16_type q ) const
666  {
667  return M_y[q](c1,c2);
668  }
669 
670  value_type
671  evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q ) const
672  {
673  return M_y[q](c1,c2);
674  }
675 
676  value_type
677  evalq( uint16_type c1, uint16_type c2, uint16_type q ) const
678  {
679  return M_y[q](c1,c2);
680  }
681 
682  GiNaC::FUNCP_CUBA M_fun;
683  gmc_ptrtype M_gmc;
684  int M_nsyms;
685  loc_type M_y;
686  vec_type M_x;
687  };
688 
689  private:
690  mutable expression_type M_fun;
691  std::vector<GiNaC::symbol> M_syms;
692  GiNaC::FUNCP_CUBA M_cfun;
693  std::string M_filename;
694  }; // GinacMatrix
696 
697  inline
698  Expr< GinacMatrix<1,1,2> >
699  expr( GiNaC::matrix const& f, std::vector<GiNaC::symbol> const& lsym, std::string filename="")
700  {
701  return Expr< GinacMatrix<1,1,2> >( GinacMatrix<1,1,2>( f, lsym, filename ) );
702  }
703 
708  template<int M, int N, int Order>
709  inline
710  Expr< GinacMatrix<M,N,Order> >
711  expr( GiNaC::matrix const& f, std::vector<GiNaC::symbol> const& lsym, std::string filename="" )
712  {
713  return Expr< GinacMatrix<M,N,Order> >( GinacMatrix<M,N,Order>( f, lsym, filename) );
714  }
715 
716  template<int M, int N, int Order>
717  inline
718  Expr< GinacMatrix<M,N,Order> >
719  expr( GiNaC::ex const& f, std::vector<GiNaC::symbol> const& lsym, std::string filename="" )
720  {
721  return Expr< GinacMatrix<M,N,Order> >( GinacMatrix<M,N,Order>( f, lsym, filename ) );
722  }
723 
724 
725  } // vf
726 } // Feel
727 #endif /* __Ginac_H */

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