yield_slip_residue.cc

The yield slip problem – residue computation.

The yield slip problem – residue computation

#include "rheolef.h"
using namespace rheolef;
using namespace std;
#include "projection.h"
int main(int argc, char**argv) {
  environment rheolef (argc, argv);
  Float tol = (argc > 1) ? atof(argv[1]) : 1e-12;
  bool sym  = (argc > 2) && (string(argv[2]) == "-symmetric") ? true : false;
  Float S, n, Cf, r;
  field uh, lambda_h;
  din >> catchmark("S")  >> S
      >> catchmark("n")  >> n
      >> catchmark("Cf") >> Cf
      >> catchmark("r")  >> r 
      >> catchmark("u")  >> uh
      >> catchmark("lambda") >> lambda_h;
  space Xh = uh.get_space();
  space Wh = lambda_h.get_space();
  geo omega = Xh.get_geo();
  geo boundary = omega["boundary"];
  trial u(Xh), lambda(Wh);
  test  v(Xh),     mu(Wh);
  form m  = integrate (u*v),
       a  = integrate (dot(grad(u),grad(v))),
       mb = integrate (lambda*mu),
       b  = integrate (boundary, u*mu);
  field lh = integrate(v);
  field r_lambda_h;
  if (!sym) {
    r_lambda_h = lazy_interpolate(Wh, 
      uh["boundary"] - compose(projection(S,n,Cf), lambda_h));
  } else {
    field beta_h = lambda_h + r*uh[boundary];
    field mr_lambda_h = b*uh - integrate(mu*compose(projection(S,n,Cf,r), beta_h));
    r_lambda_h = field(Wh);
    problem pmb (mb);
    pmb.solve (mr_lambda_h, r_lambda_h);
  }
  field mruh = a*uh - lh + b.trans_mult(lambda_h);
  field ruh(Xh);
  problem pm (m);
  pm.solve (mruh, ruh);
  ruh["boundary"] = 0;
  Float residue_u      = sqrt(dual(mruh,ruh));
  Float residue_lambda = sqrt(mb(r_lambda_h, r_lambda_h));
  dout << "residue_u      = " << residue_u      << endl
       << "residue_lambda = " << residue_lambda << endl;
  Float residue = max(residue_u, residue_lambda);
  return (residue <= tol) ? 0 : 1;
}