ClosestPointProjection.cpp

Go to the documentation of this file.

/** 
 * \file ClosestPointProjection.cpp
 * \ingroup adoc_plasticity
 * 
 * \brief Implementation of nested integration algorithm for plasticity
 *
 */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <ADOLCPlasticity.hpp>
namespace ADOLCPlasticity {
 
ClosestPointProjection::ClosestPointProjection() {
 
  if (!LogManager::checkIfChannelExist("ADOLCPlasticityWold")) {
    auto core_log = logging::core::get();
 
    core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(),
                                              "ADOLCPlasticityWold"));
    core_log->add_sink(LogManager::createSink(LogManager::getStrmSync(),
                                              "ADOLCPlasticitySync"));
    core_log->add_sink(LogManager::createSink(LogManager::getStrmSelf(),
                                              "ADOLCPlasticitySelf"));
 
    LogManager::setLog("ADOLCPlasticityWold");
    LogManager::setLog("ADOLCPlasticitySync");
    LogManager::setLog("ADOLCPlasticitySelf");
 
    MOFEM_LOG_TAG("ADOLCPlasticityWold", "ADOL-C Plasticity");
    MOFEM_LOG_TAG("ADOLCPlasticitySync", "ADOL-C Plasticity");
    MOFEM_LOG_TAG("ADOLCPlasticitySelf", "ADOL-C Plasticity");
  }
 
  tapesTags = {0, 1, 2}; //< set tapes default tags
 
  implHessianW = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "-impl_hessianW", &implHessianW, 0);
 
}
 
