yield_slip1.icc

Go to the documentation of this file.

#include "d_projection_dx.h"
yield_slip::yield_slip (Float S1, Float n1, Float Cf1, Float r1, 
    const geo& omega, const geo& boundary1, string approx)
 : S(S1), n(n1), Cf(Cf1), r(r1), boundary(boundary1), Xh(), Wh(), Yh(),
   lh(), mkh(), m(), mb(), a(), b(), c1(), pmb(), pa(), pA()
{
  Xh = space (omega,    approx);
  Wh = space (boundary, approx);
  Yh = Xh*Wh;
  trial u (Xh), lambda(Wh);
  test  v (Xh),     mu(Wh);
  m  = integrate(u*v);
  mb = integrate(lambda*mu);
  a  = integrate(dot(grad(u),grad(v))) - r*integrate(boundary, u*v);
  b  = integrate(boundary, u*mu);
  lh = integrate(v);
  pmb = problem (mb);
  pa  = problem (a);
  field vh(Xh);
  pa.solve (lh, vh);
  mkh = b*vh;
}
field yield_slip::residue (const field& beta_h) const {
  field vh (Xh);
  field rhs = b.trans_mult (beta_h);
  pa.solve (rhs, vh);
  test mu (Wh);
  field c0h = integrate(mu*compose(projection(S,n,Cf,r), beta_h));
  field mrh = b*vh + c0h - mkh;
  return mrh;
}
void yield_slip::update_derivative (const field& beta_h) const {
  trial lambda (Wh); test mu (Wh);
  c1 = integrate (lambda*mu*compose(d_projection_dx(S,n,Cf,r), beta_h));
  form A = { { a, trans(b) },
             { b,     -c1  } };
  A.set_symmetry (c1.is_symmetric());
  pA = problem(A);
}