50 #include <boost/program_options/parsers.hpp>
51 #include <boost/program_options/variables_map.hpp>
53 #include <feel/feeldiscr/functionspace.hpp>
66 uint16_type imOrder = 3*boost::fusion::result_of::value_at_c<typename Space::basis_type, 0>::basis_type::nOrder-1,
67 template<u
int16_type,u
int16_type,u
int16_type>
class Entity = Simplex>
73 typedef Space space_type;
75 static const uint16_type Dim = space_type::nDim;
76 typedef typename space_type::value_type value_type;
79 typedef Backend<value_type> backend_type;
82 typedef boost::shared_ptr<backend_type> backend_ptrtype;
83 typedef typename backend_type::sparse_matrix_type sparse_matrix_type;
84 typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
85 typedef typename backend_type::vector_type vector_type;
86 typedef typename backend_type::vector_ptrtype vector_ptrtype;
89 typedef typename space_type::mesh_type mesh_type;
90 typedef typename space_type::mesh_ptrtype mesh_ptrtype;
93 typedef bases<Lagrange<0,Scalar,Discontinuous> > basis_i_type;
96 typedef FunctionSpace<mesh_type, basis_i_type, Discontinuous,value_type > space_i_type;
97 typedef boost::shared_ptr<space_type> space_ptrtype;
98 typedef boost::shared_ptr<space_i_type> space_i_ptrtype;
99 typedef typename space_type::element_type element_type;
100 typedef typename element_type::template sub_element<0>::type element_u_type;
101 typedef typename element_type::template sub_element<1>::type element_p_type;
102 typedef typename space_i_type::element_type element_i_type;
105 static const uint16_type uOrder = boost::fusion::result_of::value_at_c<typename Space::basis_type, 0>::type::nOrder;
106 static const uint16_type pOrder = boost::fusion::result_of::value_at_c<typename Space::basis_type, 1>::type::nOrder;
110 typedef IM<Dim, Order, value_type, Entity> type;
113 Oseen(
const space_ptrtype& space,
114 const backend_ptrtype& backend,
115 const std::set<flag_type>& dirichletFlags,
116 const std::set<flag_type>& neumannFlags,
117 bool weak_dirichlet =
true,
118 bool export_matlab =
false
121 Oseen(
const space_ptrtype& space,
122 const backend_ptrtype& backend,
123 const std::set<flag_type>& dirichletFlags,
124 const std::set<flag_type>& neumannFlags,
125 po::variables_map
const& vm,
126 std::string
const& prefix =
""
130 void setBCCoeffDiff( value_type bccoeffdiff )
132 M_bcCoeffDiff = bccoeffdiff;
135 void setBCCoeffConv( value_type bccoeffconv )
137 M_bcCoeffConv = bccoeffconv;
140 void setStabCoeffDiv( value_type stabcoeffdiv )
142 M_stabCoeffDiv = stabcoeffdiv;
145 void setStabCoeffP( value_type stabcoeffp )
147 M_stabCoeffP = stabcoeffp;
150 void decouplePstab( element_p_type& pStab, value_type theta )
153 M_thetaPstab = theta;
162 void setEpsCompress( value_type epscompress )
164 M_epsCompress = epscompress;
167 void setDivDivCoeff( value_type divdivcoeff );
170 template<
typename ItRange,
typename EsigmaInc,
171 typename EnuInc,
typename EnuAbs,
172 typename Ebeta,
typename Ef,
typename Ec,
typename Eg,
174 void update(
const ItRange& itRange,
175 const EsigmaInc& sigmaInc,
182 const EnoSlip& noSlip,
183 bool updateStabilization =
true );
189 const element_u_type& velocity()
194 const element_p_type& pressure()
200 element_type& solution()
204 const element_type& solution()
const
209 value_type stabilizationEnergy()
const;
216 void solveNonSym( sparse_matrix_ptrtype
const& D, element_type& u, vector_ptrtype
const& F );
222 space_ptrtype M_space;
228 sparse_matrix_ptrtype M_matrixFull;
229 sparse_matrix_ptrtype M_matrixAu;
230 sparse_matrix_ptrtype M_matrixStab;
231 vector_ptrtype M_vectorRhsFull;
232 vector_ptrtype M_vectorSolFull;
234 value_type M_bcCoeffDiff;
235 value_type M_bcCoeffConv;
237 value_type M_stabCoeffDiv;
238 value_type M_stabCoeffP;
239 element_p_type* M_pStab;
240 value_type M_thetaPstab;
242 value_type M_epsCompress;
243 value_type M_divDivCoeff;
247 backend_ptrtype M_backend;
249 std::set<flag_type> M_dirichletFlags;
250 std::set<flag_type> M_neumannFlags;
252 space_i_ptrtype M_space_i;
253 element_i_type M_stabWeightP;
255 bool M_weak_dirichlet;
256 bool M_export_matlab;
260 template<
class Space, uint16_type imOrder,
261 template<u
int16_type,u
int16_type,u
int16_type>
class Entity>
262 Oseen<Space, imOrder, Entity>::Oseen(
const space_ptrtype& space,
263 const backend_ptrtype& backend,
264 const std::set<flag_type>& dirichletFlags,
265 const std::set<flag_type>& neumannFlags,
269 M_mesh( space->mesh() ),
271 M_sol( space,
"sol" ),
272 u( M_sol.template element<0>() ),
273 p( M_sol.template element<1>() ),
274 M_matrixFull( backend->newMatrix( M_space, M_space ) ),
275 M_matrixAu( backend->newMatrix( M_space, M_space ) ),
276 M_matrixStab( backend->newMatrix( M_space, M_space ) ),
277 M_vectorRhsFull( backend->newVector( M_space ) ),
278 M_vectorSolFull( backend->newVector( M_space ) ),
279 M_bcCoeffDiff( 100 ),
280 M_bcCoeffConv( 100 ),
281 M_stabCoeffDiv( 0.0 ),
284 M_thetaPstab( -1.0 ),
285 M_epsCompress( 0.0 ),
286 M_divDivCoeff( 0.0 ),
288 M_backend( backend ),
289 M_dirichletFlags( dirichletFlags ),
290 M_neumannFlags( neumannFlags ),
291 M_space_i( new space_i_type( space->mesh() ) ),
292 M_stabWeightP( M_space_i,
"cs" ),
293 M_weak_dirichlet( weak_dirichlet ),
294 M_export_matlab( export_matlab )
296 using namespace Feel::vf;
298 element_u_type& v = u;
299 element_p_type& q = p;
301 DVLOG(2) <<
"[Oseen::Oseen] -(p,div v) + (q,div u)\n";
302 form2( _test=M_space, _trial=M_space, _matrix=M_matrixAu, _init=
true ) =
304 - div( v ) * idt( p )
307 +
id( q ) * divt( u )
314 M_vectorRhsFull->zero();
318 template<
class Space, uint16_type imOrder,
319 template<u
int16_type,u
int16_type,u
int16_type>
class Entity>
320 Oseen<Space, imOrder, Entity>::Oseen(
const space_ptrtype& space,
321 const backend_ptrtype& backend,
322 const std::set<flag_type>& dirichletFlags,
323 const std::set<flag_type>& neumannFlags,
324 po::variables_map
const& vm,
325 std::string
const& prefix )
327 M_mesh( space->mesh() ),
329 M_sol( space,
"sol" ),
330 u( M_sol.template element<0>() ),
331 p( M_sol.template element<1>() ),
332 M_matrixFull( backend->newMatrix( M_space, M_space ) ),
333 M_matrixAu( backend->newMatrix( M_space, M_space ) ),
334 M_matrixStab( backend->newMatrix( M_space, M_space ) ),
335 M_vectorRhsFull( backend->newVector( M_space ) ),
336 M_vectorSolFull( backend->newVector( M_space ) ),
337 M_bcCoeffDiff( 0.0 ),
338 M_bcCoeffConv( 0.0 ),
339 M_stabCoeffDiv( 0.0 ),
341 M_epsCompress( 0.0 ),
342 M_divDivCoeff( 0.0 ),
344 M_backend( backend ),
345 M_dirichletFlags( dirichletFlags ),
346 M_neumannFlags( neumannFlags ),
347 M_space_i( new space_i_type( space->mesh() ) ),
348 M_stabWeightP( M_space_i,
"cs" ),
349 M_weak_dirichlet( true ),
350 M_export_matlab( true )
352 std::string _prefix = prefix;
354 if ( !_prefix.empty() )
357 M_bcCoeffDiff = vm[_prefix+
"oseen-bc-coeff-diff"].template as<double>();
358 M_bcCoeffConv = vm[_prefix+
"oseen-bc-coeff-conv"].template as<double>();
359 M_stabCoeffDiv = vm[_prefix+
"oseen-stab-coeff-div"].template as<double>();
360 M_stabCoeffP = vm[_prefix+
"oseen-stab-coeff-p"].template as<double>();
361 M_epsCompress = vm[_prefix+
"oseen-eps-compress"].template as<double>();
362 M_divDivCoeff = vm[_prefix+
"oseen-divdiv-coeff"].template as<double>();
363 M_weak_dirichlet = vm[_prefix+
"oseen-weak-dirichlet"].template as<bool>();
364 M_export_matlab = vm[_prefix+
"oseen-export-matlab"].template as<bool>();
366 using namespace Feel::vf;
368 element_u_type& v = u;
369 element_p_type& q = p;
371 DVLOG(2) <<
"[Oseen::Oseen] -(p,div v) + (q,div u)\n";
372 form2( _test=M_space, _trial=M_space, _matrix=M_matrixAu, _init=
true ) =
374 - div( v ) * idt( p )
377 +
id( q ) * divt( u )
383 M_vectorRhsFull->zero();
387 template<
class Space, uint16_type imOrder,
388 template<u
int16_type,u
int16_type,u
int16_type>
class Entity>
390 Oseen<Space, imOrder, Entity>::setDivDivCoeff( value_type divdivcoeff )
392 using namespace Feel::vf;
394 element_u_type& v = u;
396 value_type delta = divdivcoeff - M_divDivCoeff;
400 DVLOG(2) <<
"[Oseen::set_divdivcoeff] (h gamma div u,div v)\n";
401 form2( _test=M_space, _trial=M_space, _matrix=M_matrixAu ) +=
403 h()*delta*div( v )*divt( u )
405 M_divDivCoeff = divdivcoeff;
408 template<
class Space, uint16_type imOrder,
409 template<u
int16_type,u
int16_type,u
int16_type>
class Entity>
410 typename Oseen<Space, imOrder, Entity>::value_type
411 Oseen<Space, imOrder, Entity>::stabilizationEnergy()
const
413 vector_ptrtype temp( M_backend->newVector( M_space ) );
415 M_backend->prod( *M_matrixStab, M_sol.container(), *temp );
416 return M_backend->dot( *temp, M_sol.container() );
420 template<
class Space, uint16_type imOrder,
421 template<u
int16_type,u
int16_type,u
int16_type>
class Entity>
422 template<
typename ItRange,
typename EsigmaInc,
423 typename EnuInc,
typename EnuAbs,
424 typename Ebeta,
typename Ef,
typename Ec,
typename Eg,
426 void Oseen<Space, imOrder, Entity>::update(
const ItRange& itRange,
427 const EsigmaInc& sigmaInc,
434 const EnoSlip& noSlip,
435 bool updateStabilization )
437 element_u_type& v = u;
438 element_p_type& q = p;
442 using namespace Feel::vf;
444 auto bcCoeffVeloFull = ( nuAbs )*( noSlip )*M_bcCoeffDiff/hFace()-chi( ( trans( beta )*N() ) < 0 ) * trans( beta )*N();
445 auto bcCoeffVeloNorm = M_bcCoeffConv * vf::max( vf::sqrt( trans( beta )*( beta ) ), ( nuAbs )/hFace() );
450 DVLOG(2) <<
"[Oseen::update] (f,v)+(c,q)\n";
451 form1( M_space, M_vectorRhsFull, _init=
true ) =
453 DVLOG(2) <<
"[Oseen::update] (f,v)+(c,q) done\n";
455 if ( M_weak_dirichlet )
458 for (
auto diriIter = M_dirichletFlags.begin(),
459 diriEnd = M_dirichletFlags.end();
460 diriIter != diriEnd; ++diriIter )
462 DVLOG(2) <<
"[Oseen::update] <bc(g),v>_Gamma_"
463 << *diriIter <<
"\n";
464 form1( M_space, M_vectorRhsFull ) +=
466 ( -
val( noSlip )*(
val( nuAbs )*trans( N() )
467 *( grad( v )+trans( grad( v ) ) )
468 +
id( q )*trans( N() ) )
469 +
val( bcCoeffVeloFull )*trans(
id( v ) )
470 +
val( bcCoeffVeloNorm )*( trans(
id( v ) )*N() )*trans( N() ) )
473 DVLOG(2) <<
"[Oseen::update] <bc(g),v>_Gamma_"
474 << *diriIter <<
" done\n";
481 DVLOG(2) <<
"[Oseen::update] adding (nu D(u),D(v)) + (sigma u,v)\n";
482 form2( M_space, M_space, _matrix=M_matrixAu ) +=
484 val( nuInc )*trace( grad( v )*( gradt( u )+trans( gradt( u ) ) ) )
485 +
val( sigmaInc )*trans(
id( v ) )*idt( u )
487 DVLOG(2) <<
"[Oseen::update] adding (nu D(u),D(v)) + (sigma u,v) done\n";
493 DVLOG(2) <<
"[Oseen::update] ((beta*grad)u,v)\n";
494 form2( M_space, M_space, M_matrixFull, _init=
true ) =
496 trans(
id( v ) ) * ( gradt( u )*
val( beta ) )
501 if ( ( M_stabCoeffDiv != 0.0 ) ||
502 ( M_stabCoeffP != 0.0 ) ||
503 ( M_epsCompress != 0.0 ) )
505 if ( updateStabilization )
507 form2( M_space, M_space, M_matrixStab, _init=
true );
510 if ( M_stabCoeffDiv != 0.0 )
512 DVLOG(2) <<
"[Oseen::update] <gamma_div [div u],[div v]>_Gamma_int\n";
513 form( M_space, M_space, M_matrixStab ) +=
515 val( M_stabCoeffDiv * vf::pow( hFace(),2.0 ) *
516 vf::sqrt( trans( beta )*( beta ) ) )
517 * trans( jump( div( v ) ) ) * jumpt( divt( u ) )
521 if ( M_stabCoeffP != 0.0 )
523 DVLOG(2) <<
"[Oseen::update] <gamma_p [grad p],[grad q]>_Gamma_int\n";
526 M_stabCoeffP * vf::pow( h(),3.0 ) /
527 max( h() * vf::sqrt( trans( beta ) * ( beta ) ),
533 form( M_space, M_space, M_matrixStab ) +=
535 typename im<2*pOrder-2>::type(),
537 * maxface( idv( M_stabWeightP ) ) )
538 * ( leftface ( grad ( q )*N() ) *
539 leftfacet ( gradt( p )*N() ) +
540 rightface ( grad ( q )*N() ) *
541 rightfacet( gradt( p )*N() ) )
547 form( M_space, M_space, M_matrixStab ) +=
549 typename im<2*pOrder-2>::type(),
550 maxface( idv( M_stabWeightP ) )
551 * jump( grad( q ) ) * jumpt( gradt( p ) )
558 if ( M_epsCompress > 0 )
560 VLOG(1) <<
"[Oseen] adding pseudo compressibility term with coeff= " << M_epsCompress <<
"\n";
561 form( M_space, M_space, M_matrixStab ) +=
563 M_epsCompress *
id( q ) * idt( p )
565 VLOG(1) <<
"[Oseen] adding pseudo compressibility term with coeff= " << M_epsCompress <<
" done\n";
568 M_matrixStab->close();
577 if ( M_pStab && ( M_stabCoeffP != 0.0 ) )
579 DVLOG(2) <<
"[Oseen::update] added rhs stab\n";
580 const int staborder = 2*pOrder-2;
581 form( M_space, *M_vectorRhsFull ) +=
583 val( maxface( idv( M_stabWeightP ) ) )
584 * ( - ( leftface ( grad ( q )*N() ) *
585 rightfacev( gradv( *M_pStab )*N() ) +
586 rightface ( grad ( q )*N() ) *
587 leftfacev ( gradv( *M_pStab )*N() ) )
588 + ( M_thetaPstab - 1.0 )
589 * ( ( leftface ( grad ( q )*N() ) *
590 leftfacev ( gradv( *M_pStab )*N() ) +
591 rightface ( grad ( q )*N() ) *
592 rightfacev ( gradv( *M_pStab )*N() ) )
596 DVLOG(2) <<
"[Oseen::update] added rhs stab done\n";
602 for (
auto diriIter = M_dirichletFlags.begin(),
603 diriEnd = M_dirichletFlags.end();
604 diriIter != diriEnd; ++diriIter )
606 if ( M_weak_dirichlet )
608 DVLOG(2) <<
"[Oseen::update] <bc(u),v>_Gamma_"
609 << *diriIter <<
"\n";
610 form2( M_space, M_space, M_matrixFull ) +=
612 -
val( noSlip )*(
val( ( nuAbs )*trans( N() ) )
613 *( ( grad( v )+trans( grad( v ) ) ) *idt( u ) +
614 ( gradt( u )+trans( gradt( u ) ) )*
id( v ) )
615 +
id( q )*trans( N() )*idt( u )
616 - idt( p )*trans( N() )*
id( v ) )
617 +
val( bcCoeffVeloFull )*( trans(
id( v ) )*idt( u ) )
618 +
val( bcCoeffVeloNorm )
619 * ( trans(
id( v ) )*N() ) * ( trans( N() )*idt( u ) )
621 DVLOG(2) <<
"[Oseen::update] <bc(u),v>_Gamma_"
622 << *diriIter <<
" done\n";
627 M_matrixFull->close();
628 M_vectorRhsFull->close();
629 DVLOG(2) <<
"[Oseen::update] on(u)_Gamma_"
630 << *diriIter <<
"\n";
631 form2( M_space, M_space, M_matrixFull ) +=
633 u, *M_vectorRhsFull, g );
634 DVLOG(2) <<
"[Oseen::update] on(u)_Gamma_"
635 << *diriIter <<
"done \n";
643 M_matrixFull->close();
644 M_vectorRhsFull->close();
649 VLOG(1) <<
"[Oseen] added stabilisation matrix\n";
650 M_matrixFull->addMatrix( 1.0, M_matrixAu );
651 VLOG(1) <<
"[Oseen] added standard oseen matrix\n";
653 if ( M_export_matlab )
655 M_matrixFull->printMatlab(
"M_oseen.m" );
656 M_vectorRhsFull->printMatlab(
"F_oseen.m" );
661 template<
class Space, uint16_type imOrder,
662 template<u
int16_type,u
int16_type,u
int16_type>
class Entity>
663 void Oseen<Space, imOrder, Entity>::solve()
675 solveNonSym( M_matrixFull, M_vectorSolFull, M_vectorRhsFull );
678 std::copy( M_vectorSolFull.begin(),
679 M_vectorSolFull.end(),
684 solveNonSym( M_matrixFull, M_sol, M_vectorRhsFull );
692 template<
class Space, uint16_type imOrder,
693 template<u
int16_type,u
int16_type,u
int16_type>
class Entity>
695 Oseen<Space, imOrder, Entity>::solveNonSym( sparse_matrix_ptrtype
const& D,
697 vector_ptrtype
const& F )
699 M_backend->solve( _matrix=D, _solution=u, _rhs=F );
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::marker_face_const_iterator, typename MeshTraits< MeshType >::marker_face_const_iterator > markedfaces(MeshType const &mesh)
Definition: filters.hpp:969
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
Polynomial< Pset, Scalar > project(Pset const &pset, Func const &f, IM const &im)
Definition: operations.hpp:436
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
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