BCs_RVELagrange_Periodic.hpp

Go to the documentation of this file.

/* 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/>. */
 
 
#ifndef __BCS_RVELAGRANGE_PERIIODIC_HPP
#define __BCS_RVELAGRANGE_PERIIODIC_HPP
 
struct BCs_RVELagrange_Periodic: public BCs_RVELagrange_Disp {
 
  boost::ptr_vector<MethodForForceScaling> methodsOp;  //This is to apply strain vector in steps (for given steps) for nonlinear problem
 
  struct RVEBC_Data_Periodic {
    Range pRisms; // All boundary surfaces
  };
  map<int,RVEBC_Data_Periodic> setOfRVEBCPrisms; ///< maps side set id with appropriate FluxData
 
  struct MyPrismFE: public FlatPrismElementForcesAndSurcesCore {
    MyPrismFE(MoFEM::Interface &_mField): FlatPrismElementForcesAndSurcesCore(_mField) {}
    int getRule(int order) { return order; };
  };
 
  MyPrismFE feRVEBCRhs; //To calculate the Rhs or RVE BCs
  MyPrismFE feRVEBCLhs; //To calculate the Lhs or RVE BCs
  MyPrismFE feRVEBCStress; //To calculate the Rhs or RVE BCs
  MyPrismFE feRVEBCRhsResidual; //used to calculate the residual
 
  MyPrismFE& getLoopFeRVEBCLhs() { return feRVEBCLhs; }
  MyPrismFE& getLoopFeRVEBCRhs() { return feRVEBCRhs; }
 
  MyPrismFE& getLoopFeRVEBCRhsResidual() { return feRVEBCRhsResidual; }
  MyPrismFE& getLoopFeRVEBCStress() { return feRVEBCStress; }
 
  BCs_RVELagrange_Periodic(MoFEM::Interface &m_field):
  BCs_RVELagrange_Disp(m_field),
  feRVEBCRhs(m_field),
  feRVEBCLhs(m_field),
  feRVEBCRhsResidual(m_field),
  feRVEBCStress(m_field) {
  }
 
  PetscErrorCode addLagrangiangElement(
    const string element_name,
    const string field_name,
    const string lagrang_field_name,
    const string mesh_nodals_positions,
    Range &periodic_prisms
  ) {
    PetscFunctionBegin;
    
    ierr = mField.add_finite_element(element_name,MF_ZERO); CHKERRQ(ierr);
    //============================================================================================================
    //C row as Lagrange_mul_disp and col as DISPLACEMENT
    ierr = mField.modify_finite_element_add_field_row(element_name,lagrang_field_name); CHKERRQ(ierr);
    ierr = mField.modify_finite_element_add_field_col(element_name,field_name); CHKERRQ(ierr);
    //CT col as Lagrange_mul_disp and row as DISPLACEMENT
    ierr = mField.modify_finite_element_add_field_col(element_name,lagrang_field_name); CHKERRQ(ierr);
    ierr = mField.modify_finite_element_add_field_row(element_name,field_name); CHKERRQ(ierr);
    //data
    ierr = mField.modify_finite_element_add_field_data(element_name,lagrang_field_name); CHKERRQ(ierr);
    ierr = mField.modify_finite_element_add_field_data(element_name,field_name); CHKERRQ(ierr);
    //============================================================================================================
 
    if(mField.check_field(mesh_nodals_positions)) { //for high order geometry
      ierr = mField.modify_finite_element_add_field_data(element_name,mesh_nodals_positions); CHKERRQ(ierr);
    }
 
    //Adding Prisims to Element Lagrange_elem (to loop over these prisims)
    setOfRVEBCPrisms[1].pRisms = periodic_prisms;
    ierr = mField.add_ents_to_finite_element_by_type(periodic_prisms,MBPRISM,element_name); CHKERRQ(ierr);
 
    
    Range tris;
    rval = mField.get_moab().get_adjacencies(periodic_prisms,2,false,tris,moab::Interface::UNION); CHKERRQ_MOAB(rval);
    setOfRVEBC[1].tRis = tris.subset_by_type(MBTRI);
 
    PetscFunctionReturn(0);
  }
 
  struct CommonFunctionsPeriodic {
    PetscErrorCode shapeMat(
      int rank, unsigned int gg, DataForcesAndSourcesCore::EntData &col_data, MatrixDouble &N_mat, int div=1
    ) {
      PetscFunctionBegin;
 
      //we need it only for one face for periodic BC  (col_data.getN() are shape functions inclding both faces, i.e. 6 shape functions for nodes)
      int shape_size=col_data.getN().size2()/div;
      //      cout<<"shape_size = "<<shape_size<<endl;
      //      cout<<"rank = "<<rank<<endl;
      N_mat.resize(rank,shape_size*rank);    N_mat.clear();
      int gg1=0;       int kk=0;
      for(int ii=0; ii<shape_size; ii++){ //number of shape functions
        for(int jj=0; jj<rank; jj++){
          //              cout<<"ii "<<ii<<endl;
          N_mat(jj,kk)=col_data.getN()(gg,gg1); kk++;
        }
        gg1++;
      }
      PetscFunctionReturn(0);
    }
 
 
 
    PetscErrorCode shapeMatNew(
      int rank, unsigned int gg, MatrixDouble &RowN, MatrixDouble &N_mat, int div=1
    ) {
      PetscFunctionBegin;
 
      //we need it only for one face for periodic BC  (col_data.getN() are shape functions inclding both faces, i.e. 6 shape functions for nodes)
      int shape_size=RowN.size2()/div;
      //      cout<<"shape_size = "<<shape_size<<endl;
      //      cout<<"rank = "<<rank<<endl;
      N_mat.resize(rank,shape_size*rank);    N_mat.clear();
      int gg1=0;       int kk=0;
      for(int ii=0; ii<shape_size; ii++){ //number of shape functions
        for(int jj=0; jj<rank; jj++){
          //              cout<<"ii "<<ii<<endl;
          N_mat(jj,kk)=RowN(gg,gg1); kk++;
        }
        gg1++;
      }
      PetscFunctionReturn(0);
    }
 
  };
  CommonFunctionsPeriodic commonFunctionsPeriodic;
 
 
  struct CommonDataPeriodic {
    ublas::vector<MatrixDouble > D_mat;
    ublas::vector<map<EntityType, map<int, map<EntityType,  map<int, MatrixDouble > > > > > Cmat;  //Cmat[ff][row_type][row_side][col_type][col_side]
    ublas::vector<map<EntityType, map<int, ublas::vector<int> > > > RowInd;  //RowInd[ff][row_sdie][row_side]
    ublas::vector<map<EntityType, map<int, ublas::vector<int> > > > ColInd;  //ColInd[ff][col_type][col_side]
 
 
 
