← Back to branch summary

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

« back to all changes in this revision

Viewing changes to nfem/tst/ic0_tst.cc

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2010-06-12 09:08:59 UTC
  • Revision ID: james.westby@ubuntu.com-20100612090859-8gpm2gc7j3ab43et
Tags: upstream-5.89
Import upstream version 5.89

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/tst/ic0_tst.cc
 
1
// usage: prog mesh Pk precond max_iter
 
2
#include "rheolef.h"
 
3
using namespace std;
 
4
 
 
5
const Float a = acos(Float(-1));
 
6
struct f : unary_function<point,Float> {
 
7
  Float operator() (const point& x) const {
 
8
   switch (d) {
 
9
    case 1 : return    sqr(a)*sin(a*x[0]);
 
10
    case 2 : return  2*sqr(a)*sin(a*x[0])*sin(a*x[1]); 
 
11
    default: return  3*sqr(a)*sin(a*x[0])*sin(a*x[1])*sin(a*x[2]); 
 
12
   }
 
13
  }
 
14
  f(size_t d1) : d(d1) {}
 
15
  size_t d;
 
16
};
 
17
struct u : unary_function<point,Float> {
 
18
  Float operator() (const point& x) const {
 
19
   switch (d) {
 
20
    case 1 : return  sin(a*x[0]); 
 
21
    case 2 : return  sin(a*x[0])*sin(a*x[1]); 
 
22
    default: return  sin(a*x[0])*sin(a*x[1])*sin(a*x[2]); 
 
23
   }
 
24
  }
 
25
  u(size_t d1) : d(d1) {}
 
26
  size_t d;
 
27
};
 
28
struct grad_u : unary_function<point,point> {
 
29
  point operator() (const point& x) const {
 
30
   switch (d) {
 
31
    case 1 : return  point(a*cos(a*x[0]));
 
32
    case 2 : return  point(a*cos(a*x[0])*sin(a*x[1]),
 
33
                           a*sin(a*x[0])*cos(a*x[1]));
 
34
    default: return  point(a*cos(a*x[0])*sin(a*x[1])*sin(a*x[2]),
 
35
                           a*sin(a*x[0])*cos(a*x[1])*sin(a*x[2]), 
 
36
                           a*sin(a*x[0])*sin(a*x[1])*cos(a*x[2])); 
 
37
   }
 
38
  }
 
39
  grad_u(size_t d1) : d(d1) {}
 
40
  size_t d;
 
41
};
 
42
int main(int argc, char**argv) {
 
43
    geo omega (argv[1]);
 
44
    string approx   = (argc > 2) ?      argv[2]  : "P1";
 
45
    string precond  = (argc > 3) ?      argv[3]  : "ic0";
 
46
    size_t max_iter = (argc > 4) ? atoi(argv[4]) : 10000;
 
47
    Float tol = 1e-15;
 
48
    size_t d = omega.dimension();
 
49
    space Vh (omega, argv[2]);
 
50
    Vh.block ("boundary");
 
51
    form a(Vh, Vh, "grad_grad");
 
52
    form m(Vh, Vh, "mass");
 
53
    field uh (Vh, 0);
 
54
    field fh = interpolate(Vh, f(d));
 
55
    if (precond != "ic0") {
 
56
      int status = pcg (a.uu, uh.u, (m*fh).u - a.ub*uh.b, EYE, max_iter, tol, &cout, "eye");
 
57
    } else {
 
58
      int status = pcg (a.uu, uh.u, (m*fh).u - a.ub*uh.b, ic0(a.uu), max_iter, tol, &cout, "ic0");
 
59
    }
 
60
    cerr << "n_el     = " << omega.size() << endl
 
61
         << "max_iter = " << max_iter << endl
 
62
         << "residue  = " << tol << endl
 
63
        ;
 
64
    return (tol < 1e-10) ? 0 : 1;
 
65
}