/usr/share/doc/rheolef-doc/examples/stokes_taylor_error_dg.cc is in rheolef-doc 6.7-6.
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 | #include "rheolef.h"
using namespace rheolef;
using namespace std;
#include "taylor_exact.icc"
int main(int argc, char**argv) {
environment rheolef(argc, argv);
Float err_linf_expected = (argc > 1) ? atof(argv[1]) : 1e+38;
bool dump = (argc > 2);
field uh, ph; din >> uh >> ph;
space Xh = uh.get_space(),
Qh = ph.get_space();
geo omega = Xh.get_geo();
size_t k = Xh.degree();
size_t d = omega.dimension();
quadrature_option_type qopt;
qopt.set_family(quadrature_option_type::gauss);
qopt.set_order(2*k+1);
Float meas_omega = integrate (omega);
Float moy_ph = integrate (omega, ph, qopt);
ph = ph - moy_ph/meas_omega; // ph with zero average value
string high_approx = "P"+itos(k+1)+"d";
space Xh1 (omega, high_approx, "vector"),
Qh1 (omega, high_approx);
field euh = interpolate (Xh1, uh-u_exact());
field eph = interpolate (Qh1, ph-p_exact());
Float err_u_l2 = sqrt(integrate (omega, norm2(uh-u_exact()), qopt));
Float err_u_linf = euh.max_abs();
Float err_u_h1 = sqrt(integrate (omega, norm2(grad_h(euh)), qopt)
+ integrate (omega.sides(), (1/h_local())*norm2(jump(euh)), qopt));
Float err_p_l2 = sqrt(integrate (omega, sqr(ph-p_exact()), qopt));
Float err_p_linf = eph.max_abs();
derr << "err_u_l2 = " << err_u_l2 << endl
<< "err_u_linf = " << err_u_linf << endl
<< "err_u_h1 = " << err_u_h1 << endl
<< "err_p_l2 = " << err_p_l2 << endl
<< "err_p_linf = " << err_p_linf << endl;
if (dump) {
dout << catchmark("uh") << uh
<< catchmark("u") << interpolate (Xh, u_exact())
<< catchmark("eu") << euh
<< catchmark("ph") << ph
<< catchmark("p") << interpolate (Qh, p_exact())
<< catchmark("ep") << eph;
}
return (max(err_u_linf,err_p_linf) <= err_linf_expected) ? 0 : 1;
}
|