    ublas::vector<map<EntityType, map<int, MatrixDouble > > > RowN;  //RowN[ff][col_type][col_side]
    ublas::vector<map<EntityType, map<int, MatrixDouble > > > ColN;  //ColN[ff][col_type][col_side]
 
 
    //these are used in nonlineaer problem to caluclate residual vector
    ublas::vector<ublas::vector<VectorDouble > > DispAtGaussPts;   //DispAtGaussPts[ff][gg]   -- different on both faces of prism
    ublas::vector<VectorDouble >                 LagMulAtGaussPts; //LagMulAtGaussPts[gg] -- same on both faces of prism
   };
  CommonDataPeriodic commonDataPeriodic;
 
 
 
 
  //Operator to calculate column indices and colum shapue functions (N)
  struct OpRVEBCsPeriodicColInd:public FlatPrismElementForcesAndSurcesCore::UserDataOperator {
 
    RVEBC_Data_Periodic &dAta;
    bool hoGeometry;
    CommonDataPeriodic &commonDataPeriodic;
    CommonFunctionsPeriodic &commonFunctionsPeriodic;
 
    OpRVEBCsPeriodicColInd(
      const string field_name,
      RVEBC_Data_Periodic &data,
      CommonDataPeriodic &common_data_periodic,
      CommonFunctionsPeriodic &common_functions_periodic,
      bool ho_geometry = false
    ):
    FlatPrismElementForcesAndSurcesCore::UserDataOperator(
      field_name,UserDataOperator::OPCOL
    ),
    dAta(data),
    hoGeometry(ho_geometry),
    commonDataPeriodic(common_data_periodic),
    commonFunctionsPeriodic(common_functions_periodic) {
    }
 
 
    ublas::vector<VectorDouble > field_data_nodes;  //field_data_nodes[ff]
    MatrixDouble N_mat;
    
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
 
        if(type==MBVERTEX){
          commonDataPeriodic.ColInd.resize(2);
          commonDataPeriodic.ColN.resize(2);
        }
        int nb;
        //Global indices Columns
        if(data.getIndices().size()==0) PetscFunctionReturn(0);
 
