p_laplacian_fixed_point.cc

The p-Laplacian problem by the fixed-point method.

The p-Laplacian problem by the fixed-point method

#include "rheolef.h"
using namespace rheolef;
using namespace std;
#include "eta.h"
#include "dirichlet.icc"
int main(int argc, char**argv) {
  environment rheolef (argc,argv);
  geo omega (argv[1]);
  Float  eps    = std::numeric_limits<Float>::epsilon();
  string approx = (argc > 2) ?      argv[2]  : "P1";
  Float  p      = (argc > 3) ? atof(argv[3]) : 1.5;
  Float  w      = (argc > 4) ? (is_float(argv[4]) ? atof(argv[4]) :2/p) :1;
  Float  tol    = (argc > 5) ? atof(argv[5]) : 1e5*eps;
  size_t max_it = (argc > 6) ? atoi(argv[6]) : 500;
  derr << "# P-Laplacian problem by fixed-point:" << endl
       << "# geo    = " << omega.name() << endl
       << "# approx = " << approx << endl
       << "# p      = " << p << endl
       << "# w      = " << w << endl
       << "# tol    = " << tol << endl;
  space Xh (omega, approx);
  Xh.block ("boundary");
  trial u (Xh); test v (Xh);
  form m = integrate (u*v);
  problem pm (m);
  field uh (Xh), uh_star (Xh, 0.);
  uh["boundary"] = uh_star["boundary"] = 0;
  field lh = integrate (v);
  dirichlet (lh, uh);
  derr << "# n r v" << endl;
  Float r = 1, r0 = 1;
  size_t n = 0;
  do {
    form a = integrate(compose(eta(p),norm2(grad(uh)))*dot(grad(u),grad(v)));
    field mrh = a*uh - lh;
    field rh (Xh, 0);
    pm.solve (mrh, rh);
    r = rh.max_abs();
    if (n == 0) { r0 = r; }
    Float v = (n == 0) ? 0 : log10(r0/r)/n;
    derr << n << " " << r << " " << v << endl;
    if (r <= tol || n++ >= max_it) break;
    problem p (a);
    p.solve (lh, uh_star);
    uh = w*uh_star + (1-w)*uh;
  } while (true);
  dout << catchmark("p") << p << endl
       << catchmark("u") << uh;
  return (r <= tol) ? 0 : 1;
}