1 #ifndef PENALISATION_HPP
2 #define PENALISATION_HPP
7 #include <feel/feeldiscr/functionspace.hpp>
11 #include <feel/feelfilters/exporterensight.hpp>
19 #define VELOCITY_ORDER 3
20 #define PRESSURE_ORDER 2
29 po::options_description mesOptions(
"Mes options pour l'appli" );
30 mesOptions.add_options()
31 (
"OutFolder", po::value<std::string>()->default_value(
"" ),
"Dossier Sortie" )
32 (
"hsize", po::value<double>()->default_value( 0.1 ),
"h size" )
33 (
"x0", po::value<double>()->default_value( 0 ),
"x0" )
34 (
"y0", po::value<double>()->default_value( 0.333 ),
"y0" )
35 (
"z0", po::value<double>()->default_value( 0 ),
"z0" )
36 (
"nu", po::value<double>()->default_value( 0.001002 ),
"viscosity" )
37 (
"radius", po::value<double>()->default_value( 0.35 ),
"Radius" )
38 (
"Q", po::value<double>()->default_value( 2. ),
"In Flow rate" )
39 (
"RQ", po::value<double>()->default_value( 0.5 ),
"Flow rate ratio Q1/Q0" )
40 (
"epsilonpress", po::value<double>()->default_value( 0.00001 ),
"penalisation term : div u = epsilonpress p " )
41 (
"epsilon", po::value<double>()->default_value( 0.8 ),
"penalisation term for particle" )
42 (
"Tfinal", po::value<double>()->default_value( 3 ),
"Final time" )
43 (
"Ylimit", po::value<double>()->default_value( 0 ),
"Stop the simulation if y reaches Ylimit (never stop if 0)" )
44 (
"DT", po::value<double>()->default_value( 0.01 ),
"time step" )
45 (
"test", po::value<int>()->default_value( 1 ),
"test application" );
47 .add( Feel::feel_options() )
54 AboutData about(
"Penalisation",
"penalisation",
"0.1",
"desc" );
59 class Penalisation :
public Application
62 typedef Simplex<Dim,1,Dim> entity_type;
64 typedef Mesh<entity_type> mesh_type;
65 typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
68 typedef Lagrange<VELOCITY_ORDER, Vectorial> basis_veloc_type;
69 typedef Lagrange<PRESSURE_ORDER, Scalar> basis_pressure_type;
70 typedef Lagrange<0, Scalar> basis_lag_type;
72 typedef bases<Lagrange<CARAC_ORDER,Scalar> > basis_carac_type;
73 typedef bases<Lagrange<0,Scalar,Discontinuous> > basis_p0_type;
74 typedef bases<basis_veloc_type, basis_pressure_type, basis_lag_type> basis_type;
79 typedef FunctionSpace<mesh_type, basis_type> space_type;
80 typedef boost::shared_ptr<space_type> space_ptrtype;
81 typedef FunctionSpace<mesh_type, basis_carac_type,Discontinuous> space_carac_type;
82 typedef boost::shared_ptr<space_carac_type> space_carac_ptrtype;
83 typedef FunctionSpace<mesh_type, basis_p0_type,Discontinuous> space_p0_type;
84 typedef boost::shared_ptr<space_p0_type> space_p0_ptrtype;
87 typedef typename space_type::element_type element_type;
88 typedef boost::shared_ptr<element_type> element_ptrtype;
89 typedef typename element_type::template sub_element<0>::type element_veloc_type;
90 typedef typename element_type::template sub_element<1>::type element_pressure_type;
91 typedef typename element_type::template sub_element<2>::type element_lag_type;
92 typedef typename space_carac_type::element_type element_carac_type;
93 typedef boost::shared_ptr<element_carac_type> element_carac_ptrtype;
94 typedef typename space_p0_type::element_type element_p0_type;
95 typedef boost::shared_ptr<element_p0_type> element_p0_ptrtype;
98 typedef Backend<double> backend_type;
99 typedef boost::shared_ptr<backend_type> backend_ptrtype;
102 typedef typename backend_type::sparse_matrix_type sparse_matrix_type;
103 typedef typename backend_type::vector_type vector_type;
104 typedef boost::shared_ptr<sparse_matrix_type> sparse_matrix_ptrtype;
105 typedef boost::shared_ptr<vector_type> vector_ptrtype;
108 typedef Exporter<mesh_type> export_type;
109 typedef boost::shared_ptr<export_type> export_ptrtype;
116 void exportResults();
117 void updatePosition();
136 space_carac_ptrtype Yh;
137 space_p0_ptrtype P0h;
139 mesh_ptrtype mesh,mesh_visu;
140 export_ptrtype exporter;
141 backend_ptrtype backend;
143 sparse_matrix_ptrtype C;
144 sparse_matrix_ptrtype D;
148 element_carac_type carac;
153 const double epsilonpress, epsilon;
165 #endif //PENALISATION_HPP