Rheolef  7.1
an efficient C++ finite element environment
p_laplacian_fixed_point.cc
Go to the documentation of this file.
1 #include "rheolef.h"
26 using namespace rheolef;
27 using namespace std;
28 #include "eta.h"
29 #include "dirichlet.icc"
30 int main(int argc, char**argv) {
31  environment rheolef (argc,argv);
32  geo omega (argv[1]);
34  string approx = (argc > 2) ? argv[2] : "P1";
35  Float p = (argc > 3) ? atof(argv[3]) : 1.5;
36  Float w = (argc > 4) ? (is_float(argv[4]) ? atof(argv[4]) :2/p) :1;
37  Float tol = (argc > 5) ? atof(argv[5]) : 1e5*eps;
38  size_t max_it = (argc > 6) ? atoi(argv[6]) : 500;
39  derr << "# P-Laplacian problem by fixed-point:" << endl
40  << "# geo = " << omega.name() << endl
41  << "# approx = " << approx << endl
42  << "# p = " << p << endl
43  << "# w = " << w << endl
44  << "# tol = " << tol << endl;
45  space Xh (omega, approx);
46  Xh.block ("boundary");
47  trial u (Xh); test v (Xh);
48  form m = integrate (u*v);
49  problem pm (m);
50  field uh (Xh), uh_star (Xh, 0.);
51  uh["boundary"] = uh_star["boundary"] = 0;
52  field lh = integrate (v);
53  dirichlet (lh, uh);
54  derr << "# n r v" << endl;
55  Float r = 1, r0 = 1;
56  size_t n = 0;
57  do {
58  form a = integrate(compose(eta(p),norm2(grad(uh)))*dot(grad(u),grad(v)));
59  field mrh = a*uh - lh;
60  field rh (Xh, 0);
61  pm.solve (mrh, rh);
62  r = rh.max_abs();
63  if (n == 0) { r0 = r; }
64  Float v = (n == 0) ? 0 : log10(r0/r)/n;
65  derr << n << " " << r << " " << v << endl;
66  if (r <= tol || n++ >= max_it) break;
67  problem p (a);
68  p.solve (lh, uh_star);
69  uh = w*uh_star + (1-w)*uh;
70  } while (true);
71  dout << catchmark("p") << p << endl
72  << catchmark("u") << uh;
73  return (r <= tol) ? 0 : 1;
74 }
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
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
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:402
geo
see the geo page for the full documentation