/usr/share/doc/rheolef-doc/examples/p_laplacian_fixed_point.cc is in rheolef-doc 6.5-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 | #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_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);
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 = 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)),
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_it) 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;
}
|