SmallStrainPlasticity.cpp

Go to the documentation of this file.

/** \fiele SmallStrainPlasticity.cpp
 * \brief Operators and data structures for small strain plasticity
 * \ingroup small_strain_plasticity
 *
 * Implementation of small strain plasticity
 *
 * This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM 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 Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
*/
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
 
#include <adolc/adolc.h>
#include <SmallStrainPlasticity.hpp>
 
PetscErrorCode SmallStrainPlasticity::createInternalVariablesTag() {
  PetscFunctionBegin;
  int def_length = 0;
  rval = mField.get_moab().tag_get_handle(
    "PLASTIC_STRAIN",def_length,MB_TYPE_DOUBLE,
    thPlasticStrain,MB_TAG_CREAT|MB_TAG_SPARSE|MB_TAG_VARLEN,NULL
  ); CHKERRQ_MOAB(rval);
  rval = mField.get_moab().tag_get_handle(
    "INTERNAL_VARIABLES",def_length,MB_TYPE_DOUBLE,
    thInternalVariables,MB_TAG_CREAT|MB_TAG_SPARSE|MB_TAG_VARLEN,NULL
  ); CHKERRQ_MOAB(rval);
  PetscFunctionReturn(0);
}
 
SmallStrainPlasticity::OpGetDataAtGaussPts::OpGetDataAtGaussPts(const string field_name,
  vector<VectorDouble > &values_at_gauss_pts,
  vector<MatrixDouble > &gardient_at_gauss_pts
):
VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
valuesAtGaussPts(values_at_gauss_pts),
gradientAtGaussPts(gardient_at_gauss_pts),
zeroAtType(MBVERTEX) {
}
 
PetscErrorCode SmallStrainPlasticity::OpGetDataAtGaussPts::doWork(
  int side,EntityType type,DataForcesAndSourcesCore::EntData &data
) {
  PetscFunctionBegin;
  try {
    int nb_dofs = data.getFieldData().size();
    if(nb_dofs == 0) {
      PetscFunctionReturn(0);
    }
    int nb_gauss_pts = data.getN().size1();
    int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
    //initialize
    if(type == zeroAtType) {
      valuesAtGaussPts.resize(nb_gauss_pts);
      gradientAtGaussPts.resize(nb_gauss_pts);
      for(int gg = 0;gg<nb_gauss_pts;gg++) {
        valuesAtGaussPts[gg].resize(rank,false);
        gradientAtGaussPts[gg].resize(rank,3,false);
        valuesAtGaussPts[gg].clear();
        gradientAtGaussPts[gg].clear();
      }
    }
    VectorDouble& values = data.getFieldData();
    for(int gg = 0;gg<nb_gauss_pts;gg++) {
      VectorDouble &values_at_gauss_pts = valuesAtGaussPts[gg];
      MatrixDouble &gradient_at_gauss_pts = gradientAtGaussPts[gg];
      VectorAdaptor N = data.getN(gg,nb_dofs/rank);
      MatrixAdaptor diffN = data.getDiffN(gg,nb_dofs/rank);
      for(int dd = 0;dd<nb_dofs/rank;dd++) {
        for(int rr1 = 0;rr1<rank;rr1++) {
          const double v = values[rank*dd+rr1];
          values_at_gauss_pts[rr1] += N[dd]*v;
          for(int rr2 = 0;rr2<3;rr2++) {
            gradient_at_gauss_pts(rr1,rr2) += diffN(dd,rr2)*v;
          }
        }
      }
    }
    // cerr << rowFieldName << endl;
    // cerr << side << " " << type << endl;
    // cerr << values << endl;
    // cerr << valuesAtGaussPts[0] << endl;
    // cerr << gradientAtGaussPts[0] << endl;
    // cerr << endl;
  } catch (const std::exception& ex) {
    ostringstream ss;
    ss << "throw in method: " << ex.what() << endl;
    SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
  }
  PetscFunctionReturn(0);
}
 
SmallStrainPlasticity::OpGetCommonDataAtGaussPts::OpGetCommonDataAtGaussPts(
  const string field_name,CommonData &common_data
):
OpGetDataAtGaussPts(
  field_name,common_data.dataAtGaussPts[field_name],common_data.gradAtGaussPts[field_name]
) {
}
 
SmallStrainPlasticity::OpUpdate::OpUpdate(
  string field_name,
  CommonData &common_data
):
VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
commonData(common_data) {
}
 
