2
using namespace rheolef;
4
int main (int argc, char **argv) {
5
environment rheolef (argc, argv);
7
size_t n_max = (argc > 2) ? atoi(argv[2]) : 100;
8
Float delta_t = 0.5/n_max;
9
space Xh (omega, "P1");
10
Xh.block ("boundary");
11
form m (Xh, Xh, "mass");
12
form a (Xh, Xh, "grad_grad");
13
form c = m + delta_t*a;
14
solver sc = ldlt (c.uu());
15
field lh = riesz (Xh, 1);
17
branch event ("t","u");
18
dout << event (0, uh);
19
for (size_t n = 1; n <= n_max; n++) {
20
field kh = m*uh + delta_t*lh;
21
uh.set_u() = sc.solve (kh.u() - c.ub()*uh.b());
22
dout << event (Float(n)*delta_t, uh);