← Back to branch summary

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

« back to all changes in this revision

Viewing changes to doc/usrman/p-laplacian-fixed-point.h

  • 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:
doc/usrman/p-laplacian-fixed-point.h
 
1
int p_laplacian_fixed_point (Float p, field fh, field& uh, Float& r, size_t& n) {
 
2
    Float tol = r;
 
3
    Float r0 = 0;
 
4
    size_t max_iter =  n;
 
5
    const geo& omega_h = uh.get_geo();
 
6
    const space& Vh = uh.get_space();
 
7
    string grad_approx = (Vh.get_approx() == "P2") ? "P1d" : "P0";
 
8
    space Th (omega_h, grad_approx, "vector");
 
9
    form m (Vh, Vh, "mass");
 
10
    form inv_mt (Th, Th, "inv_mass");
 
11
    form grad (Vh, Th, "grad");
 
12
    cerr << "# Fixed-point algorithm on p-laplancian: p = " << p << endl
 
13
         << "# n r v" << endl;
 
14
    n = 0;
 
15
    do {
 
16
      field grad_uh = inv_mt*(grad*uh); 
 
17
      field wh = pow(sqr(grad_uh[0]) + sqr(grad_uh[1]), p/2.-1);
 
18
      form a (Vh, Vh, "grad_grad", wh);
 
19
      field rh = a*uh - m*fh;
 
20
      r = rh.u.max_abs();
 
21
      if (n == 0) r0 = r;
 
22
      Float v = (n == 0) ? 0 : log10(r0/r)/n;
 
23
      cerr << n << " " << r << " " << v << endl;
 
24
      if (r <= tol || n++ >= max_iter) break;
 
25
      ssk<Float> fact_a = ldlt(a.uu);
 
26
      uh.u = fact_a.solve (m.uu*fh.u + m.ub*fh.b - a.ub*uh.b);
 
27
    } while (true);
 
28
    return r <= tol;
 
29
}