PetscErrorCode SmallStrainPlasticity::OpUpdate::doWork(
  int side,EntityType type,DataForcesAndSourcesCore::EntData &data
) {
  PetscFunctionBegin;
 
  try {
    //do it only once, no need to repeat this for edges,faces or tets
    if(type != MBVERTEX) PetscFunctionReturn(0);
 
    int nb_dofs = data.getFieldData().size();
    if(nb_dofs==0) PetscFunctionReturn(0);
    int nb_gauss_pts = data.getN().size1();
 
    for(int gg = 0;gg<nb_gauss_pts;gg++) {
 
      if(commonData.deltaGamma[gg]>0) {
 
        VectorAdaptor plastic_strain = VectorAdaptor(
          6,ublas::shallow_array_adaptor<double>(6,&commonData.plasticStrainPtr[gg*6])
        );
        int nb_internal_variables = commonData.internalVariables[gg].size();
        VectorAdaptor internal_variables = VectorAdaptor(
          nb_internal_variables,
          ublas::shallow_array_adaptor<double>(
            nb_internal_variables,&commonData.internalVariablesPtr[gg*nb_internal_variables]
          )
        );
 
        // cerr << plastic_strain << " " << commonData.plasticStrain[gg] << endl;
        // cerr << internal_variables << " " << commonData.internalVariables[gg] << endl;
 
        plastic_strain = commonData.plasticStrain[gg];
        internal_variables = commonData.internalVariables[gg];
      }
 
    }
 
  } catch (const std::exception& ex) {
    ostringstream ss;
    ss << "throw in method: " << ex.what() << endl;
    SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
  }
 
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::MakeB::makeB(const MatrixAdaptor &diffN,MatrixDouble &B) {
  PetscFunctionBegin;
  unsigned int nb_dofs = diffN.size1();
  if(B.size2()!=3*nb_dofs) {
    B.resize(6,3*nb_dofs,false);
  }
  B.clear();
  for(int dd = 0;dd<nb_dofs;dd++) {
    const double diff[] = {
      diffN(dd,0), diffN(dd,1), diffN(dd,2)
    };
    const int dd3 = 3*dd;
    for(int rr = 0;rr<3;rr++) {
      B(rr,dd3+rr) = diff[rr];
    }
    // gamma_xy
    B(3,dd3+0) = diff[1];
    B(3,dd3+1) = diff[0];
    // gamma_yz
    B(4,dd3+1) = diff[2];
    B(4,dd3+2) = diff[1];
    // gamma_xz
    B(5,dd3+0) = diff[2];
    B(5,dd3+2) = diff[0];
  }
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::MakeB::addVolumetricB(
  const MatrixAdaptor &diffN,MatrixDouble &volB,double alpha
) {
  PetscFunctionBegin;
  int nb_dofs = diffN.size1();
  for(int dd = 0;dd<nb_dofs;dd++) {
    const double diff[] = {
      diffN(dd,0), diffN(dd,1), diffN(dd,2)
    };
    const int dd3 = 3*dd;
    for(int rr1 = 0;rr1<3;rr1++) {
      for(int rr2 = 0;rr2<3;rr2++) {
        volB(rr1,dd3+rr2) += alpha*diff[rr2]/3.;
      }
    }
  }
  PetscFunctionReturn(0);
}
 
SmallStrainPlasticity::OpAssembleRhs::OpAssembleRhs(
  string field_name,
  CommonData &common_data
):
VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,UserDataOperator::OPROW),
commonData(common_data) {
}
 
PetscErrorCode SmallStrainPlasticity::OpAssembleRhs::doWork(
  int side,EntityType type,DataForcesAndSourcesCore::EntData &data
) {
  PetscFunctionBegin;
 
  try {
 
    int nb_dofs = data.getFieldData().size();
    if(nb_dofs==0) PetscFunctionReturn(0);
    int nb_gauss_pts = data.getN().size1();
    nF.resize(nb_dofs,false);
    nF.clear();
 
    if(commonData.bBar) {
      double v = 0;
      volB.resize(6,nb_dofs,false);
      volB.clear();
      for(int gg = 0;gg<nb_gauss_pts;gg++) {
        double val = getVolume()*getGaussPts()(3,gg);
        if(getHoGaussPtsDetJac().size()>0) {
          val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
        }
        v += val;
        const MatrixAdaptor &diffN = data.getDiffN(gg,nb_dofs/3);
        ierr = addVolumetricB(diffN,volB,val); CHKERRQ(ierr);
      }
      volB /= v;
    }
 
    for(int gg = 0;gg<nb_gauss_pts;gg++) {
 
      double val = getVolume()*getGaussPts()(3,gg);
      if(getHoGaussPtsDetJac().size()>0) {
        val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
      }
      const MatrixAdaptor &diffN = data.getDiffN(gg,nb_dofs/3);
      ierr = makeB(diffN,B); CHKERRQ(ierr);
      if(commonData.bBar) {
        ierr = addVolumetricB(diffN,B,-1); CHKERRQ(ierr);
        B += volB;
      }
      // cerr << gg << endl;
      // cerr << diffN << endl;
      // cerr << B << endl;
      // cerr << commonData.sTress[gg] << endl;
      // cerr << prod(trans(B),commonData.sTress[gg]) << endl;
 
      noalias(nF) += val*prod(trans(B),commonData.sTress[gg]);
      // cerr << "nF " << nF << endl;
 
    }
 
    // cerr << side << " " << type << " " << data.getN().size1() << " "<< data.getN().size2() << endl;
 
    int *indices_ptr;
    indices_ptr = &data.getIndices()[0];
 
    ierr = VecSetValues(
      getFEMethod()->snes_f,
      nb_dofs,
      indices_ptr,
      &nF[0],
      ADD_VALUES
    ); CHKERRQ(ierr);
 
 
  } catch (const std::exception& ex) {
    ostringstream ss;
    ss << "throw in method: " << ex.what() << endl;
    SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
  }
 
  PetscFunctionReturn(0);
}
 
SmallStrainPlasticity::OpAssembleLhs::OpAssembleLhs(
  string field_name,
  CommonData &common_data
):
VolumeElementForcesAndSourcesCore::UserDataOperator(field_name,field_name,UserDataOperator::OPROWCOL),
commonData(common_data) {
  //  sYmm = false;
}
 
PetscErrorCode SmallStrainPlasticity::OpAssembleLhs::doWork(
  int row_side,int col_side,
  EntityType row_type,EntityType col_type,
  DataForcesAndSourcesCore::EntData &row_data,
  DataForcesAndSourcesCore::EntData &col_data
) {
  PetscFunctionBegin;
  try {
 
    int nb_dofs_row = row_data.getFieldData().size();
    if(nb_dofs_row==0) PetscFunctionReturn(0);
    int nb_dofs_col = col_data.getFieldData().size();
    if(nb_dofs_col==0) PetscFunctionReturn(0);
 
    K.resize(nb_dofs_row,nb_dofs_col,false);
    K.clear();
 
    int nb_gauss_pts = row_data.getN().size1();
 
    if(commonData.bBar) {
      double v = 0;
      rowVolB.resize(6,nb_dofs_row,false);
      rowVolB.clear();
      colVolB.resize(6,nb_dofs_col,false);
      colVolB.clear();
      for(int gg = 0;gg<nb_gauss_pts;gg++) {
        double val = getVolume()*getGaussPts()(3,gg);
        if(getHoGaussPtsDetJac().size()>0) {
          val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
        }
        v += val;
        const MatrixAdaptor &diffN_row = row_data.getDiffN(gg,nb_dofs_row/3);
        const MatrixAdaptor &diffN_col = col_data.getDiffN(gg,nb_dofs_col/3);
        ierr = addVolumetricB(diffN_row,rowVolB,val); CHKERRQ(ierr);
        ierr = addVolumetricB(diffN_col,colVolB,val); CHKERRQ(ierr);
      }
      rowVolB /= v;
      colVolB /= v;
    }
 
    for(int gg = 0;gg<nb_gauss_pts;gg++) {
 
      double val = getVolume()*getGaussPts()(3,gg);
      if(getHoGaussPtsDetJac().size()>0) {
        val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
      }
 
      const MatrixAdaptor &diffN_row = row_data.getDiffN(gg,nb_dofs_row/3);
      const MatrixAdaptor &diffN_col = col_data.getDiffN(gg,nb_dofs_col/3);
      ierr = makeB(diffN_row,rowB); CHKERRQ(ierr);
      ierr = makeB(diffN_col,colB); CHKERRQ(ierr);
      if(commonData.bBar) {
        ierr = addVolumetricB(diffN_row,rowB,-1); CHKERRQ(ierr); // B-bar
        ierr = addVolumetricB(diffN_col,colB,-1); CHKERRQ(ierr);
        rowB += rowVolB;
        colB += colVolB;
      }
 
      CB.resize(6,nb_dofs_col,false);
      cblas_dgemm(
        CblasRowMajor,CblasNoTrans,CblasNoTrans,
        6,nb_dofs_col,6,
        1.,
        &commonData.materialTangentOperator[gg](0,0),6,
        &colB(0,0),nb_dofs_col,
        0,
        &CB(0,0),nb_dofs_col
      );
      cblas_dgemm(
        CblasRowMajor,CblasTrans,CblasNoTrans,
        nb_dofs_row,nb_dofs_col,6,
        val,
        &rowB(0,0),nb_dofs_row,
        &CB(0,0),nb_dofs_col,
        1,
        &K(0,0),nb_dofs_col
      );
      //noalias(CB) = prod(commonData.materialTangentOperator[gg],colB);
      //noalias(K) += val*prod(trans(rowB),CB);
 
    }
 
    // cerr << K << endl;
    // cerr << row_data.getIndices() << endl;
    // cerr << col_data.getIndices() << endl;
    // cerr << getFEMethod()->snes_A << endl;
    // cerr << getFEMethod()->snes_B << endl;
 
    int *row_indices_ptr,*col_indices_ptr;
    row_indices_ptr = &row_data.getIndices()[0];
    col_indices_ptr = &col_data.getIndices()[0];
 
    ierr = MatSetValues(
      getFEMethod()->snes_B,
      nb_dofs_row,row_indices_ptr,
      nb_dofs_col,col_indices_ptr,
      &K(0,0),ADD_VALUES
    ); CHKERRQ(ierr);
 
    //is symmetric
    if(row_side != col_side || row_type != col_type) {
 
      transK.resize(nb_dofs_col,nb_dofs_row,false);
      noalias(transK) = trans(K);
      ierr = MatSetValues(
        getFEMethod()->snes_B,
        nb_dofs_col,col_indices_ptr,
        nb_dofs_row,row_indices_ptr,
        &transK(0,0),ADD_VALUES
      ); CHKERRQ(ierr);
 
    }
 
  } catch (const std::exception& ex) {
    ostringstream ss;
    ss << "throw in method: " << ex.what() << endl;
    SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
  }
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::OpCalculateStress::getTags() {
  PetscFunctionBegin;
  
  rval = mField.get_moab().tag_get_handle(
    "PLASTIC_STRAIN",thPlasticStrain
  ); CHKERRQ_MOAB(rval);
  rval = mField.get_moab().tag_get_handle(
    "INTERNAL_VARIABLES",thInternalVariables
  ); CHKERRQ_MOAB(rval);
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::OpCalculateStress::setTagsData(
  const EntityHandle tet,const int nb_gauss_pts,const int nb_internal_variables
) {
  PetscFunctionBegin;
  
  VectorDouble v;
  {
    rval = mField.get_moab().tag_get_by_ptr(
      thPlasticStrain,&tet,1,(const void **)&commonData.plasticStrainPtr,&commonData.plasticStrainSize
    );
    if(rval != MB_SUCCESS || commonData.plasticStrainSize != 6*nb_gauss_pts) {
      v.resize(6*nb_gauss_pts,false);
      v.clear();
      int tag_size[1];
      tag_size[0] = v.size();
      void const* tag_data[] = { &v[0] };
      rval = mField.get_moab().tag_set_by_ptr(
        thPlasticStrain,&tet,1,tag_data,tag_size
      ); CHKERRQ_MOAB(rval);
      rval = mField.get_moab().tag_get_by_ptr(
        thPlasticStrain,&tet,1,(const void **)&commonData.plasticStrainPtr,&commonData.plasticStrainSize
      ); CHKERRQ_MOAB(rval);
    }
  }
  if(nb_internal_variables>0) {
    rval = mField.get_moab().tag_get_by_ptr(
      thInternalVariables,&tet,1,(const void **)&commonData.internalVariablesPtr,&commonData.internalVariablesSize
    );
    if(rval != MB_SUCCESS || commonData.internalVariablesSize != nb_internal_variables*nb_gauss_pts) {
      v.resize(nb_internal_variables*nb_gauss_pts,false);
      v.clear();
      int tag_size[1];
      tag_size[0] = v.size();
      void const* tag_data[] = { &v[0] };
      rval = mField.get_moab().tag_set_by_ptr(
        thInternalVariables,&tet,1,tag_data,tag_size
      ); CHKERRQ_MOAB(rval);
      rval = mField.get_moab().tag_get_by_ptr(
        thInternalVariables,&tet,1,(const void **)&commonData.internalVariablesPtr,&commonData.internalVariablesSize
      );
    }
  }
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::OpCalculateStress::doWork(
  int side,EntityType type,DataForcesAndSourcesCore::EntData &data
) {
  PetscFunctionBegin;
 
  try {
 
    //do it only once, no need to repeat this for edges,faces or tets
    if(type != MBVERTEX) PetscFunctionReturn(0);
 
    int nb_dofs = data.getFieldData().size();
    if(nb_dofs==0) PetscFunctionReturn(0);
 
    int nb_gauss_pts = data.getN().size1();
    int nb_internal_variables = cP.internalVariables.size();
    if(!initCp) {
      ierr = getTags(); CHKERRQ(ierr);
    }
    EntityHandle tet = getNumeredEntFiniteElementPtr()->getEnt();
    ierr = setTagsData(tet,nb_gauss_pts,nb_internal_variables); CHKERRQ(ierr);
 
    commonData.sTress.resize(nb_gauss_pts);
    commonData.materialTangentOperator.resize(nb_gauss_pts);
    commonData.plasticStrain.resize(nb_gauss_pts);
    commonData.internalVariables.resize(nb_gauss_pts);
    commonData.internalFluxes.resize(nb_gauss_pts);
    commonData.deltaGamma.resize(nb_gauss_pts);
 
    // B-bar
    double ave_tr_strain = 0;
    if(commonData.bBar) {
      double v = 0;
      for(int gg = 0;gg<nb_gauss_pts;gg++) {
        double val = getVolume()*getGaussPts()(3,gg);
        if(getHoGaussPtsDetJac().size()>0) {
          val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
        }
        v += val;
        for(int ii = 0;ii<3;ii++) {
          ave_tr_strain += val*commonData.gradAtGaussPts[rowFieldName][gg](ii,ii)/3.;
        }
        // cerr << commonData.gradAtGaussPts[rowFieldName][gg] << endl;
      }
      ave_tr_strain /= v;
      // cerr << "v " << v << " ave_tr_strain " <<  ave_tr_strain << endl;
    }
 
    for(int gg = 0;gg<nb_gauss_pts;gg++) {
 
      VectorAdaptor plastic_strain = VectorAdaptor(
        6,ublas::shallow_array_adaptor<double>(6,&commonData.plasticStrainPtr[gg*6])
      );
      VectorAdaptor internal_variables = VectorAdaptor(
        nb_internal_variables,
        ublas::shallow_array_adaptor<double>(
          nb_internal_variables,&commonData.internalVariablesPtr[gg*nb_internal_variables]
        )
      );
 
      cP.sTrain.resize(6,false);
      {
        double tr_strain = 0;
        for(int ii = 0;ii<3;ii++) {
          cP.sTrain[ii] = commonData.gradAtGaussPts[rowFieldName][gg](ii,ii);
          tr_strain += cP.sTrain[ii]/3;
        }
        if(commonData.bBar) {
          for(int ii = 0;ii<3;ii++) {
            cP.sTrain[ii] += ave_tr_strain-tr_strain;
          }
        }
        cP.sTrain[3] =
        commonData.gradAtGaussPts[rowFieldName][gg](0,1)+
        commonData.gradAtGaussPts[rowFieldName][gg](1,0);
        cP.sTrain[4] =
        commonData.gradAtGaussPts[rowFieldName][gg](1,2)+
        commonData.gradAtGaussPts[rowFieldName][gg](2,1);
        cP.sTrain[5] =
        commonData.gradAtGaussPts[rowFieldName][gg](0,2)+
        commonData.gradAtGaussPts[rowFieldName][gg](2,0);
        // cerr << "Grad  " << commonData.gradAtGaussPts[rowFieldName][gg] << endl;
        // cerr << cP.sTrain << endl;
      }
 
 
      {
        cP.plasticStrain0.resize(6,false);
        noalias(cP.plasticStrain0) = plastic_strain;
        cP.internalVariables0.resize(nb_internal_variables,false);
        noalias(cP.internalVariables0) = internal_variables;
        cP.plasticStrain.resize(6,false);
        noalias(cP.plasticStrain) = plastic_strain;
        cP.internalVariables.resize(nb_internal_variables,false);
        noalias(cP.internalVariables) = internal_variables;
        cP.deltaGamma = 0;
        commonData.sTress[gg].clear();
        commonData.internalFluxes[gg].clear();
      }
 
      if(!initCp) {
        cP.gG = 0;
        ierr = cP.evaluatePotentials(); CHKERRQ(ierr);
      } else {
        cP.gG = 1;
        ierr = cP.setActiveVariablesW(); CHKERRQ(ierr);
        if(
          getFEMethod()->snes_ctx ==
          SnesMethod::CTX_SNESSETJACOBIAN &&
          (!isLinear || !hessianWCalculated)
        ) {
          ierr = cP.pLayW(); CHKERRQ(ierr);
          hessianWCalculated = true;
        } else {
          ierr = cP.pLayW_NoHessian(); CHKERRQ(ierr);
        }
        ierr = cP.setActiveVariablesYH(); CHKERRQ(ierr);
        ierr = cP.pLayY_NoGradient(); CHKERRQ(ierr);
        // cerr << "C " << cP.C << endl;
      }
 
      double y = cP.y;
      if(y>0) {
        // ierr = PCFactorSetMatSolverPackage(cP.pC,MATSOLVERPETSC); CHKERRQ(ierr);
        // ierr = KSPMonitorCancel(cP.kSp); CHKERRQ(ierr);
        // ierr = SNESMonitorCancel(cP.sNes); CHKERRQ(ierr);
        // ierr = SNESSetTolerances(cP.sNes,cP.tOl,1e-12,0,20,1000); CHKERRQ(ierr);
        // ierr = SNESLineSearchSetType(cP.lineSearch,SNESLINESEARCHBT); CHKERRQ(ierr);
        ierr = cP.solveColasetProjection(); CHKERRQ(ierr);
      }
 
      {
        commonData.sTress[gg].resize(6,false);
        noalias(commonData.sTress[gg]) = cP.sTress;
        commonData.plasticStrain[gg].resize(6,false);
        noalias(commonData.plasticStrain[gg]) = cP.plasticStrain;
        commonData.internalVariables[gg].resize(nb_internal_variables,false);
        noalias(commonData.internalVariables[gg]) = cP.internalVariables;
        commonData.internalFluxes[gg].resize(nb_internal_variables,false);
        noalias(commonData.internalFluxes[gg]) = cP.internalFluxes;
        commonData.deltaGamma[gg] = cP.deltaGamma;
      }
 
      if(getFEMethod()->snes_ctx == SnesMethod::CTX_SNESSETJACOBIAN) {
        int iter;
        ierr = SNESGetIterationNumber(getFEMethod()->snes,&iter); CHKERRQ(ierr);
        commonData.materialTangentOperator[gg].resize(6,6,false);
        if(iter>0 && cP.deltaGamma>0) {
          ierr = cP.consistentTangent(); CHKERRQ(ierr);
          noalias(commonData.materialTangentOperator[gg]) = cP.Cep;
        } else {
          noalias(commonData.materialTangentOperator[gg]) = cP.C;
        }
      }
 
      if(getFEMethod()->snes_ctx == SnesMethod::CTX_SNESSETFUNCTION) {
        int iter;
        ierr = SNESGetIterationNumber(getFEMethod()->snes,&iter); CHKERRQ(ierr);
        if(iter>1) {
          ierr = SNESSetLagJacobian(getFEMethod()->snes,1); CHKERRQ(ierr);
        }
      }
 
 
      initCp = true;
 
    }
 
  } catch (const std::exception& ex) {
    ostringstream ss;
    ss << "throw in method: " << ex.what() << endl;
    SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
  }
 
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::setActiveVariablesW() {
  PetscFunctionBegin;
  activeVariablesW.resize(12+internalVariables.size(),false);
  for(int ii = 0;ii<6;ii++) {
    activeVariablesW[ii] = sTrain[ii];
  }
  for(int ii = 0;ii<6;ii++) {
    activeVariablesW[6+ii] = plasticStrain[ii];
  }
  for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
    activeVariablesW[12+ii] = internalVariables[ii];
  }
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::setActiveVariablesYH() {
  PetscFunctionBegin;
  activeVariablesYH.resize(6+internalVariables.size(),false);
  for(int ii = 0;ii<6;ii++) {
    activeVariablesYH[ii] = sTress[ii];
  }
  for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
    activeVariablesYH[6+ii] = internalFluxes[ii];
  }
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::rEcordW() {
  PetscFunctionBegin;
  trace_on(tAgs[0]);
  {
    //active variables
    a_sTrain.resize(6,false);
    for(int ii = 0;ii<6;ii++) {
      a_sTrain[ii] <<= sTrain[ii];
    }
    a_plasticStrain.resize(6,false);
    for(int ii = 0;ii<6;ii++) {
      a_plasticStrain[ii] <<= plasticStrain[ii];
    }
    a_internalVariables.resize(internalVariables.size(),false);
    for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
      a_internalVariables[ii] <<= internalVariables[ii];
    }
    //evaluete functions
    ierr = freeHemholtzEnergy(); CHKERRQ(ierr);
    //return variables
    a_w >>= w;
  }
  trace_off();
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::rEcordY() {
  PetscFunctionBegin;
  trace_on(tAgs[1]);
  {
    //active variables
    a_sTress.resize(6,false);
    for(int ii = 0;ii<6;ii++) {
      a_sTress[ii] <<= sTress[ii];
    }
    a_internalFluxes.resize(internalVariables.size(),false);
    for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
      a_internalFluxes[ii] <<= internalFluxes[ii];
    }
    //evaluete functions
    ierr = yieldFunction(); CHKERRQ(ierr);
    //return variables
    a_y >>= y;
  }
  trace_off();
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::rEcordH() {
  PetscFunctionBegin;
  trace_on(tAgs[2]);
  {
    //active variables
    a_sTress.resize(6,false);
    for(int ii = 0;ii<6;ii++) {
      a_sTress[ii] <<= sTress[ii];
    }
    a_internalFluxes.resize(internalVariables.size(),false);
    for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
      a_internalFluxes[ii] <<= internalFluxes[ii];
    }
    //evaluete functions
    ierr = flowPotential(); CHKERRQ(ierr);
    //return variables
    a_h >>= h;
  }
  trace_off();
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::pLayW() {
  PetscFunctionBegin;
  int r;
  int adloc_return_value = 0;
  r = ::function(tAgs[0],1,activeVariablesW.size(),&activeVariablesW[0],&w);
  if(r<adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
  }
  gradientW.resize(activeVariablesW.size(),false);
  r = gradient(tAgs[0],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);
  }
  r = hessian(tAgs[0],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");
  }
  sTress = VectorAdaptor(6,ublas::shallow_array_adaptor<double>(6,&gradientW[0]));
  internalFluxes = VectorAdaptor(
    internalVariables.size(),
    ublas::shallow_array_adaptor<double>(
      internalVariables.size(),
      &gradientW[12]
    )
  );
  C.resize(6,6,false);
  for(int ii = 0;ii<6;ii++) {
    for(int jj = 0;jj<=ii;jj++) {
      C(ii,jj) = hessianW(ii,jj);
    }
  }
  partialWStrainPlasticStrain.resize(6,6,false);
  for(int ii = 0;ii<6;ii++) {
    for(int jj = 0;jj<6;jj++) {
      partialWStrainPlasticStrain(ii,jj) = hessianW(6+jj,ii);
    }
  }
  D.resize(internalVariables.size(),internalVariables.size(),false);
  for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
    for(unsigned int jj = 0;jj<=ii;jj++) {
      D(ii,jj) = hessianW(12+ii,12+jj);
    }
  }
  PetscFunctionReturn(0);
}
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::pLayW_NoHessian() {
  PetscFunctionBegin;
  int r;
  int adloc_return_value = 0;
  r = ::function(tAgs[0],1,activeVariablesW.size(),&activeVariablesW[0],&w);
  if(r<adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
  }
  gradientW.resize(activeVariablesW.size(),false);
  r = gradient(tAgs[0],activeVariablesW.size(),&activeVariablesW[0],&gradientW[0]);
  if(r<adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
  }
  sTress = VectorAdaptor(6,ublas::shallow_array_adaptor<double>(6,&gradientW[0]));
  internalFluxes = VectorAdaptor(
    internalVariables.size(),
    ublas::shallow_array_adaptor<double>(
      internalVariables.size(),
      &gradientW[12]
    )
  );
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::pLayY() {
  PetscFunctionBegin;
  int r;
  int adloc_return_value = 0;
  r = ::function(tAgs[1],1,activeVariablesYH.size(),&activeVariablesYH[0],&y);
  if(r<adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
  }
  gradientY.resize(activeVariablesYH.size());
  r = gradient(tAgs[1],activeVariablesYH.size(),&activeVariablesYH[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(
    internalVariables.size(),
    ublas::shallow_array_adaptor<double>(
      internalVariables.size(),
      &gradientY[6]
    )
  );
  PetscFunctionReturn(0);
}
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::pLayY_NoGradient() {
  PetscFunctionBegin;
  int r;
  int adloc_return_value = 0;
  r = ::function(tAgs[1],1,activeVariablesYH.size(),&activeVariablesYH[0],&y);
  if(r<adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
  }
  PetscFunctionReturn(0);
}
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::pLayH() {
  PetscFunctionBegin;
  int r;
  int adloc_return_value = 0;
  r = ::function(tAgs[2],1,activeVariablesYH.size(),&activeVariablesYH[0],&h);
  if(r<adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
  }
  gradientH.resize(activeVariablesYH.size());
  r = gradient(tAgs[2],activeVariablesYH.size(),&activeVariablesYH[0],&gradientH[0]);
  if(r<adloc_return_value) {
    SETERRQ(PETSC_COMM_SELF,MOFEM_OPERATION_UNSUCCESSFUL,"ADOL-C function evaluation with error");
  }
  hessianH.resize(6+internalVariables.size(),6+internalVariables.size(),false);
  hessianH.clear();
  vector<double*> hessian_h(activeVariablesYH.size());
  for(int dd = 0;dd<activeVariablesYH.size();dd++) {
    hessian_h[dd] = &hessianH(dd,0);
  }
  r = hessian(tAgs[2],activeVariablesYH.size(),&activeVariablesYH[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(
    internalVariables.size(),
    ublas::shallow_array_adaptor<double>(
      internalVariables.size(),
      &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(internalVariables.size(),internalVariables.size(),false);
  for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
    for(unsigned int jj = 0;jj<=ii;jj++) {
      partial2HFlux(ii,jj) = -hessianH(6+ii,6+jj);
    }
  }
  partial2HSigmaFlux.resize(6,internalVariables.size(),false);
  partial2HSigmaFlux.clear();
  for(int ii = 0;ii<6;ii++) {
    for(unsigned int jj = 0;jj<internalVariables.size();jj++) {
      partial2HSigmaFlux(ii,jj) = -hessianH(6+jj,ii);
    }
  }
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::createMatAVecR() {
  PetscFunctionBegin;
  PetscInt n;
  n = 1+6+internalVariables.size();
  dataA.resize(n,n,false);
  ierr = MatCreateSeqDense(PETSC_COMM_SELF,n,n,&dataA(0,0),&A); CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&R); CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&Chi); CHKERRQ(ierr);
  ierr = VecCreateSeq(PETSC_COMM_SELF,n,&dChi); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::destroyMatAVecR() {
  PetscFunctionBegin;
  ierr = MatDestroy(&A); CHKERRQ(ierr);
  ierr = VecDestroy(&R); CHKERRQ(ierr);
  ierr = VecDestroy(&Chi); CHKERRQ(ierr);
  ierr = VecDestroy(&dChi); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::evaluatePotentials() {
  PetscFunctionBegin;
  if(gG == 0) {
    ierr = rEcordW(); CHKERRQ(ierr);
  }
  ierr = setActiveVariablesW(); CHKERRQ(ierr);
  ierr = pLayW(); CHKERRQ(ierr);
  if(gG == 0) {
    ierr = rEcordY(); CHKERRQ(ierr);
    ierr = rEcordH(); CHKERRQ(ierr);
  }
  ierr = setActiveVariablesYH(); CHKERRQ(ierr);
  ierr = pLayY(); CHKERRQ(ierr);
  ierr = pLayH(); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::cAlculateR(Vec R) {
  PetscFunctionBegin;
  // Residual
  double *array;
  ierr = VecGetArray(R,&array); CHKERRQ(ierr);
  for(int ii = 0;ii<6;ii++) {
    array[ii] = -plasticStrain[ii] + plasticStrain0[ii] + deltaGamma*partialHSigma[ii];
  }
  for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
    array[6+ii] = -internalVariables[ii] + internalVariables0[ii] + deltaGamma*partialHFlux[ii];
  }
  array[6+internalVariables.size()] = y;
  ierr = VecRestoreArray(R,&array); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::cAlculateA() {
  PetscFunctionBegin;
  // matrix A
  ierr = MatZeroEntries(A); CHKERRQ(ierr);
  // row 0
  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<internalVariables.size();jj++) {
      dataA(ii,6+jj) = deltaGamma*a01(ii,jj);
    }
    dataA(ii,ii) -= 1;
  }
  for(int ii = 0;ii<6;ii++) {
    dataA(ii,6+internalVariables.size()) = partialHSigma[ii];
  }
  // row 1
  MatrixDouble a11 = prod(partial2HFlux,D);
  MatrixDouble a10 = prod(trans(partial2HSigmaFlux),partialWStrainPlasticStrain);
  for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
    for(int jj = 0;jj<6;jj++) {
      dataA(6+ii,jj) += deltaGamma*a10(ii,jj);
    }
    for(unsigned int jj = 0;jj<internalVariables.size();jj++) {
      dataA(6+ii,6+jj) += deltaGamma*a11(ii,jj);
    }
    dataA(6+ii,6+ii) -= 1;
  }
  for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
    dataA(6+ii,6+internalVariables.size()) = partialHFlux[ii];
  }
  // row 3
  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+internalVariables.size(),dd) = partial_y_sigma_c[dd];
  }
  for(unsigned int dd = 0;dd<partial_y_flux_d.size();dd++) {
    dataA(6+internalVariables.size(),6+dd) = partial_y_flux_d[dd];
  }
  dataA(6+internalVariables.size(),6+internalVariables.size()) = 0;
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::snesCreate() {
  PetscFunctionBegin;
  ierr = SNESCreate(PETSC_COMM_SELF,&sNes); CHKERRQ(ierr);
  ierr = SNESSetFunction(sNes,R,SmallStrainPlasticityfRes,(void*)this); CHKERRQ(ierr);
  ierr = SNESSetJacobian(sNes,A,A,SmallStrainPlasticityfJac,(void*)this); CHKERRQ(ierr);
  ierr = SNESGetKSP(sNes,&kSp); CHKERRQ(ierr);
  ierr = KSPGetPC(kSp,&pC); CHKERRQ(ierr);
  ierr = KSPSetType(kSp,KSPPREONLY); CHKERRQ(ierr);
  ierr = PCSetType(pC,PCLU); CHKERRQ(ierr);
  ierr = SNESSetTolerances(sNes,tOl,1e-12,0,20,1000); CHKERRQ(ierr);
  ierr = SNESGetLineSearch(sNes,&lineSearch); CHKERRQ(ierr); CHKERRQ(ierr);
  //ierr = SNESLineSearchSetType(lineSearch,SNESLINESEARCHBASIC); CHKERRQ(ierr);
  ierr = SNESLineSearchSetType(lineSearch,SNESLINESEARCHBT); CHKERRQ(ierr);
 
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::snesDestroy() {
  PetscFunctionBegin;
  ierr = SNESDestroy(&sNes); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::solveColasetProjection() {
  PetscFunctionBegin;
  {
    double *array;
    ierr = VecGetArray(Chi,&array); CHKERRQ(ierr);
    for(int ii = 0;ii<6;ii++) {
      array[ii] = plasticStrain[ii];
    }
    int nb_internal_variables = internalVariables.size();
    for(unsigned int ii = 0;ii<nb_internal_variables;ii++) {
      array[6+ii] = internalVariables[ii];
    }
    deltaGamma = 0;
    ierr = VecRestoreArray(Chi,&array); CHKERRQ(ierr);
  }
  ierr = SNESSolve(sNes,PETSC_NULL,Chi); CHKERRQ(ierr);
  SNESConvergedReason reason;
  ierr = SNESGetConvergedReason(sNes,&reason); CHKERRQ(ierr);
  if(reason < 0) {
    int its;
    ierr = SNESGetIterationNumber(sNes,&its); CHKERRQ(ierr);
//    ierr = PetscPrintf(
//      PETSC_COMM_SELF,
//      "** WARNING Plasticity Closet Point Projection - number of Newton iterations = %D\n",
//      its
//    ); CHKERRQ(ierr);
    noalias(plasticStrain) = plasticStrain0;
    noalias(internalVariables) = internalVariables0;
    deltaGamma = 0;
  } else {
    double *array;
    ierr = VecGetArray(Chi,&array); CHKERRQ(ierr);
    for(int ii = 0;ii<6;ii++) {
      plasticStrain[ii] = array[ii];
    }
    int nb_internal_variables = internalVariables.size();
    for(unsigned int ii = 0;ii<nb_internal_variables;ii++) {
      internalVariables[ii] = array[6+ii];
    }
    deltaGamma = array[6+nb_internal_variables];
    ierr = VecRestoreArray(Chi,&array); CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticity::ClosestPointProjection::consistentTangent() {
  PetscFunctionBegin;
  ierr = evaluatePotentials(); CHKERRQ(ierr);
  ierr = cAlculateA(); CHKERRQ(ierr);
  Ep.resize(6,6,false);
  Lp.resize(internalVariables.size(),6,false);
  Yp.resize(6,false);
  MatrixDouble ep_row;
  ep_row = deltaGamma*prod(partial2HSigma,C);
  MatrixDouble alpha_row;
  alpha_row = deltaGamma*prod(trans(partial2HSigmaFlux),C);
  VectorDouble y_row;
  y_row = prod(trans(C),partialYSigma);
  const int n = 6+internalVariables.size();
  for(int dd = 0;dd<6;dd++) {
    ierr = VecZeroEntries(R); CHKERRQ(ierr);
    double *array;
    ierr = VecGetArray(R,&array); CHKERRQ(ierr);
    {
      for(int ii = 0;ii<6;ii++) {
        array[ii] = ep_row(ii,dd);
      }
      for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
        array[6+ii] = alpha_row(ii,dd);
      }
      array[n] = y_row[dd];
    }
    ierr = VecRestoreArray(R,&array); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(R); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(R); CHKERRQ(ierr);
    ierr = KSPSolve(kSp,R,dChi); CHKERRQ(ierr);
    ierr = VecGetArray(dChi,&array); CHKERRQ(ierr);
    {
      for(int ii = 0;ii<6;ii++) {
        Ep(ii,dd) = array[ii];
      }
      for(unsigned int ii = 0;ii<internalVariables.size();ii++) {
        Lp(ii,dd) = array[6+ii];
      }
      Yp[dd] = array[n];
    }
    ierr = VecRestoreArray(dChi,&array); CHKERRQ(ierr);
  }
  Cp.resize(6,6,false);
  noalias(Cp) = prod(C,Ep);
  Cep.resize(6,6,false);
  noalias(Cep) = C+Cp;
  // cout << endl;
  // cout << "A " << dataA << endl;
  // cout << "C " << C << endl;
  // cout << "Ep " << Ep << endl;
  // // cout << "Lp " << Lp << endl;
  // // cout << "Yp " << Yp << endl;
  // cout << "Cp " << Cp << endl;
  // cout << "Cep " << Cep << endl;
  PetscFunctionReturn(0);
}
 
 
PetscErrorCode SmallStrainPlasticityfRes(SNES snes,Vec chi,Vec r,void *ctx) {
  PetscFunctionBegin;
  SmallStrainPlasticity::ClosestPointProjection *cp;
  cp = (SmallStrainPlasticity::ClosestPointProjection*)ctx;
  
  //ierr = VecView(chi,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
  {
    #if PETSC_VERSION_GE(3,6,0)
    const double *array;
    ierr = VecGetArrayRead(chi,&array); CHKERRQ(ierr);
    #else
    const double *array;
    ierr = VecGetArray(chi,&array); CHKERRQ(ierr);
    #endif
    for(int ii = 0;ii<6;ii++) {
      cp->plasticStrain[ii] = array[ii];
    }
    for(unsigned int ii = 0;ii<cp->internalVariables.size();ii++) {
      cp->internalVariables[ii] = array[6+ii];
    }
    cp->deltaGamma = array[6+cp->internalVariables.size()];
    #if PETSC_VERSION_GE(3,6,0)
    ierr = VecRestoreArrayRead(chi,&array); CHKERRQ(ierr);
    #else
    ierr = VecRestoreArray(chi,&array); CHKERRQ(ierr);
    #endif
  }
  ierr = cp->evaluatePotentials(); CHKERRQ(ierr);
  ierr = cp->cAlculateR(r); CHKERRQ(ierr);
  // ierr = VecView(r,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
  // cout << "sTress " << cp->sTress << endl;
  // cout << "internalVariables " << cp->internalVariables << endl;
  // cout << "C " << cp->C << endl;
  // cout << "D " << cp->D << endl;
  // cout << "y " << cp->y << endl;
  // cout << "h " << cp->h << endl;
  // cout << "partialYSigma " << cp->partialYSigma << endl;
  // cout << "partialYFlux " << cp->partialYFlux << endl;
  // cout << "partialHSigma " << cp->partialHSigma << endl;
  // cout << "partialHFlux " << cp->partialHFlux << endl;
  // cout << "partial2HSigma " << cp->partial2HSigma << endl;
  // cout << "partial2HFlux " << cp->partial2HFlux << endl;
  // cout << "partial2HSigmaFlux " << cp->partial2HSigmaFlux << endl;
  // cout << "C " << cp->C << endl;
  // cout << "D " << cp->D << endl;
  // cout << "partialWStrainPlasticStrain " << cp->partialWStrainPlasticStrain << endl;
  // cout << "plasticStrain " << cp->plasticStrain << endl;
  // cout << "internalVariables " << cp->internalVariables << endl;
  // cout << "delta plasticStrain " << cp->plasticStrain-cp->plasticStrain0 << endl;
  // cout << "delta internalVariables " << cp->internalVariables-cp->internalVariables << endl;
  PetscFunctionReturn(0);
}
 
PetscErrorCode SmallStrainPlasticityfJac(SNES snes,Vec chi,Mat A,Mat,void *ctx) {
  PetscFunctionBegin;
  SmallStrainPlasticity::ClosestPointProjection *cp;
  cp = (SmallStrainPlasticity::ClosestPointProjection*)ctx;
  
  {
    #if PETSC_VERSION_GE(3,6,0)
    const double *array;
    ierr = VecGetArrayRead(chi,&array); CHKERRQ(ierr);
    #else
    double *array;
    ierr = VecGetArray(chi,&array); CHKERRQ(ierr);
    #endif
    for(int ii = 0;ii<6;ii++) {
      cp->plasticStrain[ii] = array[ii];
    }
    for(unsigned int ii = 0;ii<cp->internalVariables.size();ii++) {
      cp->internalVariables[ii] = array[6+ii];
    }
    cp->deltaGamma = array[6+cp->internalVariables.size()];
    #if PETSC_VERSION_GE(3,6,0)
    ierr = VecRestoreArrayRead(chi,&array); CHKERRQ(ierr);
    #else
    ierr = VecRestoreArray(chi,&array); CHKERRQ(ierr);
    #endif
  }
  ierr = cp->evaluatePotentials(); CHKERRQ(ierr);
  ierr = cp->cAlculateA(); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}