        switch(type){
          case MBVERTEX:
            nb=data.getIndices().size()/2;  //for one face for nodes only  (total indieces are for both faces)
            commonDataPeriodic.ColInd[0][type][side].resize(nb); //1st face of prism -- nodes
            commonDataPeriodic.ColInd[1][type][side].resize(nb); //2nd face of prism -- nodes
            //different indices for cols
            for(int ii=0; ii<nb; ii++){
              commonDataPeriodic.ColInd[0][type][side][ii]=data.getIndices()[ii];
              commonDataPeriodic.ColInd[1][type][side][ii]=data.getIndices()[ii+nb];
            }
//              cout<<"commonDataPeriodic.ColInd[0] "<<"type "<<type << " side "<< side << "     = "<<  commonDataPeriodic.ColInd[0][type][side]<<endl;
//              cout<<"commonDataPeriodic.ColInd[0] "<<"type "<<type << " side "<< side << "     = "<<  commonDataPeriodic.ColInd[1][type][side]<<endl;
//              cout<<"type "<<type << " side "<< side <<endl;
            commonDataPeriodic.ColN[0][type][side]=data.getN();
            commonDataPeriodic.ColN[1][type][side]=data.getN();
            break;
          case MBEDGE:
            if(side<3){ //1st face eges 0, 1, 2 (cananical numbering)-- same indices for both faces
              commonDataPeriodic.ColInd[0][type][side]=data.getIndices();
              commonDataPeriodic.ColN[0][type][side]=data.getN();
//              cout<<"col_type "<<col_type << " col_side "<< col_side <<endl;
//              cout<<"col_data.getIndices "<<col_data.getIndices()<<endl;
              }
            else if (side >= 6){ //2nd face edges 6,7,8
              commonDataPeriodic.ColInd[1][type][side-6]=data.getIndices();
              commonDataPeriodic.ColN[1][type][side-6]=data.getN();
//              cout<<"col_type "<<col_type << " col_side "<< col_side <<endl;
              }
            break;
          case MBTRI:
            if(side==3){//1st face face 3
              commonDataPeriodic.ColInd[0][type][side]=data.getIndices();
              commonDataPeriodic.ColN[0][type][side]=data.getN();
              }
            else{ //2nd face face 4
              commonDataPeriodic.ColInd[1][type][side-1]=data.getIndices();
              commonDataPeriodic.ColN[1][type][side-1]=data.getN();
              }
            break;
          default:
            SETERRQ(PETSC_COMM_SELF,1,"data inconsitency");
        }
      PetscFunctionReturn(0);
    }
  };
 
 
 
  //Operator to calculate row indices and row shapue functions (N)
  struct OpRVEBCsPeriodicRowInd:public FlatPrismElementForcesAndSurcesCore::UserDataOperator {
 
    RVEBC_Data_Periodic &dAta;
    bool hoGeometry;
    CommonDataPeriodic &commonDataPeriodic;
    CommonFunctionsPeriodic &commonFunctionsPeriodic;
 
    OpRVEBCsPeriodicRowInd(
      const string field_name,
      RVEBC_Data_Periodic &data,
      CommonDataPeriodic &common_data_periodic,
      CommonFunctionsPeriodic &common_functions_periodic,
      bool ho_geometry = false
    ):
    FlatPrismElementForcesAndSurcesCore::UserDataOperator(
      field_name,UserDataOperator::OPROW
    ),
    dAta(data),
    hoGeometry(ho_geometry),
    commonDataPeriodic(common_data_periodic),
    commonFunctionsPeriodic(common_functions_periodic) {
    }
 
 
    ublas::vector<VectorDouble > field_data_nodes;  //field_data_nodes[ff]
    MatrixDouble N_mat;
    
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
//        cout<<"type "<<type << " side "<< side <<endl;
        if(data.getIndices().size()==0) PetscFunctionReturn(0);
        if(type == MBEDGE && side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
        if(type == MBTRI && side == 4) PetscFunctionReturn(0);  //ignore second triangel face (cananical number 4)
 
 
        if(type==MBVERTEX){
          commonDataPeriodic.RowInd.resize(2);
          commonDataPeriodic.RowN.resize(2);
        }
 
        int nb;
        //Global indices Rows
        switch(type){
          case MBVERTEX:
            nb=data.getIndices().size()/2;
            commonDataPeriodic.RowInd[0][type][side].resize(nb); //1st face of prism -- nodes
            commonDataPeriodic.RowInd[1][type][side].resize(nb); //2nd face of prism -- nodes
 
            //same indices for rows
            for(int ii=0; ii<nb; ii++){
              commonDataPeriodic.RowInd[0][type][side][ii]=data.getIndices()[ii];
              commonDataPeriodic.RowInd[1][type][side][ii]=data.getIndices()[ii];
            }
//            cout<<"type "<<type << " side "<< side <<endl;
//            cout<<"commonDataPeriodic.RowInd[0] "<<"type "<<type << " side "<< side << "     = "<<  commonDataPeriodic.RowInd[0][type][side]<<endl;
            commonDataPeriodic.RowN[0][type][side]=data.getN();
            commonDataPeriodic.RowN[1][type][side]=data.getN();
            break;
          case MBEDGE:
            if(side<3){ //1st face eges 0, 1, 2 (cananical numbering)-- same indices for both faces
              commonDataPeriodic.RowInd[0][type][side]=data.getIndices();
              commonDataPeriodic.RowInd[1][type][side]=data.getIndices();
              commonDataPeriodic.RowN[0][type][side]=data.getN();
              commonDataPeriodic.RowN[1][type][side]=data.getN();
//              cout<<"commonDataPeriodic.RowInd[0] "<<"type "<<type << " side "<< side << "     = "<<  commonDataPeriodic.RowInd[0][type][side]<<endl;
//            cout<<"type "<<type << " side "<< side <<endl;
            }
            break;
          case MBTRI:
            if(side==3){//1st face face 3
              commonDataPeriodic.RowInd[0][type][side]=data.getIndices();
              commonDataPeriodic.RowInd[1][type][side]=data.getIndices();
              commonDataPeriodic.RowN[0][type][side]=data.getN();
              commonDataPeriodic.RowN[1][type][side]=data.getN();
              }
            break;
          default:
            SETERRQ(PETSC_COMM_SELF,1,"data inconsitency");
        }
 
//        string aaa;
//        cin>>aaa;
      PetscFunctionReturn(0);
    }
  };
 
 
  /// \biref operator to calculate and assemble Cmat
  struct OpRVEBCsPeriodicCalAssemCmat:public FlatPrismElementForcesAndSurcesCore::UserDataOperator  {
    Mat Aij;
    RVEBC_Data_Periodic &dAta;
    CommonFunctionsPeriodic commonFunctionsPeriodic;
    CommonDataPeriodic &commonDataPeriodic;
    bool hoGeometry;
 
    OpRVEBCsPeriodicCalAssemCmat(
      const string field_name,
      const string lagrang_field_name,
      Mat aij,
      RVEBC_Data_Periodic &data,
      CommonFunctionsPeriodic &common_functions_periodic,
      CommonDataPeriodic &common_data_periodic,
      bool ho_geometry = false
    ):
 
    FlatPrismElementForcesAndSurcesCore::UserDataOperator(
      lagrang_field_name, field_name, UserDataOperator::OPROWCOL
    ),Aij(aij),
    dAta(data),
    commonFunctionsPeriodic(common_functions_periodic),
    commonDataPeriodic(common_data_periodic),
    hoGeometry(ho_geometry){
      sYmm = false;  //This will make sure to loop over all intities (e.g. for order=2 it will make doWork to loop 16 time)
    }
 
    MatrixDouble Kmat, transKmat;
    ublas::vector<int> rrvec, ccvec;
    MatrixDouble NTN;
 
    PetscErrorCode doWork(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSourcesCore::EntData &row_data,
      DataForcesAndSourcesCore::EntData &col_data
    ) {
      PetscFunctionBegin;
 
      try {
//        cout<<"Hello from  OpRVEBCsPeriodicLhsAssemCmat "<<endl;
        if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
        if(col_data.getIndices().size()==0) PetscFunctionReturn(0);
        if(row_type == MBEDGE && row_side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
        if(col_type == MBEDGE && col_side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
        if(row_type == MBTRI && row_side == 4) PetscFunctionReturn(0);  //ignore second triangel face (cananical number 4)
        if(col_type == MBTRI && col_side == 4) PetscFunctionReturn(0);  //ignore second triangel face (cananical number 4)
 
 
        //Cmat calculation
        
        int rank = col_data.getFieldDofs()[0]->getNbOfCoeffs();
 
        MatrixDouble N_mat_row;
        MatrixDouble N_mat_col;
 
 
        //1st face
        for(unsigned int gg = 0;gg<col_data.getN().size1();gg++) {
          double area;
          VectorDouble normal_f3 = getNormalsAtGaussPtsF3(gg);
          area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
          double val = getGaussPts()(2,gg)*area;
 
          if(row_type == MBVERTEX) {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.RowN[0][row_type][row_side], N_mat_row, 2); CHKERRQ(ierr);
//            cerr << "N_mat_row  "<< N_mat_row <<  endl;
            }
          else {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.RowN[0][row_type][row_side], N_mat_row); CHKERRQ(ierr);}
 
          if(col_type == MBVERTEX) {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.ColN[0][col_type][col_side], N_mat_col, 2); CHKERRQ(ierr);
//            cerr << "N_mat_row  "<< N_mat_row <<  endl;
            }
          else {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.ColN[0][col_type][col_side], N_mat_col); CHKERRQ(ierr);}
          if(gg==0){
//            cerr << "row_type  "<< row_type <<" row_side "<<row_side << "    "<<  endl;
//            cerr << "col_type  "<< col_type <<" col_side "<< col_side<< "    "<<  endl;
            NTN=prod(trans(N_mat_row),N_mat_col);  //we don't need to define size of NTN
            Kmat=-1.0*val*NTN; //we don't need to defien size of K
          }else{
            NTN=prod(trans(N_mat_row),N_mat_col);
            Kmat+=-1.0*val*NTN;
          }
        }
        rrvec=commonDataPeriodic.RowInd[0][row_type][row_side];
        ccvec=commonDataPeriodic.ColInd[0][col_type][col_side];
        ierr = MatSetValues(Aij,rrvec.size(),&rrvec[0],ccvec.size(),&ccvec[0],&Kmat(0,0),ADD_VALUES); CHKERRQ(ierr);
        transKmat = trans(Kmat);
        ierr = MatSetValues(Aij,ccvec.size(),&ccvec[0],rrvec.size(),&rrvec[0],&transKmat(0,0),ADD_VALUES); CHKERRQ(ierr);
 
//        cerr << "rrvec  "<< rrvec <<  endl;
//        cerr << "ccvec  "<< ccvec <<  endl;
//        cerr << "Kmat  "<< Kmat <<  endl;
 
 
        //2nd face
        for(unsigned int gg = 0;gg<col_data.getN().size1();gg++) {
          double area;
          VectorDouble normal_f3 = getNormalsAtGaussPtsF3(gg);
          area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
          double val = getGaussPts()(2,gg)*area;
 
          if(row_type == MBVERTEX) {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.RowN[1][row_type][row_side], N_mat_row, 2); CHKERRQ(ierr);
//            cerr << "N_mat_row  "<< N_mat_row <<  endl;
            }
          else {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.RowN[1][row_type][row_side], N_mat_row); CHKERRQ(ierr);}
 
          if(col_type == MBVERTEX) {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.ColN[1][col_type][col_side], N_mat_col, 2); CHKERRQ(ierr);
