32 #include <ginac/ginac.h>
33 #include <boost/parameter/preprocessor.hpp>
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 );
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 );
47 ex diff(ex
const& f, symbol
const& l,
const int n);
48 matrix diff(matrix
const& f, symbol
const& l,
const int n);
50 ex substitute(ex
const& f, symbol
const& l,
const double val );
51 ex substitute(ex
const& f, symbol
const& l, ex
const & g );
53 matrix substitute(matrix
const& f, symbol
const& l,
const double val );
54 matrix substitute(matrix
const& f, symbol
const& l, ex
const & g );
57 ex parse( std::string
const& str, std::vector<symbol>
const& syms, std::vector<symbol>
const& params = std::vector<symbol>());
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") };}
76 symbols( std::initializer_list<std::string> l )
78 std::vector<symbol> s;
79 std::for_each( l.begin(), l.end(), [&s] ( std::string
const& sym ) { s.push_back( symbol(sym) ); } );
84 symbols( std::vector<std::string> l )
86 std::vector<symbol> s;
87 std::for_each( l.begin(), l.end(), [&s] ( std::string
const& sym ) { s.push_back( symbol(sym) ); } );
102 template<
int Order = 2>
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;
118 template<
typename Funct>
119 struct HasTestFunction
121 static const bool result =
false;
123 template<
typename Funct>
124 struct HasTrialFunction
126 static const bool result =
false;
129 typedef GiNaC::ex expression_type;
130 typedef GinacEx<Order> this_type;
131 typedef double value_type;
133 template<
typename TheExpr>
136 typedef this_type type;
138 template<
typename TheExpr>
139 typename Lambda<TheExpr>::type
140 operator()( TheExpr
const& e ) {
return *
this; }
148 explicit GinacEx( expression_type
const & fun, std::vector<GiNaC::symbol>
const& syms, std::string filename=
"")
153 M_filename(filename.empty()?filename:(fs::current_path()/filename).string())
155 DVLOG(2) <<
"Ginac constructor with expression_type \n";
156 GiNaC::lst exprs(fun);
158 std::for_each( M_syms.begin(),M_syms.end(), [&]( GiNaC::symbol
const& s ) { syml.append(s); } );
161 std::string filenameWithSuffix = M_filename +
".so";
162 if( !filename.empty() && fs::exists( filenameWithSuffix ) )
163 GiNaC::link_ex(filenameWithSuffix, M_cfun);
165 GiNaC::compile_ex(exprs, syml, M_cfun, M_filename);
168 GinacEx( GinacEx
const & fun )
173 M_filename( fun.M_filename )
175 if( !(M_fun==fun.M_fun && M_syms==fun.M_syms && M_filename==fun.M_filename) || M_filename.empty() )
177 DVLOG(2) <<
"Ginac copy constructor : compile object file \n";
178 GiNaC::lst exprs(M_fun);
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);
185 DVLOG(2) <<
"Ginac copy constructor : link with existing object file \n";
186 boost::mpi::communicator world;
189 std::string filenameWithSuffix = M_filename +
".so";
190 GiNaC::link_ex(filenameWithSuffix, M_cfun);
221 const GiNaC::FUNCP_CUBA& fun()
const
226 std::vector<GiNaC::symbol>
const& syms()
const {
return M_syms; }
231 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t>
235 typedef double value_type;
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;
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;
249 typedef Eigen::Matrix<value_type,Eigen::Dynamic,1> vec_type;
253 static const bool value =
false;
256 tensor( this_type
const& expr,
257 Geo_t
const& geom, Basis_i_t
const& , Basis_j_t
const& )
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 ) )
266 tensor( this_type
const& expr,
267 Geo_t
const& geom, Basis_i_t
const& )
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 ) )
276 tensor( this_type
const& expr, Geo_t
const& geom )
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 ) )
287 template<
typename IM>
288 void init( IM
const& im )
292 void update( Geo_t
const& geom, Basis_i_t
const& , Basis_j_t
const& )
296 void update( Geo_t
const& geom, Basis_i_t
const& )
300 void update( Geo_t
const& geom )
302 M_gmc = fusion::at_key<key_type>( geom ).
get();
307 for(
int q = 0; q < M_gmc->nPoints();++q )
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 )
313 M_fun(&ni,M_x.data(),&no,&M_y[q]);
318 void update( Geo_t
const& geom, uint16_type )
320 M_gmc = fusion::at_key<key_type>( geom ).
get();
324 for(
int q = 0; q < M_gmc->nPoints();++q )
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 )
330 M_fun(&ni,M_x.data(),&no,&M_y[q]);
336 evalij( uint16_type i, uint16_type j )
const
343 evalijq( uint16_type , uint16_type , uint16_type c1, uint16_type c2, uint16_type q )
const
351 evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q )
const
357 evalq( uint16_type c1, uint16_type c2, uint16_type q )
const
362 GiNaC::FUNCP_CUBA M_fun;
372 mutable expression_type M_fun;
373 std::vector<GiNaC::symbol> M_syms;
374 GiNaC::FUNCP_CUBA M_cfun;
375 std::string M_filename;
380 expr( GiNaC::ex
const& f, std::vector<GiNaC::symbol>
const& lsym, std::string filename=
"" )
382 return Expr< GinacEx<2> >( GinacEx<2>( f, lsym, filename ) );
387 expr( std::string
const& s, std::vector<GiNaC::symbol>
const& lsym, std::string filename=
"" )
389 return Expr< GinacEx<2> >( GinacEx<2>( parse(s,lsym), lsym, filename) );
398 Expr< GinacEx<Order> >
399 expr( GiNaC::ex
const& f, std::vector<GiNaC::symbol>
const& lsym, std::string filename=
"" )
401 return Expr< GinacEx<Order> >( GinacEx<Order>( f, lsym, filename ));
406 Expr< GinacEx<Order> >
407 expr( std::string
const& s, std::vector<GiNaC::symbol>
const& lsym, std::string filename=
"" )
409 return Expr< GinacEx<Order> >( GinacEx<Order>( parse(s,lsym), lsym, filename) );
412 template<
int M=1,
int N=1,
int Order = 2>
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;
428 template<
typename Funct>
429 struct HasTestFunction
431 static const bool result =
false;
433 template<
typename Funct>
434 struct HasTrialFunction
436 static const bool result =
false;
439 typedef GiNaC::ex expression_type;
440 typedef GinacMatrix<M,N,Order> this_type;
441 typedef double value_type;
443 template<
typename TheExpr>
446 typedef this_type type;
448 template<
typename TheExpr>
449 typename Lambda<TheExpr>::type
450 operator()( TheExpr
const& e ) {
return *
this; }
458 explicit GinacMatrix( GiNaC::matrix
const & fun, std::vector<GiNaC::symbol>
const& syms, std::string filename=
"" )
460 M_fun( fun.evalm() ),
463 M_filename(filename.empty()?filename:(fs::current_path()/filename).string())
466 for(
int i = 0; i < M_fun.nops(); ++i ) exprs.append( M_fun.op(i) );
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);
472 explicit GinacMatrix( GiNaC::ex
const & fun, std::vector<GiNaC::symbol>
const& syms, std::string filename=
"" )
477 M_filename(filename.empty()?filename:(fs::current_path()/filename).string())
480 for(
int i = 0; i < M_fun.nops(); ++i ) exprs.append( M_fun.op(i) );
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);
487 GinacMatrix( GinacMatrix
const & fun )
492 M_filename( fun.M_filename )
494 if( !(M_fun==fun.M_fun && M_syms==fun.M_syms && M_filename==fun.M_filename) || M_filename.empty() )
496 DVLOG(2) <<
"Ginac copy constructor : compile object file \n";
498 for(
int i = 0; i < fun.M_fun.nops(); ++i ) exprs.append( fun.M_fun.op(i) );
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);
506 DVLOG(2) <<
"Ginac copy constructor : link with existing object file \n";
507 boost::mpi::communicator world;
509 std::string filenameWithSuffix = M_filename +
".so";
510 GiNaC::link_ex(filenameWithSuffix, M_cfun);
542 const GiNaC::FUNCP_CUBA& fun()
const
547 std::vector<GiNaC::symbol>
const& syms()
const {
return M_syms; }
551 template<
typename Geo_t,
typename Basis_i_t,
typename Basis_j_t>
555 typedef double value_type;
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;
563 typedef typename mn_to_shape<gmc_type::nDim,M,N>::type shape;
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;
575 static const bool value =
false;
578 tensor( this_type
const& expr,
579 Geo_t
const& geom, Basis_i_t
const& , Basis_j_t
const& )
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) )
588 tensor( this_type
const& expr,
589 Geo_t
const& geom, Basis_i_t
const& )
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) )
599 tensor( this_type
const& expr, Geo_t
const& geom )
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) )
609 template<
typename IM>
610 void init( IM
const& im )
614 void update( Geo_t
const& geom, Basis_i_t
const& , Basis_j_t
const& )
618 void update( Geo_t
const& geom, Basis_i_t
const& )
622 void update( Geo_t
const& geom )
624 M_gmc = fusion::at_key<key_type>( geom ).
get();
628 for(
int q = 0; q < M_gmc->nPoints();++q )
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 )
634 M_fun(&ni,M_x.data(),&no,M_y[q].data());
640 void update( Geo_t
const& geom, uint16_type )
642 M_gmc = fusion::at_key<key_type>( geom ).
get();
646 for(
int q = 0; q < M_gmc->nPoints();++q )
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 )
652 M_fun(&ni,M_x.data(),&no,M_y[q].data());
658 evalij( uint16_type i, uint16_type j )
const
665 evalijq( uint16_type , uint16_type , uint16_type c1, uint16_type c2, uint16_type q )
const
667 return M_y[q](c1,c2);
671 evaliq( uint16_type i, uint16_type c1, uint16_type c2, uint16_type q )
const
673 return M_y[q](c1,c2);
677 evalq( uint16_type c1, uint16_type c2, uint16_type q )
const
679 return M_y[q](c1,c2);
682 GiNaC::FUNCP_CUBA M_fun;
690 mutable expression_type M_fun;
691 std::vector<GiNaC::symbol> M_syms;
692 GiNaC::FUNCP_CUBA M_cfun;
693 std::string M_filename;
698 Expr< GinacMatrix<1,1,2> >
699 expr( GiNaC::matrix
const& f, std::vector<GiNaC::symbol>
const& lsym, std::string filename=
"")
701 return Expr< GinacMatrix<1,1,2> >( GinacMatrix<1,1,2>( f, lsym, filename ) );
708 template<
int M,
int N,
int Order>
710 Expr< GinacMatrix<M,N,Order> >
711 expr( GiNaC::matrix
const& f, std::vector<GiNaC::symbol>
const& lsym, std::string filename=
"" )
713 return Expr< GinacMatrix<M,N,Order> >( GinacMatrix<M,N,Order>( f, lsym, filename) );
716 template<
int M,
int N,
int Order>
718 Expr< GinacMatrix<M,N,Order> >
719 expr( GiNaC::ex
const& f, std::vector<GiNaC::symbol>
const& lsym, std::string filename=
"" )
721 return Expr< GinacMatrix<M,N,Order> >( GinacMatrix<M,N,Order>( f, lsym, filename ) );
Expr< Val< typename mpl::if_< boost::is_arithmetic< ExprT1 >, mpl::identity< Cst< ExprT1 > >, mpl::identity< ExprT1 > >::type::type > > val(ExprT1 const &__e1)
precompute expression tensor
Definition: val.hpp:304
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319