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
eim.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 <christophe.prudhomme@feelpp.org>
6  Date: 2012-05-02
7 
8  Copyright (C) 2012 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 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 _FEELPP_EIM_HPP
30 #define _FEELPP_EIM_HPP 1
31 
32 #include <limits>
33 #include <numeric>
34 #include <string>
35 #include <iostream>
36 #include <fstream>
37 
38 #include <boost/ref.hpp>
39 #include <boost/next_prior.hpp>
40 #include <boost/type_traits.hpp>
41 #include <boost/tuple/tuple.hpp>
42 #if BOOST_VERSION >= 104700
43 #include <boost/math/special_functions/nonfinite_num_facets.hpp>
44 #endif
45 #include <boost/serialization/vector.hpp>
46 #include <boost/serialization/string.hpp>
47 #include <boost/serialization/export.hpp>
48 #include <boost/serialization/base_object.hpp>
49 
50 #include <feel/feelcrb/crbdb.hpp>
52 
53 #include <feel/feelvf/vf.hpp>
54 
55 #include <feel/feelcrb/crb.hpp>
57 
58 #include <Eigen/Core>
59 
60 namespace Feel
61 {
62 class ModelCrbBaseBase {};
63 
125 template<typename ModelType>
126 class EIM : public CRBDB
127 {
128  typedef CRBDB super;
129 
130 public:
131 
132  static const uint16_type nDim = ModelType::nDim;
136  typedef ModelType model_type;
137  typedef typename ModelType::value_type value_type;
138 
139  typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix_type;
140  typedef typename matrix_type::ColXpr column_type;
141  typedef Eigen::Matrix<double, Eigen::Dynamic, 1> vector_type;
142  typedef typename ModelType::node_type node_type;
143 
144  typedef typename ModelType::functionspace_type functionspace_type;
145  typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
146  typedef typename functionspace_type::element_type element_type;
147  typedef typename functionspace_type::element_ptrtype element_ptrtype;
148 
149  typedef typename ModelType::solution_type solution_type;
150 
151  typedef typename ModelType::parameterspace_type parameterspace_type;
152  typedef typename ModelType::parameter_type parameter_type;
153  typedef typename parameterspace_type::sampling_ptrtype sampling_ptrtype;
154 
155  typedef boost::tuple<double,Eigen::Matrix<double,nDim,1> > space_residual_type;
156  typedef boost::tuple<double,parameter_type> parameter_residual_type;
157 
159 
163 
164  EIM()
165  :
166  super(),
167  M_vm(),
168  M_is_read( false ),
169  M_is_written( false ),
170  M_name( "default" ),
171  M_M( 1 ),
172  MM_max( 1 ),
173  M_WN(),
174  M_offline_done( false ),
175  M_tol( 1e-8 ),
176  M_q(),
177  M_B(),
178  M_t(),
179  M_index_max(),
180  M_model( 0 )
181  {}
182  EIM( po::variables_map const& vm, model_type* model, sampling_ptrtype sampling, double __tol = 1e-8 )
183  :
184  super(model->modelName(), model->name(), model->name(), vm ),
185  M_vm( vm ),
186  M_is_read( false ),
187  M_is_written( false ),
188  M_name( model->name() ),
189  M_trainset( sampling ),
190  M_M( 1 ),
191  MM_max( 1 ),
192  M_WN( M_vm["eim.dimension-max"].template as<int>() ),
193  M_offline_done( false ),
194  M_tol( __tol ),
195  M_q(),
196  M_B(),
197  M_t(),
198  M_index_max(),
199  M_model( model )
200  {
201  if ( !loadDB() || M_vm["eim.rebuild-database"].template as<bool>() )
202  {
203  LOG(INFO) << "construct EIM approximation...\n";
204  M_offline_done = false;
205  }
206 
207  if ( !this->isOfflineDone() )
208  offline();
209 #if 0
210  //old version of convergence
211  if( M_vm["eim.study-cvg"].template as<bool>() )
212  {
213  M_WN=1;
214  do{
215  offline();
216  studyConvergence();
217 
218  M_offline_done = false;
219  M_WN++;
220  }while( M_WN < M_vm["eim.dimension-max"].template as<int>() );
221  }
222 #endif
223 
224  }
225 
226  EIM( EIM const & __bbf )
227  :
228  super(__bbf),
229  M_is_read( __bbf.M_is_read ),
230  M_is_written( __bbf.M_is_written ),
231  M_name( __bbf.M_name ),
232  M_M( __bbf.M_M ),
233  MM_max (__bbf.MM_max ),
234  M_WN(__bbf.M_WN ),
235  M_offline_done( __bbf.M_offline_done ),
236  M_tol( __bbf.M_tol ),
237  M_q( __bbf.M_q ),
238  M_B( __bbf.M_B ),
239  M_t( __bbf.M_t ),
240  M_index_max( __bbf.M_index_max ),
241  M_model( __bbf.M_model )
242  {}
243  ~EIM()
244  {}
245 
247 
251 
252 
254 
258 
262  size_type nDOF() const { FEELPP_ASSERT( M_model != 0 ).error( "Invalid EIM model" ); return M_model->functionSpace()->nLocalDof(); }
263 
267  std::vector<element_type> const& q() const { return M_q; }
268 
272  element_type const&
273  q( size_type __m ) const
274  {
275  FEELPP_ASSERT( __m >= 0 && __m < M_M )( __m )( M_M ).error( "out of bounds access" );
276 
277  return M_q[ __m ];
278  }
279 
280  size_type mMax() const { return MM_max; }
281 
282 
287  bool isOfflineDone() const
288  {
289  return M_offline_done;
290  }
291 
293 
294 
295 
296 
300 
301  void setTrainSet( sampling_ptrtype pset ) { M_trainset = pset; }
302 
304 
308 
312  void saveDB()
313  {
314  fs::ofstream ofs( this->dbLocalPath() / this->dbFilename() );
315 
316  if ( ofs )
317  {
318  boost::archive::text_oarchive oa( ofs );
319  // write class instance to archive
320  oa << *this;
321  // archive and stream closed when destructors are called
322  }
323  }
324 
328  bool loadDB()
329  {
330  //if ( this->rebuildDB() )
331  //return false;
332 
333  fs::path db = this->lookForDB();
334 
335  if ( db.empty() )
336  return false;
337 
338  if ( !fs::exists( db ) )
339  return false;
340 
341  //std::cout << "Loading " << db << "...\n";
342  fs::ifstream ifs( db );
343 
344  if ( ifs )
345  {
346  boost::archive::text_iarchive ia( ifs );
347  // write class instance to archive
348  ia >> *this;
349  //std::cout << "Loading " << db << " done...\n";
350  this->setIsLoaded( true );
351  // archive and stream closed when destructors are called
352  return true;
353  }
354 
355  return false;
356  }
357 
358 
359 
376  vector_type beta( parameter_type const& mu ) const { return beta( mu, this->mMax() ); }
377  vector_type beta( parameter_type const& mu, solution_type const& T ) const { return beta( mu, T, this->mMax() ); }
378  vector_type beta( parameter_type const& mu, size_type M ) const;
379  vector_type beta( parameter_type const& mu, solution_type const& T, size_type M ) const;
380 
381  std::vector<double> studyConvergence( parameter_type const & mu, solution_type & solution ) const;
382 
383  void computationalTimeStatistics( std::string appname ) { return M_model->computationalTimeStatistics(); }
384  element_type residual ( size_type M ) const;
385 
386  parameter_residual_type computeBestFit( sampling_ptrtype trainset, int __M );
387 
388  element_type operator()( parameter_type const& mu , int N) const { return expansion( M_q, beta( mu ) , N); }
389  element_type operator()( parameter_type const& mu, solution_type const& T , int N ) const { return expansion( M_q, beta( mu, T ) , N ); }
393  void orthonormalize( std::vector<element_type>& );
394 
395  friend class boost::serialization::access;
396  template<class Archive>
397  void serialize(Archive & __ar, const unsigned int __version )
398  {
399  LOG(INFO) << "serializing...\n";
400 
401  __ar & BOOST_SERIALIZATION_NVP( M_name );
402 
403  LOG(INFO) << "name saved/loaded...\n";
404 
405  __ar & BOOST_SERIALIZATION_NVP( M_offline_done );
406  LOG(INFO) << "offline status...\n";
407 
408  __ar & BOOST_SERIALIZATION_NVP( MM_max );
409  LOG(INFO) << "M saved/loaded\n";
410  M_M = MM_max;
411 
412  // save index
413  __ar & BOOST_SERIALIZATION_NVP( M_index_max );
414  LOG(INFO) << "index saved/loaded\n";
415 
416  // save t
417  __ar & BOOST_SERIALIZATION_NVP( M_t );
418  LOG(INFO) << "t saved/loaded\n";
419 
420 
421  if ( Archive::is_loading::value )
422  {
423  for( int i = 0; i < MM_max; ++ i )
424  {
425  M_q.push_back( M_model->functionSpace()->element() );
426  }
427  for( int i = 0; i < MM_max; ++ i )
428  __ar & BOOST_SERIALIZATION_NVP( M_q[i] );
429  // save q
430  LOG(INFO) << "q saved/loaded\n";
431  }
432  else
433  {
434  for( int i = 0; i < MM_max; ++ i )
435  __ar & BOOST_SERIALIZATION_NVP( M_q[i] );
436  }
437 
438 
439  // save B
440  __ar & BOOST_SERIALIZATION_NVP( M_B );
441  LOG(INFO) << "B saved/loaded\n";
442 
443  }
445 
446 
447 protected:
448 
449  po::variables_map M_vm;
450  mutable bool M_is_read;
451  mutable bool M_is_written;
452 
453  std::string M_name;
454  sampling_ptrtype M_trainset;
455  size_type M_M;
456  size_type MM_max;
457 
458  size_type M_WN;
459 
460  mutable bool M_offline_done;
461 
462  double M_tol;
463 
464  std::vector<element_type> M_g;
465 
466  std::vector<element_type> M_q;
467 
468  matrix_type M_B;
469 
470  std::vector<node_type> M_t;
471 
472  std::vector<size_type> M_index_max;
473 
474  model_type* M_model;
475 protected:
476 
477 private:
478 
516  void offline();
517 public:
518  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
519 };
520 
521 template<typename ModelType>
522 typename EIM<ModelType>::vector_type
523 EIM<ModelType>::beta( parameter_type const& mu, size_type __M ) const
524 {
525  // beta=B_M\g(Od(indx),mut(i))'
526  vector_type __beta( __M );
527  for ( size_type __m = 0;__m < __M;++__m )
528  {
529  __beta[__m] = M_model->operator()( this->M_t[__m], mu );
530  }
531  this->M_B.block(0,0,__M,__M).template triangularView<Eigen::UnitLower>().solveInPlace(__beta);
532 
533  return __beta;
534 }
535 template<typename ModelType>
536 typename EIM<ModelType>::vector_type
537 EIM<ModelType>::beta( parameter_type const& mu, solution_type const& T, size_type __M ) const
538 {
539  // beta=B_M\g(Od(indx),mut(i))'
540  vector_type __beta( __M );
541  for ( size_type __m = 0;__m < __M;++__m )
542  {
543  __beta[__m] = M_model->operator()( T, this->M_t[__m], mu );
544  }
545  this->M_B.block(0,0,__M,__M).template triangularView<Eigen::UnitLower>().solveInPlace(__beta);
546 
547  return __beta;
548 }
549 
550 template<typename ModelType>
551 typename EIM<ModelType>::element_type
552 EIM<ModelType>::residual( size_type __M ) const
553 {
554  LOG(INFO) << "compute residual for m=" << __M << "...\n";
555  vector_type rhs( __M );
556 
557  //LOG(INFO) << "g[" << __M << "]=" << M_g[__M] << "\n";
558  for ( size_type __m = 0;__m < __M;++__m )
559  {
560  LOG(INFO) << "t[" << __m << "]=" << M_t[__m] << "\n";
561  rhs[__m]= M_g[__M]( M_t[__m] )(0,0,0);
562  }
563  this->M_B.block(0,0,__M,__M).template triangularView<Eigen::UnitLower>().solveInPlace(rhs);
564  LOG(INFO) << "solve B sol = rhs with rhs = " << rhs <<"\n";
565 
566  // res(:,i)=M_g(:,i)-q(:,0:i)*sigma
567  LOG(INFO) << "compute residual..." <<"\n";
568  using namespace vf;
569  auto z = expansion( M_q, rhs, __M );
570  LOG(INFO) << "return residual..." <<"\n";
571  return vf::project( _space=M_model->functionSpace(),
572  _expr=idv(M_g[__M])-idv( z ) );
573 }
574 
575 template<typename ModelType>
576 void
577 EIM<ModelType>::orthonormalize( std::vector<element_type> & __Z )
578 {
579  size_type __M = __Z.size();
580  for ( size_type __i = 0;__i < __M-1; ++__i )
581  {
582  value_type __s = inner_product(__Z[__i],__Z[__M-1]);
583  __Z[__M-1].add(- __s,__Z[__i]);
584  }
585  __Z[__M-1].scale(inner_product(__Z[__M-1],__Z[__M-1]));
586 }
587 
588 
589 template<typename ModelType>
590 typename EIM<ModelType>::parameter_residual_type
591 EIM<ModelType>::computeBestFit( sampling_ptrtype trainset, int __M )
592 {
593  LOG(INFO) << "compute best fit for m=" << __M
594  << " and trainset of size " << trainset->size() << "...\n";
595  using namespace vf;
596  parameter_type mu = M_model->parameterSpace()->element();
597 
598  vector_type maxerr( trainset->size() );
599  maxerr.setZero();
600  int index = 0;
601  LOG(INFO) << "Compute best fit M=" << __M << "\n";
602  BOOST_FOREACH( mu, *trainset )
603  {
604  LOG(INFO) << "compute best fit check mu...\n";
605  mu.check();
606  LOG_EVERY_N(INFO, 1 ) << " (every 10 mu) compute fit at mu="<< mu <<"\n" ;
607  // evaluate model at mu
608 
609  auto Z = M_model->operator()( mu );
610 
611  vector_type rhs( __M );
612  for ( size_type __m = 0;__m < __M;++__m )
613  {
614  rhs[__m]= Z( M_t[__m] )(0,0,0);
615  }
616  this->M_B.block(0,0,__M,__M).template triangularView<Eigen::UnitLower>().solveInPlace(rhs);
617  auto res = vf::project( _space=M_model->functionSpace(),
618  _expr=idv(Z)-idv( expansion( M_q, rhs, __M ) ) );
619 
620  std::string norm_used = option(_name="eim.norm-used-for-residual").template as<std::string>();
621  bool check_name_norm = false;
622  LOG( INFO ) << "[computeBestFit] norm used : "<<norm_used;
623  if( norm_used == "Linfty" )
624  {
625  check_name_norm=true;
626  auto resmax = normLinf( _range=elements(M_model->mesh()), _pset=_Q<5>(), _expr=idv(res) );
627  LOG_ASSERT( index < trainset->size() ) << "Invalid index " << index << " should be less than trainset size = " << trainset->size() << "\n";
628  maxerr( index++ ) = resmax.template get<0>();
629  }
630  if( norm_used == "L2" )
631  {
632  check_name_norm=true;
633  double norm = math::sqrt( integrate( _range=elements( M_model->mesh() ) ,_expr=idv(res)*idv(res) ).evaluate()( 0,0 ) );
634  LOG_ASSERT( index < trainset->size() ) << "Invalid index " << index << " should be less than trainset size = " << trainset->size() << "\n";
635  maxerr( index++ ) = norm;
636  }
637  if( norm_used == "LinftyVec" )
638  {
639  check_name_norm=true;
640  double norm = res.linftyNorm();
641  LOG_ASSERT( index < trainset->size() ) << "Invalid index " << index << " should be less than trainset size = " << trainset->size() << "\n";
642  maxerr( index++ ) = norm ;
643  }
644  CHECK( check_name_norm ) <<"[EIM] The name of the norm "<<norm_used<<" is not known\n";
645 
646  int index2;
647  auto err = maxerr.array().abs().maxCoeff( &index2 );
648  LOG_EVERY_N(INFO, 1 ) << " (every 10 mu) maxerr=" << err << " at index = " << index2 << " at mu = " << trainset->at(index2) << "\n";
649  }
650  LOG_ASSERT( index == trainset->size() ) << "Invalid index " << index << " should be equal to trainset size = " << trainset->size() << "\n";
651  auto err = maxerr.array().abs().maxCoeff( &index );
652  LOG(INFO)<< "err=" << err << " reached at index " << index << " and mu=" << trainset->at(index) << "\n";
653  return boost::make_tuple( err, trainset->at(index) );
654 }
655 template<typename ModelType>
656 void
657 EIM<ModelType>::offline( )
658 {
659  using namespace vf;
660 
661  if ( this->isOfflineDone() )
662  return;
663  LOG(INFO) << "[offline] starting offline stage ...\n";
664  M_M = 1;
665  LOG(INFO) << "[offline] create mu_1...\n";
666  // min element in Dmu to start with (in // each proc have the same element)
667  auto mu = M_model->parameterSpace()->min();
668 
669  if ( !M_trainset )
670  M_trainset = M_model->parameterSpace()->sampling();
671  if ( M_trainset->empty() )
672  {
673  int sampling_size = M_vm["eim.sampling-size"].template as<int>();
674  std::string file_name = ( boost::format("eim_trainset_%1%") % sampling_size ).str();
675  std::ifstream file ( file_name );
676  if( ! file )
677  {
678  M_trainset->randomize( sampling_size );
679  M_trainset->writeOnFile(file_name);
680  }
681  else
682  {
683  M_trainset->clear();
684  M_trainset->readFromFile(file_name);
685  }
686  }
687 
688  // store residual
689  auto res = M_model->functionSpace()->element();
690 
691  //if( !M_vm["eim.study-cvg"].template as<bool>() || M_WN == 1 )
692  // {
693  LOG(INFO) << "compute finite element solution at mu_1...\n";
694  M_g.push_back( M_model->operator()( mu ) );
695 
696  LOG(INFO) << "compute T^" << 0 << "...\n";
697  // Build T^0
698  auto zmax = normLinf( _range=elements(M_model->mesh()), _pset=_Q<5>(), _expr=idv(M_g[0]) );
699  // store space coordinate where max absolute value occurs
700  M_t.push_back( zmax.template get<1>() );
701  LOG(INFO) << "norm Linf = " << zmax.template get<0>() << " at " << zmax.template get<1>() << "\n";
702  //LOG(INFO) << "g = " << M_g[0] << "\n";
703 
704  LOG(INFO) << "compute and insert q_0...\n";
705  // insert first element
706  auto q = M_g[0];
707  //q.scale( 1./zmax.template get<0>() );
708  q.scale( 1./ M_g[0]( M_t[0] )( 0, 0, 0 ) );
709  M_q.push_back( q );
710 
711  LOG(INFO) << "compute entry (0,0) of interpolation matrix...\n";
712  this->M_B.resize( 1, 1 );
713  this->M_B( 0, 0 ) = 1;
714  CHECK( math::abs( M_q[0]( M_t[0] )( 0, 0, 0 ) - 1 ) < 1e-10 )
715  << "q[0](t[0] != 1 " << "q[0] = " << M_q[0]( M_t[0] )( 0, 0, 0 )
716  << " t[0] = "<< M_t[0] << "\n";
717 
718  ++M_M;
719  //}
720  //else
721  //M_M = M_WN - 1;
722 
726  double err = 1;
727 
728  //for residual storage
729  //auto res = M_model->functionSpace()->element();
730 
731  LOG(INFO) << "start greedy algorithm...\n";
732  //for( ; M_M < M_WN; ++M_M ) //err >= this->M_tol )
733  for( ; M_M < M_vm["eim.dimension-max"].template as<int>(); ++M_M ) //err >= this->M_tol )
734  {
735 
736  LOG(INFO) << "M=" << M_M << "...\n";
737 
738  LOG(INFO) << "compute best fit error...\n";
739  // compute mu = arg max inf ||G(.;mu)-z||_infty
740  auto bestfit = computeBestFit( M_trainset, this->M_M-1 );
741 
742  auto g_bestfit = M_model->operator()( bestfit.template get<1>() );
743  auto gmax = normLinf( _range=elements(M_model->mesh()), _pset=_Q<5>(), _expr=idv(g_bestfit) );
744 
745  DVLOG(2) << "best fit max error = " << bestfit.template get<0>() << " relative error = " << bestfit.template get<0>()/gmax.template get<0>() << " at mu = "
746  << bestfit.template get<1>() << " tolerance=" << M_vm["eim.error-max"].template as<double>() << "\n";
747 
748  //if we want to impose the use of dimension-max functions, we don't want to stop here
749  if ( (bestfit.template get<0>()/gmax.template get<0>()) < M_vm["eim.error-max"].template as<double>() && ! M_vm["eim.use-dimension-max-functions"].template as<bool>() )
750  break;
751 
756  LOG(INFO) << "[offline] S(" << this->M_M-1 << ") = " << bestfit.template get<1>() << "\n";
757 
758  // update M_g(:,M-1)
759  M_g.push_back( g_bestfit );
760 
761  //orthonormalize( M_g );
762 
763  // build T^m such that T^m-1 \subset T^m
764  LOG(INFO) << "[offline] compute residual M="<< M_M << "..." <<"\n";
765  res = this->residual(M_M-1);
766  //LOG(INFO) << "residual = " << res << "\n";
767  LOG(INFO) << "[offline] compute arg sup |residual|..." <<"\n";
768  auto resmax = normLinf( _range=elements(M_model->mesh()), _pset=_Q<5>(), _expr=idv(res) );
769  LOG(INFO) << "[offline] store coordinates where max absolute value is attained " << resmax.template get<1>() << "..." <<"\n";
770  // store space coordinate where max absolute value occurs
771  M_t.push_back( resmax.template get<1>() );
772 
773  LOG(INFO) << "[offline] scale new basis function by " << 1./resmax.template get<0>() << "..." <<"\n";
774  res.scale( 1./res(resmax.template get<1>())(0,0,0) );
775  LOG(INFO) << "store new basis function..." <<"\n";
776  M_q.push_back( res );
777 
778  std::for_each( M_t.begin(), M_t.end(), []( node_type const& t ) { DVLOG(2) << "t=" << t << "\n"; } );
779  // update interpolation matrix
780  M_B.conservativeResize( M_M, M_M );
781  for( int __i = 0; __i < M_M; ++__i )
782  {
783  for( int __j = 0; __j < M_M; ++__j )
784  {
785  this->M_B( __i, __j ) = M_q[__j]( M_t[__i] )(0,0,0);
786 
787  }
788  }
789  DVLOG(2) << "[offline] Interpolation matrix: M_B = " << this->M_B <<"\n";
790 #if 0
791  for( int __i = 0; __i < M_M; ++__i )
792  {
793  for( int __j = __i+1; __j < M_M; ++__j )
794  {
795  LOG_ASSERT( math::abs( M_q[__j]( M_t[__i] )(0,0,0)) < 1e-10 )
796  << "q[j](t_i) when j > i should be 0, = "
797  << math::abs( M_q[__j]( M_t[__i] )(0,0,0)) << "\n";
798  }
799  //this->M_B(__i,__i)=1;
800  //this->M_B(__i,__i)=;
801  LOG_ASSERT( math::abs( M_q[__i]( M_t[__i] )( 0, 0, 0 ) - 1 ) < 1e-10 )
802  << "q[ " << __i << "](t[" << __i << "] != 1 "
803  << "t[" << __i << "] = "<< M_t[__i] << "\n";
804  }
805 #endif
806 
807 
808  VLOG(2) << "================================================================================\n";
809  //if we want to impose the use of dimension-max functions, we don't want to stop here
810  if ( resmax.template get<0>() < M_vm["eim.error-max"].template as<double>() && ! M_vm["eim.use-dimension-max-functions"].template as<bool>() )
811  {
812  ++M_M;
813  break;
814  }
815  }
816 
817  LOG(INFO) << "[offline] M_max = " << M_M-1 << "...\n";
818  this->MM_max = this->M_M-1;
819 
820  this->M_offline_done = true;
821  LOG(INFO) << "[offline] saving DB...\n";
822  saveDB();
823  LOG(INFO) << "[offline] done with offline stage ...\n";
824 }
825 
826 template<typename ModelType>
827 std::vector<double>
828 EIM<ModelType>::studyConvergence( parameter_type const & mu , solution_type & solution) const
829 {
830  LOG(INFO) << " Convergence study \n";
831  int proc_number = Environment::worldComm().globalRank();
832 
833  std::vector<double> l2ErrorVec(this->mMax(), 0.0);
834 
835  std::string mu_str;
836  for ( int i=0; i<mu.size(); i++ )
837  mu_str= mu_str + ( boost::format( "_%1%" ) %mu(i) ).str() ;
838 
839  std::string file_name = "cvg-eim-"+M_model->name()+"-"+mu_str+".dat";
840  if( std::ifstream( file_name ) )
841  std::remove( file_name.c_str() );
842 
843  std::ofstream conv;
844  if( proc_number == Environment::worldComm().masterRank() )
845  {
846  conv.open(file_name, std::ios::app);
847  conv << "#Nb_basis" << "\t" << "L2_error" << "\n";
848  }
849 
850  bool use_expression = option(_name="eim.compute-error-with-truth-expression").template as<bool>();
851  int max = this->mMax();
852  for(int N=1; N<=max; N++)
853  {
854 
855  double exprl2norm = 0 , diffl2norm = 0 ;
856 
857  if( use_expression )
858  {
859  exprl2norm =M_model->expressionL2Norm( solution , mu );
860  auto eim_approximation = this->operator()(mu , solution, N);
861  diffl2norm = M_model->diffL2Norm( solution , mu , eim_approximation );
862  }
863  else
864  {
865  exprl2norm =M_model->projExpressionL2Norm( solution , mu );
866  auto eim_approximation = this->operator()(mu , solution, N);
867  diffl2norm = M_model->projDiffL2Norm( solution , mu , eim_approximation );
868  }
869 
870  double l2_error = diffl2norm / exprl2norm ;
871  //interpolation error : || projection_g - g ||_L2
872  double interpolation_error = M_model->interpolationError( solution , mu );
873 #if 0
874  int size = mu.size();
875  LOG(INFO)<<" mu = [ ";
876  for ( int i=0; i<size-1; i++ ) LOG(INFO)<< mu[i] <<" , ";
877  LOG(INFO)<< mu[size-1]<<" ]\n";
878 
879  //old version
880  // Compute expression
881  auto expression = M_model->operator()(mu);
882  // Compute eim expansion
883  auto eim_approximation = this->operator()(mu , N);
884  //Compute l2error
885  LOG(INFO) << "EIM name : " << M_model->name() << "\n";
886  double norm_l2_expression = expression.l2Norm();
887  LOG(INFO) << "norm_l2 expression = " << norm_l2_expression << "\n";
888  double norm_l2_approximation = eim_approximation.l2Norm();
889  LOG(INFO) << "norm_l2_approximation = " << norm_l2_approximation << "\n";
890  auto l2_error = math::abs( norm_l2_expression - norm_l2_approximation ) / norm_l2_expression;
891  LOG(INFO) << "norm l2 error = " << l2_error << "\n";
892 #endif
893  l2ErrorVec[N-1] = l2_error; // /!\ l2ErrorVec[i] represents error with i+1 bases
894 
895  if( proc_number == Environment::worldComm().masterRank() )
896  conv << N << "\t" << l2_error << "\t" << interpolation_error << "\n";
897 
898  }//loop over basis functions
899 
900  conv.close();
901  return l2ErrorVec;
902 }
903 
904 template<typename SpaceType, typename ModelSpaceType, typename ParameterSpaceType>
905 class EIMFunctionBase
906 {
907 public:
908 
909  typedef SpaceType functionspace_type;
910  typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
911  typedef typename functionspace_type::element_type element_type;
912  typedef typename functionspace_type::element_ptrtype element_ptrtype;
913  typedef typename functionspace_type::mesh_type mesh_type;
914  typedef typename functionspace_type::mesh_ptrtype mesh_ptrtype;
915  typedef typename functionspace_type::value_type value_type;
916 
917  typedef ModelSpaceType model_functionspace_type;
918  typedef typename model_functionspace_type::element_type solution_type;
919 
920  static const uint16_type nDim = mesh_type::nDim;
921 
922  typedef ParameterSpaceType parameterspace_type;
923  typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
924  typedef typename parameterspace_type::element_type parameter_type;
925  typedef typename parameterspace_type::sampling_ptrtype sampling_ptrtype;
926  typedef Eigen::Matrix<double, nDim, 1> node_type;
927 
928  typedef EIM<EIMFunctionBase<SpaceType,model_functionspace_type, ParameterSpaceType> > eim_type;
929  typedef typename eim_type::vector_type vector_type;
930 
931  typedef boost::shared_ptr<eim_type> eim_ptrtype;
932  //typedef typename eim_type::betam_type betam_type;
933  //typedef typename eim_type::qm_type qm_type;
934 
935 
936  EIMFunctionBase( po::variables_map const& vm,
937  functionspace_ptrtype fspace,
938  parameterspace_ptrtype pspace,
939  sampling_ptrtype sampling,
940  std::string const& modelname,
941  std::string const& name )
942  :
943  M_vm( vm ),
944  M_fspace( fspace ),
945  M_pspace( pspace ),
946  M_trainset( sampling ),
947  M_modelname( modelname ),
948  M_name( name )
949  {
950  LOG(INFO)<< "EimFunctionBase constructor\n";
951  }
952  virtual ~EIMFunctionBase()
953  {}
954  std::string const& name() const { return M_name; }
955  void setName( std::string const& name ) { M_name = name; }
956  std::string modelName() const { return M_modelname; }
957  void setModelName( std::string const& name ) { M_modelname = name; }
958 
959  functionspace_ptrtype functionSpace() const { return M_fspace; }
960  functionspace_ptrtype functionSpace() { return M_fspace; }
961  parameterspace_ptrtype parameterSpace() const { return M_pspace; }
962  parameterspace_ptrtype parameterSpace() { return M_pspace; }
963  sampling_ptrtype trainSet() { return M_trainset; }
964  sampling_ptrtype trainSet() const { return M_trainset; }
965  virtual void setTrainSet( sampling_ptrtype tset ) { M_trainset = tset; }
966 
967  void addOnlineTime( const double time )
968  {
969  int size = M_online_time.size();
970  M_online_time.conservativeResize( size+1 );
971  M_online_time( size ) = time;
972  }
973  Eigen::VectorXd onlineTime() const { return M_online_time; }
974 
975  mesh_ptrtype mesh() const { return M_fspace->mesh(); }
976  mesh_ptrtype mesh() { return M_fspace->mesh(); }
977 
978  virtual element_type operator()( parameter_type const& ) = 0;
979  virtual element_type operator()( solution_type const& T, parameter_type const& ) = 0;
980  virtual element_type interpolant( parameter_type const& ) = 0;
981  value_type operator()( node_type const& x, parameter_type const& mu )
982  {
983  VLOG(2) << "calling EIMFunctionBase::operator()( x=" << x << ", mu=" << mu << ")\n";
984  element_type v = this->operator()( mu );
985  value_type res = v(x)(0,0,0);
986  VLOG(2) << "EIMFunctionBase::operator() v(x)=" << res << "\n";
987  return res;
988  }
989 
990  // evaluate eim expansion at interpolation points in space and mu in parameter where T provides the coefficient
991  //value_type operator()( vector_type const& T, parameter_type const& mu ) = 0;
992 
993  value_type operator()( solution_type const& T, node_type const& x, parameter_type const& mu )
994  {
995  VLOG(2) << "calling EIMFunctionBase::operator()( x=" << x << ", mu=" << mu << ")\n";
996  element_type v = this->operator()( T, mu );
997  value_type res = v(x)(0,0,0);
998  VLOG(2) << "EIMFunctionBase::operator() v(x)=" << res << "\n";
999  return res;
1000  }
1001 
1002  virtual element_type const& q( int m ) const = 0;
1003  virtual vector_type beta( parameter_type const& mu ) const = 0;
1004  virtual vector_type beta( parameter_type const& mu, solution_type const& T ) const = 0;
1005  virtual size_type mMax() const = 0;
1006 
1007  virtual std::vector<double> studyConvergence( parameter_type const & mu , solution_type & solution ) const = 0;
1008  virtual void computationalTimeStatistics( std::string appname ) = 0;
1009  virtual double expressionL2Norm( solution_type const& T , parameter_type const& mu ) const = 0;
1010  virtual double diffL2Norm( solution_type const& T , parameter_type const& mu , element_type const& eim_expansion ) const = 0;
1011  virtual double projExpressionL2Norm( solution_type const& T , parameter_type const& mu ) const = 0;
1012  virtual double projDiffL2Norm( solution_type const& T , parameter_type const& mu , element_type const& eim_expansion ) const = 0;
1013  virtual double interpolationError( solution_type const& T , parameter_type const& mu ) const = 0;
1014 
1015  po::variables_map M_vm;
1016  functionspace_ptrtype M_fspace;
1017  parameterspace_ptrtype M_pspace;
1018  sampling_ptrtype M_trainset;
1019  std::string M_modelname;
1020  std::string M_name;
1021  Eigen::VectorXd M_online_time;//contains online computational time
1022 };
1023 
1024 
1025 
1026 template<typename ModelType, typename SpaceType, typename ExprType>
1027 class EIMFunction
1028  : public EIMFunctionBase<SpaceType, typename ModelType::functionspace_type, typename ModelType::parameterspace_type>
1029 {
1030  typedef EIMFunctionBase<SpaceType, typename ModelType::functionspace_type, typename ModelType::parameterspace_type> super;
1031 public:
1032  typedef ModelType model_type;
1033  //typedef ModelType* model_ptrtype;
1034  typedef boost::shared_ptr<model_type> model_ptrtype;
1035 
1036  typedef SpaceType functionspace_type;
1037  typedef boost::shared_ptr<functionspace_type> functionspace_ptrtype;
1038  typedef typename super::element_type element_type;
1039  typedef typename super::element_ptrtype element_ptrtype;
1040 
1041  typedef typename ModelType::element_type solution_type;
1042 
1043  typedef typename super::parameterspace_type parameterspace_type;
1044  typedef typename super::parameter_type parameter_type;
1045  typedef typename super::sampling_ptrtype sampling_ptrtype;
1046 
1047  typedef ExprType expr_type;
1048  typedef boost::shared_ptr<expr_type> expr_ptrtype;
1049 
1050  typedef typename super::eim_type eim_type;
1051  typedef typename super::eim_ptrtype eim_ptrtype;
1052 
1053  typedef typename super::vector_type vector_type;
1054 
1055  typedef CRBModel<ModelType> crbmodel_type;
1056  typedef boost::shared_ptr<crbmodel_type> crbmodel_ptrtype;
1057  typedef CRB<crbmodel_type> crb_type;
1058  typedef boost::shared_ptr<crb_type> crb_ptrtype;
1059 
1060  EIMFunction( po::variables_map const& vm,
1061  model_ptrtype model,
1062  functionspace_ptrtype space,
1063  solution_type& u,
1064  parameter_type& mu,
1065  expr_type& expr,
1066  sampling_ptrtype sampling,
1067  std::string const& name )
1068  :
1069  super( vm, space, model->parameterSpace(), sampling, model->modelName(), name ),
1070  M_vm( vm ),
1071  M_model( model ),
1072  M_expr( expr ),
1073  M_u( u ),
1074  M_mu( mu ),
1075  M_eim( new eim_type( vm, this, sampling ) )
1076  {
1077  }
1078  element_type operator()( parameter_type const& mu )
1079  {
1080  M_mu = mu;
1081 #if !defined(NDEBUG)
1082  M_mu.check();
1083 #endif
1084  M_u = M_model->solve( mu );
1085  //LOG(INFO) << "operator() mu=" << mu << "\n" << "sol=" << M_u << "\n";
1086  return vf::project( _space=this->functionSpace(), _expr=M_expr );
1087  }
1088  element_type operator()( solution_type const& T, parameter_type const& mu )
1089  {
1090  M_mu = mu;
1091 #if !defined(NDEBUG)
1092  M_mu.check();
1093 #endif
1094  // no need to solve we have already an approximation (typically from
1095  // an nonlinear iteration procedure)
1096  M_u = T;
1097  return vf::project( _space=this->functionSpace(), _expr=M_expr );
1098  }
1099 
1100  //Let g the expression that we want to have an eim expasion
1101  //here is computed its l2 norm : || g ||_L2
1102  double expressionL2Norm( solution_type const& T , parameter_type const& mu ) const
1103  {
1104  M_mu = mu;
1105 #if !defined(NDEBUG)
1106  M_mu.check();
1107 #endif
1108  M_u = T;
1109  auto mesh = this->functionSpace()->mesh();
1110  return math::sqrt( integrate( _range=elements( mesh ), _expr=M_expr*M_expr ).evaluate()( 0,0 ) );
1111  }
1112 
1113  //Let geim the eim expansion of g
1114  //here is computed || g - geim ||_L2
1115  double diffL2Norm( solution_type const& T , parameter_type const& mu , element_type const & eim_expansion ) const
1116  {
1117  M_mu = mu;
1118 #if !defined(NDEBUG)
1119  M_mu.check();
1120 #endif
1121  M_u = T;
1122  auto mesh = this->functionSpace()->mesh();
1123  auto difference = M_expr - idv(eim_expansion);
1124  return math::sqrt( integrate( _range=elements( mesh ), _expr=difference*difference ).evaluate()( 0,0 ) );
1125  }
1126 
1127 
1128  //Let \pi_g the projection of g on the function space
1129  //here is computed || \pi_g ||_L2
1130  double projExpressionL2Norm( solution_type const& T , parameter_type const& mu ) const
1131  {
1132  M_mu = mu;
1133 #if !defined(NDEBUG)
1134  M_mu.check();
1135 #endif
1136  M_u = T;
1137  auto mesh = this->functionSpace()->mesh();
1138  auto pi_g = vf::project( _space=this->functionSpace(), _expr=M_expr );
1139  return math::sqrt( integrate( _range=elements( mesh ), _expr=idv(pi_g)*idv(pi_g) ).evaluate()( 0,0 ) );
1140  }
1141 
1142  //here is computed || \pi_g - geim ||_L2
1143  double projDiffL2Norm( solution_type const& T , parameter_type const& mu , element_type const& eim_expansion ) const
1144  {
1145  M_mu = mu;
1146 #if !defined(NDEBUG)
1147  M_mu.check();
1148 #endif
1149  M_u = T;
1150  auto mesh = this->functionSpace()->mesh();
1151  auto pi_g = vf::project( _space=this->functionSpace(), _expr=M_expr );
1152  auto diff = pi_g - eim_expansion;
1153  return math::sqrt( integrate( _range=elements( mesh ), _expr=idv(diff)*idv(diff) ).evaluate()( 0,0 ) );
1154  }
1155 
1156  double interpolationError( solution_type const& T , parameter_type const& mu ) const
1157  {
1158  M_mu = mu;
1159 #if !defined(NDEBUG)
1160  M_mu.check();
1161 #endif
1162  M_u = T;
1163  auto mesh = this->functionSpace()->mesh();
1164  auto pi_g = vf::project( _space=this->functionSpace(), _expr=M_expr );
1165  auto difference = M_expr - idv(pi_g);
1166  return math::sqrt( integrate( _range=elements( mesh ), _expr=difference*difference ).evaluate()( 0,0 ) );
1167  }
1168 
1169  void computationalTimeStatistics( std::string appname )
1170  {
1171  computationalTimeStatistics( appname, typename boost::is_base_of<ModelCrbBaseBase,model_type>::type() );
1172  }
1173  void computationalTimeStatistics( std::string appname, boost::mpl::bool_<false> )
1174  {}
1175  void computationalTimeStatistics( std::string appname, boost::mpl::bool_<true> )
1176  {
1177  //auto crbmodel = crbmodel_ptrtype( new crbmodel_type( M_vm , CRBModelMode::CRB ) );
1178  auto crbmodel = crbmodel_ptrtype( new crbmodel_type( M_model , CRBModelMode::CRB ) );
1179  //make sure that the CRB DB is already build
1180  M_crb = crb_ptrtype( new crb_type( appname,
1181  M_vm ,
1182  crbmodel ) );
1183 
1184  if ( !M_crb->isDBLoaded() || M_crb->rebuildDB() )
1185  {
1186  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1187  LOG( INFO ) << "No CRB DB available, do crb offline computations...";
1188  M_crb->setOfflineStep( true );
1189  M_crb->offline();
1190  }
1191 
1192  int n_eval = option(_name="eim.computational-time-neval").template as<int>();
1193 
1194  Eigen::Matrix<double, Eigen::Dynamic, 1> time_crb;
1195  time_crb.resize( n_eval );
1196 
1197  typename crb_type::sampling_ptrtype Sampling( new typename crb_type::sampling_type( M_model->parameterSpace() ) );
1198  Sampling->logEquidistribute( n_eval );
1199 
1200  //dimension
1201  int N = option(_name="crb.dimension").template as<int>();
1202  //reduced basis approximation space
1203  auto WN = M_crb->wn();
1204  int mu_number = 0;
1205  BOOST_FOREACH( auto mu, *Sampling )
1206  {
1207  //LOG( INFO ) << "[computational] mu = \n"<<mu;
1208  boost::mpi::timer tcrb;
1209  auto o = M_crb->run( mu, option(_name="crb.online-tolerance").template as<double>() , N);
1210  auto solutions=o.template get<2>();
1211  auto uN = solutions.template get<0>();//vector of solutions ( one solution at each time step )
1212 
1213  int size=uN.size();
1214  auto u_crb = M_crb->expansion( uN[size-1] , N , WN );
1215 
1216  boost::mpi::timer teim;
1217  this->beta( mu , u_crb );
1218  this->addOnlineTime( teim.elapsed() );
1219  time_crb( mu_number ) = tcrb.elapsed() ;
1220  mu_number++;
1221  }
1222 
1223  M_model->computeStatistics( super::onlineTime() , super::name() );
1224  M_model->computeStatistics( time_crb , super::name()+" - global crb timing" );
1225 
1226  }
1227 
1228  void setTrainSet( sampling_ptrtype tset ) { M_eim->setTrainSet( tset ); }
1229  element_type interpolant( parameter_type const& mu ) { return M_eim->operator()( mu , M_eim->mMax() ); }
1230 
1231  element_type const& q( int m ) const { return M_eim->q( m ); }
1232 
1233  vector_type beta( parameter_type const& mu ) const { return M_eim->beta( mu ); }
1234  vector_type beta( parameter_type const& mu, solution_type const& T ) const { return M_eim->beta( mu, T ); }
1235 
1236  std::vector<double> studyConvergence( parameter_type const & mu , solution_type & solution ) const { return M_eim->studyConvergence( mu , solution ) ; };
1237 
1238  size_type mMax() const { return M_eim->mMax(); }
1239 
1240 private:
1241  po::variables_map M_vm;
1242  model_ptrtype M_model;
1243  expr_type M_expr;
1244  solution_type& M_u;
1245  parameter_type& M_mu;
1246  eim_ptrtype M_eim;
1247  crb_ptrtype M_crb;
1248 };
1249 
1250 namespace detail
1251 {
1252 
1253 template<typename Args>
1254 struct compute_eim_return
1255 {
1256  typedef typename boost::remove_reference<typename boost::remove_pointer<typename parameter::binding<Args, tag::model>::type>::type>::type::element_type model1_type;
1257  typedef typename boost::remove_const<typename boost::remove_pointer<model1_type>::type>::type model_type;
1258  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::expr>::type>::type expr_type;
1259  typedef typename boost::remove_reference<typename parameter::binding<Args, tag::space>::type>::type::element_type space_type;
1260  typedef EIMFunction<model_type, space_type, expr_type> type;
1261  typedef boost::shared_ptr<EIMFunction<model_type, space_type, expr_type> > ptrtype;
1262 };
1263 }
1264 
1265 BOOST_PARAMETER_FUNCTION(
1266  ( typename Feel::detail::compute_eim_return<Args>::ptrtype ), // 1. return type
1267  eim, // 2. name of the function template
1268  tag, // 3. namespace of tag types
1269  ( required
1270  ( in_out(model), * )
1271  ( in_out(element), * )
1272  ( in_out(parameter), * )
1273  ( in_out(expr), * )
1274  ( name, * )
1275  ( space, *)
1276  ) // required
1277  ( optional
1278  ( options, *, Environment::vm())
1279  //( space, *( boost::is_convertible<mpl::_,boost::shared_ptr<FunctionSpaceBase> > ), model->functionSpace() )
1280  //( space, *, model->functionSpace() )
1281  ( sampling, *, model->parameterSpace()->sampling() )
1282  ( verbose, (int), 0 )
1283  ) // optionnal
1284 )
1285 {
1286  Feel::detail::ignore_unused_variable_warning( args );
1287  typedef typename Feel::detail::compute_eim_return<Args>::type eim_type;
1288  typedef typename Feel::detail::compute_eim_return<Args>::ptrtype eim_ptrtype;
1289  return eim_ptrtype(new eim_type( options, model, space, element, parameter, expr, sampling, name ) );
1290 } // eim
1291 
1292 template<typename ModelType>
1293 struct EimFunctionNoSolve
1294 {
1295  typedef typename ModelType::functionspace_type functionspace_type;
1296  typedef typename ModelType::functionspace_ptrtype functionspace_ptrtype;
1297  typedef typename functionspace_type::element_type element_type;
1298  typedef typename ModelType::parameterspace_type parameterspace_type;
1299  typedef typename ModelType::parameterspace_ptrtype parameterspace_ptrtype;
1300  typedef typename parameterspace_type::element_type parameter_type;
1301  typedef typename parameterspace_type::sampling_type sampling_type;
1302  typedef typename parameterspace_type::sampling_ptrtype sampling_ptrtype;
1303  typedef typename functionspace_type::value_type value_type;
1304 
1305  typedef boost::shared_ptr<ModelType> model_ptrtype;
1306 
1307  static const int nb_spaces = functionspace_type::nSpaces;
1308  typedef typename mpl::if_< boost::is_same< mpl::int_<nb_spaces> , mpl::int_<2> > , fusion::vector< mpl::int_<0>, mpl::int_<1> > ,
1309  typename mpl::if_ < boost::is_same< mpl::int_<nb_spaces> , mpl::int_<3> > , fusion::vector < mpl::int_<0> , mpl::int_<1> , mpl::int_<2> > ,
1310  typename mpl::if_< boost::is_same< mpl::int_<nb_spaces> , mpl::int_<4> >, fusion::vector< mpl::int_<0>, mpl::int_<1>, mpl::int_<2>, mpl::int_<3> >,
1311  fusion::vector< mpl::int_<0>, mpl::int_<1>, mpl::int_<2>, mpl::int_<3>, mpl::int_<4> >
1312  >::type >::type >::type index_vector_type;
1313 
1314  EimFunctionNoSolve( model_ptrtype model )
1315  :
1316  M_model( model ),
1317  M_elt( M_model->functionSpace()->element() )
1318  {}
1319 
1320  element_type solve( parameter_type const& mu )
1321  {
1322  DVLOG(2) << "no solve required\n";
1323  static const bool is_composite = functionspace_type::is_composite;
1324  return solve( mu , mpl::bool_< is_composite >() );
1325  }
1326  element_type solve( parameter_type const& mu , mpl::bool_<false> )
1327  {
1328  value_type x = boost::lexical_cast<value_type>("inf");
1329  M_elt = vf::project( _space=M_model->functionSpace(), _expr=cst(x) );
1330  return M_elt;
1331  }
1332  element_type solve( parameter_type const& mu , mpl::bool_<true> )
1333  {
1334  ProjectInfCompositeCase project_inf_composite_case( M_elt );
1335  index_vector_type index_vector;
1336  fusion::for_each( index_vector, project_inf_composite_case );
1337  return project_inf_composite_case.element();
1338  }
1339 
1340  std::string modelName() const { return M_model->modelName(); }
1341  functionspace_ptrtype functionSpace() { return M_model->functionSpace(); }
1342  parameterspace_ptrtype parameterSpace() { return M_model->parameterSpace(); }
1343 
1344  struct ProjectInfCompositeCase
1345  {
1346  ProjectInfCompositeCase( element_type & composite_element)
1347  :
1348  M_element( composite_element )
1349  {}
1350 
1351  template< typename T >
1352  void
1353  operator()( const T& t ) const
1354  {
1355  auto view = M_element.template element< T::value >();
1356  auto space = view.functionSpace();
1357  view = vf::project( _space=space, _expr=cst( boost::lexical_cast<value_type>("inf") ) );
1358  }
1359 
1360  element_type element()
1361  {
1362  return M_element;
1363  }
1364 
1365  element_type M_element;
1366 
1367  }; //struct ProjectInfOnSubspace
1368 
1369  model_ptrtype M_model;
1370  element_type M_elt;
1371 };
1372 
1373 template<typename ModelType>
1374 boost::shared_ptr<EimFunctionNoSolve<ModelType>>
1375 eim_no_solve( boost::shared_ptr<ModelType> model )
1376 {
1377  return boost::shared_ptr<EimFunctionNoSolve<ModelType>>( new EimFunctionNoSolve<ModelType>( model ) );
1378 }
1379 
1380 
1381 po::options_description eimOptions( std::string const& prefix ="");
1382 }
1383 #endif /* _FEELPP_EIM_HPP */

Generated on Sun Oct 20 2013 08:24:56 for Feel++ by doxygen 1.8.4