Rheolef  7.1
an efficient C++ finite element environment
dirichlet_hho.cc
Go to the documentation of this file.
1 #include "rheolef.h"
26 using namespace rheolef;
27 using namespace std;
28 #include "dirichlet_homogeneous.h"
29 #include "diffusion_isotropic.h"
30 int main(int argc, char**argv) {
31  environment rheolef (argc, argv);
32  geo omega (argv[1]);
33  string Pkd = (argc > 2) ? argv[2] : "P1d",
34  Pld = (argc > 3) ? argv[3] : Pkd;
35  Float beta = (argc > 4) ? atof(argv[4]) : 1;
36  space Xh (omega, Pld),
37  Mh (omega["sides"], Pkd);
38  Mh.block("boundary");
39  size_t k = Xh.degree(), l = Mh.degree(), d = omega.dimension();
40  check_macro(l == k-1 || l == k || l == k+1,
41  "invalid (k,l) = ("<<k<<","<<l<<")");
42  space Xhs(omega, "P"+itos(k+1)+"d"),
43  Zh (omega, "P0"),
44  Mht(omega, "trace(P"+itos(k)+"d)");
45  space Yh = Xh*Xhs*Xh*Mht*Zh;
46  trial x(Yh), lambda(Mh);
47  test y(Yh), mu(Mh);
48  auto u = x[0], us = x[1], ut = x[2], deltat = x[3], zeta = x[4];
49  auto v = y[0], vs = y[1], vt = y[2], gammat = y[3], xi = y[4];
50  integrate_option iopt;
51  iopt.invert = true;
52  form inv_a = integrate(dot(grad_h(us),a(d)*grad_h(vs))
53  + dot(grad_h(us)-grad_h(u),a(d)*(grad_h(vs)-grad_h(v)))
54  + (us-u)*xi + (vs-v)*zeta + (ut-us)*(vt-vs)
55  + on_local_sides(beta/h_local()*deltat*gammat
56  + u*dot(a(d)*grad_h(vs),normal())
57  + v*dot(a(d)*grad_h(us),normal())
58  + (deltat - us + ut)*(gammat - vs + vt)), iopt);
59  form b = integrate(omega,on_local_sides(-mu*dot(a(d)*grad_h(us),normal())));
60  // TODO: remove omega, autodetect it ?
61  field lh = integrate (f(d)*v);
62  // TODO: f -> -f ? check the FV...
63  field lambda_h(Mh,0);
64  form s = b*inv_a*trans(b);
65  field rh = b*(inv_a*lh);
66  problem p (s);
67  p.solve (rh, lambda_h);
68  field xh = inv_a*(lh - b.trans_mult(lambda_h));
69  dout << catchmark("beta") << beta << endl
70  << catchmark("u") << xh[0]
71  << catchmark("us") << xh[1]
72  << catchmark("ut") << xh[2]
73  << catchmark("delta") << xh[3]
74  << catchmark("zeta") << xh[4]
75  << catchmark("lambda") << lambda_h;
76 }
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::normal
details::field_expr_v2_nonlinear_terminal_function< details::normal_pseudo_function< Float > > normal()
normal: see the expression page for the full documentation
Definition: field_expr_terminal.h:439
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::h_local
details::field_expr_v2_nonlinear_terminal_function< details::h_local_pseudo_function< Float > > h_local()
h_local: see the expression page for the full documentation
Definition: field_expr_terminal.h:527
mkgeo_ball.f
f
Definition: mkgeo_ball.sh:221
space
see the space page for the full documentation
rheolef.h
rheolef - reference manual
p
Definition: sphere.icc:25
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition: integrate_option.h:125
rheolef::on_local_sides
std::enable_if< details::is_field_expr_v2_variational_arg< Expr >::value,details::field_expr_quadrature_on_sides< Expr > >::type on_local_sides(const Expr &expr)
on_local_sides(expr): see the expression page for the full documentation
Definition: field_expr_quadrature.h:321
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
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::integrate_option::invert
bool invert
Definition: integrate_option.h:168
test
see the test page for the full documentation
problem
see the problem page for the full documentation
u
Definition: leveque.h:25
Float
see the Float page for the full documentation
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
mkgeo_ball.a
a
Definition: mkgeo_ball.sh:151
diffusion_isotropic.h
Tensor diffusion – isotropic case.
mkgeo_contraction.mu
mu
Definition: mkgeo_contraction.sh:193
rheolef::grad_h
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_h(const Expr &expr)
grad_h(uh): see the expression page for the full documentation
Definition: field_expr_terminal.h:949
main
int main(int argc, char **argv)
Definition: dirichlet_hho.cc:30
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
rheolef::itos
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
rheolef::trans
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition: csr.h:455
rheolef::std
Definition: vec_expr_v2.h:391
geo
see the geo page for the full documentation
dirichlet_homogeneous.h
The Poisson problem with homogeneous Dirichlet boundary conditions – right-hand-side and boundary con...
rk::beta
Float beta[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:60
lambda
Definition: yield_slip_circle.h:34