rve_mech_plas_interface_periodic.cpp File Reference

#include <MoFEM.hpp>
#include <boost/program_options.hpp>
#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>
#include <CohesiveInterfaceElement.hpp>
#include <adolc/adolc.h>
#include <NonLinearElasticElement.hpp>
#include <SmallTransverselyIsotropic.hpp>
#include <SmallStrainPlasticity.hpp>
#include <SmallStrainPlasticityMaterialModels.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <PostProcOnRefMesh.hpp>
#include <string>

Go to the source code of this file.

Classes
struct	Face_CenPos_Handle
struct	xcoord_tag
struct	ycoord_tag
struct	zcoord_tag
struct	Tri_Hand_tag
struct	Composite_xyzcoord

Typedefs
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

Functions
int	main (int argc, char *argv[])

Variables
static char	help []

Typedef Documentation

Face_CenPos_Handle_multiIndex

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

Definition at line 115 of file rve_mech_plas_interface_periodic.cpp.

Function Documentation

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 122 of file rve_mech_plas_interface_periodic.cpp.

                                 {
  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);
  }
 
 
 
  double young_modulus = 1;
  SmallStrainParaboloidalPlasticity cp;
  {
    cp.tAgs.resize(3);
    cp.tAgs[0] = 3;
    cp.tAgs[1] = 4;
    cp.tAgs[2] = 5;
    cp.tOl = 1e-12;
 
    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_yc);
 
    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);
  }
 
 
 // double young_modulus = 1;
 // SmallStrainJ2Plasticity cp;
 // {
 //   cp.tAgs.resize(3);
 //   cp.tAgs[0] = 3;
 //   cp.tAgs[1] = 4;
 //   cp.tAgs[2] = 5;
 //   cp.tOl = 1e-12;
 //
 //   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);
 // }
 
 
  //Create MoFEM (Joseph) database
  MoFEM::Core core(moab);
  MoFEM::Interface& m_field = core;
 
  vector<BitRefLevel> bit_levels;
  Tag th_meshset_info;
  int def_meshset_info[2] = {0,0};
  rval = moab.tag_get_handle("MESHSET_INFO",2,MB_TYPE_INTEGER,th_meshset_info,MB_TAG_CREAT|MB_TAG_SPARSE,&def_meshset_info);
  int meshset_data[2];
  EntityHandle root = moab.get_root_set();
  rval = moab.tag_get_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
  if(meshset_data[0]==0) {
    meshset_data[0] = 1;
    rval = moab.tag_set_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
  }
  bit_levels.push_back(BitRefLevel().set(meshset_data[0]-1));
 
 
//  EntityHandle out_meshset;
//  rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
//  //    ierr = m_field.get_problem_finite_elements_entities("POTENTIAL_PROBLEM","POTENTIAL_ELEM",out_meshset); CHKERRQ(ierr);
//  ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
//  Range LatestRefinedTets;
//  rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
//  Range LatestRefinedPrisms;
//  rval = moab.get_entities_by_type(out_meshset, MBPRISM,LatestRefinedPrisms,true); CHKERRQ_MOAB(rval);
 
 
  BitRefLevel problem_bit_level = bit_levels.back();
 
  //    const clock_t begin_time = clock();
  ierr = m_field.build_fields(); CHKERRQ(ierr);
  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
 
 
  Range interface_prims_on_problem_bit_level;
  ierr = m_field.get_entities_by_type_and_ref_level(
          problem_bit_level,BitRefLevel().set(),MBPRISM,interface_prims_on_problem_bit_level
          ); CHKERRQ(ierr);
  Range tets_on_problem_bit_level;
  ierr = m_field.get_entities_by_type_and_ref_level(
          problem_bit_level,BitRefLevel().set(),MBTET,tets_on_problem_bit_level
          ); CHKERRQ(ierr);
 
  cout<<"interface_prims_on_problem_bit_level   "<<interface_prims_on_problem_bit_level.size()<<endl;
  cout<<"tets_on_problem_bit_level    "<<tets_on_problem_bit_level.size()<<endl;
 
  //=======================================================================================================
  //Add Periodic Prisims Between Triangles on -ve and +ve faces to implement periodic bounary conditions
  //=======================================================================================================
  Range preriodic_prisms;
  //Populating the Multi-index container with -ve triangles
 
 
  Range SurTrisNeg; //Old + new traiangles on side set 101
  ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Old+New triangles on SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
 
  Range AllTris; //All new traiangles in the new-bitlevel
  ierr = m_field.get_entities_by_type_and_ref_level(problem_bit_level,BitRefLevel().set(),MBTRI,AllTris); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"All triangles in the new bit-level = %d\n",AllTris.size()); CHKERRQ(ierr);
 
  //Only triangles in the new bit-level
  SurTrisNeg = intersect(AllTris,SurTrisNeg);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"New triangles on SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
 
