rheolef  6.6
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);
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 }
quadrature_option_type - send options to the integrate function
boost::proto::result_of::make_expr< boost::proto::tag::function, vec_domain, const vec_detail::compose_< Function >, const vec< T, M > & >::type compose(const Function &f, const vec< T, M > &x)
field - piecewise polynomial finite element field
Definition: field.h:225
const vec< T, M > & u() const
Definition: field.h:279
field_vf_expr< field_vf_expr_grad< test_basic< T, M, VfTag > >> grad(const test_basic< T, M, VfTag > &x)
Definition: operators2.h:74
int main(int argc, char **argv)
Definition: convect.cc:5
STL namespace.
vec< T, M > & set_u()
Definition: field.h:281
irheostream, orheostream - large data streams
Definition: compiler.h:7
const csr< T, M > & uu() const
Definition: form.h:140
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
const csr< T, M > & ub() const
Definition: form.h:141
space – piecewise polynomial finite element space
Definition: space.h:229
double Float
Definition: compiler.h:177
size_t d
field_nonlinear_expr< field_nonlinear_expr_uf< details::acos_,field_expr_terminal_field< T, M > > > acos(const field_basic< T, M > &x)
const vec< T, M > & b() const
Definition: field.h:280
catchmark - iostream manipulator
Definition: catchmark.h:30
Float sigma
T dot(const vec< T, M > &x, const int &y)
form - representation of a finite element bilinear form
Definition: form.h:98
solver - direct or interative solver interface
Definition: solver.h:8
characteristic - the Lagrange-Galerkin method implemented
Float u(const point &x)
odiststream dout(cout)
Definition: diststream.h:317
field_basic< T, M > integrate(const geo_basic< T, M > &domain, const test_basic< T, M, details::vf_tag_01 > &expr, const quadrature_option_type &qopt=quadrature_option_type())
integrate - integrate a function or an expression
Definition: integrate.h:90