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
advreact.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christoph Winkelmann <christoph.winkelmann@epfl.ch>
6  Date: 2006-10-04
7 
8  Copyright (C) 2006 EPFL
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 3.0 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 */
30 #ifndef _ADVREACT_HPP_
31 #define _ADVREACT_HPP_
32 
33 #include <feel/feeldiscr/functionspace.hpp>
34 #include <feel/feelpoly/im.hpp>
35 
36 #include <feel/feelvf/vf.hpp>
37 
39 
40 namespace Feel
41 {
42 
43 enum StabilizationMethods {NO, CIP, SGS, SUPG, GALS};
44 
56 template<class Space>
57 class AdvReact
58 {
59 public:
60 
61  // -- TYPEDEFS --
62  typedef Space space_type;
63  typedef typename space_type::value_type value_type;
65  typedef boost::shared_ptr<backend_type> backend_ptrtype;
66 
67  static const uint16_type Dim = space_type::nDim;
68 
69 
70  static const value_type polyOrder = space_type::basis_type::nOrder;
71 
72  /* mesh */
73  typedef typename space_type::mesh_type mesh_type;
74  typedef typename space_type::mesh_ptrtype mesh_ptrtype;
75 
76  /* space */
77  typedef boost::shared_ptr<space_type> space_ptrtype;
78  typedef typename space_type::element_type element_type;
79 
80  AdvReact( const space_ptrtype& space,
81  const backend_ptrtype& backend, bool imposeBC = true )
82  :
83  M_mesh( space->mesh() ),
84  M_space( space ),
85  M_phi( space, "phi" ),
86  M_operator( space, space, backend ),
87  M_operatorStab( space, space, backend ),
88  M_rhs( space ),
89  M_rhsStab( space ),
90  M_updated( false ),
91  M_imposeBC( imposeBC ),
92  M_stabcoeff( 0.1 * std::pow( polyOrder, -3.5 ) )
93  {
94  M_StabMethod=GALS;
95  }
96 
97  // setting of options
98  void set_stabcoeff( double stabcoeff )
99  {
100  M_stabcoeff = stabcoeff * std::pow( polyOrder, -3.5 );
101  }
102 
103  void setStabMethod( StabilizationMethods Method )
104  {
105  M_StabMethod=Method;
106  }
107 
108  // update operator and rhs with given expressions
109  template<typename Esigma, typename Ebeta,
110  typename Ef, typename Eg>
111  void update( const Esigma& sigma,
112  const Ebeta& beta,
113  const Ef& f,
114  const Eg& g,
115  bool updateStabilization = true
116  );
117 
118  // solve system, call not needed, but possible
119  void solve();
120 
121  // result
122  const element_type phi()
123  {
124  solve();
125  return M_phi;
126  }
127 
128 private:
129 
130  mesh_ptrtype M_mesh;
131 
132  space_ptrtype M_space;
133 
134  element_type M_phi;
135 
138  FsFunctionalLinear<space_type> M_rhs;
139  FsFunctionalLinear<space_type> M_rhsStab;
140 
141  bool M_updated;
142 
143  const bool M_imposeBC;
144 
145  double M_stabcoeff;
146 
147  StabilizationMethods M_StabMethod;
148 
149 }; // class AdvReact
150 
151 
152 template<class Space>
153 template<typename Esigma, typename Ebeta,
154  typename Ef, typename Eg>
155 void AdvReact<Space>::update( const Esigma& sigma,
156  const Ebeta& beta,
157  const Ef& f,
158  const Eg& g,
159  bool updateStabilization
160  )
161 {
162  M_updated = true;
163 
164  using namespace Feel::vf;
165 
166  M_operator =
167  integrate( elements( M_mesh ),
168  ( val( sigma )*idt( M_phi ) + gradt( M_phi )*val( beta ) ) * id( M_phi )
169  ) ;
170 
171  if ( M_imposeBC == true )
172  {
173  M_operator +=
174  integrate( boundaryfaces( M_mesh ),
175  - chi( trans( beta )*N() < 0 ) * /* inflow */
176  val( trans( beta )*N() ) * idt( M_phi ) * id( M_phi )
177  );
178  }
179 
180  if ( M_StabMethod != NO )
181  {
182  if ( updateStabilization )
183  {
184  //good review of stabilization methods in [Chaple 2006]
185  if ( M_StabMethod== CIP && ( M_stabcoeff != 0.0 ) )
186  {
187  /* don't work properly in 3D (because of internalfaces) */
188  M_operatorStab=
189  integrate( internalfaces( M_mesh ),
190  val( M_stabcoeff*vf::pow( hFace(),2.0 ) *
191  abs( trans( beta )*N() ) ) *
192  ( jumpt( gradt( M_phi ) ) * jump( grad( M_phi ) ) )
193  );
194  }//Continuous Interior Penalty
195 
196  else if ( M_StabMethod== SUPG )
197  {
198  AUTO( coeff, vf::h()/( 2*vf::sqrt( val( trans( beta ) )*val( beta ) ) ) );
199  AUTO( L_op, ( grad( M_phi )*val( beta ) ) );
200  AUTO( L_opt, ( gradt( M_phi )*val( beta ) + val( sigma )*idt( M_phi ) ) );
201 
202  M_operatorStab=
203  integrate( elements( M_mesh ),
204  coeff
205  * L_op
206  * L_opt );
207 
208 
209  M_operatorStab+=
210  integrate( boundaryfaces( M_mesh ),
211  coeff
212  * L_op
213  * L_opt );
214 
215  M_rhsStab=
216  integrate( elements( M_mesh ),
217  coeff*L_op*val( f ) );
218 
219 
220  M_rhsStab+=
221  integrate( boundaryfaces( M_mesh ),
222  coeff*L_op*val( f ) );
223 
224  }//Streamline Upwind Petrov Galerkin
225 
226  else if ( M_StabMethod== GALS )
227  {
228  auto coeff = val( 1.0 / ( 2*vf::sqrt( trans( beta )*beta )/vf::h()+vf::abs( sigma ) ) );
229  auto L_op = ( grad( M_phi )*val( beta ) + val( sigma )*id( M_phi ) );
230  auto L_opt = ( gradt( M_phi )*val( beta ) + val( sigma )*idt( M_phi ) );
231 
232  M_operatorStab=
233  integrate( elements( M_mesh ),
234  coeff
235  * L_op
236  * L_opt );
237 
238  M_operatorStab+=
239  integrate( boundaryfaces( M_mesh ),
240  coeff
241  * L_op
242  * L_opt );
243 
244  M_rhsStab=
245  integrate( elements( M_mesh ),
246  coeff*L_op*val( f ) );
247 
248  M_rhsStab+=
249  integrate( boundaryfaces( M_mesh ),
250  coeff*L_op*val( f ) );
251 
252  }//Galerkin Least Square
253 
254  else if ( M_StabMethod== SGS )
255  {
256  AUTO( coeff, 1.0 / ( 2*vf::sqrt( val( trans( beta ) )*val( beta ) )/vf::h()+vf::abs( val( sigma ) ) ) );
257 
258  // AUTO(coeff_bound, 1.0 / (2*vf::sqrt(val(trans(beta))*val(beta))/vf::hFace()+vf::abs(val(sigma))));
259 
260  AUTO( L_op, ( grad( M_phi )* val( beta ) - val( sigma ) * id( M_phi ) ) );
261  AUTO( L_opt, ( gradt( M_phi )*val( beta ) + val( sigma )*idt( M_phi ) ) );
262 
263  M_operatorStab=
264  integrate( elements( M_mesh ),
265  coeff
266  * L_op
267  * L_opt );
268 
269  M_operatorStab+=
270  integrate( boundaryfaces( M_mesh ),
271  coeff
272  * L_op
273  * L_opt );
274 
275  M_rhsStab=
276  integrate( elements( M_mesh ),
277  coeff*L_op*val( f ) );
278 
279  M_rhsStab+=
280  integrate( boundaryfaces( M_mesh ),
281  coeff*L_op*val( f ) );
282  }//Subgrid Scale method
283 
284  //M_operatorStab.mat().printMatlab("M_advReact_stab.m");
285  //M_rhsStab.container().printMatlab("M_rhs_stab.m");
286  }//update stabilization
287 
288  M_operator.add( 1.0, M_operatorStab );
289 
290  }//DoStabilize
291 
292  M_operator.mat().close();
293  //M_operator.mat().printMatlab("M_advReact.m");
294 
295  M_rhs =
296  integrate( elements( M_mesh ),
297  val( f ) * id( M_phi )
298  );
299 
300  if ( M_imposeBC == true )
301  {
302  M_rhs +=
303  integrate( boundaryfaces( M_mesh ),
304  val( - chi( trans( beta )*N() < 0 ) /* inflow */ *
305  ( trans( beta )*N() )*( g ) ) * id( M_phi )
306  );
307  }
308 
309  if ( M_StabMethod != NO )
310  M_rhs.add( M_rhsStab );
311 
312  //M_rhs.container().printMatlab("F_advReact.m");
313 
314 } // update
315 
316 template<class Space>
317 void AdvReact<Space>::solve()
318 {
319  // -- make sure solve is needed
320  if ( !M_updated )
321  {
322  return;
323  }
324 
325  M_updated = false;
326 
327  // -- solve
328  M_operator.applyInverse( M_phi, M_rhs );
329 
330 } // AdvReact::solve
331 
332 } // Feel
333 
334 #endif /* _ADVREACT_HPP_ */
Poly::value_type integrate(Polynomial< Poly, PsetType > const &p)
integrate a polynomial over a convex
Definition: operations.hpp:51
elements_type const & elements() const
Definition: elements.hpp:355
Expr< Val< typename mpl::if_< boost::is_arithmetic< ExprT1 >, mpl::identity< Cst< ExprT1 > >, mpl::identity< ExprT1 > >::type::type > > val(ExprT1 const &__e1)
precompute expression tensor
Definition: val.hpp:304
base class for all linear algebra backends
Definition: backend.hpp:133
Advection-Reaction solver.
Definition: advreact.hpp:57
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::location_face_const_iterator, typename MeshTraits< MeshType >::location_face_const_iterator > boundaryfaces(MeshType const &mesh)
Definition: filters.hpp:1158
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::location_face_const_iterator, typename MeshTraits< MeshType >::location_face_const_iterator > internalfaces(MeshType const &mesh)
Definition: filters.hpp:1175
#define AUTO(a, b)
allow automatic type naming of complex expression
Definition: vf.hpp:45

Generated on Sun Dec 22 2013 13:10:54 for Feel++ by doxygen 1.8.5