p_laplacian_fixed_point.cc
Go to the documentation of this file.
#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;
}
form
see the form page for the full documentation
rheolef::catchmark
see the catchmark page for the full documentation
Definition: catchmark.h:67
rheolef::dot
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
field
see the field page for the full documentation
rheolef::integrate
std::enable_if< details::is_field_expr_v2_nonlinear_arg< Expr >::value &&! is_undeterminated< Result >::value, Result >::type integrate(const geo_basic< T, M > &omega, const Expr &expr, const integrate_option &iopt, Result dummy=Result())
see the integrate page for the full documentation
Definition: integrate.h:202
rheolef::is_float
bool is_float(const string &s)
is_float: see the rheostream page for the full documentation
Definition: rheostream.cc:476
eta
Definition: eta.h:25
space
see the space page for the full documentation
rheolef::grad
std::enable_if< details::is_field_convertible< Expr >::value,details::field_expr_v2_nonlinear_terminal_field< typename Expr::scalar_type,typename Expr::memory_type,details::differentiate_option::gradient >>::type grad(const Expr &expr)
grad(uh): see the expression page for the full documentation
Definition: field_expr_terminal.h:911
eta.h
The p-Laplacian problem – the eta function.
rheolef.h
rheolef - reference manual
rheolef::norm2
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition: vec.h:379
p
Definition: sphere.icc:25
main
int main(int argc, char **argv)
Definition: p_laplacian_fixed_point.cc:30
a
Definition: diffusion_isotropic.h:25
rheolef::environment
see the environment page for the full documentation
Definition: environment.h:115
lh
field lh(Float epsilon, Float t, const test &v)
Definition: burgers_diffusion_operators.icc:25
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
test
see the test page for the full documentation
problem
see the problem page for the full documentation
u
Definition: leveque.h:25
rheolef::derr
odiststream derr(cerr)
see the diststream page for the full documentation
Definition: diststream.h:436
Float
see the Float page for the full documentation
u
Float u(const point &x)
Definition: transmission_error.cc:26
rheolef::compose
details::field_expr_v2_nonlinear_node_nary< typename details::function_traits< Function >::functor_type,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits< Exprs >::type... > ::type compose(const Function &f, const Exprs &... exprs)
see the compose page for the full documentation
Definition: compose.h:246
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
trial
see the test page for the full documentation
rheolef::dout
odiststream dout(cout)
see the diststream page for the full documentation
Definition: diststream.h:430
epsilon
Float epsilon
Definition: transmission_error.cc:25
dirichlet
void dirichlet(const field &lh, field &uh)
Definition: dirichlet.icc:25
dirichlet.icc
The Poisson problem with homogeneous Dirichlet boundary condition – solver function.
rheolef::std
Definition: vec_expr_v2.h:391
geo
see the geo page for the full documentation