homogenisation_periodic_atom_test.cpp File Reference

#include <MoFEM.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <petsctime.h>
#include <adolc/adolc.h>
#include <NonLinearElasticElement.hpp>
#include <Hooke.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/iostreams/tee.hpp>
#include <boost/iostreams/stream.hpp>
#include <fstream>
#include <iostream>
#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 <boost/ptr_container/ptr_map.hpp>
#include <PostProcOnRefMesh.hpp>
#include <PostProcStresses.hpp>

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

Macros
#define	RND_EPS   1e-6

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
double	roundn (double n)
int	main (int argc, char *argv[])

Variables
static char	help [] = "...\n\n"

Macro Definition Documentation

RND_EPS

#define RND_EPS   1e-6

Definition at line 102 of file homogenisation_periodic_atom_test.cpp.

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 96 of file homogenisation_periodic_atom_test.cpp.

Function Documentation

main()

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

Definition at line 137 of file homogenisation_periodic_atom_test.cpp.

                                 {
 
  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
 
  try {
 
  moab::Core mb_instance;
  moab::Interface& moab = mb_instance;
  int rank;
  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 
  //Reade parameters from line command
  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 = 5;
  }
 
  //Applied strain on the RVE (vector of length 6) strain=[xx, yy, zz, xy, xz, zy]^T
  double myapplied_strain[6];
  int nmax=6;
  ierr = PetscOptionsGetRealArray(PETSC_NULL,"-myapplied_strain",myapplied_strain,&nmax,&flg); CHKERRQ(ierr);
  VectorDouble applied_strain;
  applied_strain.resize(6);
  cblas_dcopy(6, &myapplied_strain[0], 1, &applied_strain(0), 1);
  //    cout<<"applied_strain ="<<applied_strain<<endl;
 
 
  //Read mesh to MOAB
  const char *option;
  option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
  if(pcomm == NULL) pcomm =  new ParallelComm(&moab,PETSC_COMM_WORLD);
 
  //We need that for code profiling
  PetscLogDouble t1,t2;
  PetscLogDouble v1,v2;
  ierr = PetscTime(&v1); CHKERRQ(ierr);
  ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
 
  //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);
 
 
//    //=======================================================================================================
//    //=======================================================================================================
 
  //*
  //Define problem
 
  //Fields
  int field_rank=3;
  ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,field_rank,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
  ierr = m_field.add_field("LAGRANGE_MUL_PERIODIC",H1,AINSWORTH_LEGENDRE_BASE,field_rank); 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 rotations)
 
  //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);
 
  //Adding only -ve surfaces to the field LAGRANGE_MUL_PERIODIC (in periodic boundary conditions size of C (3M/2 x 3N))
  //to create meshset from range  SurTrisNeg
  EntityHandle SurTrisNegMeshset;
  rval = moab.create_meshset(MESHSET_SET,SurTrisNegMeshset); CHKERRQ_MOAB(rval);
  rval = moab.add_entities(SurTrisNegMeshset,SurTrisNeg); CHKERRQ_MOAB(rval);
  ierr = m_field.add_ents_to_field_by_type(SurTrisNegMeshset,MBTRI,"LAGRANGE_MUL_PERIODIC"); CHKERRQ(ierr);
  //============================================================================================================
 
  ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",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",2); CHKERRQ(ierr);
  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
 
  //ierr = m_field.set_field_order(0,MBPRISM,"LAGRANGE_MUL_PERIODIC",order); 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);
 
  //Define FE
  //define eleatic element
  boost::shared_ptr<Hooke<adouble> > hooke_adouble(new Hooke<adouble>);
  boost::shared_ptr<Hooke<double> > hooke_double(new Hooke<double>);
  NonlinearElasticElement elastic(m_field,2);
  ierr = elastic.setBlocks(hooke_double,hooke_adouble); CHKERRQ(ierr);
  ierr = elastic.addElement("ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
  ierr = elastic.setOperators("DISPLACEMENT","MESH_NODE_POSITIONS",false,true); 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");
 
  //define problems
  ierr = m_field.add_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
 
  //set finite elements for problem
  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","ELASTIC"); CHKERRQ(ierr);
  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS"); CHKERRQ(ierr);
 
 
  //set refinement level for problem
  ierr = m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",bit_level0); CHKERRQ(ierr);
 
  /***/
  //Declare problem
  /****/
  //build database
 
  //build field
  ierr = m_field.build_fields(); CHKERRQ(ierr);
  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
 
 
  //build finite elemnts
  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
 
  //build adjacencies
  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
 
  //build problem
  ierr = m_field.build_problems(); CHKERRQ(ierr);
 
  /****/
  //mesh partitioning
 
  //partition
  ierr = m_field.partition_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
  ierr = m_field.partition_finite_elements("ELASTIC_MECHANICS"); CHKERRQ(ierr);
  //what are ghost nodes, see Petsc Manual
  ierr = m_field.partition_ghost_dofs("ELASTIC_MECHANICS"); CHKERRQ(ierr);
 
  //print bcs
  ierr = m_field.print_cubit_displacement_set(); CHKERRQ(ierr);
  ierr = m_field.print_cubit_force_set(); CHKERRQ(ierr);
  //print block sets with materials
  ierr = m_field.print_cubit_materials_set(); CHKERRQ(ierr);
 
  //create matrices (here F, D and Aij are matrices for the full problem)
  Vec D;
  vector<Vec> F(6);
  ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",ROW,&F[0]); CHKERRQ(ierr);
  for(int ii = 1;ii<6;ii++) {
    ierr = VecDuplicate(F[0],&F[ii]); CHKERRQ(ierr);
  }
 
  ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",COL,&D); CHKERRQ(ierr);
  Mat Aij;
  ierr = m_field.MatCreateMPIAIJWithArrays("ELASTIC_MECHANICS",&Aij); CHKERRQ(ierr);
 
  for(int ii = 0;ii<6;ii++) {
    ierr = VecZeroEntries(F[ii]); CHKERRQ(ierr);
    ierr = VecGhostUpdateBegin(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecGhostUpdateEnd(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
  }
 
  ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
 
 
  //Assemble F and Aij
  elastic.getLoopFeLhs().snes_B = Aij;
  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",elastic.getLoopFeLhs()); CHKERRQ(ierr);
 
  lagrangian_element_periodic.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_PERIODIC","MESH_NODE_POSITIONS",Aij,F);
  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element_periodic.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
 
  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);
  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM_TRANS",lagrangian_element_rigid_body_trans.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
 
  for(int ii = 0;ii<6;ii++) {
    ierr = VecGhostUpdateBegin(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
    ierr = VecGhostUpdateEnd(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(F[ii]); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(F[ii]); CHKERRQ(ierr);
  }
 
  ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); 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("ELASTIC_MECHANICS","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);
 
 
  //==============================================================================================================================
  //Post processing
  //==============================================================================================================================
  struct MyPostProc: public PostProcVolumeOnRefinedMesh {
    bool doPreProcess;
    bool doPostProcess;
    MyPostProc(MoFEM::Interface &m_field):
    PostProcVolumeOnRefinedMesh(m_field),
    doPreProcess(true),
    doPostProcess(true)
    {}
    void setDoPreProcess() { doPreProcess = true; }
    void unSetDoPreProcess() { doPreProcess = false; }
    void setDoPostProcess() { doPostProcess = true; }
    void unSetDoPostProcess() { doPostProcess = false; }
 
    
 
    PetscErrorCode preProcess() {
      PetscFunctionBegin;
      if(doPreProcess) {
        ierr = PostProcVolumeOnRefinedMesh::preProcess(); CHKERRQ(ierr);
      }
      PetscFunctionReturn(0);
    }
    PetscErrorCode postProcess() {
      PetscFunctionBegin;
      if(doPostProcess) {
        ierr = PostProcVolumeOnRefinedMesh::postProcess(); CHKERRQ(ierr);
      }
      PetscFunctionReturn(0);
    }
  };
 
  MyPostProc 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);
  ierr = post_proc.addFieldValuesGradientPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
 
  for(
      map<int,NonlinearElasticElement::BlockData>::iterator sit = elastic.setOfBlocks.begin();
      sit != elastic.setOfBlocks.end(); sit++
      ) {
    post_proc.getOpPtrVector().push_back(
   new PostPorcStress(
    post_proc.postProcMesh,
    post_proc.mapGaussPts,
    "DISPLACEMENT",
    sit->second,
    post_proc.commonData,
    false
    )
   );
  }
  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
  ierr = m_field.set_global_ghost_vector(
   "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
   ); CHKERRQ(ierr);
 
  //==============================================================================================================================
 
    //Solver
  KSP solver;
  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
  ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
  ierr = KSPSetUp(solver); CHKERRQ(ierr);
 
  MatrixDouble Dmat;
  Dmat.resize(6,6);
  Dmat.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);
  for(int ii = 0;ii<6;ii++) {
    ierr = VecZeroEntries(D); CHKERRQ(ierr);
    ierr = KSPSolve(solver,F[ii],D); CHKERRQ(ierr);
    ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = m_field.set_global_ghost_vector(
     "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
     ); CHKERRQ(ierr);
 
    //post-processing the resutls
    post_proc.setDoPreProcess();
    post_proc.unSetDoPostProcess();
    ierr = m_field.loop_finite_elements(
            "ELASTIC_MECHANICS","ELASTIC",post_proc
            ); CHKERRQ(ierr);
    post_proc.unSetDoPreProcess();
    post_proc.setDoPostProcess();
    {
      ostringstream sss;
      sss << "mode_" << "Disp" << "_" << ii << ".h5m";
      rval = post_proc.postProcMesh.write_file(
             sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART"
             ); CHKERRQ_MOAB(rval);
    }
 
    ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
    ierr = m_field.loop_finite_elements(
    "ELASTIC_MECHANICS","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++) {
      Dmat(jj,ii) = avec[jj];
    }
  }
 
  ierr = VecRestoreArray(stress_homo,&avec); CHKERRQ(ierr);
  PetscPrintf(
              PETSC_COMM_WORLD,"\nHomogenised Stiffens Matrix = \n\n"
              );
 
  for(int ii=0; ii<6; ii++) {
    PetscPrintf(
    PETSC_COMM_WORLD,
    "stress %d\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\n",
    ii,Dmat(ii,0),Dmat(ii,1),Dmat(ii,2),Dmat(ii,3),Dmat(ii,4),Dmat(ii,5)
    );
  }
  //==============================================================================================================================
 
 //Open mesh_file_name.txt for writing
  ofstream myfile;
  myfile.open ((string(mesh_file_name)+".txt").c_str());
 
  //Output displacements
  myfile << "<<<< Homonenised stress >>>>>" << endl;
 
  if(pcomm->rank()==0){
    for(int ii=0; ii<6; ii++){
      myfile << boost::format("%.3lf") % roundn(Dmat(ii,0)) << endl;
      avec++;
    }
  }
  //Close mesh_file_name.txt
  myfile.close();
 
 
 
  //detroy matrices
  for(int ii = 0;ii<6;ii++) {
    ierr = VecDestroy(&F[ii]); CHKERRQ(ierr);
  }
  ierr = VecDestroy(&D); CHKERRQ(ierr);
  ierr = MatDestroy(&Aij); CHKERRQ(ierr);
  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
  ierr = VecDestroy(&stress_homo); CHKERRQ(ierr);
 
  ierr = PetscTime(&v2);CHKERRQ(ierr);
  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
 
  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);
 
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
}

roundn()

double roundn

(

double

n

)

Definition at line 105 of file homogenisation_periodic_atom_test.cpp.

{
    //break n into fractional part (fract) and integral part (intp)
    double fract, intp;
    fract = modf(n,&intp);
 
//    //round up
//    if (fract>=.5)
//    {
//        n*=10;
//        ceil(n);
//        n/=10;
//    }
//  //round down
//    if (fract<=.5)
//    {
//      n*=10;
//        floor(n);
//        n/=10;
//    }
    // case where n approximates zero, set n to "positive" zero
    if (abs(intp)==0)
    {
        if(abs(fract)<=RND_EPS)
           {
               n=0.000;
           }
    }
    return n;
}

Variable Documentation

help

char help[] = "...\n\n"

static

Definition at line 54 of file homogenisation_periodic_atom_test.cpp.