rheolef  6.5
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 }
38