//            cerr << "N_mat_row  "<< N_mat_row <<  endl;
            }
          else {
            ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.ColN[1][col_type][col_side], N_mat_col); CHKERRQ(ierr);}
          if(gg==0){
//            cerr << "row_type  "<< row_type <<" row_side "<<row_side << "    "<<  endl;
//            cerr << "col_type  "<< col_type <<" col_side "<< col_side<< "    "<<  endl;
            NTN=prod(trans(N_mat_row),N_mat_col);  //we don't need to define size of NTN
            Kmat=val*NTN; //we don't need to defien size of K
          }else{
            NTN=prod(trans(N_mat_row),N_mat_col);
            Kmat+=val*NTN;
          }
        }
        rrvec=commonDataPeriodic.RowInd[1][row_type][row_side];
        ccvec=commonDataPeriodic.ColInd[1][col_type][col_side];
        ierr = MatSetValues(Aij,rrvec.size(),&rrvec[0],ccvec.size(),&ccvec[0],&Kmat(0,0),ADD_VALUES); CHKERRQ(ierr);
        transKmat = trans(Kmat);
        ierr = MatSetValues(Aij,ccvec.size(),&ccvec[0],rrvec.size(),&rrvec[0],&transKmat(0,0),ADD_VALUES); CHKERRQ(ierr);
 
 
//        cerr << "\n\n\nrrvec  "<< rrvec <<  endl;
//        cerr << "ccvec  "<< ccvec <<  endl;
//        cerr << "Kmat  "<< Kmat <<  endl;
 
 
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
//      string wait;
//      cin>>wait;
      PetscFunctionReturn(0);
    }
  };
 
 
 
 
  /// \biref operator to calculate the RHS of the constrain for the RVE boundary conditions
  struct OpDmatRhs:public FlatPrismElementForcesAndSurcesCore::UserDataOperator {
 
    RVEBC_Data_Periodic &dAta;
    bool hoGeometry;
    CommonDataPeriodic &commonDataPeriodic;
    CommonFunctionsPeriodic &commonFunctionsPeriodic;
 
    OpDmatRhs(
      const string field_name,
      const string lagrang_field_name,
      RVEBC_Data_Periodic &data,
      CommonDataPeriodic &common_data_periodic,
      CommonFunctionsPeriodic &common_functions_periodic,
      bool ho_geometry = false
    ):
    FlatPrismElementForcesAndSurcesCore::UserDataOperator(
      lagrang_field_name,UserDataOperator::OPROW
    ),
    dAta(data),
    hoGeometry(ho_geometry),
    commonDataPeriodic(common_data_periodic),
    commonFunctionsPeriodic(common_functions_periodic) {
    }
 
    ublas::vector<MatrixDouble > X_mat;  //Here X_mat are different for left and right face
    MatrixDouble N_mat;
 
    PetscErrorCode calculateDmat(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
//      if(data.getIndices().size()==0) PetscFunctionReturn(0);
//      if(type == MBEDGE && side >= 6) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
//      if(type == MBTRI && side == 4) PetscFunctionReturn(0);  //ignore second triangel face (cananical number 4)
      
      int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
 
      double x,y,z,x1,y1,z1;
      for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
        double area;
        VectorDouble normal_f3 = getNormalsAtGaussPtsF3(gg);
        area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
        double val = getGaussPts()(2,gg)*area;
 
        // Gauss point coordinates of the prism for both faces with HO geometry
        x = getHOCoordsAtGaussPtsF3()(gg,0);
        y = getHOCoordsAtGaussPtsF3()(gg,1);
        z = getHOCoordsAtGaussPtsF3()(gg,2);
 
        x1 = getHOCoordsAtGaussPtsF4()(gg,0);
        y1 = getHOCoordsAtGaussPtsF4()(gg,1);
        z1 = getHOCoordsAtGaussPtsF4()(gg,2);
 
//        cout<<"x "<<x << " y "<< y << " z   "<<  z <<endl;
//        cout<<"x1 "<<x1 << " y1 "<< y1 << " z1   "<<  z1 <<endl;
 
        commonDataPeriodic.D_mat.resize(2);
        X_mat.resize(2);
 
        switch(rank) {
          case 3: //mech problem
          X_mat[0].resize(rank,6,false);  X_mat[0].clear();
          X_mat[0](0,0)=2.0*x;  X_mat[0](0,3)=y;  X_mat[0](0,5)=z;
          X_mat[0](1,1)=2.0*y;  X_mat[0](1,3)=x;  X_mat[0](1,4)=z;
          X_mat[0](2,2)=2.0*z;  X_mat[0](2,4)=y;  X_mat[0](2,5)=x;
          X_mat[0]=0.5*X_mat[0];
 
          X_mat[1].resize(rank,6,false); X_mat[1].clear();
          X_mat[1](0,0)=2.0*x1;  X_mat[1](0,3)=y1;  X_mat[1](0,5)=z1;
          X_mat[1](1,1)=2.0*y1;  X_mat[1](1,3)=x1;  X_mat[1](1,4)=z1;
          X_mat[1](2,2)=2.0*z1;  X_mat[1](2,4)=y1;  X_mat[1](2,5)=x1;
          X_mat[1]=0.5*X_mat[1];
          break;
          case 1:  //moisture transport or thermal problem
          X_mat[0].resize(rank,3,false); X_mat[0].clear();
          X_mat[0](0,0)=x; X_mat[0](0,1)=y; X_mat[0](0,2)=z;
 
          X_mat[1].resize(rank,3,false);  X_mat[1].clear();
          X_mat[1](0,0)=x1;  X_mat[1](0,1)=y1;  X_mat[1](0,2)=z1;
          break;
          default:
          SETERRQ(PETSC_COMM_SELF,1,"not implemented");
        }
 
        //
        if(type == MBVERTEX) {
          ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat, 2); CHKERRQ(ierr);
        } else {
          ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat); CHKERRQ(ierr);
        }
 
        if(gg==0) {
          commonDataPeriodic.D_mat[0]=-1*val*prod(trans(N_mat),X_mat[0]);
          commonDataPeriodic.D_mat[1]=   val*prod(trans(N_mat),X_mat[1]);
        } else {
          commonDataPeriodic.D_mat[0]+=-1*val*prod(trans(N_mat),X_mat[0]);
          commonDataPeriodic.D_mat[1]+=   val*prod(trans(N_mat),X_mat[1]);
        }
      }
 
