dirichlet_hho.cc

The Poisson problem by the hybrid high order method

#include "rheolef.h"
using namespace rheolef;
using namespace std;
#include "dirichlet_homogeneous.h"
#include "diffusion_isotropic.h"
int main(int argc, char**argv) {
environment rheolef (argc, argv);
geo omega (argv[1]);
string Pkd = (argc > 2) ? argv[2] : "P1d",
Pld = (argc > 3) ? argv[3] : Pkd;
Float beta = (argc > 4) ? atof(argv[4]) : 1;
space Xh (omega, Pld),
Mh (omega["sides"], Pkd);
Mh.block("boundary");
size_t k = Xh.degree(), l = Mh.degree(), d = omega.dimension();
check_macro(l == k-1 || l == k || l == k+1,
"invalid (k,l) = ("<<k<<","<<l<<")");
space Xhs(omega, "P"+itos(k+1)+"d"),
Zh (omega, "P0"),
Mht(omega, "trace(P"+itos(k)+"d)");
space Yh = Xh*Xhs*Xh*Mht*Zh;
trial x(Yh), lambda(Mh);
test y(Yh), mu(Mh);
auto u = x[0], us = x[1], ut = x[2], deltat = x[3], zeta = x[4];
auto v = y[0], vs = y[1], vt = y[2], gammat = y[3], xi = y[4];
integrate_option iopt;
iopt.invert = true;
form inv_a = integrate(dot(grad_h(us),a(d)*grad_h(vs))
+ dot(grad_h(us)-grad_h(u),a(d)*(grad_h(vs)-grad_h(v)))
+ (us-u)*xi + (vs-v)*zeta + (ut-us)*(vt-vs)
+ on_local_sides(beta/h_local()*deltat*gammat
+ u*dot(a(d)*grad_h(vs),normal())
+ v*dot(a(d)*grad_h(us),normal())
+ (deltat - us + ut)*(gammat - vs + vt)), iopt);
form b = integrate(omega,on_local_sides(-mu*dot(a(d)*grad_h(us),normal())));
// TODO: remove omega, autodetect it ?
field lh = integrate (f(d)*v);
// TODO: f -> -f ? check the FV...
field lambda_h(Mh,0);
form s = b*inv_a*trans(b);
field rh = b*(inv_a*lh);
problem p (s);
p.solve (rh, lambda_h);
field xh = inv_a*(lh - b.trans_mult(lambda_h));
dout << catchmark("beta") << beta << endl
<< catchmark("u") << xh[0]
<< catchmark("us") << xh[1]
<< catchmark("ut") << xh[2]
<< catchmark("delta") << xh[3]
<< catchmark("zeta") << xh[4]
<< catchmark("lambda") << lambda_h;
}
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