/usr/share/doc/rheolef-doc/examples/neumann-laplace.cc is in rheolef-doc 5.93-2.
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 | #include "rheolef.h"
using namespace rheolef;
using namespace std;
#include "neumann-laplace-assembly.h"
size_t N;
Float f (const point& x) { return 1; }
Float g (const point& x) { return -0.5/N; }
int main(int argc, char**argv) {
geo omega (argv[1]);
N = omega.dimension();
space Xh (omega, argv[2]);
form m (Xh, Xh, "mass");
form a (Xh, Xh, "grad_grad");
field b = m*field(Xh,1.0);
csr<Float> A = neumann_laplace_assembly (a.uu, b.u);
field fh = interpolate(Xh, f);
space Wh (omega, omega["boundary"], argv[2]);
field gh = interpolate(Wh, g);
form mb (Wh, Xh, "mass");
field lh = m*fh + mb*gh;
vec<Float> L(lh.u.size()+1, 0.0);
for (size_t i = 0; i < L.size()-1; i++) L.at(i) = lh.u.at(i);
L.at(L.size()-1) = 0;
vec<Float> U (L.size());
ssk<Float> fact_A = ldlt(A);
U = fact_A.solve(L);
field uh(Xh);
for (size_t i = 0; i < U.size()-1; i++) uh.u.at(i) = U.at(i);
Float lambda = U.at(U.size()-1);
cout << setprecision(numeric_limits<Float>::digits10)
<< catchmark("u") << uh
<< catchmark("lambda") << lambda << endl;
}
|