30 #ifndef _ADVREACT_HPP_
31 #define _ADVREACT_HPP_
33 #include <feel/feeldiscr/functionspace.hpp>
43 enum StabilizationMethods {NO, CIP, SGS, SUPG, GALS};
62 typedef Space space_type;
63 typedef typename space_type::value_type value_type;
65 typedef boost::shared_ptr<backend_type> backend_ptrtype;
67 static const uint16_type Dim = space_type::nDim;
70 static const value_type polyOrder = space_type::basis_type::nOrder;
73 typedef typename space_type::mesh_type mesh_type;
74 typedef typename space_type::mesh_ptrtype mesh_ptrtype;
77 typedef boost::shared_ptr<space_type> space_ptrtype;
78 typedef typename space_type::element_type element_type;
80 AdvReact(
const space_ptrtype& space,
81 const backend_ptrtype& backend,
bool imposeBC =
true )
83 M_mesh( space->mesh() ),
85 M_phi( space,
"phi" ),
86 M_operator( space, space, backend ),
87 M_operatorStab( space, space, backend ),
91 M_imposeBC( imposeBC ),
92 M_stabcoeff( 0.1 * std::pow( polyOrder, -3.5 ) )
98 void set_stabcoeff(
double stabcoeff )
100 M_stabcoeff = stabcoeff * std::pow( polyOrder, -3.5 );
103 void setStabMethod( StabilizationMethods Method )
109 template<
typename Esigma,
typename Ebeta,
110 typename Ef,
typename Eg>
111 void update(
const Esigma& sigma,
115 bool updateStabilization =
true
122 const element_type phi()
132 space_ptrtype M_space;
138 FsFunctionalLinear<space_type> M_rhs;
139 FsFunctionalLinear<space_type> M_rhsStab;
143 const bool M_imposeBC;
147 StabilizationMethods M_StabMethod;
152 template<
class Space>
153 template<
typename Esigma,
typename Ebeta,
154 typename Ef,
typename Eg>
159 bool updateStabilization
164 using namespace Feel::vf;
168 (
val( sigma )*idt( M_phi ) + gradt( M_phi )*
val( beta ) ) *
id( M_phi )
171 if ( M_imposeBC ==
true )
175 - chi( trans( beta )*N() < 0 ) *
176 val( trans( beta )*N() ) * idt( M_phi ) *
id( M_phi )
180 if ( M_StabMethod != NO )
182 if ( updateStabilization )
185 if ( M_StabMethod== CIP && ( M_stabcoeff != 0.0 ) )
190 val( M_stabcoeff*vf::pow( hFace(),2.0 ) *
191 abs( trans( beta )*N() ) ) *
192 ( jumpt( gradt( M_phi ) ) * jump( grad( M_phi ) ) )
196 else if ( M_StabMethod== SUPG )
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 ) ) );
217 coeff*L_op*
val( f ) );
222 coeff*L_op*
val( f ) );
226 else if ( M_StabMethod== GALS )
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 ) );
246 coeff*L_op*
val( f ) );
250 coeff*L_op*
val( f ) );
254 else if ( M_StabMethod== SGS )
256 AUTO( coeff, 1.0 / ( 2*vf::sqrt(
val( trans( beta ) )*
val( beta ) )/vf::h()+vf::abs(
val( sigma ) ) ) );
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 ) ) );
277 coeff*L_op*
val( f ) );
281 coeff*L_op*
val( f ) );
288 M_operator.add( 1.0, M_operatorStab );
292 M_operator.mat().close();
297 val( f ) *
id( M_phi )
300 if ( M_imposeBC ==
true )
304 val( - chi( trans( beta )*N() < 0 ) *
305 ( trans( beta )*N() )*( g ) ) *
id( M_phi )
309 if ( M_StabMethod != NO )
310 M_rhs.add( M_rhsStab );
316 template<
class Space>
317 void AdvReact<Space>::solve()
328 M_operator.applyInverse( M_phi, M_rhs );
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