rheolef  6.6
stokes_cavity_check.cc
Go to the documentation of this file.
1 #include "rheolef.h"
2 using namespace rheolef;
3 using namespace std;
4 int main(int argc, char**argv) {
5  environment rheolef (argc, argv);
6  Float tol = (argc > 1) ? atof(argv[1]) : 1e-10;
7  field uh, ph;
8  din >> catchmark("u") >> uh
9  >> catchmark("p") >> ph;
10  const space& Xh = uh.get_space();
11  const space& Qh = ph.get_space();
12  trial u (Xh), p (Qh); test v (Xh), q (Qh);
13  form a = integrate (2*ddot(D(u),D(v)));
14  form b = integrate (-div(u)*q);
15  form mp = integrate (p*q);
16  solver fact_mp (mp.uu());
17  field m_ru = a*uh + b.trans_mult(ph);
18  field m_rp = b*uh;
19  m_ru["top"] = m_ru["bottom"] = 0;
20  if (Xh.get_geo().dimension() == 3) {
21  m_ru[1]["left"] = m_ru[1]["right"] = 0;
22  m_ru["front"] = m_ru["back"] = 0;
23  } else {
24  m_ru["left"] = m_ru["right"] = 0;
25  }
26  field rp (Qh);
27  rp.set_u() = fact_mp.solve(m_rp.u());
28  Float res_u = m_ru.u().max_abs();
29  Float res_p = sqrt(mp(rp,rp));
30  Float res = max(res_u, res_p);
31  Float p_constant = mp(ph, field(Qh,1.));
32  derr << "check: residue(uh) = " << res_u << endl
33  << "check: residue(ph) = " << res_p << endl
34  << "check: residue = " << res << endl
35  << "m(p,1) = " << p_constant << endl;
36  return (res <= tol) ? 0 : 1;
37 }
field - piecewise polynomial finite element field
Definition: field.h:225
const vec< T, M > & u() const
Definition: field.h:279
STL namespace.
field_vf_expr< field_vf_expr_grad< test_basic< T, M, VfTag > >> D(const test_basic< T, M, VfTag > &x)
Definition: operators2.h:97
field_basic< T, M > trans_mult(const field_basic< T, M > &yh) const
Definition: form.cc:101
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
field_basic< Float > field
Definition: field.h:367
space – piecewise polynomial finite element space
Definition: space.h:229
double Float
Definition: compiler.h:177
T ddot(const tensor_basic< T > &a, const tensor_basic< T > &b)
Definition: tensor.cc:200
field_nonlinear_expr< field_nonlinear_expr_uf< details::sqrt_,field_expr_terminal_field< T, M > > > sqrt(const field_basic< T, M > &x)
idiststream din(cin)
field_vf_expr< field_vf_expr_div< test_basic< T, M, VfTag > >> div(const test_basic< T, M, VfTag > &x)
Definition: operators2.h:124
Definition: rotating-hill.h:1
catchmark - iostream manipulator
Definition: catchmark.h:30
boost::proto::result_of::make_expr< boost::proto::tag::function, vec_domain, const vec_detail::max_, const int &, const vec< T, M > & >::type max(const int &x, const vec< T, M > &y)
const space_type & get_space() const
Definition: field.h:271
odiststream derr(cerr)
Definition: diststream.h:340
form - representation of a finite element bilinear form
Definition: form.h:98
solver - direct or interative solver interface
Definition: solver.h:8
int main(int argc, char **argv)
Float u(const point &x)
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