rheolef  7.0
convect.cc
Go to the documentation of this file.
1 #include "rheolef.h"
2 using namespace rheolef;
3 using namespace std;
4 #include "rotating-hill.h"
5 int main (int argc, char **argv) {
6  environment rheolef (argc,argv);
7  geo omega (argv[1]);
8  string approx = (argc > 2) ? argv[2] : "P1";
9  Float nu = (argc > 3) ? atof(argv[3]) : 1e-2;
10  size_t n_max = (argc > 4) ? atoi(argv[4]) : 50;
11  size_t d = omega.dimension();
12  Float delta_t = 2*acos(-1.)/n_max;
13  space Vh (omega, approx, "vector");
14  field uh = interpolate (Vh, u(d));
15  space Xh (omega, approx);
16  Xh.block ("boundary");
17  field phi_h = interpolate (Xh, phi(d,nu,0));
18  characteristic X (-delta_t*uh);
19  quadrature_option qopt;
21  qopt.set_order (Xh.degree());
22  trial phi (Xh); test psi (Xh);
23  branch event ("t","phi");
24  dout << catchmark("nu") << nu << endl
25  << event (0, phi_h);
26  for (size_t n = 1; n <= n_max; n++) {
27  Float t = n*delta_t;
28  Float c1 = 1 + delta_t*phi::sigma(d,nu,t);
29  Float c2 = delta_t*nu;
30  form a = integrate (c1*phi*psi + c2*dot(grad(phi),grad(psi)), qopt);
31  field lh = integrate (compose(phi_h, X)*psi, qopt);
32  solver sa (a.uu());
33  phi_h.set_u() = sa.solve (lh.u() - a.ub()*phi_h.b());
34  dout << event (t, phi_h);
35  }
36 }
field - piecewise polynomial finite element field
details::field_expr_v2_nonlinear_node_nary< typename details::function_traits< Function >::functor_type,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits< Exprs >::type... > ::type compose(const Function &f, const Exprs &... exprs)
Definition: compose.h:157
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
Definition: vec_expr_v2.h:335
quadrature_option - send options to the integrate function
int main(int argc, char **argv)
Definition: convect.cc:5
STL namespace.
vec< T, M > & set_u()
Definition: field.h:293
irheostream, orheostream - large data streams
Definition: compiler.h:7
static field_basic< T, sequential > interpolate(const space_basic< T, sequential > &Vh, const field_basic< T, sequential > &uh)
branch - a parameter-dependent sequence of field
Definition: branch.h:46
std::enable_if< details::is_field_expr_v2_nonlinear_arg< Expr >::value &&! is_undeterminated< Result >::value, Result >::type integrate(const geo_basic< T, M > &omega, const Expr &expr, const quadrature_option &qopt, Result dummy=Result())
integrate - integrate a function or an expression
Definition: integrate.h:107
space – piecewise polynomial finite element space
Definition: space.h:229
double Float
Definition: compiler.h:160
size_t d
std::enable_if< details::is_field_expr_v2_linear_arg< Expr >::value,details::field_expr_v2_nonlinear_terminal_field_grad< typename Expr::scalar_type,typename Expr::memory_type >>::type grad(const Expr &expr)
const vec< T, M > & b() const
Definition: field.h:292
const csr< T, M > & ub() const
Definition: form.h:141
catchmark - iostream manipulator
Definition: catchmark.h:30
const csr< T, M > & uu() const
Definition: form.h:140
Float sigma
form - representation of a finite element bilinear form
Definition: form.h:98
characteristic - the Lagrange-Galerkin method implemented
Float u(const point &x)
const vec< T, M > & u() const
Definition: field.h:291
void set_family(family_type type)
odiststream dout(cout)
Definition: diststream.h:313