26 #ifndef __ConvectionCrb_H
27 #define __ConvectionCrb_H 1
46 #include <feel/feeldiscr/functionspace.hpp>
50 #include <feel/feeldiscr/operatorlagrangep1.hpp>
52 #include <feel/feelfilters/exporter.hpp>
59 #include <Eigen/Dense>
61 #include <feel/feelcrb/modelcrbbase.hpp>
65 using namespace Feel::vf;
67 #if !defined( CONVECTION_DIM )
68 #define CONVECTION_DIM 2
70 #if !defined( CONVECTION_ORDER_U )
71 #define CONVECTION_ORDER_U 2
73 #if !defined( CONVECTION_ORDER_P )
74 #define CONVECTION_ORDER_P 1
76 #if !defined( CONVECTION_ORDER_T )
77 #define CONVECTION_ORDER_T 2
79 #if !defined( CRB_SOLVER )
84 class ParameterDefinition
87 static const uint16_type ParameterSpaceDimension = 2;
91 class FunctionSpaceDefinition
94 static const uint16_type Order = 1;
96 static const int Order_s = CONVECTION_ORDER_U;
97 static const int Order_p = CONVECTION_ORDER_P;
98 static const int Order_t = CONVECTION_ORDER_T;
105 typedef Lagrange<Order_s, Vectorial,Continuous,PointSetFekete> basis_u_type;
106 typedef Lagrange<Order_p, Scalar,Continuous,PointSetFekete> basis_p_type;
107 typedef Lagrange<Order_t, Scalar,Continuous,PointSetFekete> basis_t_type;
109 #if defined( FEELPP_USE_LM )
110 typedef Lagrange<0, Scalar> basis_l_type;
111 typedef bases< basis_u_type , basis_p_type , basis_t_type,basis_l_type> basis_type;
113 typedef bases< basis_u_type , basis_p_type , basis_t_type> basis_type;
129 class ConvectionCrb :
public ModelCrbBase< ParameterDefinition, FunctionSpaceDefinition >
133 typedef ModelCrbBase<ParameterDefinition,FunctionSpaceDefinition> super_type;
134 typedef typename super_type::funs_type funs_type;
135 typedef typename super_type::funsd_type funsd_type;
137 static const uint16_type Order = 1;
138 static const uint16_type ParameterSpaceDimension = 2;
139 static const bool is_time_dependent =
false;
142 static const int Order_s = CONVECTION_ORDER_U;
143 static const int Order_p = CONVECTION_ORDER_P;
144 static const int Order_t = CONVECTION_ORDER_T;
150 typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
153 typedef boost::shared_ptr<backend_type> backend_ptrtype;
159 typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
160 typedef backend_type::vector_ptrtype vector_ptrtype;
162 typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> eigen_matrix_type;
163 typedef eigen_matrix_type ematrix_type;
164 typedef boost::shared_ptr<eigen_matrix_type> eigen_matrix_ptrtype;
168 typedef Lagrange<Order_s, Vectorial,Continuous,PointSetFekete> basis_u_type;
169 typedef Lagrange<Order_p, Scalar,Continuous,PointSetFekete> basis_p_type;
170 typedef Lagrange<Order_t, Scalar,Continuous,PointSetFekete> basis_t_type;
173 typedef boost::shared_ptr<U_space_type> U_space_ptrtype;
175 typedef boost::shared_ptr<T_space_type> T_space_ptrtype;
179 typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
181 typedef parameterspace_type::element_ptrtype parameter_ptrtype;
183 typedef parameterspace_type::sampling_ptrtype sampling_ptrtype;
184 typedef std::vector< std::vector< double > > beta_vector_type;
186 #if defined( FEELPP_USE_LM )
187 typedef Lagrange<0, Scalar> basis_l_type;
188 typedef bases< basis_u_type , basis_p_type , basis_t_type,basis_l_type> basis_type;
190 typedef bases< basis_u_type , basis_p_type , basis_t_type> basis_type;
203 typedef boost::shared_ptr<space_type> space_ptrtype;
204 typedef typename space_type::element_type element_type;
205 typedef typename element_type:: sub_element<0>::type element_0_type;
207 typedef typename element_type:: sub_element<1>::type element_1_type;
208 typedef typename element_type:: sub_element<2>::type element_2_type;
209 #if defined( FEELPP_USE_LM )
210 typedef typename element_type:: sub_element<3>::type element_3_type;
214 typedef space_ptrtype functionspace_ptrtype;
216 typedef boost::shared_ptr<element_type> element_ptrtype;
219 typedef boost::shared_ptr<oplin_type> oplin_ptrtype;
220 typedef FsFunctionalLinear<space_type> funlin_type;
221 typedef boost::shared_ptr<funlin_type> funlin_ptrtype;
226 typedef boost::tuple<
227 std::vector< std::vector<sparse_matrix_ptrtype> >,
228 std::vector< std::vector<std::vector<vector_ptrtype> > >
229 > affine_decomposition_type;
231 typedef Eigen::MatrixXd matrixN_type;
238 Feel::gmsh_ptrtype createMesh();
244 std::string modelName()
246 std::ostringstream ostr;
247 ostr <<
"naturalconvection" ;
258 parameter_type refParameter()
263 affine_decomposition_type computeAffineDecomposition()
265 return boost::make_tuple( M_Aqm, M_Fqm );
285 int Ql(
int l )
const;
288 int mMaxF(
int output_index,
int q );
294 boost::tuple<beta_vector_type, std::vector<beta_vector_type> >
295 computeBetaQm( parameter_type
const& mu,
double time=0 ) ;
297 boost::tuple<beta_vector_type, std::vector<beta_vector_type> >
298 computeBetaQm(element_type
const& T, parameter_type
const& mu,
double time=0 )
300 return computeBetaQm( mu , time );
303 void update( parameter_type
const& mu );
305 void solve( sparse_matrix_ptrtype& D, element_type& u, vector_ptrtype& F );
312 void solve( parameter_type
const& mu, element_ptrtype& T );
317 element_type solve( parameter_type
const& mu );
322 void l2solve( vector_ptrtype& u, vector_ptrtype
const& f );
327 void exportResults( element_type& u );
328 void exportResults( element_ptrtype& U,
int t );
329 void exportResults( element_type& U,
double t );
335 double scalarProduct( vector_ptrtype
const& X, vector_ptrtype
const& Y );
350 void run(
const double * X,
unsigned long N,
double * Y,
unsigned long P );
356 value_type output(
int output_index, parameter_type
const& mu , element_type& unknown,
bool need_to_solve=
false);
358 sparse_matrix_ptrtype newMatrix()
const
360 return M_backend->newMatrix( Xh, Xh );
363 vector_ptrtype newVector()
const
365 return M_backend->newVector( Xh );
369 space_ptrtype functionSpace()
374 void setMeshSize(
double s )
380 po::options_description
const& optionsDescription()
const
391 po::variables_map
const&
vm()
const
398 sparse_matrix_ptrtype innerProduct()
403 void updateJacobianWithoutAffineDecomposition(
const vector_ptrtype& X, sparse_matrix_ptrtype& J );
404 void updateJacobian(
const vector_ptrtype& X, sparse_matrix_ptrtype& J );
405 void updateResidual(
const vector_ptrtype& X, vector_ptrtype& R );
406 sparse_matrix_ptrtype computeTrilinearForm(
const element_type& X );
407 sparse_matrix_ptrtype jacobian(
const element_type& X );
408 vector_ptrtype residual(
const element_type& X );
412 po::options_description M_desc;
414 po::variables_map M_vm;
415 backend_ptrtype M_backend;
418 boost::shared_ptr<OperatorLagrangeP1<typename space_type::sub_functionspace<2>::type::element_type> > P1h;
420 oplin_ptrtype M_oplin;
423 sparse_matrix_ptrtype M_L;
424 sparse_matrix_ptrtype M_D;
427 sparse_matrix_ptrtype D,M;
430 boost::shared_ptr<export_type> exporter;
433 std::map<std::string,std::pair<boost::timer,double> > timers;
435 std::vector <double> Grashofs;
436 double M_current_Grashofs;
437 double M_current_Prandtl;
445 std::vector< std::vector<sparse_matrix_ptrtype> > M_Aqm;
446 std::vector< std::vector<sparse_matrix_ptrtype> > M_Mqm;
447 std::vector< std::vector<std::vector<vector_ptrtype> > > M_Fqm;
450 parameterspace_ptrtype M_Dmu;
451 beta_vector_type M_betaAqm;
452 std::vector<beta_vector_type> M_betaFqm;
454 element_type M_unknown;
457 sparse_matrix_ptrtype M_A_tril;
Definition: solverlinear.hpp:33
simplex of dimension Dim
Definition: simplex.hpp:73
export Feel generated data to some file formatsUse the visitor and factory pattern.
Definition: exporter.hpp:82
double value_type
numerical type is double
Definition: convection_crb.hpp:194
base class for all linear algebra backends
Definition: backend.hpp:133
unifying mesh class
Definition: createsubmesh.hpp:41
Parameter space sampling class.
Definition: parameterspace.hpp:219
Definition: functionspace.hpp:1374
Parameter space class.
Definition: parameterspace.hpp:58
Definition: matrixsparse.hpp:76
parameterspace_ptrtype parameterSpace() const
return the parameter space
Definition: convection_crb.hpp:253
Definition: convection_crb.hpp:129
po::variables_map const & vm() const
Definition: convection_crb.hpp:391
element of a parameter space
Definition: parameterspace.hpp:84