neumann-laplace-lambda.cc

Go to the documentation of this file.

#include "rheolef.h"
using namespace rheolef;
using namespace std;
size_t d;
Float f (const point& x) { return 1; }
Float g (const point& x) { return -0.5/d; }
int main(int argc, char**argv) {
  environment rheolef (argc, argv);
  geo omega (argv[1]);
  d = omega.dimension();
  space Xh (omega, argv[2]);
  trial u (Xh); test v (Xh);
  form  a  = integrate (dot(grad(u),grad(v)));
  field b  = integrate(v);
  field lh = integrate (f*v) + integrate ("boundary", g*v); 
  form  A = {{   a,      b},
             { trans(b), 0}};
  field B =  {   lh,     0};
  A.set_symmetry(true);
  space IR = space::real();
  space Yh = Xh*IR;
  field Uh (Yh,0); 
  problem p (A);
  p.solve (B, Uh);
  dout << catchmark("u") << Uh[0];
  // get lambda:
  Float lambda = (Uh.comm().rank() == Uh.constraint_process_rank()) ? Uh[1].dof(0) : 0;
#ifdef _RHEOLEF_HAVE_MPI
  mpi::broadcast (Uh.comm(), lambda, Uh.constraint_process_rank()); 
#endif // _RHEOLEF_HAVE_MPI
  dout << catchmark("lambda") << lambda << endl;
  check_macro(fabs(lambda) < 1e-10, "unexpected lambda="<<lambda);
}