MoFEMErrorCode ClosestPointProjection::recordW() {
  MoFEMFunctionBegin;
 
  auto plastic_strain = getPlasticStrain();
  auto total_strain = getTotalStrain();
  auto internal_variables = getInternalVariables();
 
  trace_on(tapesTags[ClosestPointProjection::W]);
  {
    // active variables
    a_plasticStrain.resize(6, false);
    for (int ii = 0; ii < 6; ii++) {
      a_plasticStrain[ii] <<= plastic_strain[ii];
    }
    a_sTrain.resize(6, false);
    for (int ii = 0; ii < 6; ii++) {
      a_sTrain[ii] <<= total_strain[ii];
    }
    a_internalVariables.resize(nbInternalVariables, false);
    for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
      a_internalVariables[ii] <<= internal_variables[ii];
    }
    // evaluete functions
    CHKERR freeHemholtzEnergy();
    a_w >>= w;
  }
  trace_off();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ClosestPointProjection::recordY() {
  MoFEMFunctionBegin;
 
  auto stress = getStress();
  auto internal_fluxes = getInternalFluxes();
 
  trace_on(tapesTags[ClosestPointProjection::Y]);
  {
    // active variables
    a_sTress.resize(6, false);
    for (int ii = 0; ii < 6; ii++) {
      a_sTress[ii] <<= stress[ii];
    }
    a_internalFluxes.resize(nbInternalVariables, false);
    for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
      a_internalFluxes[ii] <<= internal_fluxes[ii];
    }
    // evaluete functions
    CHKERR yieldFunction();
    // return variables
    a_y >>= y;
  }
  trace_off();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ClosestPointProjection::recordH() {
  MoFEMFunctionBegin;
 
  auto stress = getStress();
  auto internal_fluxes = getInternalFluxes();
 
  trace_on(tapesTags[ClosestPointProjection::H]);
  {
    // active variables
    a_sTress.resize(6, false);
    for (int ii = 0; ii < 6; ii++) {
      a_sTress[ii] <<= stress[ii];
    }
    a_internalFluxes.resize(nbInternalVariables, false);
    for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
      a_internalFluxes[ii] <<= internal_fluxes[ii];
    }
    // evaluete functions
    CHKERR flowPotential();
    // return variables
    a_h >>= h;
  }
  trace_off();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ClosestPointProjection::playW() {
  MoFEMFunctionBegin;
  int r;
  int adloc_return_value = 0;
  r = ::function(tapesTags[ClosestPointProjection::W], 1,
                 activeVariablesW.size(), &activeVariablesW[0], &w);
  if (r < adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  r = gradient(tapesTags[ClosestPointProjection::W], activeVariablesW.size(),
               &activeVariablesW[0], &gradientW[0]);
  if (r < adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
 
  hessianW.resize(activeVariablesW.size(), activeVariablesW.size(), false);
  hessianW.clear();
  vector<double *> hessian_w(activeVariablesW.size());
  for (unsigned int dd = 0; dd < activeVariablesW.size(); dd++) {
    hessian_w[dd] = &hessianW(dd, 0);
  }
  if (!implHessianW) {
    // call adolc hessian (energy, active variables)
    r = hessian(tapesTags[ClosestPointProjection::W], activeVariablesW.size(),
                &activeVariablesW[0], &hessian_w[0]);
    if (r < adloc_return_value) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
              "ADOL-C function evaluation with error");
    }
  } else {
    // call hard coded hessian depending on material model
    CHKERR codedHessianW(hessian_w);
  }
  
  partialWStrainPlasticStrain.resize(6, 6, false);
  for (int ii = 0; ii < 6; ii++) {
    for (int jj = 0; jj < 6; jj++) {
      partialWStrainPlasticStrain(ii, jj) = hessianW(6 + ii, jj);
    }
  }
  C.resize(6, 6, false);
  for (int ii = 0; ii < 6; ii++) {
    for (int jj = 0; jj <= ii; jj++) {
      C(ii, jj) = hessianW(6 + ii, 6 + jj);
    }
  }
  D.resize(nbInternalVariables, nbInternalVariables, false);
  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
    for (unsigned int jj = 0; jj <= ii; jj++) {
      D(ii, jj) = hessianW(12 + ii, 12 + jj);
    }
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
ClosestPointProjection::playW_NoHessian() {
  MoFEMFunctionBegin;
  int r;
  int adloc_return_value = 0;
  r = ::function(tapesTags[ClosestPointProjection::W], 1,
                 activeVariablesW.size(), &activeVariablesW[0], &w);
  if (r < adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  r = gradient(tapesTags[ClosestPointProjection::W], activeVariablesW.size(),
               &activeVariablesW[0], &gradientW[0]);
  if (r < adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ClosestPointProjection::playY() {
  MoFEMFunctionBegin;
  auto active_variables_yh = getActiveVariablesYH();
  int r;
  int adloc_return_value = 0;
  r = ::function(tapesTags[ClosestPointProjection::Y], 1,
                 active_variables_yh.size(), &active_variables_yh[0], &y);
  if (r < adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  gradientY.resize(active_variables_yh.size());
  r = gradient(tapesTags[ClosestPointProjection::Y], active_variables_yh.size(),
               &active_variables_yh[0], &gradientY[0]);
  if (r < adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  partialYSigma =
      VectorAdaptor(6, ublas::shallow_array_adaptor<double>(6, &gradientY[0]));
  partialYFlux = VectorAdaptor(
      nbInternalVariables,
      ublas::shallow_array_adaptor<double>(nbInternalVariables, &gradientY[6]));
  MoFEMFunctionReturn(0);
}
MoFEMErrorCode
ClosestPointProjection::playY_NoGradient() {
  MoFEMFunctionBegin;
  auto active_variables_yh = getActiveVariablesYH();
  int r;
  int adloc_return_value = 0;
  r = ::function(tapesTags[ClosestPointProjection::Y], 1,
                 active_variables_yh.size(), &active_variables_yh[0], &y);
  if (r < adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  MoFEMFunctionReturn(0);
}
MoFEMErrorCode ClosestPointProjection::playH() {
  MoFEMFunctionBegin;
  auto active_variables_yh = getActiveVariablesYH();
  int r;
  int adloc_return_value = 0;
  r = ::function(tapesTags[ClosestPointProjection::H], 1,
                 active_variables_yh.size(), &active_variables_yh[0], &h);
  if (r < adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  gradientH.resize(active_variables_yh.size());
  r = gradient(tapesTags[ClosestPointProjection::H], active_variables_yh.size(),
               &active_variables_yh[0], &gradientH[0]);
  if (r < adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  hessianH.resize(6 + nbInternalVariables, 6 + nbInternalVariables, false);
  hessianH.clear();
  vector<double *> hessian_h(active_variables_yh.size());
  for (int dd = 0; dd < active_variables_yh.size(); dd++) {
    hessian_h[dd] = &hessianH(dd, 0);
  }
  r = hessian(tapesTags[ClosestPointProjection::H], active_variables_yh.size(),
              &active_variables_yh[0], &hessian_h[0]);
  if (r < adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  partialHSigma =
      VectorAdaptor(6, ublas::shallow_array_adaptor<double>(6, &gradientH[0]));
  partialHFlux = VectorAdaptor(
      nbInternalVariables,
      ublas::shallow_array_adaptor<double>(nbInternalVariables, &gradientH[6]));
  partialHFlux *= -1;
  partial2HSigma.resize(6, 6, false);
  for (int ii = 0; ii < 6; ii++) {
    for (int jj = 0; jj <= ii; jj++) {
      partial2HSigma(ii, jj) = hessianH(ii, jj);
    }
  }
  partial2HFlux.resize(nbInternalVariables, nbInternalVariables, false);
  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
    for (unsigned int jj = 0; jj <= ii; jj++) {
      partial2HFlux(ii, jj) = -hessianH(6 + ii, 6 + jj);
    }
  }
  partial2HSigmaFlux.resize(6, nbInternalVariables, false);
  partial2HSigmaFlux.clear();
  for (int ii = 0; ii < 6; ii++) {
    for (unsigned int jj = 0; jj < nbInternalVariables; jj++) {
      partial2HSigmaFlux(ii, jj) = -hessianH(6 + jj, ii);
    }
  }
  MoFEMFunctionReturn(0);
}
MoFEMErrorCode ClosestPointProjection::playH_NoHessian() {
  MoFEMFunctionBegin;
  auto active_variables_yh = getActiveVariablesYH();
  int r;
  int adloc_return_value = 0;
  r = ::function(tapesTags[ClosestPointProjection::H], 1,
                 active_variables_yh.size(), &active_variables_yh[0], &h);
  if (r < adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  gradientH.resize(active_variables_yh.size());
  r = gradient(tapesTags[ClosestPointProjection::H], active_variables_yh.size(),
               &active_variables_yh[0], &gradientH[0]);
  if (r < adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
            "ADOL-C function evaluation with error");
  }
  partialHSigma =
      VectorAdaptor(6, ublas::shallow_array_adaptor<double>(6, &gradientH[0]));
  partialHFlux = VectorAdaptor(
      nbInternalVariables,
      ublas::shallow_array_adaptor<double>(nbInternalVariables, &gradientH[6]));
  partialHFlux *= -1;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ClosestPointProjection::createMatAVecR() {
  MoFEMFunctionBegin;
  PetscInt n;
  n = 1 + 6 + nbInternalVariables;
  dataA.resize(n, n, false);
  Mat A_tmp;
  Vec R_tmp, Chi_tmp, dChi_tmp;
  CHKERR MatCreateSeqDense(PETSC_COMM_SELF, n, n, &dataA(0, 0), &A_tmp);
  CHKERR VecCreateSeq(PETSC_COMM_SELF, n, &R_tmp);
  CHKERR VecCreateSeq(PETSC_COMM_SELF, n, &Chi_tmp);
  CHKERR VecCreateSeq(PETSC_COMM_SELF, n, &dChi_tmp);
 
  A = SmartPetscObj<Mat>(A_tmp);
  R = SmartPetscObj<Vec>(R_tmp);
  Chi = SmartPetscObj<Vec>(Chi_tmp);
  dChi = SmartPetscObj<Vec>(dChi_tmp);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
ClosestPointProjection::evaluatePotentials() {
  MoFEMFunctionBegin;
  CHKERR recordW();
  CHKERR recordY();
  CHKERR recordH();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ClosestPointProjection::playPotentials() {
  MoFEMFunctionBegin;
  CHKERR playW();
  CHKERR playY();
  CHKERR playH();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ClosestPointProjection::playPotentials_NoHessian() {
  MoFEMFunctionBegin;
  CHKERR playW_NoHessian();
  CHKERR playY_NoGradient();
  CHKERR playH_NoHessian();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
ClosestPointProjection::calculateR(Vec R) {
  MoFEMFunctionBegin;
 
  auto plastic_strain = getPlasticStrain();
  auto internal_variables = getInternalVariables();
 
  // Residual
  double *array;
  CHKERR VecGetArray(R, &array);
  for (int ii = 0; ii < 6; ii++) {
    array[ii] = -plastic_strain[ii] + plasticStrain0[ii] +
                deltaGamma * partialHSigma[ii];
  }
  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
    array[6 + ii] = -internal_variables[ii] + internalVariables0[ii] +
                    deltaGamma * partialHFlux[ii];
  }
  array[6 + nbInternalVariables] = y;
  CHKERR VecRestoreArray(R, &array);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ClosestPointProjection::calculateA() {
  MoFEMFunctionBegin;
  // matrix A
  CHKERR MatZeroEntries(A);
  // row 0 (Strain)
  MatrixDouble a00 = prod(partial2HSigma, partialWStrainPlasticStrain);
  MatrixDouble a01 = prod(partial2HSigmaFlux, D);
  for (int ii = 0; ii < 6; ii++) {
    for (int jj = 0; jj < 6; jj++) {
      dataA(ii, jj) = deltaGamma * a00(ii, jj);
    }
    for (unsigned int jj = 0; jj < nbInternalVariables; jj++) {
      dataA(ii, 6 + jj) = deltaGamma * a01(ii, jj);
    }
    dataA(ii, ii) -= 1;
  }
  for (int ii = 0; ii < 6; ii++) {
    dataA(ii, 6 + nbInternalVariables) = partialHSigma[ii];
  }
  // row 1 (Fluxes)
  MatrixDouble a11 = prod(partial2HFlux, D);
  MatrixDouble a10 =
      prod(trans(partial2HSigmaFlux), partialWStrainPlasticStrain);
  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
    for (int jj = 0; jj < 6; jj++) {
      dataA(6 + ii, jj) += deltaGamma * a10(ii, jj);
    }
    for (unsigned int jj = 0; jj < nbInternalVariables; jj++) {
      dataA(6 + ii, 6 + jj) += deltaGamma * a11(ii, jj);
    }
    dataA(6 + ii, 6 + ii) -= 1;
  }
  for (unsigned int ii = 0; ii < nbInternalVariables; ii++) {
    dataA(6 + ii, 6 + nbInternalVariables) = partialHFlux[ii];
  }
  // row 3 (Plastic multiplier)
  VectorDouble partial_y_sigma_c =
      prod(trans(partialWStrainPlasticStrain), partialYSigma);
  VectorDouble partial_y_flux_d = prod(trans(D), partialYFlux);
  for (unsigned int dd = 0; dd < partial_y_sigma_c.size(); dd++) {
    dataA(6 + nbInternalVariables, dd) = partial_y_sigma_c[dd];
  }
  for (unsigned int dd = 0; dd < partial_y_flux_d.size(); dd++) {
    dataA(6 + nbInternalVariables, 6 + dd) = partial_y_flux_d[dd];
  }
  dataA(6 + nbInternalVariables, 6 + nbInternalVariables) = 0;
  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ClosestPointProjection::snesCreate() {
  MoFEMFunctionBegin;
  sNes = createSNES(PETSC_COMM_SELF);
  CHKERR SNESSetFunction(sNes, R, ADOLCPlasticityRes, (void *)this);
  CHKERR SNESSetJacobian(sNes, A, A, ADOLCPlasticityJac, (void *)this);
 
  KSP ksp;
  PC pc;
  SNESLineSearch line_search;
  CHKERR SNESGetKSP(sNes, &ksp);
  CHKERR KSPGetPC(ksp, &pc);
  CHKERR KSPSetType(ksp, KSPPREONLY);
  CHKERR PCSetType(pc, PCLU);
  CHKERR SNESSetTolerances(sNes, 1e-8, 1e-6, 1e-12, 20, 1000);
  CHKERR SNESGetLineSearch(sNes, &line_search);
  // CHKERR SNESLineSearchSetType(line_search, SNESLINESEARCHBASIC);
  CHKERR SNESLineSearchSetType(line_search, SNESLINESEARCHBT);
  // CHKERR SNESLineSearchSetType(line_search, SNESLINESEARCHL2);
  
  CHKERR SNESAppendOptionsPrefix(sNes, "cp_");
  CHKERR SNESSetFromOptions(sNes);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
ClosestPointProjection::solveClosestProjection() {
  MoFEMFunctionBegin;
 
  auto plastic_strain = getPlasticStrain();
  auto internal_variables = getInternalVariables();
 
  {
    double *array;
    CHKERR VecGetArray(Chi, &array);
    for (int ii = 0; ii < 6; ii++) {
      array[ii] = plastic_strain[ii];
    }
    int nb_internal_variables = nbInternalVariables;
    for (unsigned int ii = 0; ii < nb_internal_variables; ii++) {
      array[6 + ii] = internal_variables[ii];
    }
    deltaGamma = 0;
    CHKERR VecRestoreArray(Chi, &array);
  }
 
  CHKERR SNESSolve(sNes, PETSC_NULLPTR, Chi);
  SNESConvergedReason reason;
  CHKERR SNESGetConvergedReason(sNes, &reason);
  if (reason < 0) {
    int its;
    CHKERR SNESGetIterationNumber(sNes, &its);
    MOFEM_LOG_C(
        "ADOLCPlasticitySelf", Sev::warning,
        "Plasticity Closest Point Projection - number of Newton iterations = %d",
        its);
    plastic_strain = plasticStrain0;
    internal_variables = internalVariables0;
    deltaGamma = 0;
  } else {
    double *array;
    CHKERR VecGetArray(Chi, &array);
    for (int ii = 0; ii < 6; ii++) {
      plastic_strain[ii] = array[ii];
    }
    int nb_internal_variables = nbInternalVariables;
    for (unsigned int ii = 0; ii < nb_internal_variables; ii++) {
      internal_variables[ii] = array[6 + ii];
    }
    deltaGamma = array[6 + nb_internal_variables];
    CHKERR VecRestoreArray(Chi, &array);
  }
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode
ClosestPointProjection::consistentTangent() {
  MoFEMFunctionBegin;
  CHKERR playPotentials();
  CHKERR calculateA();
  Ep.resize(6, 6, false);
  MatrixDouble ep_row; // inrement internal varaibles as result of incremen of
                       // total strain
  ep_row = deltaGamma * prod(partial2HSigma, C);
  MatrixDouble alpha_row;
  alpha_row = deltaGamma * prod(trans(partial2HSigmaFlux), C);
  VectorDouble y_row;
  y_row = prod(C, partialYSigma);
  const int n = 6 + nbInternalVariables;
  // Iterate strains
  for (int dd = 0; dd < 6; dd++) {
    CHKERR VecZeroEntries(R);
    double *array;
    CHKERR VecGetArray(R, &array);
    {
      for (auto ii = 0; ii < 6; ii++) {
        array[ii] = ep_row(ii, dd);
      }
      for (auto ii = 0; ii < nbInternalVariables; ii++) {
        array[6 + ii] = alpha_row(ii, dd);
      }
      array[n] = y_row[dd];
    }
    CHKERR VecRestoreArray(R, &array);
    CHKERR VecAssemblyBegin(R);
    CHKERR VecAssemblyEnd(R);
 
    KSP ksp;
    CHKERR SNESGetKSP(sNes, &ksp);
    CHKERR KSPSolve(ksp, R, dChi);
    CHKERR VecGetArray(dChi, &array);
    {
      for (auto ii = 0; ii < 6; ii++) {
        Ep(ii, dd) = array[ii];
      }
    }
    CHKERR VecRestoreArray(dChi, &array);
  }
  Cp.resize(6, 6, false);
  noalias(Cp) = prod(C, Ep);
  Cep.resize(6, 6, false);
  noalias(Cep) = C + Cp;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ADOLCPlasticityRes(SNES snes, Vec chi, Vec r, void *ctx) {
  MoFEMFunctionBegin;
  ClosestPointProjection *cp;
  cp = (ClosestPointProjection *)ctx;
 
  auto plastic_strain = cp->getPlasticStrain();
  auto internal_variables = cp->getInternalVariables();
 
  {
    const double *array;
    CHKERR VecGetArrayRead(chi, &array);
    for (auto ii = 0; ii < 6; ii++) {
      plastic_strain[ii] = array[ii];
    }
    for (auto ii = 0; ii < cp->nbInternalVariables; ii++) {
      internal_variables[ii] = array[6 + ii];
    }
    cp->deltaGamma = array[6 + cp->nbInternalVariables];
    CHKERR VecRestoreArrayRead(chi, &array);
  }
  CHKERR cp->playPotentials_NoHessian();
  CHKERR cp->calculateR(r);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ADOLCPlasticityJac(SNES snes, Vec chi, Mat A, Mat,
                                         void *ctx) {
  MoFEMFunctionBegin;
  ClosestPointProjection *cp;
  cp = (ClosestPointProjection *)ctx;
 
  auto plastic_strain = cp->getPlasticStrain();
  auto internal_variables = cp->getInternalVariables();
 
  {
    const double *array;
    CHKERR VecGetArrayRead(chi, &array);
    for (auto ii = 0; ii < 6; ii++) {
      plastic_strain[ii] = array[ii];
    }
    for (auto ii = 0; ii < cp->nbInternalVariables; ii++) {
      internal_variables[ii] = array[6 + ii];
    }
    cp->deltaGamma = array[6 + cp->nbInternalVariables];
    CHKERR VecRestoreArrayRead(chi, &array);
  }
  CHKERR cp->playPotentials();
  CHKERR cp->calculateA();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode ClosestPointProjection::recordTapes() {
  MoFEMFunctionBegin;
  VectorDouble active(12 + nbInternalVariables);
  VectorDouble gradient(12 + nbInternalVariables);
  auto tmp_active = getVectorAdaptor(&(active[0]), 12 + nbInternalVariables);
  auto tmp_gradient =
      getVectorAdaptor(&(gradient[0]), 12 + nbInternalVariables);
  activeVariablesW.swap(tmp_active);
  gradientW.swap(tmp_gradient);
  CHKERR evaluatePotentials();
  MoFEMFunctionReturn(0);
}
 
} // namespace ADOLCPlasticity