rve_mech_plas_periodic.cpp

Go to the documentation of this file.

/** \file small_strain_plasticity.cpp
 * \ingroup small_strain_plasticity
 * \brief Small strain plasticity example
 *
 */
 
/* 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/program_options.hpp>
using namespace std;
namespace po = boost::program_options;
 
#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 <MethodForForceScaling.hpp>
#include <TimeForceScale.hpp>
#include <VolumeCalculation.hpp>
#include <BCs_RVELagrange_Disp.hpp>
#include <BCs_RVELagrange_Trac.hpp>
#include "BCs_RVELagrange_Periodic.hpp"
#include <BCs_RVELagrange_Trac_Rigid_Trans.hpp>
 
#ifndef WITH_MODULE_SMALL_STRAIN_PLASTICITY
  #error "Install module small strain plasticity https://bitbucket.org/likask/mofem_um_small_strain_plasticity"
#endif
 
#include <adolc/adolc.h>
#include <SmallStrainPlasticity.hpp>
#include <SmallStrainPlasticityMaterialModels.hpp>
 
#include <Projection10NodeCoordsOnField.hpp>
#include <PostProcOnRefMesh.hpp>
 
#include <string>
 
using namespace boost::numeric;
 
 
 
 
static char help[] =
"...\n"
"\n";
 
 
 
 
//=================================================================================================================================
//Define class and multindex container to store data for traiangles on the boundary of the RVE (it cannot be defined within main)
//=================================================================================================================================
 
struct Face_CenPos_Handle {
  double xcoord, ycoord, zcoord;
  const EntityHandle  Tri_Hand;
  Face_CenPos_Handle(double _xcoord, double _ycoord,  double _zcoord,  const EntityHandle _Tri_Hand):xcoord(_xcoord),
  ycoord(_ycoord), zcoord(_zcoord), Tri_Hand(_Tri_Hand) {}
};
 
struct xcoord_tag {}; //tags to used in the multindex container
struct ycoord_tag {};
struct zcoord_tag {};
struct Tri_Hand_tag {};
struct Composite_xyzcoord {};
 
typedef multi_index_container<
Face_CenPos_Handle,
indexed_by<
ordered_non_unique<
tag<xcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord> >,
 
ordered_non_unique<
tag<ycoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord> >,
 
ordered_non_unique<
tag<zcoord_tag>, member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> >,
 
ordered_unique<
tag<Tri_Hand_tag>, member<Face_CenPos_Handle,const EntityHandle,&Face_CenPos_Handle::Tri_Hand> >,
 
ordered_unique<
tag<Composite_xyzcoord>,
composite_key<
Face_CenPos_Handle,
member<Face_CenPos_Handle,double,&Face_CenPos_Handle::xcoord>,
member<Face_CenPos_Handle,double,&Face_CenPos_Handle::ycoord>,
member<Face_CenPos_Handle,double,&Face_CenPos_Handle::zcoord> > >
> > Face_CenPos_Handle_multiIndex;
 
//============================================================================================
//============================================================================================
 
 
 
 
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
 
  try {
 
  moab::Core mb_instance;
  moab::Interface& moab = mb_instance;
  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
  if(pcomm == NULL) pcomm =  new ParallelComm(&moab,PETSC_COMM_WORLD);
 
  PetscBool flg = PETSC_TRUE;
  char mesh_file_name[255];
  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
  if(flg != PETSC_TRUE) {
    SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
  }
 
 PetscInt order;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
  if(flg != PETSC_TRUE) {
    order = 2;
  }
  PetscInt bubble_order;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_bubble_order",&bubble_order,&flg); CHKERRQ(ierr);
  if(flg != PETSC_TRUE) {
    bubble_order = order;
  }
 
  // use this if your mesh is partitioned and you run code on parts,
  // you can solve very big problems
  PetscBool is_partitioned = PETSC_FALSE;
  ierr = PetscOptionsGetBool(PETSC_NULL,"-my_is_partitioned",&is_partitioned,&flg); CHKERRQ(ierr);
 
  //Applied strain on the RVE (vector of length 6) strain=[xx, yy, zz, xy, xz, zy]^T
  double mygiven_strain[6];
  int nmax=6;
  ierr = PetscOptionsGetRealArray(PETSC_NULL,"-my_given_strain",mygiven_strain,&nmax,&flg); CHKERRQ(ierr);
  VectorDouble given_strain;
  given_strain.resize(6);
  cblas_dcopy(6, &mygiven_strain[0], 1, &given_strain(0), 1);
  cout<<"given_strain ="<<given_strain<<endl;
 
  //Read mesh to MOAB
  if(is_partitioned == PETSC_TRUE) {
    const char *option;
    option = "PARALLEL=BCAST_DELETE;PARALLEL_RESOLVE_SHARED_ENTS;PARTITION=PARALLEL_PARTITION;";
    rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
  } else {
    const char *option;
    option = "";
    rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
  }
 
 
 //Create MoFEM (Joseph) database
  MoFEM::Core core(moab);
  MoFEM::Interface& m_field = core;
 
  //ref meshset ref level 0
  ierr = m_field.seed_ref_level_3D(0,0); CHKERRQ(ierr);
 
  // stl::bitset see for more details
  BitRefLevel bit_level0;
  bit_level0.set(0);
  EntityHandle meshset_level0;
  rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERRQ_MOAB(rval);
  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
  ierr = m_field.get_entities_by_ref_level(bit_level0,BitRefLevel().set(),meshset_level0); CHKERRQ(ierr);
 
 
 
    //=======================================================================================================
    //Add Periodic Prisims Between Triangles on -ve and +ve faces to implement periodic bounary conditions
    //=======================================================================================================
 
    //Populating the Multi-index container with -ve triangles
    Range SurTrisNeg;
    ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
    Face_CenPos_Handle_multiIndex Face_CenPos_Handle_varNeg, Face_CenPos_Handle_varPos;
    double TriCen[3], coords_Tri[9];
    for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end();  it++) {
      const EntityHandle* conn_face;  int num_nodes_Tri;
 
      //get nodes attached to the face
      rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
      //get nodal coordinates
      rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
 
      //Find triangle centriod
      TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
      TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
      TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
 
      //round values to 3 disimal places
      if(TriCen[0]>=0) TriCen[0]=double(int(TriCen[0]*1000.0+0.5))/1000.0;  else TriCen[0]=double(int(TriCen[0]*1000.0-0.5))/1000.0; //-ve and +ve value
      if(TriCen[1]>=0) TriCen[1]=double(int(TriCen[1]*1000.0+0.5))/1000.0;  else TriCen[1]=double(int(TriCen[1]*1000.0-0.5))/1000.0;
      if(TriCen[2]>=0) TriCen[2]=double(int(TriCen[2]*1000.0+0.5))/1000.0;  else TriCen[2]=double(int(TriCen[2]*1000.0-0.5))/1000.0;
 
      //        cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << "   TriCen[1]= "<< TriCen[1] << "   TriCen[2]= "<< TriCen[2] <<endl;
      //fill the multi-index container with centriod coordinates and triangle handles
      Face_CenPos_Handle_varNeg.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
      //        for(int ii=0; ii<3; ii++) cout<<"TriCen "<<TriCen[ii]<<endl;
      //        cout<<endl<<endl;
    }
 
    //    double aaa;
    //    aaa=0.5011;
    //    cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*1000.0+0.5))/1000.0<<endl<<endl<<endl<<endl;
    //    aaa=-0.5011;
    //    cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*1000.0-0.5))/1000.0<<endl<<endl<<endl<<endl;
 
 
    //Populating the Multi-index container with +ve triangles
    Range SurTrisPos;
    ierr = m_field.get_cubit_msId_entities_by_dimension(102,SIDESET,2,SurTrisPos,true); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
    for(Range::iterator it = SurTrisPos.begin(); it!=SurTrisPos.end();  it++) {
      const EntityHandle* conn_face;  int num_nodes_Tri;
 
      //get nodes attached to the face
      rval = moab.get_connectivity(*it,conn_face,num_nodes_Tri,true); CHKERRQ_MOAB(rval);
      //get nodal coordinates
      rval = moab.get_coords(conn_face,num_nodes_Tri,coords_Tri); CHKERRQ_MOAB(rval);
 
      //Find triangle centriod
      TriCen[0]= (coords_Tri[0]+coords_Tri[3]+coords_Tri[6])/3.0;
      TriCen[1]= (coords_Tri[1]+coords_Tri[4]+coords_Tri[7])/3.0;
      TriCen[2]= (coords_Tri[2]+coords_Tri[5]+coords_Tri[8])/3.0;
 
      //round values to 3 disimal places
      if(TriCen[0]>=0) TriCen[0]=double(int(TriCen[0]*1000.0+0.5))/1000.0;  else TriCen[0]=double(int(TriCen[0]*1000.0-0.5))/1000.0;
      if(TriCen[1]>=0) TriCen[1]=double(int(TriCen[1]*1000.0+0.5))/1000.0;  else TriCen[1]=double(int(TriCen[1]*1000.0-0.5))/1000.0;
      if(TriCen[2]>=0) TriCen[2]=double(int(TriCen[2]*1000.0+0.5))/1000.0;  else TriCen[2]=double(int(TriCen[2]*1000.0-0.5))/1000.0;
      //        cout<<"\n\n\nTriCen[0]= "<<TriCen[0] << "   TriCen[1]= "<< TriCen[1] << "   TriCen[2]= "<< TriCen[2] <<endl;
 
      //fill the multi-index container with centriod coordinates and triangle handles
      Face_CenPos_Handle_varPos.insert(Face_CenPos_Handle(TriCen[0], TriCen[1], TriCen[2], *it));
    }
 
    //Find minimum and maximum X, Y and Z coordinates of the RVE (using multi-index container)
    double XcoordMin, YcoordMin, ZcoordMin, XcoordMax, YcoordMax, ZcoordMax;
    typedef Face_CenPos_Handle_multiIndex::index<xcoord_tag>::type::iterator Tri_Xcoord_iterator;
    typedef Face_CenPos_Handle_multiIndex::index<ycoord_tag>::type::iterator Tri_Ycoord_iterator;
    typedef Face_CenPos_Handle_multiIndex::index<zcoord_tag>::type::iterator Tri_Zcoord_iterator;
    Tri_Xcoord_iterator XcoordMin_it, XcoordMax_it;
    Tri_Ycoord_iterator YcoordMin_it, YcoordMax_it;
    Tri_Zcoord_iterator ZcoordMin_it, ZcoordMax_it;
 
    //XcoordMax_it-- because .end() will point iterator after the data range but .begin() will point the iteratore to the first value of range
    XcoordMin_it=Face_CenPos_Handle_varNeg.get<xcoord_tag>().begin();                  XcoordMin=XcoordMin_it->xcoord;
    XcoordMax_it=Face_CenPos_Handle_varPos.get<xcoord_tag>().end();    XcoordMax_it--; XcoordMax=XcoordMax_it->xcoord;
    YcoordMin_it=Face_CenPos_Handle_varNeg.get<ycoord_tag>().begin();                  YcoordMin=YcoordMin_it->ycoord;
    YcoordMax_it=Face_CenPos_Handle_varPos.get<ycoord_tag>().end();    YcoordMax_it--; YcoordMax=YcoordMax_it->ycoord;
    ZcoordMin_it=Face_CenPos_Handle_varNeg.get<zcoord_tag>().begin();                  ZcoordMin=ZcoordMin_it->zcoord;
    ZcoordMax_it=Face_CenPos_Handle_varPos.get<zcoord_tag>().end();    ZcoordMax_it--; ZcoordMax=ZcoordMax_it->zcoord;
 
//    cout<<"XcoordMin "<<XcoordMin << "      XcoordMax "<<XcoordMax <<endl;
//    cout<<"YcoordMin "<<YcoordMin << "      YcoordMax "<<YcoordMax <<endl;
//    cout<<"ZcoordMin "<<ZcoordMin << "      ZcoordMax "<<ZcoordMax <<endl;
    double LxRVE, LyRVE, LzRVE;
    LxRVE=XcoordMax-XcoordMin;
    LyRVE=YcoordMax-YcoordMin;
    LzRVE=ZcoordMax-ZcoordMin;
 
 
    //Creating Prisims between triangles on -ve and +ve faces
    typedef Face_CenPos_Handle_multiIndex::index<Tri_Hand_tag>::type::iterator Tri_Hand_iterator;
    Tri_Hand_iterator Tri_Neg;
    typedef Face_CenPos_Handle_multiIndex::index<Composite_xyzcoord>::type::iterator xyzcoord_iterator;
    xyzcoord_iterator Tri_Pos;
    Range PrismRange;
    double XPos, YPos, ZPos;
    //int count=0;
 
    //loop over -ve triangles to create prisims elemenet between +ve and -ve triangles
    for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end();  it++) {
 
      Tri_Neg=Face_CenPos_Handle_varNeg.get<Tri_Hand_tag>().find(*it);
      //cout<<"xcoord= "<<Tri_iit->xcoord << "   ycoord= "<< Tri_iit->ycoord << "   ycoord= "<< Tri_iit->zcoord <<endl;
 
      //corresponding +ve triangle
      if(Tri_Neg->xcoord==XcoordMin){XPos=XcoordMax;              YPos=Tri_Neg->ycoord;  ZPos=Tri_Neg->zcoord;};
      if(Tri_Neg->ycoord==YcoordMin){XPos=YPos=Tri_Neg->xcoord;   YPos=YcoordMax;        ZPos=Tri_Neg->zcoord;};
      if(Tri_Neg->zcoord==ZcoordMin){XPos=YPos=Tri_Neg->xcoord;   YPos=Tri_Neg->ycoord;  ZPos=ZcoordMax;      };
      Tri_Pos=Face_CenPos_Handle_varPos.get<Composite_xyzcoord>().find(boost::make_tuple(XPos, YPos, ZPos));
      //        cout<<"Tri_Neg->xcoord= "<<Tri_Neg->xcoord << "   Tri_Neg->ycoord "<< Tri_Neg->ycoord << "   Tri_Neg->zcoord= "<< Tri_Neg->zcoord <<endl;
      //        cout<<"Tri_Pos->xcoord= "<<Tri_Pos->xcoord << "   Tri_Pos->ycoord "<< Tri_Pos->ycoord << "   Tri_Pos->zcoord= "<< Tri_Pos->zcoord <<endl;
 
 
      //+ve and -ve nodes and their coords (+ve and -ve tiangles nodes can have matching problems, which can produce twisted prism)
      EntityHandle PrismNodes[6];
      vector<EntityHandle> TriNodesNeg, TriNodesPos;
      double CoordNodeNeg[9], CoordNodePos[9];
      rval = moab.get_connectivity(&(Tri_Neg->Tri_Hand),1,TriNodesNeg,true); CHKERRQ_MOAB(rval);
      rval = moab.get_connectivity(&(Tri_Pos->Tri_Hand),1,TriNodesPos,true); CHKERRQ_MOAB(rval);
      rval = moab.get_coords(&TriNodesNeg[0],3,CoordNodeNeg);  MOAB_THROW(rval);
      rval = moab.get_coords(&TriNodesPos[0],3,CoordNodePos);  MOAB_THROW(rval);
      for(int ii=0; ii<3; ii++){
        PrismNodes[ii]=TriNodesNeg[ii];
      }
      //        for(int ii=0; ii<3; ii++){
      //            cout<<"xcoord= "<<CoordNodeNeg[3*ii] << "   ycoord= "<< CoordNodeNeg[3*ii+1] << "   zcoord= "<< CoordNodeNeg[3*ii+2] <<endl;
      //        }
      //        for(int ii=0; ii<3; ii++){
      //            cout<<"xcoord= "<<CoordNodePos[3*ii] << "   ycoord= "<< CoordNodePos[3*ii+1] << "   zcoord= "<< CoordNodePos[3*ii+2] <<endl;
      //        }
 
      //Match exact nodes to each other to avoide the problem of twisted prisms
      double XNodeNeg, YNodeNeg, ZNodeNeg, XNodePos, YNodePos, ZNodePos;
      for(int ii=0; ii<3; ii++){
        if(Tri_Neg->xcoord==XcoordMin){XNodeNeg=XcoordMax;          YNodeNeg=CoordNodeNeg[3*ii+1];   ZNodeNeg=CoordNodeNeg[3*ii+2];};
        if(Tri_Neg->ycoord==YcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=YcoordMax;              ZNodeNeg=CoordNodeNeg[3*ii+2];};
        if(Tri_Neg->zcoord==ZcoordMin){XNodeNeg=CoordNodeNeg[3*ii]; YNodeNeg=CoordNodeNeg[3*ii+1];   ZNodeNeg=ZcoordMax;};
        for(int jj=0; jj<3; jj++){
          //Round nodal coordinates to 3 dicimal places only for comparison
          //round values to 3 disimal places
          if(XNodeNeg>=0) XNodeNeg=double(int(XNodeNeg*1000.0+0.5))/1000.0;   else XNodeNeg=double(int(XNodeNeg*1000.0-0.5))/1000.0;
          if(YNodeNeg>=0) YNodeNeg=double(int(YNodeNeg*1000.0+0.5))/1000.0;   else YNodeNeg=double(int(YNodeNeg*1000.0-0.5))/1000.0;
          if(ZNodeNeg>=0) ZNodeNeg=double(int(ZNodeNeg*1000.0+0.5))/1000.0;   else ZNodeNeg=double(int(ZNodeNeg*1000.0-0.5))/1000.0;
 
          XNodePos=CoordNodePos[3*jj]; YNodePos=CoordNodePos[3*jj+1]; ZNodePos=CoordNodePos[3*jj+2];
          if(XNodePos>=0) XNodePos=double(int(XNodePos*1000.0+0.5))/1000.0;   else XNodePos=double(int(XNodePos*1000.0-0.5))/1000.0;
          if(YNodePos>=0) YNodePos=double(int(YNodePos*1000.0+0.5))/1000.0;   else YNodePos=double(int(YNodePos*1000.0-0.5))/1000.0;
          if(ZNodePos>=0) ZNodePos=double(int(ZNodePos*1000.0+0.5))/1000.0;   else ZNodePos=double(int(ZNodePos*1000.0-0.5))/1000.0;
 
          if(XNodeNeg==XNodePos  &&  YNodeNeg==YNodePos  &&  ZNodeNeg==ZNodePos){
            PrismNodes[3+ii]=TriNodesPos[jj];
            break;
          }
        }
      }
      //prism nodes and their coordinates
      double CoordNodesPrisms[18];
      rval = moab.get_coords(PrismNodes,6,CoordNodesPrisms);  MOAB_THROW(rval);
      //        for(int ii=0; ii<6; ii++){
      //            cout<<"xcoord= "<<CoordNodesPrisms[3*ii] << "   ycoord= "<< CoordNodesPrisms[3*ii+1] << "   zcoord= "<< CoordNodesPrisms[3*ii+2] <<endl;
      //        }
      //        cout<<endl<<endl;
      //insertion of individula prism element and its addition to range PrismRange
      EntityHandle PeriodicPrism;
      rval = moab.create_element(MBPRISM,PrismNodes,6,PeriodicPrism); CHKERRQ_MOAB(rval);
      PrismRange.insert(PeriodicPrism);
 
      //Adding Prisims to Element LAGRANGE_ELEM (to loop over these prisims)
      EntityHandle PrismRangeMeshset;
      rval = moab.create_meshset(MESHSET_SET,PrismRangeMeshset); CHKERRQ_MOAB(rval);
      rval = moab.add_entities(PrismRangeMeshset,PrismRange); CHKERRQ_MOAB(rval);
      ierr = m_field.seed_ref_level_3D(PrismRangeMeshset,bit_level0); CHKERRQ(ierr);
      ierr = m_field.get_entities_by_ref_level(bit_level0,BitRefLevel().set(),meshset_level0); CHKERRQ(ierr);
 
      //        //to see individual prisms
      //        Range Prism1;
      //        Prism1.insert(PeriodicPrism);
      //        EntityHandle out_meshset1;
      //        rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
      //        rval = moab.add_entities(out_meshset1,Prism1); CHKERRQ_MOAB(rval);
      //        ostringstream sss;
      //        sss << "Prism" << count << ".vtk"; count++;
      //        rval = moab.write_file(sss.str().c_str(),"VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
 
    }
 
 
    //cout<<"PrismRange "<<PrismRange<<endl;
    //    //Saving prisms in interface.vtk
    //      EntityHandle out_meshset1;
    //      rval = moab.create_meshset(MESHSET_SET,out_meshset1); CHKERRQ_MOAB(rval);
    //      rval = moab.add_entities(out_meshset1,PrismRange); CHKERRQ_MOAB(rval);
    //      rval = moab.write_file("Prisms.vtk","VTK","",&out_meshset1,1); CHKERRQ_MOAB(rval);
 
 
//    //=======================================================================================================
//    //=======================================================================================================
 
 
 
  SmallStrainParaboloidalPlasticity cp;
  {
    cp.tAgs.resize(3);
    cp.tAgs[0] = 0;
    cp.tAgs[1] = 1;
    cp.tAgs[2] = 2;
    cp.tOl = 1e-12;
 
    double young_modulus = 1;
    double poisson_ratio = 0.25;
    cp.sIgma_yt = 1;
    cp.sIgma_yc = 1;
 
    cp.Ht = 0.1;
    cp.Hc = 0.1;
 
    cp.nup = 0.3;
 
 
    {
      ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
      ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
      cp.mu = MU(young_modulus,poisson_ratio);
      cp.lambda = LAMBDA(young_modulus,poisson_ratio);
      ierr = PetscOptionsGetReal(0,"-my_sigma_yt",&cp.sIgma_yt,0); CHKERRQ(ierr);
      ierr = PetscOptionsGetReal(0,"-my_sigma_yc",&cp.sIgma_yc,0); CHKERRQ(ierr);
 
      ierr = PetscOptionsGetReal(0,"-my_Ht",&cp.Ht,0); CHKERRQ(ierr);
      ierr = PetscOptionsGetReal(0,"-my_Hc",&cp.Hc,0); CHKERRQ(ierr);
 
      ierr = PetscOptionsGetReal(0,"-my_nt",&cp.nt,0); CHKERRQ(ierr);
      ierr = PetscOptionsGetReal(0,"-my_nc",&cp.nc,0); CHKERRQ(ierr);
 
 
      ierr = PetscOptionsGetReal(0,"-my_nup",&cp.nup,0); CHKERRQ(ierr);
    }
 
    PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
    PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
 
    PetscPrintf(PETSC_COMM_WORLD,"sIgma_yt = %4.2e \n",cp.sIgma_yt);
    PetscPrintf(PETSC_COMM_WORLD,"sIgma_yc = %4.2e \n",cp.sIgma_yt);
 
    PetscPrintf(PETSC_COMM_WORLD,"nup = %4.2e \n",cp.nup);
 
    cp.sTrain.resize(6,false);
    cp.sTrain.clear();
    cp.plasticStrain.resize(6,false);
    cp.plasticStrain.clear();
    cp.internalVariables.resize(2,false);
    cp.internalVariables.clear();
    cp.createMatAVecR();
    cp.snesCreate();
    // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
  }
 
 
 
//  SmallStrainJ2Plasticity cp;
//  {
//    cp.tAgs.resize(3);
//    cp.tAgs[0] = 0;
//    cp.tAgs[1] = 1;
//    cp.tAgs[2] = 2;
//    cp.tOl = 1e-12;
//
//    double young_modulus = 1;
//    double poisson_ratio = 0.25;
//    cp.sIgma_y = 1;
//    cp.H = 0.1;
//    cp.K = 0;
//    cp.phi = 1;
//    {
//      ierr = PetscOptionsGetReal(0,"-my_young_modulus",&young_modulus,0); CHKERRQ(ierr);
//      ierr = PetscOptionsGetReal(0,"-my_poisson_ratio",&poisson_ratio,0); CHKERRQ(ierr);
//      cp.mu = MU(young_modulus,poisson_ratio);
//      cp.lambda = LAMBDA(young_modulus,poisson_ratio);
//      ierr = PetscOptionsGetReal(0,"-my_sigma_y",&cp.sIgma_y,0); CHKERRQ(ierr);
//      ierr = PetscOptionsGetReal(0,"-my_H",&cp.H,0); CHKERRQ(ierr);
//      ierr = PetscOptionsGetReal(0,"-my_K",&cp.K,0); CHKERRQ(ierr);
//      ierr = PetscOptionsGetReal(0,"-my_phi",&cp.phi,0); CHKERRQ(ierr);
//    }
//
//    PetscPrintf(PETSC_COMM_WORLD,"young_modulus = %4.2e \n",young_modulus);
//    PetscPrintf(PETSC_COMM_WORLD,"poisson_ratio = %4.2e \n",poisson_ratio);
//    PetscPrintf(PETSC_COMM_WORLD,"sIgma_y = %4.2e \n",cp.sIgma_y);
//    PetscPrintf(PETSC_COMM_WORLD,"H = %4.2e \n",cp.H);
//    PetscPrintf(PETSC_COMM_WORLD,"K = %4.2e \n",cp.K);
//    PetscPrintf(PETSC_COMM_WORLD,"phi = %4.2e \n",cp.phi);
//
//
//    cp.sTrain.resize(6,false);
//    cp.sTrain.clear();
//    cp.plasticStrain.resize(6,false);
//    cp.plasticStrain.clear();
//    cp.internalVariables.resize(7,false);
//    cp.internalVariables.clear();
//    cp.createMatAVecR();
//    cp.snesCreate();
//    // ierr = SNESSetFromOptions(cp.sNes); CHKERRQ(ierr);
//
//  }
 
 
//    BitRefLevel problem_bit_level = bit_levels.back();
//    Range CubitSideSets_meshsets;
//    ierr = m_field.get_cubit_meshsets(SIDESET,CubitSideSets_meshsets); CHKERRQ(ierr);
 
    //Fields
    ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
    ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
    ierr = m_field.add_field("LAGRANGE_MUL_PERIODIC",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);  //For lagrange multipliers to control the periodic motion
    ierr = m_field.add_field("LAGRANGE_MUL_RIGID_TRANS",NOFIELD,NOBASE,3); CHKERRQ(ierr);  //To control the rigid body motion (3 Traslations and 3
 
 
    //add entitities (by tets) to the field
    ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
    ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
 
 
    Range surface_negative;
    ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,surface_negative,true); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",surface_negative.size()); CHKERRQ(ierr);
    ierr = m_field.add_ents_to_field_by_type(surface_negative,MBTRI,"LAGRANGE_MUL_PERIODIC"); CHKERRQ(ierr);
 
 
    //set app. order
    ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",bubble_order); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
 
    ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
 
    ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_PERIODIC",order); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_PERIODIC",1); CHKERRQ(ierr);
 
 
    //build field
    ierr = m_field.build_fields(); CHKERRQ(ierr);
    //10 node tets
    Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
    ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material,0); CHKERRQ(ierr);
 
    //FE
    ierr = m_field.add_finite_element("PLASTIC"); CHKERRQ(ierr);
    //Define rows/cols and element data
    ierr = m_field.modify_finite_element_add_field_row("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_col("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("PLASTIC","DISPLACEMENT"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("PLASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
    ierr = m_field.add_ents_to_finite_element_by_type(0,MBTET,"PLASTIC"); CHKERRQ(ierr);
 
    BCs_RVELagrange_Periodic lagrangian_element_periodic(m_field);
    BCs_RVELagrange_Trac lagrangian_element_trac(m_field);
 
    lagrangian_element_periodic.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",PrismRange);
    lagrangian_element_trac.addLagrangiangElement("LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS");
 
    //build finite elements
    ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
    //build adjacencies
    ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
 
 
    DMType dm_name = "PLASTIC_PROB";
    ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
 
    DM dm;
    ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
    ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
 
    //set dm datastruture which created mofem datastructures
    ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
    ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
    ierr = DMMoFEMSetIsPartitioned(dm,is_partitioned); CHKERRQ(ierr);
    //add elements to dm
    ierr = DMMoFEMAddElement(dm,"PLASTIC"); CHKERRQ(ierr);
    ierr = DMMoFEMAddElement(dm,"LAGRANGE_ELEM"); CHKERRQ(ierr);
    ierr = DMMoFEMAddElement(dm,"LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
    ierr = DMSetUp(dm); CHKERRQ(ierr);
 
    //create matrices
    Vec F,D;
    Mat Aij;
    ierr = DMCreateGlobalVector_MoFEM(dm,&D); CHKERRQ(ierr);
    ierr = VecDuplicate(D,&F); CHKERRQ(ierr);
    ierr = DMCreateMatrix_MoFEM(dm,&Aij); CHKERRQ(ierr);
 
    ierr = VecZeroEntries(D); CHKERRQ(ierr);
    ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = DMoFEMMeshToLocalVector(dm,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
 
    ierr = VecZeroEntries(F); CHKERRQ(ierr);
    ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = DMoFEMMeshToLocalVector(dm,F,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
 
    ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
 
    vector<Vec> Fvec(6); //jthis will be used to caluclate the homogenised stiffness matix
    for(int ii = 0;ii<6;ii++) {
      ierr = VecDuplicate(D,&Fvec[ii]); CHKERRQ(ierr);
      ierr = VecZeroEntries(Fvec[ii]); CHKERRQ(ierr);
      ierr = VecGhostUpdateBegin(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
      ierr = VecGhostUpdateEnd(Fvec[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    }
 
    //assemble Aij and F
    SmallStrainPlasticity small_strain_plasticity(m_field);
    {
      PetscBool bbar = PETSC_TRUE;
      ierr = PetscOptionsGetBool(0,"-my_bbar",&bbar,0); CHKERRQ(ierr);
      small_strain_plasticity.commonData.bBar = bbar;
    }
    {
      small_strain_plasticity.feRhs.getOpPtrVector().push_back(
        new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
          "DISPLACEMENT",small_strain_plasticity.commonData
        )
      );
      small_strain_plasticity.feRhs.getOpPtrVector().push_back(
        new SmallStrainPlasticity::OpCalculateStress(
          m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
        )
      );
      small_strain_plasticity.feRhs.getOpPtrVector().push_back(
        new SmallStrainPlasticity::OpAssembleRhs(
          "DISPLACEMENT",small_strain_plasticity.commonData
        )
      );
      small_strain_plasticity.feLhs.getOpPtrVector().push_back(
        new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
          "DISPLACEMENT",small_strain_plasticity.commonData
        )
      );
      small_strain_plasticity.feLhs.getOpPtrVector().push_back(
        new SmallStrainPlasticity::OpCalculateStress(
          m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
        )
      );
      small_strain_plasticity.feLhs.getOpPtrVector().push_back(
        new SmallStrainPlasticity::OpAssembleLhs(
          "DISPLACEMENT",small_strain_plasticity.commonData
        )
      );
      small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
        new SmallStrainPlasticity::OpGetCommonDataAtGaussPts(
          "DISPLACEMENT",small_strain_plasticity.commonData
        )
      );
      small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
        new SmallStrainPlasticity::OpCalculateStress(
          m_field,"DISPLACEMENT",small_strain_plasticity.commonData,cp,false
        )
      );
      small_strain_plasticity.feUpdate.getOpPtrVector().push_back(
        new SmallStrainPlasticity::OpUpdate(
          "DISPLACEMENT",small_strain_plasticity.commonData
        )
      );
      ierr = small_strain_plasticity.createInternalVariablesTag(); CHKERRQ(ierr);
    }
 
    lagrangian_element_periodic.setRVEBCsOperatorsNonlinear("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",Aij,Fvec,F,given_strain);
 
    BCs_RVELagrange_Trac_Rigid_Trans lagrangian_element_rigid_body_trans(m_field);
    lagrangian_element_rigid_body_trans.setRVEBCsRigidBodyTranOperators("DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS",Aij, lagrangian_element_periodic.setOfRVEBC);
 
 
    TimeForceScale time_force_scale("-my_macro_strian_history",false);
    lagrangian_element_periodic.methodsOp.push_back(new TimeForceScale("-my_macro_strian_history",false));
 
 
 
 
    //Adding elements to DMSnes
    //Rhs
    ierr = DMMoFEMSNESSetFunction(dm,"PLASTIC",&small_strain_plasticity.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
    ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
    ierr = DMMoFEMSNESSetFunction(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCRhsResidual,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
 
    //Lhs
    ierr = DMMoFEMSNESSetJacobian(dm,"PLASTIC",&small_strain_plasticity.feLhs,NULL,NULL); CHKERRQ(ierr);
    ierr = DMMoFEMSNESSetJacobian(dm,"LAGRANGE_ELEM",&lagrangian_element_periodic.feRVEBCLhs,NULL,NULL); CHKERRQ(ierr);
    ierr = DMMoFEMSNESSetJacobian(dm,"LAGRANGE_ELEM_TRANS",&lagrangian_element_rigid_body_trans.feRVEBCLhs,NULL,NULL); CHKERRQ(ierr);
 
 
    // Create Time Stepping solver
    SNES snes;
    SnesCtx *snes_ctx;
    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    //ierr = SNESSetDM(snes,dm); CHKERRQ(ierr);
    ierr = DMMoFEMGetSnesCtx(dm,&snes_ctx); CHKERRQ(ierr);
    ierr = SNESSetFunction(snes,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
    ierr = SNESSetJacobian(snes,Aij,Aij,SnesMat,snes_ctx); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
 
    PostProcVolumeOnRefinedMesh post_proc(m_field);
    ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
    ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
    ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
    ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
 
       // Volume calculation
    //==============================================================================================================================
    VolumeElementForcesAndSourcesCore vol_elem(m_field);
    Vec volume_vec;
    int volume_vec_ghost[] = { 0 };
    ierr = VecCreateGhost(
            PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
            );  CHKERRQ(ierr);
    ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
    vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
 
    ierr = m_field.loop_finite_elements("PLASTIC_PROB","PLASTIC", vol_elem);  CHKERRQ(ierr);
 
    ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
    double rve_volume;
    ierr = VecSum(volume_vec,&rve_volume);  CHKERRQ(ierr);
 
    //            ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
    ierr = VecDestroy(&volume_vec);
    //=============================================================================================================================
 
    double final_time = 1,delta_time = 0.1;
    ierr = PetscOptionsGetReal(0,"-my_final_time",&final_time,0); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(0,"-my_delta_time",&delta_time,0); CHKERRQ(ierr);
    double delta_time0 = delta_time;
 
    // ierr = MatView(Aij,PETSC_VIEWER_DRAW_SELF); CHKERRQ(ierr);
    // std::string wait;
    // std::cin >> wait;
 
    Vec D0;
    ierr = VecDuplicate(D,&D0); CHKERRQ(ierr);
 
    int step = 0;
    double t = 0;
    SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
    for(;t<final_time;step++) {
      t += delta_time;
 
      if(t>final_time){
      break;
      }
      PetscPrintf(
        PETSC_COMM_WORLD,"Step %d Time %6.4g final time %3.2g\n",step,t,final_time
      );
      //set time
      lagrangian_element_periodic.getLoopFeRVEBCRhs().ts_t = t;
      //solve problem at step
      ierr = VecAssemblyBegin(D);
      ierr = VecAssemblyEnd(D);
      ierr = VecCopy(D,D0); CHKERRQ(ierr);
 
      if(step == 0 || reason < 0) {
        ierr = SNESSetLagJacobian(snes,-2); CHKERRQ(ierr);
      } else {
        ierr = SNESSetLagJacobian(snes,-1); CHKERRQ(ierr);
      }
 
 
      ierr = SNESSolve(snes,PETSC_NULL,D); CHKERRQ(ierr);
      int its;
      ierr = SNESGetIterationNumber(snes,&its); CHKERRQ(ierr);
 
      ierr = PetscPrintf(
        PETSC_COMM_WORLD,"number of Newton iterations = %D\n",its
      ); CHKERRQ(ierr);
 
      ierr = SNESGetConvergedReason(snes,&reason); CHKERRQ(ierr);
//      cout<<"reason ====== "<<reason<<endl;
 
          if(reason<0) {
            t -= delta_time;
            delta_time *= 0.1;
            ierr = VecCopy(D0,D); CHKERRQ(ierr);
          } else {const int its_d = 20;
            const double gamma = 0.5;
            const double max_reudction = 1;
            const double min_reduction = 1e-1;
            double reduction;
            reduction = pow((double)its_d/(double)(its+1),gamma);
            if(delta_time >= max_reudction*delta_time0 && reduction > 1) {
              reduction = 1;
            } else if(delta_time <= min_reduction*delta_time0 && reduction < 1) {
              reduction = 1;
            }
 
            ierr = PetscPrintf(
              PETSC_COMM_WORLD,
              "reduction delta_time = %6.4e delta_time = %6.4e\n",
              reduction,delta_time
            ); CHKERRQ(ierr);
            delta_time *= reduction;
            if(reduction>1 && delta_time < min_reduction*delta_time0) {
              delta_time = min_reduction*delta_time0;
            }
 
           ierr = DMoFEMMeshToGlobalVector(
              dm,D,INSERT_VALUES,SCATTER_REVERSE
            ); CHKERRQ(ierr);
            ierr = DMoFEMLoopFiniteElements(
              dm,"PLASTIC",&small_strain_plasticity.feUpdate
            ); CHKERRQ(ierr);
 
 
           //Save data on mesh
            ierr = DMoFEMLoopFiniteElements(
              dm,"PLASTIC",&post_proc
            ); CHKERRQ(ierr);
            string out_file_name;
            {
              std::ostringstream stm;
              stm << "out_" << step << ".h5m";
              out_file_name = stm.str();
            }
            ierr = PetscPrintf(
              PETSC_COMM_WORLD,"out file %s\n",out_file_name.c_str()
            ); CHKERRQ(ierr);
            rval = post_proc.postProcMesh.write_file(
              out_file_name.c_str(),"MOAB","PARALLEL=WRITE_PART"
            ); CHKERRQ_MOAB(rval);
 
            //===================================== Homgenised stress (for given strain) ====================================================
            VectorDouble StressHomo;
            StressHomo.resize(6);
            StressHomo.clear();
 
            //calculate honmogenised stress
            //create a vector for 6 components of homogenized stress
            Vec stress_homo;
            int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
            ierr = VecCreateGhost(
                                  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
                                  );  CHKERRQ(ierr);
 
            lagrangian_element_periodic.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",stress_homo);
 
            PetscScalar *avec;
            ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
            ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
            ierr = m_field.loop_finite_elements(
                    "PLASTIC_PROB","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCStress()
                    ); CHKERRQ(ierr);
            ierr = PetscOptionsGetReal(
                   PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
                   ); CHKERRQ(ierr);
            ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
            ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
            ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
            ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
            ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
 
            for(int jj=0; jj<6; jj++) {
              StressHomo(jj) = avec[jj];
            }
 
            double scale;
            ierr = time_force_scale.getForceScale(t,scale); CHKERRQ(ierr);
 
            PetscPrintf(PETSC_COMM_WORLD,
              "Macro_Strain Homo_Stress Path "
              );
 
            //cout<<"Macro-strain ";
            for(int ii=0; ii<6; ii++) {
              PetscPrintf(
              PETSC_COMM_WORLD,
              "%8.5e\t",
              t*given_strain(ii)
              );
            }
 
            //cout<<"Homo stress ";
            for(int ii=0; ii<6; ii++) {
              PetscPrintf(
              PETSC_COMM_WORLD,
              "%8.5e\t",
              StressHomo(ii)
              );
            }
 
            PetscPrintf(PETSC_COMM_WORLD,
              "\n"
            );
 
            //==============================================================================================================================
 
        }
      }
 
  ierr = VecDestroy(&D0); CHKERRQ(ierr);
  ierr = MatDestroy(&Aij); CHKERRQ(ierr);
  ierr = VecDestroy(&F); CHKERRQ(ierr);
  ierr = VecDestroy(&D); CHKERRQ(ierr);
 
    }
    CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
}