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
crbmodeltrilinear.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 __CRBModelTrilinear_H
30 #define __CRBModelTrilinear_H 1
31 
32 #include <boost/shared_ptr.hpp>
33 
34 #include <vector>
35 
36 
38 #include <feel/feelvf/vf.hpp>
39 
41 
42 
43 namespace Feel
44 {
45 
61 template<typename ModelType>
63 {
64 public:
65 
66 
70 
71  static const uint16_type ParameterSpaceDimension = ModelType::ParameterSpaceDimension;
72  static const bool is_time_dependent = ModelType::is_time_dependent;
73 
75 
79 
81  typedef ModelType model_type;
82  typedef boost::shared_ptr<ModelType> model_ptrtype;
83 
85  typedef typename model_type::value_type value_type;
87  typedef typename ModelType::mesh_type mesh_type;
88 
90  typedef typename ModelType::mesh_ptrtype mesh_ptrtype;
91 
93  typedef typename ModelType::space_type space_type;
94 
96  typedef typename model_type::functionspace_type functionspace_type;
97  typedef typename model_type::functionspace_ptrtype functionspace_ptrtype;
98 
100  typedef typename model_type::element_type element_type;
101  typedef typename model_type::element_ptrtype element_ptrtype;
102 
103  typedef typename model_type::backend_type backend_type;
104  typedef boost::shared_ptr<backend_type> backend_ptrtype;
105  typedef typename model_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
106  typedef typename model_type::vector_ptrtype vector_ptrtype;
107  typedef typename model_type::vector_type vector_type;
108 
109  typedef typename model_type::eigen_matrix_type eigen_matrix_type;
110 
111  typedef typename model_type::parameterspace_type parameterspace_type;
112  typedef typename model_type::parameterspace_ptrtype parameterspace_ptrtype;
113  typedef typename model_type::parameter_type parameter_type;
114  typedef typename model_type::parameter_ptrtype parameter_ptrtype;
115  typedef typename model_type::sampling_type sampling_type;
116  typedef typename model_type::sampling_ptrtype sampling_ptrtype;
117 
118 
119  typedef Eigen::VectorXd vectorN_type;
120  typedef std::vector< std::vector< double > > beta_vector_type;
121 
122  typedef typename boost::tuple<sparse_matrix_ptrtype,
123  sparse_matrix_ptrtype,
124  std::vector<vector_ptrtype>
125  > offline_merge_type;
126 
127 
128 
129  typedef typename boost::tuple<std::vector< std::vector<sparse_matrix_ptrtype> >,
130  std::vector< std::vector<sparse_matrix_ptrtype> >,
131  std::vector< std::vector< std::vector<vector_ptrtype> > >
132  > affine_decomposition_type;
133 
134 
135  typedef typename boost::tuple< beta_vector_type,
136  beta_vector_type,
137  std::vector<beta_vector_type>
138  > betaqm_type;
139 
140 
141  static const int nb_spaces = functionspace_type::nSpaces;
142  typedef typename mpl::if_< boost::is_same< mpl::int_<nb_spaces> , mpl::int_<2> > , fusion::vector< mpl::int_<0>, mpl::int_<1> > ,
143  typename mpl::if_ < boost::is_same< mpl::int_<nb_spaces> , mpl::int_<3> > , fusion::vector < mpl::int_<0> , mpl::int_<1> , mpl::int_<2> > ,
144  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> >,
145  fusion::vector< mpl::int_<0>, mpl::int_<1>, mpl::int_<2>, mpl::int_<3>, mpl::int_<4> >
146  >::type >::type >::type index_vector_type;
147 
148 
150 
154 
156  :
157  M_Aqm(),
158  M_Fqm(),
159  M_is_initialized( false ),
160  M_mode( CRBModelMode::PFEM ),
161  M_model( new model_type() ),
162  M_backend( backend_type::build( BACKEND_PETSC ) )
163  {
164  this->init();
165  }
166 
167  CRBModelTrilinear( po::variables_map const& vm, CRBModelMode mode = CRBModelMode::PFEM )
168  :
169  M_Aqm(),
170  M_Fqm(),
171  M_is_initialized( false ),
172  M_vm( vm ),
173  M_mode( mode ),
174  M_model( new model_type( vm ) ),
175  M_backend( backend_type::build( vm ) )
176  {
177  this->init();
178  }
179 
183  CRBModelTrilinear( model_ptrtype & model )
184  :
185  M_Aqm(),
186  M_Fqm(),
187  M_is_initialized( false ),
188  M_vm(),
189  M_mode( CRBModelMode::PFEM ),
190  M_model( model ),
191  M_backend( backend_type::build( model->vm ) )
192  {
193  this->init();
194  }
195 
200  :
201  M_Aqm( o.M_Aqm ),
202  M_Fqm( o.M_Fqm ),
203  M_is_initialized( o.M_is_initialized ),
204  M_vm( o.M_vm ),
205  M_mode( o.M_mode ),
206  M_model( o.M_model ),
207  M_backend( o.M_backend )
208  {
209  this->init();
210  }
211 
214  {}
215 
217  FEELPP_DONT_INLINE void init()
218  {
219  if ( M_is_initialized )
220  return;
221 
222  M_is_initialized=true;
223 
224  if( ! M_model->isInitialized() )
225  {
226  LOG( INFO ) << "CRBModel Model is not initialized";
227  M_model->initModel();
228  M_model->setInitialized( true );
229  }
230 
231  M_is_initialized=true;
232  }
233 
235 
239 
242  {
243  if ( this != &o )
244  {
245  M_Aqm = o.M_Aqm;
246  M_Fqm = o.M_Fqm;
247  M_model = o.M_model;
248  M_backend = o.M_backend;
249  }
250 
251  return *this;
252  }
254 
258 
262  po::variables_map vm() const
263  {
264  return M_vm;
265  }
266 
271  virtual sparse_matrix_ptrtype newMatrix() const
272  {
273  return M_model->newMatrix();
274  }
275 
280  virtual vector_ptrtype newVector() const
281  {
282  return M_model->newVector();
283  }
284 
286  functionspace_ptrtype functionSpace() const
287  {
288  return M_model->functionSpace();
289  }
290 
292  size_type Qa() const
293  {
294  return M_model->Qa()-M_model->QaTri();
295  }
296 
297 
299  size_type QaTri() const
300  {
301  return M_model->QaTri();
302  }
303 
304  sparse_matrix_ptrtype computeTrilinearForm( element_type const& xi )
305  {
306  return M_model->computeTrilinearForm( xi );
307  }
308 
309  sparse_matrix_ptrtype jacobian( element_type const& xi )
310  {
311  return M_model->jacobian( xi );
312  }
313 
314  vector_ptrtype residual( element_type const& xi )
315  {
316  return M_model->residual( xi );
317  }
318 
319 
320  vectorN_type computeStatistics ( Eigen::VectorXd vector , std::string name )
321  {
322  return M_model->computeStatistics( vector , name );
323  }
324 
325 
327  size_type Nl() const
328  {
329  return M_model->Nl();
330  }
331 
333  size_type Ql( int l ) const
334  {
335  return M_model->Ql( l );
336  }
337 
339  parameterspace_ptrtype parameterSpace() const
340  {
341  return M_model->parameterSpace();
342  }
343 
344  parameter_type refParameter()
345  {
346  return M_model->refParameter();
347  }
348 
349 
351 
355 
359  void setMeshSize( double s )
360  {
361  M_model->setMeshSize( s );
362  }
363 
364 
366 
370 
374  betaqm_type computeBetaQm( parameter_type const& mu , double time=0 )
375  {
376  return computeBetaQm( mu , mpl::bool_<model_type::is_time_dependent>(), time );
377  }
378  betaqm_type computeBetaQm( parameter_type const& mu , mpl::bool_<true>, double time=0 )
379  {
380  return M_model->computeBetaQm( mu , time );
381  }
382  betaqm_type computeBetaQm( parameter_type const& mu , mpl::bool_<false>, double time=0)
383  {
384  beta_vector_type betaAqm;
385  beta_vector_type betaMqm;
386  std::vector<beta_vector_type> betaFqm;
387  boost::tuple<
388  beta_vector_type,
389  std::vector<beta_vector_type> >
390  model_beta;
391 
392  model_beta = M_model->computeBetaQm( mu );
393  betaAqm = model_beta.get<0>();
394  betaFqm = model_beta.get<1>();
395 
396  return boost::make_tuple( betaMqm, betaAqm, betaFqm );
397  }
398 
399 
400 
408  affine_decomposition_type computeAffineDecomposition()
409  {
410  std::vector< std::vector<sparse_matrix_ptrtype> > mass;
411  boost::tie( M_Aqm, M_Fqm ) = M_model->computeAffineDecomposition();
412  // to have compatibility with SCM, we need to provide the same interface than CRBModel
413  return boost::make_tuple( mass, M_Aqm, M_Fqm );
414  }
415 
416 
417 
426  sparse_matrix_ptrtype Aqm( uint16_type q, uint16_type m, bool transpose = false ) const
427  {
428  if ( transpose )
429  return M_Aqm[q][m]->transpose();
430 
431  return M_Aqm[q][m];
432  }
433 
434 
445  value_type Aqm( uint16_type q, uint16_type m, element_type const& xi_i, element_type const& xi_j, bool transpose = false )
446  {
447  return M_Aqm[q][m]->energy( xi_j, xi_i, transpose );
448  }
449 
455  vector_ptrtype Fqm( uint16_type l, uint16_type q, int m ) const
456  {
457  return M_Fqm[l][q][m];
458  }
459 
471  value_type Fqm( uint16_type l, uint16_type q, uint16_type m, element_type const& xi )
472  {
473  element_ptrtype eltF( new element_type( M_model->functionSpace() ) );
474  for(int i=0; i<eltF->localSize();i++)
475  eltF->operator()(i)=M_Fqm[l][q][m]->operator()(i);
476  return inner_product( *eltF , xi );
477  }
478 
490  value_type Fqm( uint16_type l, uint16_type q, uint16_type m, element_ptrtype const& xi )
491  {
492  return inner_product( M_Fqm[l][q][m], xi );
493  }
494 
498  double scalarProduct( vector_type const& X, vector_type const& Y )
499  {
500  return M_model->scalarProduct( X, Y );
501  }
505  double scalarProduct( vector_ptrtype const& X, vector_ptrtype const& Y )
506  {
507  return M_model->scalarProduct( X, Y );
508  }
509 
513  double scalarProductForMassMatrix( vector_type const& X, vector_type const& Y )
514  {
515  return M_model->scalarProductForMassMatrix( X, Y );
516  }
520  double scalarProductForMassMatrix( vector_ptrtype const& X, vector_ptrtype const& Y )
521  {
522  return M_model->scalarProductForMassMatrix( X, Y );
523  }
524 
525 
526 
530  double scalarProductForPod( vector_ptrtype const& X, vector_ptrtype const& Y )
531  {
532  return scalarProductForPod( X, Y , mpl::bool_<model_type::is_time_dependent>() );
533  }
534  double scalarProductForPod( vector_ptrtype const& X, vector_ptrtype const& Y , mpl::bool_<true> )
535  {
536  return M_model->scalarProductForPod( X, Y );
537  }
538  double scalarProductForPod( vector_ptrtype const& X, vector_ptrtype const& Y , mpl::bool_<false> )
539  {
540  return 0;
541  }
542 
543 
547  element_type solve( parameter_type const& mu )
548  {
549  return M_model->solve( mu );
550  }
551 
552 
557  void l2solve( vector_ptrtype& u, vector_ptrtype const& f )
558  {
559  return M_model->l2solve( u, f );
560  }
561 
562 
566  void run( const double * X, unsigned long N, double * Y, unsigned long P );
567 
572  value_type output( int output_index, parameter_type const& mu, element_type &u, bool need_to_solve=false )
573  {
574  return M_model->output( output_index, mu , u , need_to_solve);
575  }
576 
577 
578 
579  double timeStep()
580  {
581  return timeStep( mpl::bool_<model_type::is_time_dependent>() );
582  }
583  double timeStep( mpl::bool_<true> )
584  {
585  double timestep;
586 
587  if ( M_model->isSteady() )
588  timestep=1e30;
589  else
590  timestep = M_model->timeStep();
591  return timestep;
592  }
593  double timeStep( mpl::bool_<false> )
594  {
595  return 1e30;
596  }
597 
598  double timeInitial()
599  {
600  return timeInitial( mpl::bool_<model_type::is_time_dependent>() );
601  }
602  double timeInitial( mpl::bool_<true> )
603  {
604  return M_model->timeInitial();
605  }
606  double timeInitial( mpl::bool_<false> )
607  {
608  return 0;
609  }
610 
611  double timeFinal()
612  {
613  return timeFinal( mpl::bool_<model_type::is_time_dependent>() );
614  }
615  double timeFinal( mpl::bool_<true> )
616  {
617  double timefinal;
618  if ( M_model->isSteady() )
619  timefinal=1e30;
620  else
621  timefinal = M_model->timeFinal();
622  return timefinal;
623  }
624  double timeFinal( mpl::bool_<false> )
625  {
626  return 1e30;
627  }
628 
629 
630  //only to compile
631  element_type solveFemUsingAffineDecompositionFixedPoint( parameter_type const& mu )
632  {
633  auto zero=this->functionSpace()->element();
634  return zero;
635  }
636  element_type solveFemDualUsingAffineDecompositionFixedPoint( parameter_type const& mu )
637  {
638  auto zero=this->functionSpace()->element();
639  return zero;
640  }
641 
642  element_type solveFemUsingOfflineEim( parameter_type const& mu ){};
643  offline_merge_type result_offline_merge_type;
644  offline_merge_type update( parameter_type const& mu, double time=0 ) // for scm
645  {
646  return result_offline_merge_type ;
647  }
648  sparse_matrix_ptrtype const& innerProduct() const // for the scm
649  {
650  return sparse_matrix ;
651  }
652  sparse_matrix_ptrtype const& innerProductForMassMatrix() const // for the scm
653  {
654  return sparse_matrix ;
655  }
656 
657  sparse_matrix_ptrtype Mqm( uint16_type q, uint16_type m, bool transpose = false ) const
658  {
659  return sparse_matrix;
660  }
661  value_type Mqm( uint16_type q, uint16_type m, element_type const& xi_i, element_type const& xi_j, bool transpose = false ) const
662  {
663  return 0;
664  }
665  size_type Qm() const
666  {
667  return 0;
668  }
669  int mMaxA( int q )
670  {
671  return 1;
672  }
673  int mMaxM( int q )
674  {
675  return 1;
676  }
677  int mMaxMF( int q )
678  {
679  return 1;
680  }
681  int mMaxF(int output_index, int q )
682  {
683  return 1;
684  }
685 
687  bool isSteady()
688  {
689  return isSteady( mpl::bool_<model_type::is_time_dependent>() );
690  }
691  bool isSteady( mpl::bool_<true> )
692  {
693  return M_model->isSteady();
694  }
695  bool isSteady( mpl::bool_<false> )
696  {
697  return true;
698  }
699 
703  typename model_type::funs_type scalarContinuousEim()
704  {
705  return M_model->scalarContinuousEim();
706  }
707 
711  typename model_type::funsd_type scalarDiscontinuousEim()
712  {
713  return M_model->scalarDiscontinuousEim();
714  }
715 
716 
717 protected:
718 
720  std::vector< std::vector<sparse_matrix_ptrtype> > M_Aqm;
721 
723  std::vector< std::vector<std::vector<vector_ptrtype> > > M_Fqm;
724 
725 private:
726 
727  bool M_is_initialized;
728 
730  po::variables_map M_vm;
731 
733  CRBModelMode M_mode;
734 
736  model_ptrtype M_model;
737 
738  backend_ptrtype M_backend;
739 
740  sparse_matrix_ptrtype sparse_matrix; // only to compile
741 
742 };
743 
744 
745 
746 
747 template<typename TruthModelType>
748 void
749 CRBModelTrilinear<TruthModelType>::run( const double * X, unsigned long N, double * Y, unsigned long P )
750 {
751  M_model->run( X, N, Y, P );
752 }
753 
754 
755 
756 }
757 #endif /* __CRBModelTrilinear_H */

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