29 #ifndef __CRBModelTrilinear_H
30 #define __CRBModelTrilinear_H 1
32 #include <boost/shared_ptr.hpp>
61 template<
typename ModelType>
71 static const uint16_type ParameterSpaceDimension = ModelType::ParameterSpaceDimension;
72 static const bool is_time_dependent = ModelType::is_time_dependent;
82 typedef boost::shared_ptr<ModelType> model_ptrtype;
97 typedef typename model_type::functionspace_ptrtype functionspace_ptrtype;
101 typedef typename model_type::element_ptrtype element_ptrtype;
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;
109 typedef typename model_type::eigen_matrix_type eigen_matrix_type;
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;
119 typedef Eigen::VectorXd vectorN_type;
120 typedef std::vector< std::vector< double > > beta_vector_type;
122 typedef typename boost::tuple<sparse_matrix_ptrtype,
123 sparse_matrix_ptrtype,
124 std::vector<vector_ptrtype>
125 > offline_merge_type;
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;
135 typedef typename boost::tuple< beta_vector_type,
137 std::vector<beta_vector_type>
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;
159 M_is_initialized( false ),
160 M_mode( CRBModelMode::PFEM ),
162 M_backend( backend_type::build( BACKEND_PETSC ) )
171 M_is_initialized( false ),
187 M_is_initialized( false ),
189 M_mode( CRBModelMode::PFEM ),
191 M_backend( backend_type::build( model->
vm ) )
203 M_is_initialized( o.M_is_initialized ),
206 M_model( o.M_model ),
207 M_backend( o.M_backend )
219 if ( M_is_initialized )
222 M_is_initialized=
true;
224 if( ! M_model->isInitialized() )
226 LOG( INFO ) <<
"CRBModel Model is not initialized";
227 M_model->initModel();
228 M_model->setInitialized(
true );
231 M_is_initialized=
true;
248 M_backend = o.M_backend;
262 po::variables_map
vm()
const
273 return M_model->newMatrix();
282 return M_model->newVector();
288 return M_model->functionSpace();
294 return M_model->Qa()-M_model->QaTri();
301 return M_model->QaTri();
304 sparse_matrix_ptrtype computeTrilinearForm(
element_type const& xi )
306 return M_model->computeTrilinearForm( xi );
309 sparse_matrix_ptrtype jacobian(
element_type const& xi )
311 return M_model->jacobian( xi );
316 return M_model->residual( xi );
320 vectorN_type computeStatistics ( Eigen::VectorXd vector , std::string name )
322 return M_model->computeStatistics( vector , name );
329 return M_model->Nl();
335 return M_model->Ql( l );
341 return M_model->parameterSpace();
344 parameter_type refParameter()
346 return M_model->refParameter();
361 M_model->setMeshSize( s );
376 return computeBetaQm( mu , mpl::bool_<model_type::is_time_dependent>(), time );
378 betaqm_type
computeBetaQm( parameter_type
const& mu , mpl::bool_<true>,
double time=0 )
380 return M_model->computeBetaQm( mu , time );
382 betaqm_type
computeBetaQm( parameter_type
const& mu , mpl::bool_<false>,
double time=0)
384 beta_vector_type betaAqm;
385 beta_vector_type betaMqm;
386 std::vector<beta_vector_type> betaFqm;
389 std::vector<beta_vector_type> >
392 model_beta = M_model->computeBetaQm( mu );
393 betaAqm = model_beta.get<0>();
394 betaFqm = model_beta.get<1>();
396 return boost::make_tuple( betaMqm, betaAqm, betaFqm );
410 std::vector< std::vector<sparse_matrix_ptrtype> > mass;
411 boost::tie(
M_Aqm,
M_Fqm ) = M_model->computeAffineDecomposition();
413 return boost::make_tuple( mass,
M_Aqm,
M_Fqm );
426 sparse_matrix_ptrtype
Aqm( uint16_type q, uint16_type m,
bool transpose =
false )
const
429 return M_Aqm[q][m]->transpose();
447 return M_Aqm[q][m]->energy( xi_j, xi_i, transpose );
455 vector_ptrtype
Fqm( uint16_type l, uint16_type q,
int m )
const
457 return M_Fqm[l][q][m];
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);
490 value_type Fqm( uint16_type l, uint16_type q, uint16_type m, element_ptrtype
const& xi )
500 return M_model->scalarProduct( X, Y );
507 return M_model->scalarProduct( X, Y );
515 return M_model->scalarProductForMassMatrix( X, Y );
522 return M_model->scalarProductForMassMatrix( X, Y );
534 double scalarProductForPod( vector_ptrtype
const& X, vector_ptrtype
const& Y , mpl::bool_<true> )
536 return M_model->scalarProductForPod( X, Y );
538 double scalarProductForPod( vector_ptrtype
const& X, vector_ptrtype
const& Y , mpl::bool_<false> )
549 return M_model->solve( mu );
557 void l2solve( vector_ptrtype& u, vector_ptrtype
const& f )
559 return M_model->l2solve( u, f );
566 void run(
const double * X,
unsigned long N,
double * Y,
unsigned long P );
574 return M_model->output( output_index, mu , u , need_to_solve);
581 return timeStep( mpl::bool_<model_type::is_time_dependent>() );
583 double timeStep( mpl::bool_<true> )
587 if ( M_model->isSteady() )
590 timestep = M_model->timeStep();
593 double timeStep( mpl::bool_<false> )
600 return timeInitial( mpl::bool_<model_type::is_time_dependent>() );
602 double timeInitial( mpl::bool_<true> )
604 return M_model->timeInitial();
606 double timeInitial( mpl::bool_<false> )
613 return timeFinal( mpl::bool_<model_type::is_time_dependent>() );
615 double timeFinal( mpl::bool_<true> )
618 if ( M_model->isSteady() )
621 timefinal = M_model->timeFinal();
624 double timeFinal( mpl::bool_<false> )
631 element_type solveFemUsingAffineDecompositionFixedPoint( parameter_type
const& mu )
636 element_type solveFemDualUsingAffineDecompositionFixedPoint( parameter_type
const& mu )
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 )
646 return result_offline_merge_type ;
648 sparse_matrix_ptrtype
const& innerProduct() const
650 return sparse_matrix ;
652 sparse_matrix_ptrtype
const& innerProductForMassMatrix() const
654 return sparse_matrix ;
657 sparse_matrix_ptrtype Mqm( uint16_type q, uint16_type m,
bool transpose =
false )
const
659 return sparse_matrix;
681 int mMaxF(
int output_index,
int q )
689 return isSteady( mpl::bool_<model_type::is_time_dependent>() );
691 bool isSteady( mpl::bool_<true> )
693 return M_model->isSteady();
695 bool isSteady( mpl::bool_<false> )
705 return M_model->scalarContinuousEim();
713 return M_model->scalarDiscontinuousEim();
720 std::vector< std::vector<sparse_matrix_ptrtype> >
M_Aqm;
723 std::vector< std::vector<std::vector<vector_ptrtype> > >
M_Fqm;
727 bool M_is_initialized;
730 po::variables_map M_vm;
736 model_ptrtype M_model;
738 backend_ptrtype M_backend;
740 sparse_matrix_ptrtype sparse_matrix;
747 template<
typename TruthModelType>
751 M_model->run( X, N, Y, P );