//      cout<<endl<<endl;
 
      PetscFunctionReturn(0);
    }
  };
 
  struct OpRVEBCsPeriodicRhs:public OpDmatRhs {
 
    vector<Vec> &F;
 
    OpRVEBCsPeriodicRhs(
      const string field_name,
      const string lagrang_field_name,
      vector<Vec> &f,
      RVEBC_Data_Periodic &data,
      CommonDataPeriodic &common_data_periodic,
      CommonFunctionsPeriodic &common_functions_periodic,
      bool ho_geometry = false
    ):
    OpDmatRhs(
      field_name,lagrang_field_name,data,common_data_periodic,common_functions_periodic,ho_geometry
    ),
    F(f) {
    }
 
    VectorDouble appliedStrain;
    VectorDouble nF;
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
      if(data.getIndices().size()==0) PetscFunctionReturn(0);
      if(type == MBEDGE && side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
      if(type == MBTRI && side == 4) PetscFunctionReturn(0);  //ignore second triangel face (cananical number 4)
 
      
      ierr = calculateDmat(side,type,data); CHKERRQ(ierr);
 
//      if(type == MBVERTEX) {
//      cout<<"\n\nNew element "<<endl;
//      }
//      cout<<"type "<<type << " side "<< side <<endl;
//      cout<< "commonDataPeriodic.D_mat[0] "<<commonDataPeriodic.D_mat[0]<<endl;
//      cout<< "commonDataPeriodic.D_mat[1] "<<commonDataPeriodic.D_mat[1]<<endl<<endl;
//      string aaa;
//      cin >> aaa;
 
 
      ublas::vector<int> rowvec;
      if(type == MBVERTEX) {
        int nb=data.getIndices().size()/2;  //for one face
        rowvec.resize(nb);
        for(int ii=0; ii<nb; ii++){
          rowvec[ii]=data.getIndices()[ii];
        }
      } else{
        rowvec=data.getIndices();
      }
 
      appliedStrain.resize(commonDataPeriodic.D_mat[0].size2(),false);
 
      int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
      int size = (rank == 3) ? 6 : 3;
      for(int ii = 0;ii<size;ii++) {
        appliedStrain.clear();
        appliedStrain[ii] = 1.0;
        nF = prod(commonDataPeriodic.D_mat[0], appliedStrain);
        ierr = VecSetValues(F[ii],rowvec.size(),&rowvec[0],&nF[0],ADD_VALUES); CHKERRQ(ierr);
        nF = prod(commonDataPeriodic.D_mat[1], appliedStrain);
        ierr = VecSetValues(F[ii],rowvec.size(),&rowvec[0],&nF[0],ADD_VALUES); CHKERRQ(ierr);
      }
 
      PetscFunctionReturn(0);
    }
  };
 
 
 
 
  struct OpRVEBCsPeriodicRhs_givenStrain:public OpDmatRhs {
 
    Vec F;
    VectorDouble givenStrain;
    boost::ptr_vector<MethodForForceScaling> &methodsOp;
 
    OpRVEBCsPeriodicRhs_givenStrain(
      const string field_name,
      const string lagrang_field_name,
      Vec f,
      VectorDouble given_strain,
      boost::ptr_vector<MethodForForceScaling> &methods_op,
      RVEBC_Data_Periodic &data,
      CommonDataPeriodic &common_data_periodic,
      CommonFunctionsPeriodic &common_functions_periodic,
      bool ho_geometry = false
    ):
    OpDmatRhs(
      field_name,lagrang_field_name,data,common_data_periodic,common_functions_periodic,ho_geometry
    ),
    F(f), givenStrain(given_strain),methodsOp(methods_op){
    }
 
    VectorDouble nF;
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
      if(data.getIndices().size()==0) PetscFunctionReturn(0);
      if(type == MBEDGE && side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
      if(type == MBTRI && side == 4) PetscFunctionReturn(0);  //ignore second triangel face (cananical number 4)
 
      
      ierr = calculateDmat(side,type,data); CHKERRQ(ierr);
 
      ublas::vector<int> rowvec;
      if(type == MBVERTEX) {
        int nb=data.getIndices().size()/2;  //for one face
        rowvec.resize(nb);
        for(int ii=0; ii<nb; ii++){
          rowvec[ii]=data.getIndices()[ii];
        }
      } else{
        rowvec=data.getIndices();
      }
 
      //This is to scale givenStrain vector to apply the strian in steps
      VectorDouble scaled_given_strain;
      scaled_given_strain.resize(6);
      scaled_given_strain=givenStrain;
//      cout <<"givenStrain Before = "<<scaled_given_strain<<endl;
      ierr = MethodForForceScaling::applyScale(getFEMethod(),methodsOp,scaled_given_strain); CHKERRQ(ierr);
//      cout <<"givenStrain  After = "<<scaled_given_strain<<endl;
 
//      SNES snes = getFEMethod()->snes;
//      int iter;
//      SNESGetIterationNumber(snes,&iter);
//      if(iter>0)  PetscFunctionReturn(0);
 
      nF = prod(commonDataPeriodic.D_mat[0], scaled_given_strain);
      nF*=-1;
      ierr = VecSetValues(F,rowvec.size(),&rowvec[0],&nF[0],ADD_VALUES); CHKERRQ(ierr);
 
      nF = prod(commonDataPeriodic.D_mat[1], scaled_given_strain);
      nF*=-1;
      ierr = VecSetValues(F,rowvec.size(),&rowvec[0],&nF[0],ADD_VALUES); CHKERRQ(ierr);
 
      PetscFunctionReturn(0);
    }
  };
 
 
  struct OpRVEBCsPeriodicCalDispAtGaussPts:public FlatPrismElementForcesAndSurcesCore::UserDataOperator {
 
    RVEBC_Data_Periodic &dAta;
    bool hoGeometry;
    CommonDataPeriodic &commonDataPeriodic;
    CommonFunctionsPeriodic &commonFunctionsPeriodic;
 
    OpRVEBCsPeriodicCalDispAtGaussPts(
      const string field_name,
      RVEBC_Data_Periodic &data,
      CommonDataPeriodic &common_data_periodic,
      CommonFunctionsPeriodic &common_functions_periodic,
      bool ho_geometry = false
    ):
    FlatPrismElementForcesAndSurcesCore::UserDataOperator(
      field_name,UserDataOperator::OPCOL
    ),
    dAta(data),
    hoGeometry(ho_geometry),
    commonDataPeriodic(common_data_periodic),
    commonFunctionsPeriodic(common_functions_periodic) {
    }
 
 
    ublas::vector<VectorDouble > field_data_nodes;  //field_data_nodes[ff]
    MatrixDouble N_mat;
    
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
//      cout<<" OpRVEBCsPeriodicCalDispAtGaussPts "<<endl;
      if(data.getIndices().size()==0) PetscFunctionReturn(0);
      int nb_gauss_pts = data.getN().size1();
      int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
 
       //initialize
      if(type == MBVERTEX) {
        commonDataPeriodic.DispAtGaussPts.resize(2);
        field_data_nodes.resize(2);
 
        commonDataPeriodic.DispAtGaussPts[0].resize(nb_gauss_pts);
        commonDataPeriodic.DispAtGaussPts[1].resize(nb_gauss_pts);
        for(int gg = 0;gg<nb_gauss_pts;gg++) {
          commonDataPeriodic.DispAtGaussPts[0][gg].resize(rank);  commonDataPeriodic.DispAtGaussPts[0][gg].clear();
          commonDataPeriodic.DispAtGaussPts[1][gg].resize(rank);  commonDataPeriodic.DispAtGaussPts[1][gg].clear();
        }
 
        int nb=data.getFieldData().size()/2;  //for one face for nodes only  (total indieces are for both faces)
        field_data_nodes[0].resize(nb); //1st face of prism -- nodes
        field_data_nodes[1].resize(nb); //2nd face of prism -- nodes
 
        for(int ii=0; ii<nb; ii++){
          field_data_nodes[0][ii]=data.getFieldData()[ii];
          field_data_nodes[1][ii]=data.getFieldData()[ii+nb];
        }
//        cout<<"field_data_nodes[0] "<<field_data_nodes[0]<<endl;
//        cout<<"field_data_nodes[1] "<<field_data_nodes[1]<<endl;
      }
 
      for(unsigned int gg = 0; gg<nb_gauss_pts;gg++) {
          if(type == MBVERTEX) {
            ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat, 2); CHKERRQ(ierr);}
          else {
            ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat); CHKERRQ(ierr);}
 
          if(type == MBVERTEX) {
            commonDataPeriodic.DispAtGaussPts[0][gg]=prod(N_mat, field_data_nodes[0]);
            commonDataPeriodic.DispAtGaussPts[1][gg]=prod(N_mat, field_data_nodes[1]);}
          else if(type == MBEDGE && side < 3){
            commonDataPeriodic.DispAtGaussPts[0][gg]+=prod(N_mat, data.getFieldData());}
          else if(type == MBEDGE && side  >= 6){
            commonDataPeriodic.DispAtGaussPts[1][gg]+=prod(N_mat, data.getFieldData());}
          else if(type == MBTRI && side  == 3){
            commonDataPeriodic.DispAtGaussPts[0][gg]+=prod(N_mat, data.getFieldData());}
          else if(type == MBTRI && side  == 4){
            commonDataPeriodic.DispAtGaussPts[1][gg]+=prod(N_mat, data.getFieldData());}
      }
 
    PetscFunctionReturn(0);
    }
  };
 
 
 
 
  struct OpRVEBCsPeriodicCalCU:public FlatPrismElementForcesAndSurcesCore::UserDataOperator {
 
    Vec F;
    RVEBC_Data_Periodic &dAta;
    bool hoGeometry;
    CommonDataPeriodic &commonDataPeriodic;
    CommonFunctionsPeriodic &commonFunctionsPeriodic;
 
    OpRVEBCsPeriodicCalCU(
      const string field_name,
      RVEBC_Data_Periodic &data,
      Vec f,
      CommonDataPeriodic &common_data_periodic,
      CommonFunctionsPeriodic &common_functions_periodic,
      bool ho_geometry = false
    ):
    FlatPrismElementForcesAndSurcesCore::UserDataOperator(
      field_name,UserDataOperator::OPROW
    ),
    F(f),
    dAta(data),
    hoGeometry(ho_geometry),
    commonDataPeriodic(common_data_periodic),
    commonFunctionsPeriodic(common_functions_periodic) {
    }
 
 
    ublas::vector<int> row_ind;
    VectorDouble fe1, fe2;
    MatrixDouble N_mat;
    
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
//      cout<<" OpRVEBCsPeriodicCalCU "<<endl;
      if(data.getIndices().size()==0) PetscFunctionReturn(0);
      if(type == MBQUAD) PetscFunctionReturn(0);
      if(type == MBEDGE && side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
      if(type == MBTRI && side == 4) PetscFunctionReturn(0);  //ignore second triangel face (cananical number 4)
 
      int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
 
      if(type == MBVERTEX) {
        int nb=data.getIndices().size()/2;  //for one face for nodes only  (total indieces are for both faces)
        row_ind.resize(nb); //1st face of prism -- nodes
        for(int ii=0; ii<nb; ii++){
          row_ind[ii]=data.getIndices()[ii];
        }
      }
 
      for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
          double area;
          VectorDouble normal_f3 = getNormalsAtGaussPtsF3(gg);
          area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
          double val = getGaussPts()(2,gg)*area;
 
          if(type == MBVERTEX) {
            ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat, 2); CHKERRQ(ierr);}
          else {
            ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat); CHKERRQ(ierr);}
 
          if(gg==0){
            fe1=-val*prod(trans(N_mat), commonDataPeriodic.DispAtGaussPts[0][gg]);
            fe2= val*prod(trans(N_mat), commonDataPeriodic.DispAtGaussPts[1][gg]);}
          else{
            fe1+=-val*prod(trans(N_mat), commonDataPeriodic.DispAtGaussPts[0][gg]);
            fe2+= val*prod(trans(N_mat), commonDataPeriodic.DispAtGaussPts[1][gg]);}
      }
 
      if(type == MBVERTEX) {
          ierr = VecSetValues(F,row_ind.size(),&row_ind[0],&fe1[0],ADD_VALUES); CHKERRQ(ierr);
          ierr = VecSetValues(F,row_ind.size(),&row_ind[0],&fe2[0],ADD_VALUES); CHKERRQ(ierr);}
      else{
          ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&fe1[0],ADD_VALUES); CHKERRQ(ierr);
          ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&fe2[0],ADD_VALUES); CHKERRQ(ierr);}
 
    PetscFunctionReturn(0);
    }
  };
 
 
 
    struct OpRVEBCsPeriodicCalLagMulAtGaussPts:public FlatPrismElementForcesAndSurcesCore::UserDataOperator {
 
    RVEBC_Data_Periodic &dAta;
    bool hoGeometry;
    CommonDataPeriodic &commonDataPeriodic;
    CommonFunctionsPeriodic &commonFunctionsPeriodic;
 
    OpRVEBCsPeriodicCalLagMulAtGaussPts(
      const string field_name,
      RVEBC_Data_Periodic &data,
      CommonDataPeriodic &common_data_periodic,
      CommonFunctionsPeriodic &common_functions_periodic,
      bool ho_geometry = false
    ):
    FlatPrismElementForcesAndSurcesCore::UserDataOperator(
      field_name,UserDataOperator::OPROW
    ),
    dAta(data),
    hoGeometry(ho_geometry),
    commonDataPeriodic(common_data_periodic),
    commonFunctionsPeriodic(common_functions_periodic) {
    }
 
 
    VectorDouble field_data_nodes;
    MatrixDouble N_mat;
    
 
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
 
 
      if(data.getIndices().size()==0) PetscFunctionReturn(0);
      if(type == MBEDGE && side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
      if(type == MBTRI && side == 4) PetscFunctionReturn(0);  //ignore second triangel face (cananical number 4)
 
      int nb_gauss_pts = data.getN().size1();
      int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
 
       //initialize
      if(type == MBVERTEX) {
        commonDataPeriodic.LagMulAtGaussPts.resize(nb_gauss_pts);
        for(int gg = 0;gg<nb_gauss_pts;gg++) {
          commonDataPeriodic.LagMulAtGaussPts[gg].resize(rank);  commonDataPeriodic.LagMulAtGaussPts[gg].clear();
        }
 
        int nb=data.getFieldData().size()/2;  //for one face for nodes only  (total indieces are for both faces)
        field_data_nodes.resize(nb); //1st face of prism -- nodes
 
        for(int ii=0; ii<nb; ii++){
          field_data_nodes[ii]=data.getFieldData()[ii];
        }
      }
 
 
      for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
          if(type == MBVERTEX) {
            ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat, 2); CHKERRQ(ierr);}
          else {
            ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat); CHKERRQ(ierr);}
 
          if(type == MBVERTEX) {
            commonDataPeriodic.LagMulAtGaussPts[gg]=prod(N_mat, field_data_nodes);}
          else {
            commonDataPeriodic.LagMulAtGaussPts[gg]+=prod(N_mat, data.getFieldData());}
      }
 
    PetscFunctionReturn(0);
    }
  };
 
 
 
    struct OpRVEBCsPeriodicCalCTLam:public FlatPrismElementForcesAndSurcesCore::UserDataOperator {
 
    Vec F;
    RVEBC_Data_Periodic &dAta;
    bool hoGeometry;
    CommonDataPeriodic &commonDataPeriodic;
    CommonFunctionsPeriodic &commonFunctionsPeriodic;
 
    OpRVEBCsPeriodicCalCTLam(
      const string field_name,
      RVEBC_Data_Periodic &data,
      Vec f,
      CommonDataPeriodic &common_data_periodic,
      CommonFunctionsPeriodic &common_functions_periodic,
      bool ho_geometry = false
    ):
    FlatPrismElementForcesAndSurcesCore::UserDataOperator(
      field_name,UserDataOperator::OPCOL
    ),
    F(f),
    dAta(data),
    hoGeometry(ho_geometry),
    commonDataPeriodic(common_data_periodic),
    commonFunctionsPeriodic(common_functions_periodic) {
    }
 
 
    ublas::vector<ublas::vector<int> > col_ind;  //indices for two faces
    VectorDouble fe1, fe2;
    
    MatrixDouble N_mat;
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
//      cout<<" OpRVEBCsPeriodicCalCTLam "<<endl;
      if(data.getIndices().size()==0) PetscFunctionReturn(0);
      int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
      if(type == MBVERTEX) {
        col_ind.resize(2);
        int nb=data.getIndices().size()/2;  //for one face for nodes only  (total indieces are for both faces)
        col_ind[0].resize(nb); //1st face of prism -- nodes
        col_ind[1].resize(nb); //2nd face of prism -- nodes
        for(int ii=0; ii<nb; ii++){
          col_ind[0][ii]=data.getIndices()[ii];
          col_ind[1][ii]=data.getIndices()[ii+nb];
        }
      }
 
 
//      cout<<"data.getN().size1() "<<data.getN().size1()<<endl;
      for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
//         cout<<"gg "<<gg<<endl;
          double area;
          VectorDouble normal_f3 = getNormalsAtGaussPtsF3(gg);
          area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
          double val = getGaussPts()(2,gg)*area;
 
          if(type == MBVERTEX) {
            ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat, 2); CHKERRQ(ierr);}
          else {
            ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat); CHKERRQ(ierr);}
 
 
          if(gg==0){
              if(type == MBVERTEX) {
                fe1=-val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);
                fe2= val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
              else if(type == MBEDGE && side < 3){
                fe1=-val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
              else if(type == MBEDGE && side  >= 6){
                fe2= val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
              else if(type == MBTRI && side  == 3){
                fe1=-val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
              else if(type == MBTRI && side  == 4){
                fe2= val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
          }else{
          if(type == MBVERTEX) {
            fe1+=-val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);
            fe2+= val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
          else if(type == MBEDGE && side < 3){
            fe1+=-val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
          else if(type == MBEDGE && side  >= 6){
            fe2+= val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
          else if(type == MBTRI && side  == 3){
            fe1+=-val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
          else if(type == MBTRI && side  == 4){
            fe2+= val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
        }
 
      }
//      cout<<" Before assembley of F"<<endl;
 
      if(type == MBVERTEX) {
          ierr = VecSetValues(F,col_ind[0].size(),&col_ind[0][0],&fe1[0],ADD_VALUES); CHKERRQ(ierr);
          ierr = VecSetValues(F,col_ind[1].size(),&col_ind[1][0],&fe2[0],ADD_VALUES); CHKERRQ(ierr);}
      else if(type == MBEDGE && side < 3){
          ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&fe1[0],ADD_VALUES); CHKERRQ(ierr);}
      else if(type == MBEDGE && side  >= 6){
          ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&fe2[0],ADD_VALUES); CHKERRQ(ierr);}
      else if(type == MBTRI && side  == 3){
          ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&fe1[0],ADD_VALUES); CHKERRQ(ierr);}
      else if(type == MBTRI && side  == 4){
          ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&fe2[0],ADD_VALUES); CHKERRQ(ierr);}
    PetscFunctionReturn(0);
    }
  };
 
 
 
  //for nonlinear problems
  PetscErrorCode setRVEBCsOperatorsNonlinear(
    string field_name,
    string lagrang_field_name,
    string mesh_nodals_positions,
    Mat aij,vector<Vec> &fvec,Vec f,VectorDouble given_strain
  ) {
    PetscFunctionBegin;
 
    bool ho_geometry = false;
    if(mField.check_field(mesh_nodals_positions)) {
      ho_geometry = true;
    }
 
    map<int,RVEBC_Data_Periodic>::iterator sit = setOfRVEBCPrisms.begin();
    for(;sit!=setOfRVEBCPrisms.end();sit++) {
 
      //Col indices and Col N
      feRVEBCLhs.getOpPtrVector().push_back(
        new OpRVEBCsPeriodicColInd(
          field_name,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
        )
      );
 
 
      //Row indices and Row N
      feRVEBCLhs.getOpPtrVector().push_back(
        new OpRVEBCsPeriodicRowInd(
          lagrang_field_name,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
        )
      );
 
      //Calculate and assemble Cmat
      feRVEBCLhs.getOpPtrVector().push_back(
        new OpRVEBCsPeriodicCalAssemCmat(
          field_name,lagrang_field_name,aij,sit->second,commonFunctionsPeriodic, commonDataPeriodic,ho_geometry
        )
      );
 
      //RHS
      feRVEBCRhs.getOpPtrVector().push_back(
        new OpRVEBCsPeriodicRhs_givenStrain(
          field_name,lagrang_field_name,f,given_strain, methodsOp, sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
        )
      );
 
 
      //Calculate displacement at gauss point
      feRVEBCRhsResidual.getOpPtrVector().push_back(
        new OpRVEBCsPeriodicCalDispAtGaussPts(
          field_name,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
        )
      );
 
      //Calculate and assemble CU = NT * u
      feRVEBCRhsResidual.getOpPtrVector().push_back(
        new OpRVEBCsPeriodicCalCU(
          lagrang_field_name,sit->second,f,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
        )
      );
 
 
      //Calculate lagrange multipliers at gauss point
      feRVEBCRhsResidual.getOpPtrVector().push_back(
        new OpRVEBCsPeriodicCalLagMulAtGaussPts(
          lagrang_field_name,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
        )
      );
 
      //Calculate and assemble CT Lam = NT * Lam
      feRVEBCRhsResidual.getOpPtrVector().push_back(
        new OpRVEBCsPeriodicCalCTLam(
          field_name,sit->second,f,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
        )
      );
 
//      cout<<"Hi 1 from setRVEBCsOperators periodic "<<endl;
//      string aaa;
//      cin>>aaa;
    }
    PetscFunctionReturn(0);
  }
 
 
 
  //for linear problems
  PetscErrorCode setRVEBCsOperators(
    string field_name,string lagrang_field_name,string mesh_nodals_positions,Mat aij,vector<Vec> &f
  ) {
    PetscFunctionBegin;
 
    bool ho_geometry = false;
    if(mField.check_field(mesh_nodals_positions)) {
      ho_geometry = true;
    }
    //      cout<<"Hi 1 from setRVEBCsOperators periodic "<<endl;
    map<int,RVEBC_Data_Periodic>::iterator sit = setOfRVEBCPrisms.begin();
    for(;sit!=setOfRVEBCPrisms.end();sit++) {
 
      //        LHS
      //Col indices and Col N
      feRVEBCLhs.getOpPtrVector().push_back(
        new OpRVEBCsPeriodicColInd(
          field_name,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
        )
      );
 
 
      //Row indices and Row N
      feRVEBCLhs.getOpPtrVector().push_back(
        new OpRVEBCsPeriodicRowInd(
          lagrang_field_name,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
        )
      );
 
      //Calculate and assemble Cmat
      feRVEBCLhs.getOpPtrVector().push_back(
        new OpRVEBCsPeriodicCalAssemCmat(
          field_name,lagrang_field_name,aij,sit->second,commonFunctionsPeriodic, commonDataPeriodic,ho_geometry
        )
      );
 
 
 
      //RHS
      feRVEBCRhs.getOpPtrVector().push_back(
        new OpRVEBCsPeriodicRhs(
          field_name,lagrang_field_name,f,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
        )
      );
 
    }
    PetscFunctionReturn(0);
  }
 
 
  /// \biref operator to calculate the RVE homogenised stress
  struct OpRVEHomoStress: public OpDmatRhs  {
 
    Vec Stress_Homo;
 
    OpRVEHomoStress(
      const string field_name,
      const string lagrang_field_name,
      Vec stress_homo,
      RVEBC_Data_Periodic &data,
      CommonDataPeriodic &common_data_periodic,
      CommonFunctionsPeriodic &common_functions_periodic,
      bool ho_geometry = false
    ):
    OpDmatRhs(
      field_name,lagrang_field_name,data,common_data_periodic,common_functions_periodic,ho_geometry
    ),
    Stress_Homo(stress_homo) {
    }
 
    ublas::vector<VectorDouble > Stress_Homo_elem;
    VectorDouble field_data;
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
 
      if(data.getIndices().size()==0) PetscFunctionReturn(0);
      if(type == MBEDGE && side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
      if(type == MBTRI && side == 4) PetscFunctionReturn(0);  //ignore second triangel face (cananical number 4)
 
      
      ierr = calculateDmat(side,type,data); CHKERRQ(ierr);
 
//      cerr << commonDataPeriodic.D_mat[0] << endl;
//      cerr << data.getFieldData() << endl;
//      cerr << rowFieldName << endl;
//      cerr << data.getFieldDofs()[0]->getName() << endl;
//      cerr << data.getFieldDofs()[0]->getNbOfCoeffs() << endl;
//      cerr << data.getFieldDofs().size() << endl;
//      cerr << data.getIndices() << endl;
//      */
 
 
      Stress_Homo_elem.resize(2);
      if(type == MBVERTEX) {
        int nb=data.getFieldData().size()/2;  //for one face for nodes only  (total indieces are for both faces)
        field_data.resize(nb);
 
        for(int ii=0; ii<nb; ii++){
          field_data[ii]=data.getFieldData()[ii];
        }
        Stress_Homo_elem[0] = prod(trans(commonDataPeriodic.D_mat[0]), -1*field_data);   //Lamda=data.getFieldData() is reaction force (so multiply for -1 to get the force)
        Stress_Homo_elem[1] = prod(trans(commonDataPeriodic.D_mat[1]), -1*field_data);   //Lamda=data.getFieldData() is reaction force (so multiply for -1 to get the force)
      }
      else{
        Stress_Homo_elem[0] = prod(trans(commonDataPeriodic.D_mat[0]), -1*data.getFieldData());   //Lamda=data.getFieldData() is reaction force (so multiply for -1 to get the force)
        Stress_Homo_elem[1] = prod(trans(commonDataPeriodic.D_mat[1]), -1*data.getFieldData());   //Lamda=data.getFieldData() is reaction force (so multiply for -1 to get the force)
      }
 
 
//      cerr << "field_data "<<field_data << endl;
//      cerr << "Stress_Homo_elem[0] "<<Stress_Homo_elem[0] << endl;
//      cerr << "Stress_Homo_elem[1] "<<Stress_Homo_elem[1] << endl;
//      cerr << "commonDataPeriodic.D_mat[0] "<<commonDataPeriodic.D_mat[0] << endl;
 
      int Indices6[6]={0, 1, 2, 3, 4, 5};
      int Indices3[3]={0, 1, 2};
 
      int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
      switch(rank) {
        case 3:
        ierr = VecSetValues(Stress_Homo,6,Indices6,&Stress_Homo_elem[0][0],ADD_VALUES); CHKERRQ(ierr);
        ierr = VecSetValues(Stress_Homo,6,Indices6,&Stress_Homo_elem[1][0],ADD_VALUES); CHKERRQ(ierr);
        break;
        case 1:
        ierr = VecSetValues(Stress_Homo,3,Indices3,&Stress_Homo_elem[0][0],ADD_VALUES); CHKERRQ(ierr);
        ierr = VecSetValues(Stress_Homo,3,Indices3,&Stress_Homo_elem[1][0],ADD_VALUES); CHKERRQ(ierr);
        break;
        default:
        SETERRQ(PETSC_COMM_SELF,1,"not implemented");
      }
 
      PetscFunctionReturn(0);
    }
  };
 
  PetscErrorCode setRVEBCsHomoStressOperators(
    string field_name,
    string lagrang_field_name,
    string mesh_nodals_positions,
    Vec stress_homo
  ) {
    PetscFunctionBegin;
    bool ho_geometry = false;
    if(mField.check_field(mesh_nodals_positions)) {
      ho_geometry = true;
    }
    map<int,RVEBC_Data_Periodic>::iterator sit = setOfRVEBCPrisms.begin();
    for(;sit!=setOfRVEBCPrisms.end();sit++) {
      feRVEBCStress.getOpPtrVector().push_back(
        new OpRVEHomoStress(
          field_name,lagrang_field_name,stress_homo,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
        )
      );
    }
    PetscFunctionReturn(0);
  }
 
};
 
#endif