p_laplacian_fixed_point.cc

Go to the documentation of this file.

#include "rheolef.h"
using namespace rheolef;
using namespace std;
#include "eta.icc"
#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_iter = (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");
  string grad_approx = "P" + itos(Xh.degree()-1) + "d";
  space Th (omega, grad_approx, "vector");
  form inv_mt (Th, Th, "inv_mass");
  form m (Xh, Xh, "mass");
  form b (Xh, Th, "grad");
  solver sm (m.uu());
  quadrature_option_type qopt;
  qopt.set_family (quadrature_option_type::gauss);
  qopt.set_order  (2*Xh.degree()-1);
  field uh (Xh);
  uh ["boundary"] = 0;
  field lh = riesz (Xh, 1);
  dirichlet (lh, uh);
  derr << "# n r v" << endl;
  Float r = 1, r0 = 1;
  size_t n = 0;
  do {
    field grad_uh = inv_mt*(b*uh);
    form a (Xh, Xh, "grad_grad", compose(eta(p), norm2(grad_uh)), qopt);
    field mrh = a*uh - lh;
    field rh (Xh, 0);
    rh.set_u() = sm.solve (mrh.u());
    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_iter) break;
    solver sa (a.uu());
    vec<Float> u_star = sa.solve (lh.u() - a.ub()*uh.b());
    uh.set_u() = w*u_star + (1-w)*uh.u();
  } while (true);
  dout << catchmark("p") << p << endl
       << catchmark("u") << uh;
  return (r <= tol) ? 0 : 1;
}