//  cout<<"SurTrisNeg   "<<SurTrisNeg<<endl;
  Face_CenPos_Handle_multiIndex Face_CenPos_Handle_varNeg, Face_CenPos_Handle_varPos;
  double TriCen[3], coords_Tri[9];
 
 
 
  double roundfact=1e3;
 
  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;
 
//    cout<<"\nTriCen[0]= "<<TriCen[0] << "   TriCen[1]= "<< TriCen[1] << "   TriCen[2]= "<< TriCen[2] <<endl;
    TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
    TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
    TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
 
//    cout<<"TriCen[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*roundfact+0.5))/roundfact<<endl<<endl<<endl<<endl;
  //    aaa=-0.5011;
  //    cout<<"\n\n\n\n\nfloor(aaa+0.5) = "<<double(int(aaa*roundfact-0.5))/roundfact<<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,"Old+New triangles on SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
 
  //Only triangles in the new bit-level
  SurTrisPos = intersect(AllTris,SurTrisPos);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"New triangles on 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
    TriCen[0]=round(TriCen[0]*roundfact)/roundfact;
    TriCen[1]=round(TriCen[1]*roundfact)/roundfact;
    TriCen[2]=round(TriCen[2]*roundfact)/roundfact;
 
//    cout<<"TriCen[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;
  double XPos, YPos, ZPos;
  //int count=0;
 
  //loop over -ve triangles to create prisims elemenet between +ve and -ve triangles
  int count=1;
  int count1=0;
  for(Range::iterator it = SurTrisNeg.begin(); it!=SurTrisNeg.end();  it++) {
 
//    count1++;
//    cout<<"count1 = "<<count1<< endl;
//    if(count1 == 24) {
 
 
    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
        XNodeNeg=round(XNodeNeg*roundfact)/roundfact;
        YNodeNeg=round(YNodeNeg*roundfact)/roundfact;
        ZNodeNeg=round(ZNodeNeg*roundfact)/roundfact;
 
        XNodePos=CoordNodePos[3*jj]; YNodePos=CoordNodePos[3*jj+1]; ZNodePos=CoordNodePos[3*jj+2];
        XNodePos=round(XNodePos*roundfact)/roundfact;
        YNodePos=round(YNodePos*roundfact)/roundfact;
        ZNodePos=round(ZNodePos*roundfact)/roundfact;
 
        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 preriodic_prisms
    EntityHandle PeriodicPrism;
    rval = moab.create_element(MBPRISM,PrismNodes,6,PeriodicPrism); CHKERRQ_MOAB(rval);
    preriodic_prisms.insert(PeriodicPrism);
 
//    //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);
//
//    string aaa;
//    cin>>aaa;
//
//    }
 
  }
 
 
 
 
  //Adding the newely added prisms to the bit-level (which is known as seeding)
  /*{
    Range ents;
    moab.get_adjacencies(preriodic_prisms,1,true,ents,moab::Interface::UNION);
    moab.get_adjacencies(preriodic_prisms,2,true,ents,moab::Interface::UNION);
    ierr = m_field.seed_ref_level(ents.subset_by_type(MBEDGE),problem_bit_level); CHKERRQ(ierr);
    ierr = m_field.seed_ref_level(ents.subset_by_type(MBQUAD),problem_bit_level); CHKERRQ(ierr);
  }*/
  ierr = m_field.seed_ref_level(preriodic_prisms,problem_bit_level); CHKERRQ(ierr);
  ierr = m_field.seed_finite_elements(preriodic_prisms); CHKERRQ(ierr);
 
 
  //=======================================================================================================
  //=======================================================================================================
  //=======================================================================================================
 
  EntityHandle out_meshset;
  rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
  ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
 
  Range LatestRefinedTets;
  rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
 
  Range Interface_Periodic_Prisms;
  rval = moab.get_entities_by_type(out_meshset, MBPRISM,Interface_Periodic_Prisms,true); CHKERRQ_MOAB(rval);
  cout<<"Interface_Periodic_Prisms    " <<Interface_Periodic_Prisms.size()<<endl;
 
 
  //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,MB_TAG_SPARSE,MF_ZERO); 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);
 
