← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty

« back to all changes in this revision

Viewing changes to doc/pexamples/stokes_cavity_check.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

Show diffs side-by-side

added added

removed removed

Lines of Context:
doc/pexamples/stokes_cavity_check.cc
 
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
  form a  (Xh, Xh, "2D_D");
 
13
  form b  (Xh, Qh, "div"); b = -b;
 
14
  form mp (Qh, Qh, "mass");
 
15
  solver fact_mp (mp.uu());
 
16
  field m_ru = a*uh + b.trans_mult(ph);
 
17
  field m_rp = b*uh;
 
18
  for (size_t i_comp = 0; i_comp < m_ru.size(); i_comp++) {
 
19
    m_ru[i_comp]["top"] = 0;
 
20
    m_ru[i_comp]["bottom"] = 0;
 
21
  }
 
22
  if (Xh.get_geo().dimension() == 3) {
 
23
    m_ru[1]["left"] = m_ru[1]["right"] = 0;
 
24
    for (size_t i_comp = 0; i_comp < m_ru.size(); i_comp++) {
 
25
      m_ru[i_comp]["front"] = 0;
 
26
      m_ru[i_comp]["back"] = 0;
 
27
    }
 
28
  } else {
 
29
    for (size_t i_comp = 0; i_comp < m_ru.size(); i_comp++) {
 
30
      m_ru[i_comp]["left"] = 0;
 
31
      m_ru[i_comp]["right"] = 0;
 
32
    }
 
33
  }
 
34
  field rp (Qh);
 
35
  rp.set_u() = fact_mp.solve(m_rp.u());
 
36
  Float res_u = m_ru.u().max_abs();
 
37
  Float res_p = sqrt(mp(rp,rp));
 
38
  Float res = max(res_u, res_p);
 
39
  Float p_constant = mp(ph, field(Qh,1.));
 
40
  derr << "check: residue(uh) = " << res_u << endl
 
41
       << "check: residue(ph) = " << res_p << endl
 
42
       << "check: residue     = " << res << endl
 
43
       << "m(p,1)             = " << p_constant << endl;
 
44
  return (res <= tol) ? 0 : 1;
 
45
}