incompressible-elasticity.cc

The incompressible elasticity problem for the embankment geometry.

The incompressible elasticity problem for the embankment geometry

#include "rheolef.h"
using namespace rheolef;
using namespace std;
#include "embankment.icc"
int main(int argc, char**argv) {
  environment rheolef (argc, argv);
  geo omega (argv[1]);
  Float inv_lambda = (argc > 2 ? atof(argv[2]) : 0);
  size_t d = omega.dimension();
  space Xh = embankment_space(omega, "P2");
  space Qh (omega, "P1");
  point f (0,0,0);
  f [d-1] = -1;
  trial u (Xh), p (Qh); test v (Xh), q (Qh);
  field lh = integrate (dot(f,v));
  form  a  = integrate (2*ddot(D(u),D(v)));
  form  b  = integrate (-div(u)*q);
  form  c  = integrate (inv_lambda*p*q);
  field uh (Xh, 0), ph (Qh, 0);
  problem_mixed elasticity (a, b, c);
  elasticity.solve (lh, field(Qh,0), uh, ph);
  dout << catchmark("inv_lambda")  << inv_lambda << endl
       << catchmark("u")  << uh
       << catchmark("p")  << ph;
}