/usr/share/doc/rheolef-doc/examples/convect.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 | #include "rheolef.h"
using namespace rheolef;
using namespace std;
#include "rotating-hill.h"
int main (int argc, char **argv) {
environment rheolef (argc,argv);
geo omega (argv[1]);
string approx = (argc > 2) ? argv[2] : "P1";
Float nu = (argc > 3) ? atof(argv[3]) : 1e-2;
size_t n_max = (argc > 4) ? atoi(argv[4]) : 50;
size_t d = omega.dimension();
Float delta_t = 2*acos(-1.)/n_max;
space Vh (omega, approx, "vector");
field uh = interpolate (Vh, u(d));
space Xh (omega, approx);
Xh.block ("boundary");
field phi_h = interpolate (Xh, phi(d,nu,0));
characteristic X (-delta_t*uh);
quadrature_option_type qopt;
qopt.set_family (quadrature_option_type::gauss_lobatto);
qopt.set_order (Xh.degree());
trial phi (Xh); test psi (Xh);
branch event ("t","phi");
dout << catchmark("nu") << nu << endl
<< event (0, phi_h);
for (size_t n = 1; n <= n_max; n++) {
Float t = n*delta_t;
Float c1 = 1 + delta_t*phi::sigma(d,nu,t);
Float c2 = delta_t*nu;
form a = integrate (c1*phi*psi + c2*dot(grad(phi),grad(psi)), qopt);
field lh = integrate (compose(phi_h, X)*psi, qopt);
solver sa (a.uu());
phi_h.set_u() = sa.solve (lh.u() - a.ub()*phi_h.b());
dout << event (t, phi_h);
}
}
|