///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
//! @examplefile yield_slip1.icc The yield slip problem -- class body
#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);
}
