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
crbmodel.hpp
Go to the documentation of this file.
1 /* -*- mode: c++ -*-
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2009-08-09
7 
8  Copyright (C) 2009 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 __CRBModel_H
30 #define __CRBModel_H 1
31 
32 #include <boost/shared_ptr.hpp>
33 
34 #include <vector>
35 
36 
38 #include <feel/feeldiscr/bdf2.hpp>
39 #include <feel/feelvf/vf.hpp>
40 
42 #include <feel/feeldiscr/operatorlinearcomposite.hpp>
44 #include <feel/feeldiscr/fsfunctionallinearcomposite.hpp>
45 
47 #include <feel/feelcore/pslogger.hpp>
48 
49 
50 namespace Feel
51 {
52 enum class CRBModelMode
53 {
54  PFEM = 0, SCM = 1, CRB = 2, SCM_ONLINE=3, CRB_ONLINE=4
55 };
56 
72 template<typename ModelType>
73 class CRBModel : public boost::enable_shared_from_this<CRBModel<ModelType> >
74 {
75 public:
76 
77 
81 
82  static const uint16_type ParameterSpaceDimension = ModelType::ParameterSpaceDimension;
83  static const bool is_time_dependent = ModelType::is_time_dependent;
84 
86 
90 
92  typedef ModelType model_type;
93  typedef boost::shared_ptr<ModelType> model_ptrtype;
94 
96  typedef typename model_type::value_type value_type;
98  typedef typename ModelType::mesh_type mesh_type;
99 
101  typedef typename ModelType::mesh_ptrtype mesh_ptrtype;
102 
104  typedef typename ModelType::space_type space_type;
105 
107  typedef typename model_type::functionspace_type functionspace_type;
108  typedef typename model_type::functionspace_ptrtype functionspace_ptrtype;
109 
111  typedef typename model_type::element_type element_type;
112  typedef typename model_type::element_ptrtype element_ptrtype;
113 
114  typedef typename model_type::backend_type backend_type;
115  typedef boost::shared_ptr<backend_type> backend_ptrtype;
116  typedef typename model_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
117  typedef typename model_type::vector_ptrtype vector_ptrtype;
118  typedef typename model_type::vector_type vector_type;
119 
120  typedef typename model_type::eigen_matrix_type eigen_matrix_type;
121 
122  typedef typename model_type::parameterspace_type parameterspace_type;
123  typedef typename model_type::parameterspace_ptrtype parameterspace_ptrtype;
124  typedef typename model_type::parameter_type parameter_type;
125  typedef typename model_type::parameter_ptrtype parameter_ptrtype;
126  typedef typename model_type::sampling_type sampling_type;
127  typedef typename model_type::sampling_ptrtype sampling_ptrtype;
128 
129 
130  typedef typename std::vector< std::vector < element_ptrtype > > initial_guess_type;
131  //typedef Eigen::VectorXd theta_vector_type;
132  typedef Eigen::VectorXd vectorN_type;
133  typedef std::vector< std::vector< double > > beta_vector_type;
134 
135  typedef typename boost::tuple<sparse_matrix_ptrtype,
136  sparse_matrix_ptrtype,
137  std::vector<vector_ptrtype>
138  > offline_merge_type;
139 
140 
141  typedef typename boost::tuple<std::vector< std::vector<sparse_matrix_ptrtype> >,
142  std::vector< std::vector<sparse_matrix_ptrtype> >,
143  std::vector< std::vector< std::vector<vector_ptrtype> > >
144  > affine_decomposition_type;
145 
146 
147  typedef typename boost::tuple< beta_vector_type,
148  beta_vector_type,
149  std::vector<beta_vector_type>
150  > betaqm_type;
151 
154  typedef boost::shared_ptr<bdf_type> bdf_ptrtype;
155 
156  typedef OperatorLinearComposite< space_type , space_type > operatorcomposite_type;
157  typedef boost::shared_ptr<operatorcomposite_type> operatorcomposite_ptrtype;
158 
159  typedef FsFunctionalLinearComposite< space_type > functionalcomposite_type;
160  typedef boost::shared_ptr<functionalcomposite_type> functionalcomposite_ptrtype;
161 
163  typedef boost::shared_ptr<preconditioner_type> preconditioner_ptrtype;
164 
165  static const int nb_spaces = functionspace_type::nSpaces;
166  typedef typename mpl::if_< boost::is_same< mpl::int_<nb_spaces> , mpl::int_<2> > , fusion::vector< mpl::int_<0>, mpl::int_<1> > ,
167  typename mpl::if_ < boost::is_same< mpl::int_<nb_spaces> , mpl::int_<3> > , fusion::vector < mpl::int_<0> , mpl::int_<1> , mpl::int_<2> > ,
168  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> >,
169  fusion::vector< mpl::int_<0>, mpl::int_<1>, mpl::int_<2>, mpl::int_<3>, mpl::int_<4> >
170  >::type >::type >::type index_vector_type;
171 
172 
174 
178 
179  CRBModel()
180  :
181  M_Aqm(),
182  M_Mqm(),
183  M_Fqm(),
184  M_is_initialized( false ),
185  M_mode( CRBModelMode::PFEM ),
186  M_model( new model_type() ),
187  M_backend( backend_type::build( BACKEND_PETSC ) ),
188  M_B()
189  {
190  this->init();
191  }
192 
193  CRBModel( po::variables_map const& vm, CRBModelMode mode = CRBModelMode::PFEM )
194  :
195  M_Aqm(),
196  M_InitialGuessV(),
197  M_InitialGuessVector(),
198  M_Mqm(),
199  M_Fqm(),
200  M_is_initialized( false ),
201  M_vm( vm ),
202  M_mode( mode ),
203  M_model( new model_type( vm ) ),
204  M_backend( backend_type::build( vm ) ),
205  M_B()
206  {
207  this->init();
208  }
209 
213  CRBModel( model_ptrtype & model )
214  :
215  M_Aqm(),
216  M_InitialGuessV(),
217  M_InitialGuessVector(),
218  M_Mqm(),
219  M_Fqm(),
220  M_is_initialized( false ),
221  M_vm(),
222  M_mode( CRBModelMode::PFEM ),
223  M_model( model ),
224  M_backend( backend_type::build( model->vm ) ),
225  M_B()
226  {
227  this->init();
228  }
229 
230  CRBModel( model_ptrtype & model , CRBModelMode mode )
231  :
232  M_Aqm(),
233  M_InitialGuessV(),
234  M_InitialGuessVector(),
235  M_Mqm(),
236  M_Fqm(),
237  M_is_initialized( false ),
238  M_vm(),
239  M_mode( mode ),
240  M_model( model ),
241  M_backend( backend_type::build( Environment::vm() ) ),
242  M_B()
243  {
244  this->init();
245  }
246 
250  CRBModel( CRBModel const & o )
251  :
252  M_Aqm( o.M_Aqm ),
253  M_InitialGuessV( o.M_InitialGuessV ),
254  M_InitialGuessVector( o.M_InitialGuessVector ),
255  M_Mqm( o.M_Mqm ),
256  M_Fqm( o.M_Fqm ),
257  M_is_initialized( o.M_is_initialized ),
258  M_vm( o.M_vm ),
259  M_mode( o.M_mode ),
260  M_model( o.M_model ),
261  M_backend( o.M_backend ),
262  M_B( o.M_B )
263  {
264  this->init();
265  }
266 
268  virtual ~CRBModel()
269  {}
270 
272  FEELPP_DONT_INLINE void init()
273  {
274  if ( M_is_initialized )
275  return;
276 
277  M_preconditioner = preconditioner(_pc=(PreconditionerType) M_backend->pcEnumType(), // by default : lu in seq or wirh mumps, else gasm in parallel
278  _backend= M_backend,
279  _pcfactormatsolverpackage=(MatSolverPackageType) M_backend->matSolverPackageEnumType(),// mumps if is installed ( by defaut )
280  _worldcomm=M_backend->comm(),
281  _prefix=M_backend->prefix() ,
282  _rebuild=true);
283  M_is_initialized=true;
284 
285  if( ! M_model->isInitialized() )
286  {
287  LOG( INFO ) << "CRBModel Model is not initialized";
288  M_model->initModel();
289  M_model->setInitialized( true );
290  }
291 
292  initB();
293 
294  if ( M_mode != CRBModelMode::CRB_ONLINE &&
295  M_mode != CRBModelMode::SCM_ONLINE )
296  {
297  //the model is already initialized
298  //std::cout << " -- init FEM model\n";
299  //M_model->init();
300  }
301 
302  auto Xh = M_model->functionSpace();
303  u = Xh->element();
304  v = Xh->element();
305 
306  }
307 
309 
313 
316  {
317  if ( this != &o )
318  {
319  M_Aqm = o.M_Aqm;
320  M_Fqm = o.M_Fqm;
321  M_Mqm = o.M_Mqm;
322  M_model = o.M_model;
323  M_backend = o.M_backend;
324  M_B = o.M_B;
325  }
326 
327  return *this;
328  }
330 
334 
338  po::variables_map vm() const
339  {
340  return M_vm;
341  }
342 
347  virtual sparse_matrix_ptrtype newMatrix() const
348  {
349  return M_model->newMatrix();
350  }
351 
356  virtual vector_ptrtype newVector() const
357  {
358  return M_model->newVector();
359  }
360 
364  sparse_matrix_ptrtype const& innerProduct() const
365  {
366  return M_B;
367  }
368 
372  sparse_matrix_ptrtype innerProduct()
373  {
374  return M_B;
375  }
376 
380  sparse_matrix_ptrtype const& innerProductForMassMatrix() const
381  {
382  return M_model->innerProductForMassMatrix();
383  }
384 
388  sparse_matrix_ptrtype innerProductForMassMatrix()
389  {
390  return M_model->innerProductForMassMatrix();
391  }
392 
393 
394 
395 
399  sparse_matrix_ptrtype const& h1() const
400  {
401  return M_B;
402  }
403 
404  sparse_matrix_ptrtype h1()
405  {
406  return M_B;
407  }
408 
410  functionspace_ptrtype functionSpace() const
411  {
412  return M_model->functionSpace();
413  }
414 
416  size_type Qa() const
417  {
418  return M_model->Qa();
419  }
420 
422  //int Qm() const { return 1; }//return M_model->Qm(); }
423 
424  size_type Qm() const
425  {
426  return Qm( mpl::bool_<model_type::is_time_dependent>() );
427  }
428  size_type Qm( mpl::bool_<true> ) const
429  {
430  return M_model->Qm();
431  }
432  size_type Qm( mpl::bool_<false> ) const
433  {
434  if( M_model->constructOperatorCompositeM() )
435  return functionspace_type::nSpaces;
436  else
437  return 1;
438  }
439 
440  int QInitialGuess() const
441  {
442  return M_model->QInitialGuess();
443  }
444 
445  int mMaxA(int q )
446  {
447  return M_model->mMaxA( q );
448  }
449 
450 
451  int mMaxM( int q )
452  {
453  return mMaxM(q, mpl::bool_<model_type::is_time_dependent>() );
454  }
455  int mMaxM( int q, mpl::bool_<true> )
456  {
457  return M_model->mMaxM( q );
458  }
459  int mMaxM( int q , mpl::bool_<false> )
460  {
461  return 1;
462  }
463 
464  int mMaxInitialGuess( int q ) const
465  {
466  return M_model->mMaxInitialGuess( q );
467  }
468 
469  int mMaxF(int output_index, int q )
470  {
471  return M_model->mMaxF( output_index, q );
472  }
473 
475  size_type Nl() const
476  {
477  return M_model->Nl();
478  }
479 
481  size_type Ql( int l ) const
482  {
483  return M_model->Ql( l );
484  }
485 
487  parameterspace_ptrtype parameterSpace() const
488  {
489  return M_model->parameterSpace();
490  }
491 
492 
493  parameter_type refParameter()
494  {
495  return M_model->refParameter();
496  }
498 
502 
506 #if 0
507  void setMeshSize( double s )
508  {
509  M_model->setMeshSize( s );
510  }
511 #endif
512 
514 
518 
519  beta_vector_type computeBetaInitialGuess( parameter_type const& mu ) const
520  {
521  return M_model->computeBetaInitialGuess( mu );
522  }
523 
527  betaqm_type computeBetaQm( parameter_type const& mu , double time=0 )
528  {
529  return computeBetaQm( mu , mpl::bool_<model_type::is_time_dependent>(), time );
530  }
531  betaqm_type computeBetaQm( parameter_type const& mu , mpl::bool_<true>, double time=0 )
532  {
533  return M_model->computeBetaQm( mu , time );
534  }
535  betaqm_type computeBetaQm( parameter_type const& mu , mpl::bool_<false>, double time=0 )
536  {
537  beta_vector_type betaAqm;
538  beta_vector_type betaMqm;
539  std::vector<beta_vector_type> betaFqm;
540  boost::tuple<
541  beta_vector_type,
542  std::vector<beta_vector_type> >
543  steady_beta;
544 
545  steady_beta = M_model->computeBetaQm( mu , time );
546  betaAqm = steady_beta.get<0>();
547  betaFqm = steady_beta.get<1>();
548 
549  int nspace = functionspace_type::nSpaces;
550  //if model provides implementation of operator composite M
551  if ( M_model->constructOperatorCompositeM() )
552  {
553  betaMqm.resize( nspace );
554  for(int q=0; q<nspace; q++)
555  {
556  betaMqm[q].resize(1);
557  betaMqm[q][0] = 1 ;
558  }
559  }
560  else
561  {
562  betaMqm.resize( 1 );
563  betaMqm[0].resize(1);
564  betaMqm[0][0] = 1 ;
565  }
566 
567  return boost::make_tuple( betaMqm, betaAqm, betaFqm );
568  }
569 
570 
571  betaqm_type computeBetaQm( element_type const& T, parameter_type const& mu , double time=0 )
572  {
573  return computeBetaQm( T , mu , mpl::bool_<model_type::is_time_dependent>(), time );
574  }
575  betaqm_type computeBetaQm( element_type const& T, parameter_type const& mu , mpl::bool_<true>, double time=0 )
576  {
577  return M_model->computeBetaQm( T, mu, time );
578  }
579  betaqm_type computeBetaQm( element_type const& T, parameter_type const& mu , mpl::bool_<false>, double time=0 )
580  {
581  beta_vector_type betaAqm, betaMqm ;
582  std::vector<beta_vector_type> betaFqm;
583  boost::tuple<
584  beta_vector_type,
585  std::vector<beta_vector_type> >
586  steady_beta;
587 
588  steady_beta = M_model->computeBetaQm(T, mu , time );
589  betaAqm = steady_beta.get<0>();
590  betaFqm = steady_beta.get<1>();
591 
592  int nspace = functionspace_type::nSpaces;
593  if ( M_model->constructOperatorCompositeM() )
594  {
595  betaMqm.resize( nspace );
596  for(int q=0; q<nspace; q++)
597  {
598  betaMqm[q].resize(1);
599  betaMqm[q][0] = 1 ;
600  }
601  }
602  else
603  {
604  betaMqm.resize( 1 );
605  betaMqm[0].resize(1);
606  betaMqm[0][0] = 1 ;
607  }
608 
609 
610  return boost::make_tuple( betaMqm, betaAqm, betaFqm );
611  }
612 
613 
614  element_ptrtype assembleInitialGuess( parameter_type const& mu );
615 
619  offline_merge_type update( parameter_type const& mu, double time=0 )
620  {
621  auto all_beta = this->computeBetaQm( mu , time );
622  offline_merge_type offline_merge;
623 
624  if( option(_name="crb.stock-matrices").template as<bool>() )
625  offline_merge = offlineMerge( all_beta , mu );
626  else
627  offline_merge = offlineMergeOnFly( all_beta, mu );
628 
629  return offline_merge;
630  }
631  offline_merge_type update( parameter_type const& mu, element_type const& T, double time=0 )
632  {
633  auto all_beta = this->computeBetaQm( T , mu , time );
634 
635  offline_merge_type offline_merge;
636 
637  if( option(_name="crb.stock-matrices").template as<bool>() )
638  offline_merge = offlineMerge( all_beta , mu );
639  else
640  offline_merge = offlineMergeOnFly( all_beta, mu );
641 
642  return offline_merge;
643 
644  }
645 
646  element_type solveFemUsingAffineDecompositionFixedPoint( parameter_type const& mu );
647  element_type solveFemDualUsingAffineDecompositionFixedPoint( parameter_type const& mu );
648  element_type solveFemUsingOfflineEim( parameter_type const& mu );
649 
650 
654  void initModel()
655  {
656  return M_model->initModel();
657  }
658  void setInitialized( const bool & b)
659  {
660  return M_model->setInitialized( b );
661  }
662  bool isInitialized()
663  {
664  return M_model->isInitialized();
665  }
666 
667 
668 
672  typename model_type::funs_type scalarContinuousEim()
673  {
674  return M_model->scalarContinuousEim();
675  }
676 
680  typename model_type::funsd_type scalarDiscontinuousEim()
681  {
682  return M_model->scalarDiscontinuousEim();
683  }
684 
685 
686  struct ComputeNormL2InCompositeCase
687  {
688 
689  ComputeNormL2InCompositeCase( element_type const composite_u1 ,
690  element_type const composite_u2 )
691  :
692  M_composite_u1 ( composite_u1 ),
693  M_composite_u2 ( composite_u2 )
694  {}
695 
696  template< typename T >
697  void
698  operator()( const T& t ) const
699  {
700  int i = T::value;
701  if( i == 0 )
702  M_vec.resize( 1 );
703  else
704  M_vec.conservativeResize( i+1 );
705 
706  auto u1 = M_composite_u1.template element< T::value >();
707  auto u2 = M_composite_u2.template element< T::value >();
708  mesh_ptrtype mesh = u1.functionSpace()->mesh();
709  double norm = normL2(_range=elements( mesh ),_expr=( idv(u1)-idv(u2) ) );
710  M_vec(i)= norm ;
711  }
712 
713  double norm()
714  {
715  return M_vec.sum();
716  }
717 
718  mutable vectorN_type M_vec;
719  element_type M_composite_u1;
720  element_type M_composite_u2;
721  };
722 
723 
724  double computeNormL2( element_type u1 , element_type u2 )
725  {
726  static const bool is_composite = functionspace_type::is_composite;
727  return computeNormL2( u1, u2 , mpl::bool_< is_composite >() );
728  }
729  double computeNormL2( element_type u1 , element_type u2 , mpl::bool_<false>)
730  {
731  double norm=normL2(_range=elements( u2.mesh() ),_expr=( idv(u1)-idv(u2) ) );
732  return norm;
733  }
734  double computeNormL2( element_type u1 , element_type u2 , mpl::bool_<true>)
735  {
736  //in that case we accumulate the norm of every components
737  ComputeNormL2InCompositeCase compute_normL2_in_composite_case( u1, u2 );
738  index_vector_type index_vector;
739  fusion::for_each( index_vector, compute_normL2_in_composite_case );
740  return compute_normL2_in_composite_case.norm();
741  }
742 
743 
753  //bilinear form
754  operatorcomposite_ptrtype operatorCompositeA() const
755  {
756  return M_model->operatorCompositeA();
757  }
758  //linear form
759  std::vector< functionalcomposite_ptrtype > functionalCompositeF() const
760  {
761  return M_model->functionalCompositeF();
762  }
763  //mass matrix
764  operatorcomposite_ptrtype operatorCompositeM() const
765  {
766  return operatorCompositeM( mpl::bool_<model_type::is_time_dependent>() );
767  }
768  operatorcomposite_ptrtype operatorCompositeM( mpl::bool_<true> ) const
769  {
770  return M_model->operatorCompositeM();
771  }
772  operatorcomposite_ptrtype operatorCompositeM( mpl::bool_<false> ) const
773  {
774  bool constructed_by_model = M_model->constructOperatorCompositeM();
775  if( constructed_by_model )
776  return M_model->operatorCompositeM();
777  else
778  return preAssembleMassMatrix();
779  }
780 
781 
791  affine_decomposition_type computeAffineDecomposition()
792  {
793  return computeAffineDecomposition( mpl::bool_<model_type::is_time_dependent>() );
794  }
795  affine_decomposition_type computeAffineDecomposition( mpl::bool_<true> )
796  {
797  boost::tie( M_Mqm, M_Aqm, M_Fqm ) = M_model->computeAffineDecomposition();
798 
799  if( M_Aqm.size() == 0 )
800  {
801  auto compositeM = operatorCompositeM();
802  int q_max = this->Qm();
803  M_Mqm.resize( q_max);
804  for(int q=0; q<q_max; q++)
805  {
806  int m_max = this->mMaxM(q);
807  M_Mqm[q].resize(m_max);
808  for(int m=0; m<m_max;m++)
809  {
810  auto operatorfree = compositeM->operatorlinear(q,m);
811  size_type pattern = operatorfree->pattern();
812  auto trial = operatorfree->domainSpace();
813  auto test=operatorfree->dualImageSpace();
814  M_Mqm[q][m]= M_backend->newMatrix( _test=test , _trial=trial , _pattern=pattern );
815  operatorfree->matPtr(M_Mqm[q][m]);//fill the matrix
816  }//m
817  }//q
818 
819  auto compositeA = operatorCompositeA();
820  q_max = this->Qa();
821  M_Aqm.resize( q_max);
822  for(int q=0; q<q_max; q++)
823  {
824  int m_max = this->mMaxA(q);
825  M_Aqm[q].resize(m_max);
826  for(int m=0; m<m_max;m++)
827  {
828  auto operatorfree = compositeA->operatorlinear(q,m);
829  size_type pattern = operatorfree->pattern();
830  auto trial = operatorfree->domainSpace();
831  auto test=operatorfree->dualImageSpace();
832  M_Aqm[q][m]= M_backend->newMatrix( _test=test , _trial=trial , _pattern=pattern );
833  operatorfree->matPtr(M_Aqm[q][m]);//fill the matrix
834  }//m
835  }//q
836 
837  auto vector_compositeF = functionalCompositeF();
838  int number_outputs = vector_compositeF.size();
839  M_Fqm.resize(number_outputs);
840  for(int output=0; output<number_outputs; output++)
841  {
842  auto composite_f = vector_compositeF[output];
843  int q_max = this->Ql(output);
844  M_Fqm[output].resize( q_max);
845  for(int q=0; q<q_max; q++)
846  {
847  int m_max = this->mMaxF(output,q);
848  M_Fqm[output][q].resize(m_max);
849  for(int m=0; m<m_max;m++)
850  {
851  auto operatorfree = composite_f->functionallinear(q,m);
852  auto space = operatorfree->space();
853  M_Fqm[output][q][m]= M_backend->newVector( space );
854  operatorfree->containerPtr(M_Fqm[output][q][m]);//fill the vector
855  }//m
856  }//q
857  }//output
858  }
859 
860  return boost::make_tuple( M_Mqm, M_Aqm, M_Fqm );
861  }
862  affine_decomposition_type computeAffineDecomposition( mpl::bool_<false> )
863  {
864  initial_guess_type initial_guess;
865  boost::tie( M_Aqm, M_Fqm ) = M_model->computeAffineDecomposition();
866 
867  if ( M_Aqm.size() > 0 )
868  {
869  assembleMassMatrix();
870  }
871  else
872  {
873  auto compositeM = operatorCompositeM();
874  int q_max = this->Qm();
875  M_Mqm.resize( q_max);
876  for(int q=0; q<q_max; q++)
877  {
878  int m_max = this->mMaxM(q);
879  M_Mqm[q].resize(m_max);
880  for(int m=0; m<m_max;m++)
881  {
882  auto operatorfree = compositeM->operatorlinear(q,m);
883  size_type pattern = operatorfree->pattern();
884  auto trial = operatorfree->domainSpace();
885  auto test=operatorfree->dualImageSpace();
886  M_Mqm[q][m]= M_backend->newMatrix( _test=test , _trial=trial , _pattern=pattern );
887  operatorfree->matPtr(M_Mqm[q][m]);//fill the matrix
888  }//m
889  }//q
890 
891  auto compositeA = operatorCompositeA();
892  q_max = this->Qa();
893  M_Aqm.resize( q_max);
894  for(int q=0; q<q_max; q++)
895  {
896  int m_max = this->mMaxA(q);
897  M_Aqm[q].resize(m_max);
898  for(int m=0; m<m_max;m++)
899  {
900  auto operatorfree = compositeA->operatorlinear(q,m);
901  size_type pattern = operatorfree->pattern();
902  auto trial = operatorfree->domainSpace();
903  auto test=operatorfree->dualImageSpace();
904  M_Aqm[q][m]= M_backend->newMatrix( _test=test , _trial=trial , _pattern=pattern );
905  operatorfree->matPtr(M_Aqm[q][m]);//fill the matrix
906  }//m
907  }//q
908 
909  auto vector_compositeF = functionalCompositeF();
910  int number_outputs = M_model->Nl();
911  M_Fqm.resize(number_outputs);
912  for(int output=0; output<number_outputs; output++)
913  {
914  auto composite_f = vector_compositeF[output];
915  int q_max = this->Ql(output);
916  M_Fqm[output].resize(q_max);
917  for(int q=0; q<q_max; q++)
918  {
919  int m_max = this->mMaxF(output,q);
920  M_Fqm[output][q].resize(m_max);
921  for(int m=0; m<m_max;m++)
922  {
923  auto functionalfree = composite_f->functionallinear(q,m);
924  auto space = functionalfree->space();
925  M_Fqm[output][q][m]= M_backend->newVector( space );
926  functionalfree->containerPtr(M_Fqm[output][q][m]);//fill the vector
927  }//m
928  }//q
929  }//output
930  }
931 
932  return boost::make_tuple( M_Mqm, M_Aqm, M_Fqm );
933  }
934 
935 
936  std::vector< std::vector<element_ptrtype> > computeInitialGuessAffineDecomposition( )
937  {
938  return M_model->computeInitialGuessAffineDecomposition( );
939  }
940 
941  std::vector< std::vector<element_ptrtype> > computeInitialGuessVAffineDecomposition( )
942  {
943  initial_guess_type initial_guess_v;
944  initial_guess_v = M_model->computeInitialGuessAffineDecomposition();
945  this->assembleInitialGuessV( initial_guess_v);
946  for(int q=0; q<M_InitialGuessVector.size(); q++)
947  {
948  for(int m=0; m<M_InitialGuessVector[q].size(); m++)
949  *M_InitialGuessV[q][m] = *M_InitialGuessVector[q][m];
950  }
951  return M_InitialGuessV;
952  }
953 
963  value_type h1( element_type const& xi_i, element_type const& xi_j ) const
964  {
965  return M_B->energy( xi_j, xi_i );
966  }
976  value_type h1( element_type const& xi_i ) const
977  {
978  return M_B->energy( xi_i, xi_i );
979  }
980 
981 
982 
991  sparse_matrix_ptrtype Aqm( uint16_type q, uint16_type m, bool transpose = false )
992  {
993  if( option(_name="crb.stock-matrices").template as<bool>() )
994  {
995  //in this case matrices have already been stocked
996  if ( transpose )
997  return M_Aqm[q][m]->transpose();
998 
999  return M_Aqm[q][m];
1000  }
1001  else
1002  {
1003  auto composite = operatorCompositeA();
1004  auto opfree = composite->operatorlinear(q,m);
1005  size_type pattern = opfree->pattern();
1006  auto trial = opfree->domainSpace();
1007  auto test = opfree->dualImageSpace();
1008  auto matrix = M_backend->newMatrix( _trial=trial , _test=test , _pattern=pattern );
1009  opfree->matPtr( matrix ); //assemble the matrix
1010 
1011  if ( transpose )
1012  return matrix->transpose();
1013 
1014  return matrix;
1015  }
1016  }
1017 
1018 
1027  const sparse_matrix_ptrtype Mqm( uint16_type q, uint16_type m, bool transpose = false ) const
1028  {
1029  if( option(_name="crb.stock-matrices").template as<bool>() )
1030  {
1031  //in this case matrices have already been stocked
1032  if ( transpose )
1033  return M_Mqm[q][m]->transpose();
1034 
1035  return M_Mqm[q][m];
1036  }
1037  else
1038  {
1039  auto composite = operatorCompositeM();
1040  auto opfree = composite->operatorlinear(q,m);
1041  size_type pattern = opfree->pattern();
1042  auto trial = opfree->domainSpace();
1043  auto test = opfree->dualImageSpace();
1044  auto matrix = M_backend->newMatrix( _trial=trial , _test=test , _pattern=pattern );
1045  opfree->matPtr( matrix ); //assemble the matrix
1046 
1047  if ( transpose )
1048  return matrix->transpose();
1049 
1050  return matrix;
1051  }
1052  }
1053 
1062  sparse_matrix_ptrtype Mqm( uint16_type q, uint16_type m, bool transpose = false )
1063  {
1064  if( option(_name="crb.stock-matrices").template as<bool>() )
1065  {
1066  //in this case matrices have already been stocked
1067  if ( transpose )
1068  return M_Mqm[q][m]->transpose();
1069 
1070  return M_Mqm[q][m];
1071  }
1072  else
1073  {
1074  auto composite = operatorCompositeM();
1075  auto opfree = composite->operatorlinear(q,m);
1076  size_type pattern = opfree->pattern();
1077  auto trial = opfree->domainSpace();
1078  auto test = opfree->dualImageSpace();
1079  auto matrix = M_backend->newMatrix( _trial=trial , _test=test , _pattern=pattern );
1080  opfree->matPtr( matrix ); //assemble the matrix
1081 
1082  if ( transpose )
1083  return matrix->transpose();
1084 
1085  return matrix;
1086  }
1087  }
1088 
1089  const vector_ptrtype InitialGuessVector( uint16_type q, uint16_type m ) const
1090  {
1091  return M_InitialGuessVector[q][m];
1092  }
1093 
1094  vector_ptrtype InitialGuessVector( uint16_type q, uint16_type m )
1095  {
1096  return M_InitialGuessVector[q][m];
1097  }
1098 
1099 
1110  value_type Aqm( uint16_type q, uint16_type m, element_type const& xi_i, element_type const& xi_j, bool transpose = false ) const
1111  {
1112  if( option(_name="crb.stock-matrices").template as<bool>() )
1113  {
1114  //in this case matrices have already been stocked
1115  return M_Aqm[q][m]->energy( xi_j, xi_i, transpose );
1116  }
1117  else
1118  {
1119  auto composite = operatorCompositeA();
1120  auto opfree = composite->operatorlinear(q,m);
1121  size_type pattern = opfree->pattern();
1122  auto trial = opfree->domainSpace();
1123  auto test = opfree->dualImageSpace();
1124  auto matrix = M_backend->newMatrix( _trial=trial , _test=test , _pattern=pattern );
1125  opfree->matPtr( matrix ); //assemble the matrix
1126 
1127  return matrix->energy( xi_j, xi_i, transpose );
1128  }
1129  }
1130 
1141  value_type Mqm( uint16_type q, uint16_type m, element_type const& xi_i, element_type const& xi_j, bool transpose = false ) const
1142  {
1143  if( option(_name="crb.stock-matrices").template as<bool>() )
1144  {
1145  //in this case matrices have already been stocked
1146  return M_Mqm[q][m]->energy( xi_j, xi_i, transpose );
1147  }
1148  else
1149  {
1150  auto composite = operatorCompositeM();
1151  auto opfree = composite->operatorlinear(q,m);
1152  size_type pattern = opfree->pattern();
1153  auto trial = opfree->domainSpace();
1154  auto test = opfree->dualImageSpace();
1155  auto matrix = M_backend->newMatrix( _trial=trial , _test=test , _pattern=pattern );
1156  opfree->matPtr( matrix ); //assemble the matrix
1157 
1158  return matrix->energy( xi_j, xi_i, transpose );
1159  }
1160  }
1161 
1162 
1163 
1164 
1165  beta_vector_type const& betaInitialGuessQm( mpl::bool_<true> ) const
1166  {
1167  return M_model->betaInitialGuessQm();
1168  }
1169 
1170 
1176  vector_ptrtype Fqm( uint16_type l, uint16_type q, int m ) const
1177  {
1178  if( option(_name="crb.stock-matrices").template as<bool>() )
1179  {
1180  return M_Fqm[l][q][m];
1181  }
1182  else
1183  {
1184  auto vector_composite = functionalCompositeF();
1185  auto composite = vector_composite[l];
1186  auto functional = composite->functionallinear(q,m);
1187  auto space = functional->space();
1188  auto vector = M_backend->newVector( space );
1189  functional->containerPtr( vector );
1190  return vector;
1191  }
1192  }
1193 
1194  element_ptrtype InitialGuessQm( uint16_type q, int m ) const
1195  {
1196  return M_InitialGuessV[q][m];
1197  }
1198 
1210  value_type Fqm( uint16_type l, uint16_type q, uint16_type m, element_type const& xi )
1211  {
1212 
1213  value_type result=0;
1214 
1215  if( option(_name="crb.stock-matrices").template as<bool>() )
1216  {
1217  result = inner_product( *M_Fqm[l][q][m] , xi );
1218  }
1219  else
1220  {
1221  auto vector_composite = functionalCompositeF();
1222  auto composite = vector_composite[l];
1223  auto functional = composite->functionallinear(q,m);
1224  auto space = functional->space();
1225  auto vector = M_backend->newVector( space );
1226  functional->containerPtr( vector );
1227  result = inner_product( *vector, xi ) ;
1228  }
1229 
1230  return result;
1231  }
1232 
1233 
1245  value_type Fqm( uint16_type l, uint16_type q, uint16_type m, element_ptrtype const& xi )
1246  {
1247  value_type result=0;
1248 
1249  if( option(_name="crb.stock-matrices").template as<bool>() )
1250  {
1251  result = inner_product( *M_Fqm[l][q][m] , xi );
1252  }
1253  else
1254  {
1255  auto vector_composite = functionalCompositeF();
1256  auto composite = vector_composite[l];
1257  auto functional = composite->functionallinear(q,m);
1258  auto space = functional->space();
1259  auto vector = M_backend->newVector( space );
1260  functional->containerPtr( vector );
1261  result = inner_product( *vector, xi ) ;
1262  }
1263 
1264  return result;
1265  }
1266 
1267 
1271  double scalarProduct( vector_type const& X, vector_type const& Y )
1272  {
1273  return M_model->scalarProduct( X, Y );
1274  }
1278  double scalarProduct( vector_ptrtype const& X, vector_ptrtype const& Y )
1279  {
1280  return M_model->scalarProduct( X, Y );
1281  }
1282 
1283 
1287  double scalarProductForMassMatrix( vector_type const& X, vector_type const& Y )
1288  {
1289  return M_model->scalarProductForMassMatrix( X, Y );
1290  }
1294  double scalarProductForMassMatrix( vector_ptrtype const& X, vector_ptrtype const& Y )
1295  {
1296  return M_model->scalarProductForMassMatrix( X, Y );
1297  }
1298 
1299 
1303  double scalarProductForPod( vector_type const& X, vector_type const& Y )
1304  {
1305  return scalarProductForPod( X, Y ,mpl::bool_<model_type::is_time_dependent>() );
1306  }
1307  double scalarProductForPod( vector_type const& X, vector_type const& Y , mpl::bool_<true> )
1308  {
1309  return M_model->scalarProductForPod( X, Y );
1310  }
1311  double scalarProductForPod( vector_type const& X, vector_type const& Y , mpl::bool_<false> )
1312  {
1313  return 0;
1314  }
1315 
1316 
1320  double scalarProductForPod( vector_ptrtype const& X, vector_ptrtype const& Y )
1321  {
1322  return scalarProductForPod( X, Y , mpl::bool_<model_type::is_time_dependent>() );
1323  }
1324  double scalarProductForPod( vector_ptrtype const& X, vector_ptrtype const& Y , mpl::bool_<true> )
1325  {
1326  return M_model->scalarProductForPod( X, Y );
1327  }
1328  double scalarProductForPod( vector_ptrtype const& X, vector_ptrtype const& Y , mpl::bool_<false> )
1329  {
1330  return 0;
1331  }
1332 
1333 
1337  element_type solve( parameter_type const& mu )
1338  {
1339  //return this->solveFemUsingAffineDecompositionFixedPoint( mu );
1340  return M_model->solve( mu );
1341  }
1342 
1343 
1347  void solve( parameter_type const& mu, element_ptrtype& u )
1348  {
1349  return M_model->solve( mu, u );
1350  }
1351 
1356  void solve( parameter_type const& mu, element_ptrtype& u, vector_ptrtype const& L, bool transpose = false )
1357  {
1358  return M_model->solve( mu, u, L, transpose );
1359  }
1360 
1365  void l2solve( vector_ptrtype& u, vector_ptrtype const& f )
1366  {
1367  return M_model->l2solve( u, f );
1368  }
1369 
1370 
1383  boost::tuple<int, value_type>
1384  solve( sparse_matrix_ptrtype const& A,
1385  vector_ptrtype & x,
1386  vector_ptrtype const& f,
1387  bool tranpose = false
1388  )
1389  {
1390  //return M_model->solve( A, x, f );
1391  }
1392 
1398 #if 0
1399  bool
1400  exportResults( double time, std::vector<element_ptrtype> const& v )
1401  {
1402  //return model->export( time, v );
1403  }
1404 #endif
1405 
1408  void run( const double * X, unsigned long N, double * Y, unsigned long P );
1409 
1414  value_type output( int output_index, parameter_type const& mu , element_type &u , bool need_to_solve=false)
1415  {
1416  return M_model->output( output_index, mu , u , need_to_solve );
1417  }
1418 
1419  int computeNumberOfSnapshots()
1420  {
1421  return computeNumberOfSnapshots( mpl::bool_<model_type::is_time_dependent>() );
1422  }
1423  int computeNumberOfSnapshots( mpl::bool_<true> )
1424  {
1425  return M_model->computeNumberOfSnapshots();
1426  }
1427  int computeNumberOfSnapshots( mpl::bool_<false> )
1428  {
1429  return 1;
1430  }
1431 
1432  vectorN_type computeStatistics ( Eigen::VectorXd vector , std::string name )
1433  {
1434  return M_model->computeStatistics( vector , name );
1435  }
1436 
1437  double timeStep()
1438  {
1439  return timeStep( mpl::bool_<model_type::is_time_dependent>() );
1440  }
1441  double timeStep( mpl::bool_<true> )
1442  {
1443  double timestep;
1444 
1445  if ( M_model->isSteady() ) timestep=1e30;
1446 
1447  else timestep = M_model->timeStep();
1448 
1449  return timestep;
1450  }
1451  double timeStep( mpl::bool_<false> )
1452  {
1453  return 1e30;
1454  }
1455 
1456  double timeInitial()
1457  {
1458  return timeInitial( mpl::bool_<model_type::is_time_dependent>() );
1459  }
1460  double timeInitial( mpl::bool_<true> )
1461  {
1462  return M_model->timeInitial();
1463  }
1464  double timeInitial( mpl::bool_<false> )
1465  {
1466  return 0;
1467  }
1468 
1469  double timeFinal()
1470  {
1471  return timeFinal( mpl::bool_<model_type::is_time_dependent>() );
1472  }
1473  double timeFinal( mpl::bool_<true> )
1474  {
1475  double timefinal;
1476 
1477  if ( M_model->isSteady() ) timefinal=1e30;
1478 
1479  else timefinal = M_model->timeFinal();
1480 
1481  return timefinal;
1482  }
1483  double timeFinal( mpl::bool_<false> )
1484  {
1485  return 1e30;
1486  }
1487 
1488  int timeOrder()
1489  {
1490  return timeOrder( mpl::bool_<model_type::is_time_dependent>() );
1491  }
1492  int timeOrder( mpl::bool_<true> )
1493  {
1494  return M_model->timeOrder();
1495  }
1496  int timeOrder( mpl::bool_<false> )
1497  {
1498  return 0;
1499  }
1500 
1501 
1502  bool isSteady()
1503  {
1504  return isSteady( mpl::bool_<model_type::is_time_dependent>() );
1505  }
1506  bool isSteady( mpl::bool_<true> )
1507  {
1508  return M_model->isSteady();
1509  }
1510  bool isSteady( mpl::bool_<false> )
1511  {
1512  return true;
1513  }
1514 
1515 
1516  void initializationField( element_ptrtype& initial_field,parameter_type const& mu )
1517  {
1518  return initializationField( initial_field,mu,mpl::bool_<model_type::is_time_dependent>() );
1519  }
1520  void initializationField( element_ptrtype& initial_field,parameter_type const& mu,mpl::bool_<true> )
1521  {
1522  return M_model->initializationField( initial_field,mu );
1523  }
1524  void initializationField( element_ptrtype& initial_field,parameter_type const& mu,mpl::bool_<false> ) {};
1525 
1526 
1528 
1529 
1530 
1531 protected:
1532 
1533 
1535  std::vector< std::vector<sparse_matrix_ptrtype> > M_Aqm;
1536 
1537  mutable std::vector< std::vector<element_ptrtype> > M_InitialGuessV;
1538  mutable std::vector< std::vector<vector_ptrtype> > M_InitialGuessVector;
1539 
1541  mutable std::vector< std::vector<sparse_matrix_ptrtype> > M_Mqm;
1542 
1544  std::vector< std::vector<std::vector<vector_ptrtype> > > M_Fqm;
1545 
1546 private:
1547 
1548  bool M_is_initialized;
1549 
1551  po::variables_map M_vm;
1552 
1554  CRBModelMode M_mode;
1555 
1557  model_ptrtype M_model;
1558 
1559  backend_ptrtype M_backend;
1560 
1561  // ! matrix associated with inner product
1562  sparse_matrix_ptrtype M_B;
1563  sparse_matrix_ptrtype M_H1;
1564 
1565  beta_vector_type M_dummy_betaMqm;
1566 
1568  void initB();
1569 
1576  offline_merge_type offlineMerge( betaqm_type const& all_beta, parameter_type const& mu );
1577  offline_merge_type offlineMergeOnFly( betaqm_type const& all_beta , parameter_type const& mu );
1578 
1579  void assembleMassMatrix( );
1580  void assembleMassMatrix( mpl::bool_<true> );
1581  void assembleMassMatrix( mpl::bool_<false> );
1582 
1583  operatorcomposite_ptrtype preAssembleMassMatrix( ) const ;
1584  operatorcomposite_ptrtype preAssembleMassMatrix( mpl::bool_<true> ) const ;
1585  operatorcomposite_ptrtype preAssembleMassMatrix( mpl::bool_<false> ) const ;
1586 
1587  void assembleInitialGuessV( initial_guess_type & initial_guess );
1588  void assembleInitialGuessV( initial_guess_type & initial_guess, mpl::bool_<true> );
1589  void assembleInitialGuessV( initial_guess_type & initial_guess, mpl::bool_<false> );
1590  element_type u,v;
1591 
1592  preconditioner_ptrtype M_preconditioner;
1593 
1594 };
1595 
1596 
1597 template <typename ModelType>
1598 struct PreAssembleMassMatrixInCompositeCase
1599 {
1600 
1601 
1603  typedef typename ModelType::mesh_type mesh_type;
1604 
1606  typedef typename ModelType::mesh_ptrtype mesh_ptrtype;
1607 
1609  typedef typename ModelType::space_type space_type;
1610 
1612  typedef typename ModelType::functionspace_type functionspace_type;
1613  typedef typename ModelType::functionspace_ptrtype functionspace_ptrtype;
1614 
1616  typedef typename ModelType::element_type element_type;
1617  typedef typename ModelType::element_ptrtype element_ptrtype;
1618 
1619 
1620  typedef OperatorLinearComposite<space_type, space_type> operatorcomposite_type;
1621  typedef boost::shared_ptr<operatorcomposite_type> operatorcomposite_ptrtype;
1622 
1623  PreAssembleMassMatrixInCompositeCase( element_type const u ,
1624  element_type const v )
1625  :
1626  M_composite_u ( u ),
1627  M_composite_v ( v )
1628  {
1629  auto Xh = M_composite_u.functionSpace();
1630  op_mass = opLinearComposite( _imageSpace=Xh, _domainSpace=Xh );
1631  }
1632 
1633  template< typename T >
1634  void
1635  operator()( const T& t ) const
1636  {
1637 
1638  using namespace Feel::vf;
1639 
1640  auto u = M_composite_u.template element< T::value >();
1641  auto v = M_composite_v.template element< T::value >();
1642  auto Xh = M_composite_u.functionSpace();
1643  mesh_ptrtype mesh = Xh->mesh();
1644 
1645  auto expr = integrate( _range=elements( mesh ) , _expr=trans( idt( u ) )*id( v ) ) ;
1646  auto opfree = opLinearFree( _imageSpace=Xh, _domainSpace=Xh, _expr=expr );
1647 
1648  //each composant of the affine decomposition
1649  //is the subspace contribution
1650  int q=T::value;
1651  int m=0;
1652  auto tuple = boost::make_tuple( q , m );
1653  op_mass->addElement( tuple , opfree );
1654  }
1655 
1656 
1657  operatorcomposite_ptrtype opmass()
1658  {
1659  return op_mass;
1660  }
1661 
1662  element_type M_composite_u;
1663  element_type M_composite_v;
1664  operatorcomposite_ptrtype op_mass;
1665 
1666 };
1667 
1668 
1669 template <typename ModelType>
1670 struct AssembleMassMatrixInCompositeCase
1671 {
1672 
1674  typedef typename ModelType::mesh_type mesh_type;
1675 
1677  typedef typename ModelType::mesh_ptrtype mesh_ptrtype;
1678 
1680  typedef typename ModelType::space_type space_type;
1681 
1683  typedef typename ModelType::functionspace_type functionspace_type;
1684  typedef typename ModelType::functionspace_ptrtype functionspace_ptrtype;
1685 
1687  typedef typename ModelType::element_type element_type;
1688  typedef typename ModelType::element_ptrtype element_ptrtype;
1689 
1690 
1691  AssembleMassMatrixInCompositeCase( element_type const u ,
1692  element_type const v ,
1693  boost::shared_ptr<CRBModel<ModelType> > crb_model)
1694  :
1695  M_composite_u ( u ),
1696  M_composite_v ( v ),
1697  M_crb_model( crb_model )
1698  {}
1699 
1700  template< typename T >
1701  void
1702  operator()( const T& t ) const
1703  {
1704 
1705  using namespace Feel::vf;
1706 
1707  auto u = M_composite_u.template element< T::value >();
1708  auto v = M_composite_v.template element< T::value >();
1709  auto Xh = M_composite_u.functionSpace();
1710  mesh_ptrtype mesh = Xh->mesh();
1711 
1712  form2( _test=Xh, _trial=Xh, _matrix=M_crb_model->Mqm(0,0) ) +=
1713  integrate( _range=elements( mesh ), _expr=trans( idt( u ) )*id( v ) );
1714 
1715  }
1716 
1717  element_type M_composite_u;
1718  element_type M_composite_v;
1719  mutable boost::shared_ptr<CRBModel<ModelType> > M_crb_model;
1720 };
1721 
1722 template <typename ModelType>
1723 struct AssembleInitialGuessVInCompositeCase
1724 {
1725 
1727  typedef typename ModelType::mesh_type mesh_type;
1728 
1730  typedef typename ModelType::mesh_ptrtype mesh_ptrtype;
1731 
1733  typedef typename ModelType::space_type space_type;
1734 
1736  typedef typename ModelType::functionspace_type functionspace_type;
1737  typedef typename ModelType::functionspace_ptrtype functionspace_ptrtype;
1738 
1740  typedef typename ModelType::element_type element_type;
1741  typedef typename ModelType::element_ptrtype element_ptrtype;
1742 
1743  typedef typename std::vector< std::vector < element_ptrtype > > initial_guess_type;
1744 
1745  AssembleInitialGuessVInCompositeCase( element_type const v ,
1746  initial_guess_type const initial_guess ,
1747  boost::shared_ptr<CRBModel<ModelType> > crb_model)
1748  :
1749  M_composite_v ( v ),
1750  M_composite_initial_guess ( initial_guess ),
1751  M_crb_model ( crb_model )
1752  {}
1753 
1754  template< typename T >
1755  void
1756  operator()( const T& t ) const
1757  {
1758  auto v = M_composite_v.template element< T::value >();
1759  auto Xh = M_composite_v.functionSpace();
1760  mesh_ptrtype mesh = Xh->mesh();
1761  int q_max = M_crb_model->QInitialGuess();
1762  for(int q = 0; q < q_max; q++)
1763  {
1764  int m_max = M_crb_model->mMaxInitialGuess(q);
1765  for( int m = 0; m < m_max ; m++)
1766  {
1767  auto initial_guess_qm = M_crb_model->InitialGuessVector(q,m);
1768  auto view = M_composite_initial_guess[q][m]->template element< T::value >();
1769  form1( _test=Xh, _vector=initial_guess_qm ) +=
1770  integrate ( _range=elements( mesh ), _expr=trans( Feel::vf::idv( view ) )*Feel::vf::id( v ) );
1771  }
1772  }
1773 
1774  }
1775 
1776  element_type M_composite_v;
1777  initial_guess_type M_composite_initial_guess;
1778  mutable boost::shared_ptr<CRBModel<ModelType> > M_crb_model;
1779 };
1780 
1781 
1782 
1783 
1784 template<typename TruthModelType>
1785 void
1786 CRBModel<TruthModelType>::initB()
1787 {
1788 
1789  //the matrix associated with H1 scalar product is now given by the model
1790  M_B = M_model->innerProduct();
1791 #if 0
1792  LOG(INFO) << "[CRBModel::initB] initialize scalar product\n";
1793  M_B = M_backend->newMatrix( M_model->functionSpace(), M_model->functionSpace() );
1794  using namespace Feel::vf;
1795  typename functionspace_type::element_type u( M_model->functionSpace() );
1796  form2( M_model->functionSpace(), M_model->functionSpace(), M_B, _init=true ) =
1797  integrate( elements( M_model->functionSpace()->mesh() ),
1798  gradt( u )*trans( grad( u ) ) );
1799 
1800  M_B->close();
1801 
1802  auto M = M_backend->newMatrix( M_model->functionSpace(), M_model->functionSpace() );
1803  form2( M_model->functionSpace(), M_model->functionSpace(), M, _init=true ) =
1804  integrate( elements( M_model->functionSpace()->mesh() ),
1805  idt( u )*id( u ) );
1806  M_B->printMatlab( "ipB.m" );
1807  M->printMatlab( "ipM.m" );
1808  M->close();
1809  LOG(INFO) << "[CRBModel::initB] starting eigen solve\n";
1810 #if 0
1811  SolverEigen<double>::eigenmodes_type modesmin=
1812  eigs( _matrixA=M_B,
1813  _matrixB=M,
1814  _problem=( EigenProblemType )GHEP,
1815  _solver=( EigenSolverType )M_vm["solvereigen-solver-type"].as<int>(),
1816  //_spectrum=LARGEST_MAGNITUDE,
1817  _spectrum=SMALLEST_MAGNITUDE,
1818  //_transform=SINVERT,
1819  _ncv=M_vm["solvereigen-ncv"].as<int>(),
1820  _nev=M_vm["solvereigen-nev"].as<int>(),
1821  _tolerance=M_vm["solvereigen-tol"].as<double>(),
1822  _maxit=M_vm["solvereigen-maxiter"].as<int>()
1823  );
1824  double eigmin = 1;
1825 
1826  if ( modesmin.empty() || modesmin.begin()->second.get<0>()<1e-6 )
1827  {
1828  LOG(INFO) << "coercivity constant not computed, taking 1\n";
1829  }
1830 
1831  else
1832  {
1833  eigmin = modesmin.begin()->second.get<0>();
1834  }
1835 
1836  LOG(INFO) << "[CRBModel::initB] coercivity constant (tau) = " << eigmin << "\n";
1837 #else
1838  double eigmin = 1;
1839 #endif
1840  M_B->addMatrix( eigmin, M );
1841 
1842 #endif
1843 }
1844 
1845 
1846 //create a vector of preassemble objects
1847 //for the mass matrix
1848 
1849 template<typename TruthModelType>
1850 typename CRBModel<TruthModelType>::operatorcomposite_ptrtype
1851 CRBModel<TruthModelType>::preAssembleMassMatrix() const
1852 {
1853  static const bool is_composite = functionspace_type::is_composite;
1854  return preAssembleMassMatrix( mpl::bool_< is_composite >() );
1855 }
1856 
1857 template<typename TruthModelType>
1858 typename CRBModel<TruthModelType>::operatorcomposite_ptrtype
1859 CRBModel<TruthModelType>::preAssembleMassMatrix( mpl::bool_<false> ) const
1860 {
1861 
1862  auto Xh = M_model->functionSpace();
1863  auto mesh = Xh->mesh();
1864 
1865  auto expr=integrate( _range=elements( mesh ) , _expr=idt( u )*id( v ) );
1866  auto op_mass = opLinearComposite( _domainSpace=Xh , _imageSpace=Xh );
1867  auto opfree = opLinearFree( _domainSpace=Xh , _imageSpace=Xh , _expr=expr );
1868  opfree->setName("mass operator (automatically created)");
1869  //in this case, the affine decompositon has only one element
1870  op_mass->addElement( boost::make_tuple(0,0) , opfree );
1871  return op_mass;
1872 }
1873 
1874 template<typename TruthModelType>
1875 typename CRBModel<TruthModelType>::operatorcomposite_ptrtype
1876 CRBModel<TruthModelType>::preAssembleMassMatrix( mpl::bool_<true> ) const
1877 {
1878  auto Xh = M_model->functionSpace();
1879 
1880  index_vector_type index_vector;
1881  PreAssembleMassMatrixInCompositeCase<TruthModelType> preassemble_mass_matrix_in_composite_case ( u , v );
1882  fusion::for_each( index_vector, preassemble_mass_matrix_in_composite_case );
1883 
1884  auto op_mass = preassemble_mass_matrix_in_composite_case.opmass();
1885  return op_mass;
1886 }
1887 
1888 
1889 template<typename TruthModelType>
1890 void
1891 CRBModel<TruthModelType>::assembleMassMatrix( void )
1892 {
1893  static const bool is_composite = functionspace_type::is_composite;
1894  return assembleMassMatrix( mpl::bool_< is_composite >() );
1895 }
1896 
1897 template<typename TruthModelType>
1898 void
1899 CRBModel<TruthModelType>::assembleMassMatrix( mpl::bool_<false> )
1900 {
1901  using namespace Feel::vf;
1902  auto Xh = M_model->functionSpace();
1903  M_Mqm.resize( 1 );
1904  M_Mqm[0].resize( 1 );
1905  M_Mqm[0][0] = M_backend->newMatrix( _test=Xh , _trial=Xh );
1906  auto mesh = Xh->mesh();
1907  form2( _test=Xh, _trial=Xh, _matrix=M_Mqm[0][0] ) =
1908  integrate( _range=elements( mesh ), _expr=idt( u )*id( v ) );
1909  M_Mqm[0][0]->close();
1910 }
1911 
1912 template<typename TruthModelType>
1913 void
1914 CRBModel<TruthModelType>::assembleMassMatrix( mpl::bool_<true> )
1915 {
1916  auto Xh = M_model->functionSpace();
1917 
1918  index_vector_type index_vector;
1919  int size = functionspace_type::nSpaces;
1920  M_Mqm.resize( 1 );
1921  M_Mqm[0].resize(1);
1922  M_Mqm[0][0]=M_backend->newMatrix( _test=Xh , _trial=Xh );
1923 
1924  AssembleMassMatrixInCompositeCase<TruthModelType> assemble_mass_matrix_in_composite_case ( u , v , this->shared_from_this());
1925  fusion::for_each( index_vector, assemble_mass_matrix_in_composite_case );
1926 
1927  M_Mqm[0][0]->close();
1928 }
1929 
1930 
1931 
1932 template<typename TruthModelType>
1933 void
1934 CRBModel<TruthModelType>::assembleInitialGuessV( initial_guess_type & initial_guess )
1935 {
1936  static const bool is_composite = functionspace_type::is_composite;
1937  return assembleInitialGuessV( initial_guess, mpl::bool_< is_composite >() );
1938 }
1939 
1940 template<typename TruthModelType>
1941 void
1942 CRBModel<TruthModelType>::assembleInitialGuessV( initial_guess_type & initial_guess, mpl::bool_<true> )
1943 {
1944  auto Xh = M_model->functionSpace();
1945  auto mesh = Xh->mesh();
1946 
1947  int q_max= this->QInitialGuess();
1948  M_InitialGuessV.resize( q_max );
1949  M_InitialGuessVector.resize( q_max );
1950  for(int q = 0; q < q_max; q++ )
1951  {
1952  int m_max= this->mMaxInitialGuess(q);
1953  M_InitialGuessV[q].resize( m_max );
1954  M_InitialGuessVector[q].resize( m_max );
1955  for(int m = 0; m < m_max; m++ )
1956  {
1957  M_InitialGuessV[q][m] = Xh->elementPtr();
1958  M_InitialGuessVector[q][m] = M_model->newVector();
1959  }
1960  }
1961 
1962  index_vector_type index_vector;
1963  AssembleInitialGuessVInCompositeCase<TruthModelType> assemble_initial_guess_v_in_composite_case ( v , initial_guess , this->shared_from_this());
1964  fusion::for_each( index_vector, assemble_initial_guess_v_in_composite_case );
1965 
1966  for(int q = 0; q < q_max; q++ )
1967  {
1968  int m_max = this->mMaxInitialGuess(q) ;
1969  for(int m = 0; m < m_max; m++ )
1970  M_InitialGuessVector[q][m]->close();
1971  }
1972 
1973 }
1974 
1975 template<typename TruthModelType>
1976 void
1977 CRBModel<TruthModelType>::assembleInitialGuessV( initial_guess_type & initial_guess, mpl::bool_<false> )
1978 {
1979  using namespace Feel::vf;
1980  auto Xh = M_model->functionSpace();
1981  auto mesh = Xh->mesh();
1982 
1983  int q_max= this->QInitialGuess();
1984  M_InitialGuessV.resize( q_max );
1985  M_InitialGuessVector.resize( q_max );
1986  for(int q = 0; q < q_max; q++ )
1987  {
1988  int m_max= this->mMaxInitialGuess(q);
1989  M_InitialGuessV[q].resize( m_max );
1990  M_InitialGuessVector[q].resize( m_max );
1991  for(int m = 0; m < m_max; m++ )
1992  {
1993  M_InitialGuessV[q][m] = Xh->elementPtr();
1994  M_InitialGuessVector[q][m] = M_model->newVector();
1995  form1( _test=Xh, _vector=M_InitialGuessVector[q][m]) =
1996  integrate( _range=elements( mesh ), _expr=idv( initial_guess[q][m] )*id( v ) );
1997  M_InitialGuessVector[q][m]->close();
1998  }
1999  }
2000 }
2001 
2002 
2003 template<typename TruthModelType>
2004 typename CRBModel<TruthModelType>::element_ptrtype
2005 CRBModel<TruthModelType>::assembleInitialGuess( parameter_type const& mu )
2006 {
2007  auto Xh = M_model->functionSpace();
2008  element_ptrtype initial_guess = Xh->elementPtr();
2009  initial_guess_type vector_initial_guess;
2010  beta_vector_type beta;
2011  vector_initial_guess = M_model->computeInitialGuessAffineDecomposition();
2012  beta = M_model->computeBetaInitialGuess( mu );
2013 
2014  int q_max = vector_initial_guess.size();
2015  for ( size_type q = 0; q < q_max; ++q )
2016  {
2017  int m_max=vector_initial_guess[q].size();
2018  for ( size_type m = 0; m < m_max ; ++m )
2019  {
2020  element_type temp = Xh->element();
2021  temp = *vector_initial_guess[q][m];
2022  temp.scale( beta[q][m] );
2023  *initial_guess += temp;
2024  }
2025  }
2026  return initial_guess;
2027 }
2028 
2029 
2030 template<typename TruthModelType>
2031 typename CRBModel<TruthModelType>::offline_merge_type
2032 CRBModel<TruthModelType>::offlineMergeOnFly(betaqm_type const& all_beta, parameter_type const& mu )
2033 {
2034  //recovery from the model of free operators Aqm in A = \sum_q \sum_m \beta_qm( u ; mu ) A_qm( u, v )
2035  auto compositeA = this->operatorCompositeA();
2036  //idem with other operator / functional
2037  auto compositeM = this->operatorCompositeM();
2038  auto vector_compositeF = this->functionalCompositeF();
2039 
2040  //acces to beta coefficients
2041  auto beta_M = all_beta.template get<0>();
2042  auto beta_A = all_beta.template get<1>();
2043  //warning : beta_F is a vector of beta_coefficients
2044  auto beta_F = all_beta.template get<2>();
2045 
2046  //associate beta coefficients to operators
2047  compositeA->setScalars( beta_A );
2048  compositeM->setScalars( beta_M );
2049 
2050  //merge
2051  auto A = M_model->newMatrix();
2052  auto M = M_model->newMatrix();
2053  compositeA->sumAllMatrices( A );
2054  //auto A = compositeA->sumAllMatrices();
2055  //auto M = compositeM->sumAllMatrices();
2056  compositeM->sumAllMatrices( M );
2057 
2058  std::vector<vector_ptrtype> F( Nl() );
2059 
2060  for(int output=0; output<Nl(); output++)
2061  {
2062  auto compositeF = vector_compositeF[output];
2063  compositeF->setScalars( beta_F[output] );
2064  F[output] = M_model->newVector();
2065  compositeF->sumAllVectors( F[output] );
2066  }
2067 
2068  return boost::make_tuple( M, A, F );
2069 }
2070 
2071 
2072 template<typename TruthModelType>
2073 typename CRBModel<TruthModelType>::offline_merge_type
2074 CRBModel<TruthModelType>::offlineMerge( betaqm_type const& all_beta , parameter_type const& mu )
2075 {
2076 
2077 #if 0
2078  sparse_matrix_ptrtype A( M_backend->newMatrix(
2079  _test=M_model->functionSpace(),
2080  _trial=M_model->functionSpace()
2081  ) );
2082 
2083  sparse_matrix_ptrtype M( M_backend->newMatrix(
2084  _test=M_model->functionSpace(),
2085  _trial=M_model->functionSpace()
2086  ) );
2087 
2088 #else
2089 
2090  auto A = this->newMatrix();
2091  auto M = this->newMatrix();
2092  //size_type pattern = operatorCompositeA()->operatorlinear(0,0)->pattern();
2093 
2094 #endif
2095  std::vector<vector_ptrtype> F( Nl() );
2096 
2097  //acces to beta coefficients
2098  auto beta_M = all_beta.template get<0>();
2099  auto beta_A = all_beta.template get<1>();
2100  //warning : beta_F is a vector of beta_coefficients
2101  auto beta_F = all_beta.template get<2>();
2102 
2103  A->zero();
2104  for ( size_type q = 0; q < Qa(); ++q )
2105  {
2106  for(size_type m = 0; m < mMaxA(q); ++m )
2107  A->addMatrix( beta_A[q][m], M_Aqm[q][m] );
2108  }
2109 
2110  if( Qm() > 0 )
2111  {
2112  for ( size_type q = 0; q < Qm(); ++q )
2113  {
2114  for(size_type m = 0; m < mMaxM(q) ; ++m )
2115  M->addMatrix( beta_M[q][m] , M_Mqm[q][m] );
2116  }
2117  }
2118 
2119  for ( size_type l = 0; l < Nl(); ++l )
2120  {
2121  F[l] = M_backend->newVector( M_model->functionSpace() );
2122  F[l]->zero();
2123 
2124  for ( size_type q = 0; q < Ql( l ); ++q )
2125  {
2126  for ( size_type m = 0; m < mMaxF(l,q); ++m )
2127  F[l]->add( beta_F[l][q][m] , M_Fqm[l][q][m] );
2128  }
2129  F[l]->close();
2130  }
2131 
2132  return boost::make_tuple( M, A, F );
2133 }
2134 
2135 
2136 template<typename TruthModelType>
2137 typename CRBModel<TruthModelType>::element_type
2138 CRBModel<TruthModelType>::solveFemUsingOfflineEim( parameter_type const& mu )
2139 {
2140  auto Xh= this->functionSpace();
2141 
2142  bdf_ptrtype mybdf;
2143  mybdf = bdf( _space=Xh, _vm=this->vm() , _name="mybdf" );
2144  sparse_matrix_ptrtype A;
2145  sparse_matrix_ptrtype M;
2146  std::vector<vector_ptrtype> F;
2147  element_ptrtype InitialGuess = Xh->elementPtr();
2148  vector_ptrtype Rhs( M_backend->newVector( Xh ) );
2149  auto u = Xh->element();
2150 
2151  double time_initial;
2152  double time_step;
2153  double time_final;
2154 
2155  if ( this->isSteady() )
2156  {
2157  time_initial=0;
2158  time_step = 1e30;
2159  time_final = 1e30;
2160  }
2161  else
2162  {
2163  time_initial=this->timeInitial();
2164  time_step=this->timeStep();
2165  time_final=this->timeFinal();
2166  this->initializationField( InitialGuess, mu );
2167  }
2168 
2169  mybdf->setTimeInitial( time_initial );
2170  mybdf->setTimeStep( time_step );
2171  mybdf->setTimeFinal( time_final );
2172 
2173  double bdf_coeff ;
2174  auto vec_bdf_poly = M_backend->newVector( Xh );
2175 
2176  for( mybdf->start(*InitialGuess); !mybdf->isFinished(); mybdf->next() )
2177  {
2178  bdf_coeff = mybdf->polyDerivCoefficient( 0 );
2179  auto bdf_poly = mybdf->polyDeriv();
2180  *vec_bdf_poly = bdf_poly;
2181  boost::tie(M, A, F) = this->update( mu , mybdf->time() );
2182  *Rhs = *F[0];
2183  if( !isSteady() )
2184  {
2185  A->addMatrix( bdf_coeff, M );
2186  Rhs->addVector( *vec_bdf_poly, *M );
2187  }
2188  M_backend->solve( _matrix=A , _solution=u, _rhs=Rhs );
2189  mybdf->shiftRight(u);
2190  }
2191 
2192  return u;
2193 
2194 }
2195 
2196 template<typename TruthModelType>
2197 typename CRBModel<TruthModelType>::element_type
2198 CRBModel<TruthModelType>::solveFemUsingAffineDecompositionFixedPoint( parameter_type const& mu )
2199 {
2200  auto Xh= this->functionSpace();
2201 
2202  bdf_ptrtype mybdf;
2203  mybdf = bdf( _space=Xh, _vm=this->vm() , _name="mybdf" );
2204  sparse_matrix_ptrtype A;
2205  sparse_matrix_ptrtype M;
2206  std::vector<vector_ptrtype> F;
2207  element_ptrtype InitialGuess = Xh->elementPtr();
2208  auto u = Xh->element();
2209  auto uold = Xh->element();
2210  vector_ptrtype Rhs( M_backend->newVector( Xh ) );
2211 
2212  double time_initial;
2213  double time_step;
2214  double time_final;
2215 
2216  if ( this->isSteady() )
2217  {
2218  time_initial=0;
2219  time_step = 1e30;
2220  time_final = 1e30;
2221  //we want to have the initial guess given by function update
2222  InitialGuess = this->assembleInitialGuess( mu ) ;
2223  }
2224  else
2225  {
2226  time_initial=this->timeInitial();
2227  time_step=this->timeStep();
2228  time_final=this->timeFinal();
2229  this->initializationField( InitialGuess, mu );
2230  }
2231 
2232  mybdf->setTimeInitial( time_initial );
2233  mybdf->setTimeStep( time_step );
2234  mybdf->setTimeFinal( time_final );
2235 
2236  u=*InitialGuess;
2237  double norm=0;
2238  int iter=0;
2239 
2240  double bdf_coeff ;
2241  auto vec_bdf_poly = M_backend->newVector( Xh );
2242 
2243  int max_fixedpoint_iterations = this->vm()["crb.max-fixedpoint-iterations"].template as<int>();
2244  double increment_fixedpoint_tol = this->vm()["crb.increment-fixedpoint-tol"].template as<double>();
2245  for( mybdf->start(*InitialGuess); !mybdf->isFinished(); mybdf->next() )
2246  {
2247  iter=0;
2248  bdf_coeff = mybdf->polyDerivCoefficient( 0 );
2249  auto bdf_poly = mybdf->polyDeriv();
2250  *vec_bdf_poly = bdf_poly;
2251  do {
2252  boost::tie(M, A, F) = this->update( mu , u , mybdf->time() );
2253  *Rhs = *F[0];
2254 
2255  if( !isSteady() )
2256  {
2257  A->addMatrix( bdf_coeff, M );
2258  Rhs->addVector( *vec_bdf_poly, *M );
2259  }
2260  uold = u;
2261  M_preconditioner->setMatrix( A );
2262  M_backend->solve( _matrix=A , _solution=u, _rhs=Rhs , _prec=M_preconditioner);
2263  norm = this->computeNormL2( uold , u );
2264  iter++;
2265  } while( norm > increment_fixedpoint_tol && iter<max_fixedpoint_iterations );
2266  mybdf->shiftRight(u);
2267  }
2268  return u;
2269 }
2270 
2271 template<typename TruthModelType>
2272 typename CRBModel<TruthModelType>::element_type
2273 CRBModel<TruthModelType>::solveFemDualUsingAffineDecompositionFixedPoint( parameter_type const& mu )
2274 {
2275  int output_index = option(_name="crb.output-index").template as<int>();
2276 
2277  auto Xh= this->functionSpace();
2278 
2279  bdf_ptrtype mybdf;
2280  mybdf = bdf( _space=Xh, _vm=this->vm() , _name="mybdf" );
2281  sparse_matrix_ptrtype A,Adu;
2282  sparse_matrix_ptrtype M;
2283  std::vector<vector_ptrtype> F;
2284  auto udu = Xh->element();
2285  auto uold = Xh->element();
2286  vector_ptrtype Rhs( M_backend->newVector( Xh ) );
2287  auto dual_initial_field = Xh->elementPtr();
2288 
2289  double time_initial;
2290  double time_step;
2291  double time_final;
2292 
2293  if ( this->isSteady() )
2294  {
2295  time_initial=0;
2296  time_step = 1e30;
2297  time_final = 1e30;
2298  //InitialGuess = this->assembleInitialGuess( mu ) ;
2299  }
2300  else
2301  {
2302  time_initial=this->timeFinal()+this->timeStep();
2303  time_step=-this->timeStep();
2304  time_final=this->timeInitial()+this->timeStep();
2305  }
2306 
2307  mybdf->setTimeInitial( time_initial );
2308  mybdf->setTimeStep( time_step );
2309  mybdf->setTimeFinal( time_final );
2310 
2311  double norm=0;
2312  int iter=0;
2313 
2314  double bdf_coeff ;
2315  auto vec_bdf_poly = M_backend->newVector( Xh );
2316 
2317  if ( this->isSteady() )
2318  udu.zero() ;
2319  else
2320  {
2321  boost::tie( M, A, F) = this->update( mu , mybdf->timeInitial() );
2322  *Rhs=*F[output_index];
2323  M_preconditioner->setMatrix( M );
2324  M_backend->solve( _matrix=M, _solution=dual_initial_field, _rhs=Rhs, _prec=M_preconditioner );
2325  udu=*dual_initial_field;
2326  }
2327 
2328 
2329  int max_fixedpoint_iterations = this->vm()["crb.max-fixedpoint-iterations"].template as<int>();
2330  double increment_fixedpoint_tol = this->vm()["crb.increment-fixedpoint-tol"].template as<double>();
2331  for( mybdf->start(udu); !mybdf->isFinished(); mybdf->next() )
2332  {
2333  iter=0;
2334  bdf_coeff = mybdf->polyDerivCoefficient( 0 );
2335  auto bdf_poly = mybdf->polyDeriv();
2336  *vec_bdf_poly = bdf_poly;
2337  do {
2338  boost::tie(M, A, F) = this->update( mu , udu , mybdf->time() );
2339 
2340  if( ! isSteady() )
2341  {
2342  A->addMatrix( bdf_coeff, M );
2343  Rhs->zero();
2344  *vec_bdf_poly = bdf_poly;
2345  Rhs->addVector( *vec_bdf_poly, *M );
2346  }
2347  else
2348  {
2349  *Rhs = *F[output_index];
2350  Rhs->scale( -1 );
2351  }
2352 
2353  if( option("crb.use-symmetric-matrix").template as<bool>() )
2354  Adu = A;
2355  else
2356  A->transpose( Adu );
2357 
2358  uold = udu;
2359  M_preconditioner->setMatrix( Adu );
2360  M_backend->solve( _matrix=Adu , _solution=udu, _rhs=Rhs , _prec=M_preconditioner);
2361  norm = this->computeNormL2( uold , udu );
2362  iter++;
2363  } while( norm > increment_fixedpoint_tol && iter<max_fixedpoint_iterations );
2364  mybdf->shiftRight(udu);
2365  }
2366  return udu;
2367 }
2368 
2369 
2370 template<typename TruthModelType>
2371 void
2372 CRBModel<TruthModelType>::run( const double * X, unsigned long N, double * Y, unsigned long P )
2373 {
2374  M_model->run( X, N, Y, P );
2375 }
2376 
2377 
2378 
2379 }
2380 #endif /* __CRBModel_H */

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