rheolef  6.3
p_laplacian_fixed_point.cc
Go to the documentation of this file.
1 #include "rheolef.h"
2 using namespace rheolef;
3 using namespace std;
4 #include "eta.icc"
5 #include "dirichlet.icc"
6 int main(int argc, char**argv) {
7  environment rheolef (argc,argv);
8  geo omega (argv[1]);
10  string approx = (argc > 2) ? argv[2] : "P1";
11  Float p = (argc > 3) ? atof(argv[3]) : 1.5;
12  Float w = (argc > 4) ? (is_float(argv[4]) ? atof(argv[4]) : 2/p) : 1;
13  Float tol = (argc > 5) ? atof(argv[5]) : 1e5*eps;
14  size_t max_iter = (argc > 6) ? atoi(argv[6]) : 500;
15  derr << "# P-Laplacian problem by fixed-point:" << endl
16  << "# geo = " << omega.name() << endl
17  << "# approx = " << approx << endl
18  << "# p = " << p << endl
19  << "# w = " << w << endl
20  << "# tol = " << tol << endl;
21  space Xh (omega, approx);
22  Xh.block ("boundary");
23  string grad_approx = "P" + itos(Xh.degree()-1) + "d";
24  space Th (omega, grad_approx, "vector");
25  form inv_mt (Th, Th, "inv_mass");
26  form m (Xh, Xh, "mass");
27  form b (Xh, Th, "grad");
28  solver sm (m.uu());
31  qopt.set_order (2*Xh.degree()-1);
32  field uh (Xh);
33  uh ["boundary"] = 0;
34  field lh = riesz (Xh, 1);
35  dirichlet (lh, uh);
36  derr << "# n r v" << endl;
37  Float r = 1, r0 = 1;
38  size_t n = 0;
39  do {
42  field mrh = a*uh - lh;
43  field rh (Xh, 0);
44  rh.set_u() = sm.solve (mrh.u());
45  r = rh.max_abs();
46  if (n == 0) { r0 = r; }
47  Float v = (n == 0) ? 0 : log10(r0/r)/n;
48  derr << n << " " << r << " " << v << endl;
49  if (r <= tol || n++ >= max_iter) break;
50  solver sa (a.uu());
51  vec<Float> u_star = sa.solve (lh.u() - a.ub()*uh.b());
52  uh.set_u() = w*u_star + (1-w)*uh.u();
53  } while (true);
54  dout << catchmark("p") << p << endl
55  << catchmark("u") << uh;
56  return (r <= tol) ? 0 : 1;
57 }
58