Logo
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ on Travis-CI Feel++ on Twitter Feel++ on YouTube Feel++ community
 All Classes Files Functions Variables Typedefs Pages
penalisation.hpp
1 #ifndef PENALISATION_HPP
2 #define PENALISATION_HPP
3 
4 #include <feel/feelcore/application.hpp>
5 #include <feel/options.hpp>
6 
7 #include <feel/feeldiscr/functionspace.hpp>
8 #include <feel/feelpoly/im.hpp>
9 
10 #include <feel/feelfilters/loadmesh.hpp>
11 #include <feel/feelfilters/exporterensight.hpp>
12 #include <feel/feelvf/vf.hpp>
13 
14 #include <feel/feelalg/backend.hpp>
15 
16 #include <iostream>
17 #include <fstream>
18 
19 #define VELOCITY_ORDER 3
20 #define PRESSURE_ORDER 2
21 #define CARAC_ORDER 2
22 
23 
24 namespace Feel
25 {
26 
27 inline po::options_description makeOptions()
28 {
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" );
46  return mesOptions
47  .add( Feel::feel_options() )
48  .add( backend_options( "stokes_backend" ) )
49  ;
50 }//makeoptions()
51 
52 inline AboutData makeAbout()
53 {
54  AboutData about( "Penalisation", "penalisation", "0.1", "desc" );
55  return about;
56 }
57 
58 template <int Dim>
59 class Penalisation : public Application
60 {
61  // mesh
62  typedef Simplex<Dim,1,Dim> entity_type;
63 
64  typedef Mesh<entity_type> mesh_type;
65  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
66 
67  // bases
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;
71  // typedef bases<Lagrange<CARAC_ORDER,Scalar,Discontinuous> > basis_carac_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;
75  //typedef bases<basis_veloc_type, basis_pressure_type> basis_type;
76 
77 
78  // space
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;
85 
86  // elements
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;
96 
97  // backend
98  typedef Backend<double> backend_type;
99  typedef boost::shared_ptr<backend_type> backend_ptrtype;
100 
101  //algebra
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;
106 
107  //export
108  typedef Exporter<mesh_type> export_type;
109  typedef boost::shared_ptr<export_type> export_ptrtype;
110 
111 public :
112  Penalisation();
113  void run();
114 
115 private :
116  void exportResults();
117  void updatePosition();
118  void updateChi();
119  void stokes();
120  void addCL();
121  void initStokes();
122 
123  const double nu;
124  const double radius;
125 
126  const double H1;
127  const double H2;
128  const double L1;
129  const double E;
130 
131  const double Q;
132  const double Qtop;
133  const double Qbot;
134 
135  space_ptrtype Xh;
136  space_carac_ptrtype Yh;
137  space_p0_ptrtype P0h;
138 
139  mesh_ptrtype mesh,mesh_visu;
140  export_ptrtype exporter;
141  backend_ptrtype backend;
142 
143  sparse_matrix_ptrtype C;
144  sparse_matrix_ptrtype D;
145  vector_ptrtype F;
146 
147  element_ptrtype U;
148  element_carac_type carac;
149  element_p0_type p0;
150 
151  double xp,yp,zp;
152  double Vx, Vy, Vz;
153  const double epsilonpress, epsilon;
154  double Tfinal,t;
155  const double dt;
156  const double Ylimit;
157  int iter;
158 
159  boost::timer chrono;
160 
161 }; //Penalisation
162 
163 }//namespace Feel
164 
165 #endif //PENALISATION_HPP
Definition: penalisation.hpp:59
Definition: qs_laplacian.doc:2