//  ierr = PetscPrintf(PETSC_COMM_WORLD,"number of triangle in SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
  ierr = m_field.add_ents_to_field_by_type(SurTrisNeg,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);
 
    Range newtets;
    for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|UNKNOWNNAME,it)) {
      if(it->getName() != "MAT_PLASTIC") continue;
      EntityHandle meshset = it->getMeshset();
      int id = it->getMeshsetId();
 
      rval = m_field.get_moab().get_entities_by_type(
                meshset,MBTET,newtets,true
             ); CHKERRQ_MOAB(rval);
 
      cout<< "========================== newtets "<<newtets.size()<<endl;
 
      //intersection of new and old tets for plastic
      newtets = intersect(newtets,LatestRefinedTets);
      ierr = m_field.seed_finite_elements(newtets); CHKERRQ(ierr);
    }
 
    ierr = m_field.add_ents_to_finite_element_by_type(newtets,MBTET,"PLASTIC"); CHKERRQ(ierr);
    cout<< "========================== newtets "<<newtets.size()<<endl;
 
 
  //================================================================================================================================
  //          Trans-Isotropic Yarns
  //================================================================================================================================
 
  NonlinearElasticElement trans_elastic(m_field,2);
  trans_elastic.commonData.spatialPositions = "DISPLACEMENT";
  trans_elastic.commonData.meshPositions = "MESH_NODE_POSITIONS";
  std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicADouble> > tranversly_isotropic_adouble_ptr_map;
  std::map<int,boost::shared_ptr<SmallStrainTranverslyIsotropicDouble> > tranversly_isotropic_double_ptr_map;
  bool trans_iso_blocks = false;
  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
    //Get block name
    string name = it->getName();
    if (name.compare(0,20,"MAT_ELASTIC_TRANSISO") == 0) {
      cout<<"================================ it is MAT_ELASTIC_TRANSISO "<<endl;
      trans_iso_blocks = true;
      int id = it->getMeshsetId();
      Mat_Elastic_TransIso mydata;
      ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
      tranversly_isotropic_adouble_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicADouble>();
      tranversly_isotropic_double_ptr_map[id] = boost::make_shared<SmallStrainTranverslyIsotropicDouble>();
      //nu_p, nu_pz, E_p, E_z, G_zp
      tranversly_isotropic_adouble_ptr_map.at(id)->E_p = mydata.data.Youngp;
      tranversly_isotropic_double_ptr_map.at(id)->E_p = mydata.data.Youngp;
      tranversly_isotropic_adouble_ptr_map.at(id)->E_z = mydata.data.Youngz;
      tranversly_isotropic_double_ptr_map.at(id)->E_z = mydata.data.Youngz;
      tranversly_isotropic_adouble_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
      tranversly_isotropic_double_ptr_map.at(id)->nu_p = mydata.data.Poissonp;
      tranversly_isotropic_adouble_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
      tranversly_isotropic_double_ptr_map.at(id)->nu_pz = mydata.data.Poissonpz;
      double shear_zp;
      if(mydata.data.Shearzp!=0) {
        shear_zp = mydata.data.Shearzp;
      } else {
        shear_zp = mydata.data.Youngz/(2*(1+mydata.data.Poissonpz));
      }
      tranversly_isotropic_adouble_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
      tranversly_isotropic_double_ptr_map.at(it->getMeshsetId())->G_zp = shear_zp;
      //get tets from block where material is defined
      EntityHandle meshset = it->getMeshset();
      rval = m_field.get_moab().get_entities_by_type(
        meshset,MBTET,trans_elastic.setOfBlocks[id].tEts,true
      ); CHKERRQ_MOAB(rval);
 
//      cout<<"================================ trans_elastic.setOfBlocks[id].tEts "<<trans_elastic.setOfBlocks[id].tEts.size()<<endl;
      //intersection of new and old tets for tran-iso part
      trans_elastic.setOfBlocks[id].tEts = intersect(trans_elastic.setOfBlocks[id].tEts,LatestRefinedTets);
//      cout<<"================================ trans_elastic.setOfBlocks[id].tEts "<<trans_elastic.setOfBlocks[id].tEts.size()<<endl;
 
      //adding material to nonlinear class
      trans_elastic.setOfBlocks[id].iD = id;
      //note that material parameters are defined internally in material model
      trans_elastic.setOfBlocks[id].E = 0; // this is not working for this material
      trans_elastic.setOfBlocks[id].PoissonRatio = 0; // this is not working for this material
      trans_elastic.setOfBlocks[id].materialDoublePtr = tranversly_isotropic_double_ptr_map.at(id);
      trans_elastic.setOfBlocks[id].materialAdoublePtr = tranversly_isotropic_adouble_ptr_map.at(id);
      ierr = m_field.seed_finite_elements(trans_elastic.setOfBlocks[id].tEts); CHKERRQ(ierr);
    }
  }
 
 
    ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
    ierr = m_field.add_finite_element("TRAN_ISOTROPIC_ELASTIC",MF_ZERO); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_row("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_col("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","POTENTIAL_FIELD"); CHKERRQ(ierr);
 
    if(m_field.check_field("MESH_NODE_POSITIONS")) {
      ierr = m_field.modify_finite_element_add_field_data("TRAN_ISOTROPIC_ELASTIC","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
    }
    for(
        map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
        sit!=trans_elastic.setOfBlocks.end();sit++
        ) {
 
 
      cout<<" sit->second.tEts "<<sit->second.tEts.size()<<endl;
      ierr = m_field.add_ents_to_finite_element_by_type(sit->second.tEts,MBTET,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
    }
 
 
      //Rhs
    trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData));
    trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData));
    if(m_field.check_field("MESH_NODE_POSITIONS")) {
      trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData));
    }
    map<int,NonlinearElasticElement::BlockData>::iterator sit = trans_elastic.setOfBlocks.begin();
    for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
      trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress("DISPLACEMENT",sit->second,trans_elastic.commonData,2,false,false,true));
      trans_elastic.feRhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpRhsPiolaKirchhoff("DISPLACEMENT",sit->second,trans_elastic.commonData));
    }
 
    //Lhs
    trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("DISPLACEMENT",trans_elastic.commonData));
    trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("POTENTIAL_FIELD",trans_elastic.commonData));
    if(m_field.check_field("MESH_NODE_POSITIONS")) {
      trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpGetCommonDataAtGaussPts("MESH_NODE_POSITIONS",trans_elastic.commonData));
    }
    sit = trans_elastic.setOfBlocks.begin();
    for(;sit!=trans_elastic.setOfBlocks.end();sit++) {
      trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress("DISPLACEMENT",sit->second,trans_elastic.commonData,2,true,false,true));
      trans_elastic.feLhs.getOpPtrVector().push_back(new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx("DISPLACEMENT","DISPLACEMENT",sit->second,trans_elastic.commonData));
    }
 
 
 
 
  //================================================================================================================================
  //          Interface Cohesive elemetns
  //================================================================================================================================
  //FE Interface
 
    ierr = m_field.add_finite_element("INTERFACE"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_row("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_col("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("INTERFACE","DISPLACEMENT"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("INTERFACE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
 
    boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation> interface_materials;
 
    //FIXME this in fact allow only for one type of interface,
    //problem is Young Moduls in interface mayoung_modulusterial
 
 
    double interface_beta, interface_ft, interface_Gf, interface_h;
    ierr = PetscOptionsGetReal(0,"-interface_beta",&interface_beta,0); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(0,"-interface_ft",&interface_ft,0); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(0,"-interface_Gf",&interface_Gf,0); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(0,"-interface_h",&interface_h,0); CHKERRQ(ierr);
 
 
    for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
      cout << endl << *it << endl;
      //Get block name
      string name = it->getName();
      if (name.compare(0,10,"MAT_INTERF") == 0) {
        Mat_Interf mydata;
        ierr = it->getAttributeDataStructure(mydata); CHKERRQ(ierr);
        cout << mydata;
        interface_materials.push_back(new CohesiveInterfaceElement::PhysicalEquation(m_field));
//        interface_materials.back().h = 0.134;  // we used this because E0=35Gpa and Em=4.7Gpa
        interface_materials.back().youngModulus = mydata.data.alpha;
//        interface_materials.back().beta = mydata.data.beta;
//        interface_materials.back().ft = mydata.data.ft;
//        interface_materials.back().Gf = mydata.data.Gf;
 
        interface_materials.back().h    = interface_h;
        interface_materials.back().beta = interface_beta;
        interface_materials.back().ft   = interface_ft;
        interface_materials.back().Gf   = interface_Gf;
 
        EntityHandle meshset = it->getMeshset();
        Range tris;  //All the triangles after splitting
        rval = moab.get_entities_by_type(meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
//        cout<<"tris.size()  "<< tris.size() <<endl;
        Range ents3d;
        rval = moab.get_adjacencies(tris,3,false,ents3d,moab::Interface::UNION); CHKERRQ_MOAB(rval);
        interface_materials.back().pRisms = ents3d.subset_by_type(MBPRISM);
//        cout<<"unite( interface_prisms, ents3d.subset_by_type(MBPRISM) ) "<<unite( interface_prisms, ents3d.subset_by_type(MBPRISM) )<<endl;
      }
    }
 
      ierr = m_field.seed_finite_elements(interface_materials.back().pRisms); CHKERRQ(ierr);
//      cout<<"interface_materials.back().pRisms.size()  "<< interface_materials.back().pRisms <<endl;
      ierr = m_field.add_ents_to_finite_element_by_type(interface_materials.back().pRisms,MBPRISM,"INTERFACE"); CHKERRQ(ierr);
 
 
   //  cout<<"young_modulus  "<<young_modulus<<endl;
    { //FIXME
      boost::ptr_vector<CohesiveInterfaceElement::PhysicalEquation>::iterator pit = interface_materials.begin();
      for(; pit != interface_materials.end();pit++) {
        pit->youngModulus = young_modulus;  //Young's modulus of bulk elastic material (then for interface E0=E/h)
      }
    }
 
    CohesiveInterfaceElement cohesive_elements(m_field); //make it golbal, so that its other uses can see it
    ierr = cohesive_elements.addOps("DISPLACEMENT",interface_materials); 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",preriodic_prisms);
 
 
    //Adding element LAGRANGE_ELEM_TRANS
//    lagrangian_element_trac.addLagrangiangElement("LAGRANGE_ELEM_TRANS","DISPLACEMENT","LAGRANGE_MUL_RIGID_TRANS","MESH_NODE_POSITIONS");
    ierr = m_field.add_finite_element("LAGRANGE_ELEM_TRANS",MF_ZERO); CHKERRQ(ierr);
    //============================================================================================================
    //C row as Lagrange_mul_disp and col as DISPLACEMENT
    ierr = m_field.modify_finite_element_add_field_row("LAGRANGE_ELEM_TRANS","LAGRANGE_MUL_RIGID_TRANS"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_col("LAGRANGE_ELEM_TRANS","DISPLACEMENT"); CHKERRQ(ierr);
    //CT col as Lagrange_mul_disp and row as DISPLACEMENT
    ierr = m_field.modify_finite_element_add_field_col("LAGRANGE_ELEM_TRANS","LAGRANGE_MUL_RIGID_TRANS"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_row("LAGRANGE_ELEM_TRANS","DISPLACEMENT"); CHKERRQ(ierr);
    //data
    ierr = m_field.modify_finite_element_add_field_data("LAGRANGE_ELEM_TRANS","LAGRANGE_MUL_RIGID_TRANS"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("LAGRANGE_ELEM_TRANS","DISPLACEMENT"); CHKERRQ(ierr);
    //============================================================================================================
 
 
    Range Tri103;
    for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,it)) {
      if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
 
        //total tri including old + new
        rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,Tri103,true); CHKERRQ_MOAB(rval);
        ierr = PetscPrintf(PETSC_COMM_WORLD,"Old triangles on SideSet 103 = %d\n",Tri103.size()); CHKERRQ(ierr);
        //Only triangles in the new bit-level
        Tri103 = intersect(AllTris, Tri103);
        ierr = PetscPrintf(PETSC_COMM_WORLD,"New triangles on SideSet 103 = %d\n",Tri103.size()); CHKERRQ(ierr);
        ierr = m_field.add_ents_to_finite_element_by_type(Tri103,MBTRI,"LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
      }
    }
 
//    build finite elements
    ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
    //build adjacencies
    ierr = m_field.build_adjacencies(problem_bit_level); 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,problem_bit_level); 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,"TRAN_ISOTROPIC_ELASTIC"); CHKERRQ(ierr);
    ierr = DMMoFEMAddElement(dm,"INTERFACE"); 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 = 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_load_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,"TRAN_ISOTROPIC_ELASTIC",&trans_elastic.feRhs,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
    ierr = DMMoFEMSNESSetFunction(dm,"INTERFACE",&cohesive_elements.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,"TRAN_ISOTROPIC_ELASTIC",&trans_elastic.feLhs,NULL,NULL); CHKERRQ(ierr);
    ierr = DMMoFEMSNESSetJacobian(dm,"INTERFACE",&cohesive_elements.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;
    SNES snes_one_step_only;
 
 
    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);
 
    ierr = SNESCreate(PETSC_COMM_WORLD,&snes_one_step_only); CHKERRQ(ierr);
    ierr = SNESSetFunction(snes_one_step_only,F,SnesRhs,snes_ctx); CHKERRQ(ierr);
    ierr = SNESSetJacobian(snes_one_step_only,Aij,Aij,SnesMat,snes_ctx); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes_one_step_only); CHKERRQ(ierr);
    double atol,rtol,stol;
    int maxit,maxf;
    SNESGetTolerances(snes_one_step_only,&atol,&rtol,&stol,&maxit,&maxf);
    maxit = 1;
    atol *= 0;
    SNESSetTolerances(snes_one_step_only,atol,rtol,stol,maxit,maxf);
    SNESLineSearch linesearch;
    SNESGetLineSearch(snes_one_step_only,&linesearch);
    SNESLineSearchSetType(linesearch,SNESLINESEARCHBASIC);
 
 
 
    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 = m_field.loop_finite_elements("PLASTIC_PROB","TRAN_ISOTROPIC_ELASTIC", 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 = SNESSetLagJacobian(snes_one_step_only,-2);
        ierr = SNESSolve(snes_one_step_only,PETSC_NULL,D); CHKERRQ(ierr);
        ierr = SNESGetConvergedReason(snes_one_step_only,&reason); CHKERRQ(ierr);
        if(reason<0) {
          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);
 
        if(reason<0) {
          t -= delta_time;
          delta_time *= 0.1;
          ierr = VecCopy(D0,D); CHKERRQ(ierr);
        } else {const int its_d = 60;
          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);
 
 
           ierr = DMoFEMLoopFiniteElements(
              dm,"INTERFACE",&cohesive_elements.feHistory
            ); CHKERRQ(ierr);
 
         //Save data on mesh
          {
            ierr = DMoFEMLoopFiniteElements(
              dm,"PLASTIC",&post_proc
            ); CHKERRQ(ierr);
 
            string out_file_name;
            {
              std::ostringstream stm;
              stm << "out_matrix_" << 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);
 
 
 
            ierr = DMoFEMLoopFiniteElements(
              dm,"TRAN_ISOTROPIC_ELASTIC",&post_proc
            ); CHKERRQ(ierr);
                       {
              std::ostringstream stm;
              stm << "out_fibres_" << 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();
}

Variable Documentation

help

char help[]

static

Initial value:

=
"...\n"
"\n"

Definition at line 71 of file rve_mech_plas_interface_periodic.cpp.