BCs_RVELagrange_Disp.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_DISP_HPP
#define __BCS_RVELAGRANGE_DISP_HPP
 
struct BCs_RVELagrange_Disp {
 
  struct MyTriFE: public FaceElementForcesAndSourcesCore {
    MyTriFE(MoFEM::Interface &_mField): FaceElementForcesAndSourcesCore(_mField) {}
    int getRule(int order) { return order; };
  };
  
  boost::ptr_vector<MethodForForceScaling> methodsOp;  //This is to apply strain vector in steps (for given steps) for nonlinear problem
 
  
  MoFEM::Interface &mField;
  MyTriFE feRVEBCRhs; //To calculate the Rhs or RVE BCs
  MyTriFE feRVEBCLhs; //To calculate the Lhs or RVE BCs
  MyTriFE feRVEBCStress; //To calculate the Rhs or RVE BCs
  MyTriFE feRVEBCRhsResidual; //To calculate the Rhs or RVE BCs
  MyTriFE feRVEBCRhsHomoC; //To be used for the calculation of homgenised stiffness matrix
 
  MyTriFE& getLoopFeRVEBCLhs() { return feRVEBCLhs; }
  MyTriFE& getLoopFeRVEBCRhs() { return feRVEBCRhs; }
  
  MyTriFE& getLoopFeRVEBCRhsResidual() { return feRVEBCRhsResidual; }
  MyTriFE& getLoopFeRVEBCStress() { return feRVEBCStress; }
  
  MyTriFE& getLoopFeRVEBCRhsHomoC() { return feRVEBCRhsHomoC; }
 
 
  BCs_RVELagrange_Disp(MoFEM::Interface &m_field):
  mField(m_field),
  feRVEBCRhs(m_field),
  feRVEBCLhs(m_field),
  feRVEBCRhsResidual(m_field),
  feRVEBCStress(m_field),
  feRVEBCRhsHomoC(m_field){
  }
 
  struct RVEBC_Data {
    Range tRis; // All boundary surfaces
  };
  map<int,RVEBC_Data> setOfRVEBC; ///< maps side set id with appropriate FluxData
 
  PetscErrorCode addLagrangiangElement(
    const string element_name,
    const string field_name,
    const string lagrang_field_name,
    const string mesh_nodals_positions
  ) {
    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);
    }
 
    //this is alternative method for setting boundary conditions, to bypass bu in cubit file reader.
    //not elegant, but good enough
    for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField,SIDESET,it)) {
      if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
        //          cout<<"Hi from addLagrangiangElement "<<endl;
        rval = mField.get_moab().get_entities_by_type(it->meshset,MBTRI,setOfRVEBC[it->getMeshsetId()].tRis,true); CHKERRQ_MOAB(rval);
        ierr = mField.add_ents_to_finite_element_by_type(setOfRVEBC[it->getMeshsetId()].tRis,MBTRI,element_name); CHKERRQ(ierr);
      }
    }
 
    PetscFunctionReturn(0);
  }
 
  /// \biref operator to calculate the LHS for the RVE bounary conditions
  struct OpRVEBCsLhs:public FaceElementForcesAndSourcesCore::UserDataOperator {
 
    RVEBC_Data &dAta;
    bool hoGeometry;
    Mat Aij;
    OpRVEBCsLhs(
      const string field_name, const string lagrang_field_name, Mat aij,
      RVEBC_Data &data,bool ho_geometry = false
    ):
    FaceElementForcesAndSourcesCore::UserDataOperator(
      lagrang_field_name, field_name, UserDataOperator::OPROWCOL
    ),
    dAta(data),hoGeometry(ho_geometry),Aij(aij) {
      // This will make sure to loop over all entities
      // (e.g. for order=2 it will make doWork to loop 16 time)
      sYmm = false;
    }
 
    MatrixDouble K,transK;
 
    /** \brief calculate thermal convection term in the lhs of equations
    *
    * C = intS N^T  N dS
    */
    PetscErrorCode doWork(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSurcesCore::EntData &row_data,
      DataForcesAndSurcesCore::EntData &col_data
    ) {
      PetscFunctionBegin;
 
      try {
        int nb_row = row_data.getIndices().size();
        int nb_col = col_data.getIndices().size();
 
        if(nb_row == 0) PetscFunctionReturn(0);
        if(nb_col == 0) PetscFunctionReturn(0);
 
        
        const FENumeredDofEntity *dof_ptr;
        ierr = getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(row_data.getIndices()[0],&dof_ptr); CHKERRQ(ierr);
        int rank = dof_ptr->getNbOfCoeffs();
 
        K.resize(nb_row/rank,nb_col/rank);
        K.clear();
 
        for(unsigned int gg = 0;gg<row_data.getN().size1();gg++) {
          double area;
          if(hoGeometry) {
            area = norm_2(getNormalsAtGaussPt(gg))*0.5;
          }   else {
            area = getArea();
          }
          double val = getGaussPts()(2,gg)*area;
          noalias(K) += val*outer_prod(
            row_data.getN(gg,nb_row/rank),col_data.getN(gg,nb_col/rank)
          );
        }
 
        VectorInt row_indices,col_indices;
        row_indices.resize(nb_row/rank);
        col_indices.resize(nb_col/rank);
 
        for(int rr = 0;rr < rank; rr++) {
          unsigned int nb_rows;
          unsigned int nb_cols;
          int *rows;
          int *cols;
          if(rank > 1) {
 
            ublas::noalias(row_indices) = ublas::vector_slice<VectorInt >
            (row_data.getIndices(), ublas::slice(rr, rank, row_data.getIndices().size()/rank));
 
            ublas::noalias(col_indices) = ublas::vector_slice<VectorInt >
            (col_data.getIndices(), ublas::slice(rr, rank, col_data.getIndices().size()/rank));
 
            nb_rows = row_indices.size();
            nb_cols = col_indices.size();
            rows = &*row_indices.data().begin();
            cols = &*col_indices.data().begin();
 
          } else {
 
            nb_rows = row_data.getIndices().size();
            nb_cols = col_data.getIndices().size();
            rows = &*row_data.getIndices().data().begin();
            cols = &*col_data.getIndices().data().begin();
 
          }
 
          // Matrix C
          ierr = MatSetValues(Aij,nb_rows,rows,nb_cols,cols,&K(0,0),ADD_VALUES); CHKERRQ(ierr);
 
          // Matrix CT
          transK.resize(nb_col/rank,nb_row/rank);
          noalias(transK) = trans(K);
          ierr = MatSetValues(Aij,nb_cols,cols,nb_rows,rows,&transK(0,0),ADD_VALUES); CHKERRQ(ierr);
 
        }
 
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
 
      PetscFunctionReturn(0);
    }
  };
 
  struct OpDmatRhs: public FaceElementForcesAndSourcesCore::UserDataOperator {
 
    RVEBC_Data &dAta;
    bool hoGeometry;
 
    OpDmatRhs(
      const string field_name, const string lagrang_field_name,
      RVEBC_Data &data,bool ho_geometry = false
    ):
    FaceElementForcesAndSourcesCore::UserDataOperator(
      lagrang_field_name,UserDataOperator::OPROW
    ),
    dAta(data),hoGeometry(ho_geometry) {
    }
 
    VectorDouble f,applied_strain;
    MatrixDouble X_mat,N_mat,D_mat;
 
    int rank,nb_row_dofs;
 
    PetscErrorCode calculateDmat(int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
      PetscFunctionBegin;
 
      double x,y,z;
      for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
 
        double area;
        if(hoGeometry) {
          area = norm_2(getNormalsAtGaussPt(gg))*0.5;
        } else {
          area = getArea();
        }
        double val = getGaussPts()(2,gg)*area;
 
        x = getHoCoordsAtGaussPts()(gg,0);
        y = getHoCoordsAtGaussPts()(gg,1);
        z = getHoCoordsAtGaussPts()(gg,2);
 
        switch(rank) {
          case 3: //mech problem
          X_mat.resize(3,6,false);
          X_mat.clear();
          X_mat(0,0) = 2.0*x;  X_mat(0,3)=y;  X_mat(0,5)=z;
          X_mat(1,1) = 2.0*y;  X_mat(1,3)=x;  X_mat(1,4)=z;
          X_mat(2,2) = 2.0*z;  X_mat(2,4)=y;  X_mat(2,5)=x;
          X_mat = 0.5*X_mat;
          break;
          case 1:  //moisture transport or thermal problem
          X_mat.resize(3,3,false);
          X_mat.clear();
          X_mat(0,0)= x;  X_mat(0,1)= y;  X_mat(0,2) = z;
          break;
          default:
          SETERRQ(PETSC_COMM_SELF,1,"not implemented");
        }
 
        int shape_size = data.getN().size2();
        N_mat.resize(rank,shape_size*rank);
        N_mat.clear();
 
        // FIXME: this is inefficient implementation
        {
          int kk=0;
          for(int ii=0; ii<shape_size; ii++){
            //number of shape functions
            double val = data.getN()(gg,ii);
            for(int jj=0; jj<rank; jj++) {
              N_mat(jj,kk) = val;
              kk++;
            }
          }
          if(gg==0) {
            D_mat = val*prod(trans(N_mat),X_mat);
          } else{
            D_mat += val*prod(trans(N_mat),X_mat);
          }
        }
 
      }
 
      PetscFunctionReturn(0);
    }
 
  };
  
  
  /// \biref operator to calculate the RHS of the constrain for the RVE boundary conditions
  struct OpRVEBCsRhs:public OpDmatRhs {
    Vec F;
    VectorDouble givenStrain;
    boost::ptr_vector<MethodForForceScaling> &methodsOp;
    
    OpRVEBCsRhs(
      const string field_name, const string lagrang_field_name,
      Vec f, VectorDouble given_strain, boost::ptr_vector<MethodForForceScaling> &methods_op, RVEBC_Data &data,bool ho_geometry = false
    ):
    OpDmatRhs(field_name,lagrang_field_name,data,ho_geometry),F(f),givenStrain(given_strain),methodsOp(methods_op) {
    }
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
      PetscFunctionBegin;
 
      nb_row_dofs = data.getIndices().size();
 
      if(nb_row_dofs ==0) PetscFunctionReturn(0);
      if(dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt())==dAta.tRis.end()) PetscFunctionReturn(0);
 
      
      const FENumeredDofEntity *dof_ptr;
      ierr = getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(data.getIndices()[0],&dof_ptr); CHKERRQ(ierr);
      rank = dof_ptr->getNbOfCoeffs();
      nb_row_dofs /= rank;
 
      ierr = calculateDmat(side,type,data); CHKERRQ(ierr);
      
      //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;
      f.resize(D_mat.size1(),false);
      noalias(f) = prod(D_mat, scaled_given_strain);
      f*=-1;
 
//      cout <<"f   = "<<f<<endl;
//      string wait;
//      cin>>wait;
 
      Vec myF = F;
      if(F == PETSC_NULL) {
        myF = getFEMethod()->snes_f;
      }
 
      ierr = VecSetValues(myF,data.getIndices().size(),&data.getIndices()[0],&f[0],ADD_VALUES); CHKERRQ(ierr);
      PetscFunctionReturn(0);
    }
  };
  
  /// \biref operator to calculate the RHS for the calculation of the homgoensied stiffness matrix
  struct OpRVEBCsRhsHomoC:public OpDmatRhs {
    
    vector<Vec> &F;
    
    OpRVEBCsRhsHomoC(
                const string field_name, const string lagrang_field_name,
                vector<Vec> &f,RVEBC_Data &data,bool ho_geometry = false
                ):
    OpDmatRhs(field_name,lagrang_field_name,data,ho_geometry),F(f) {
    }
    
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
      PetscFunctionBegin;
      
      nb_row_dofs = data.getIndices().size();
      
      if(nb_row_dofs ==0) PetscFunctionReturn(0);
      if(dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt())==dAta.tRis.end()) PetscFunctionReturn(0);
      
      
      const FENumeredDofEntity *dof_ptr;
      ierr = getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(data.getIndices()[0],&dof_ptr); CHKERRQ(ierr);
      rank = dof_ptr->getNbOfCoeffs();
      nb_row_dofs /= rank;
      
      ierr = calculateDmat(side,type,data); CHKERRQ(ierr);
      
      applied_strain.resize(D_mat.size2(),false);
      f.resize(D_mat.size1(),false);
      
      int size = (rank == 3) ? 6 : 3;
      for(int ii = 0;ii<size;ii++) {
        applied_strain.clear();
        applied_strain[ii] = 1.0;
        noalias(f) = prod(D_mat, applied_strain);
        ierr = VecSetValues(
                            F[ii],data.getIndices().size(),&data.getIndices()[0],&f[0],ADD_VALUES
                            ); CHKERRQ(ierr);
      }
      
      PetscFunctionReturn(0);
      
    }
  };
 
  struct CommonFunctions {
    PetscErrorCode shapeMat(int rank, unsigned int gg, DataForcesAndSurcesCore::EntData &col_data, MatrixDouble &N_mat) {
      PetscFunctionBegin;
      int shape_size=col_data.getN().size2();
      //      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 funcitons
        for(int jj=0; jj<rank; jj++){
          //              cout<<"ii "<<ii<<endl;
          N_mat(jj,kk)=col_data.getN()(gg,gg1); kk++;
        }
        gg1++;
      }
      PetscFunctionReturn(0);
    }
  };
  CommonFunctions common_functions;
  
  
  /// \biref operator to calculate and assemble CU
  struct OpRVEBCsRhsForceExternal_CU:public FaceElementForcesAndSourcesCore::UserDataOperator {
 
    bool hoGeometry;
    Vec F;
    RVEBC_Data &dAta;
    CommonFunctions commonFunctions;
 
 
    OpRVEBCsRhsForceExternal_CU(
      const string field_name, const string lagrang_field_name, Vec f,
      CommonFunctions &_common_functions,  RVEBC_Data &data, bool ho_geometry = false
    ):
    FaceElementForcesAndSourcesCore::UserDataOperator(
      lagrang_field_name, field_name, UserDataOperator::OPROWCOL
    ),
    dAta(data),hoGeometry(ho_geometry),F(f), commonFunctions(_common_functions){
      // This will make sure to loop over all entities
      // (e.g. for order=2 it will make doWork to loop 16 time)
      sYmm = false;
    }
 
    MatrixDouble K,transK, N_mat_row, N_mat_col;
    VectorDouble CU, CTLam;
 
    PetscErrorCode doWork(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSurcesCore::EntData &row_data,
      DataForcesAndSurcesCore::EntData &col_data
    ) {
      PetscFunctionBegin;
 
      try {
        int nb_row = row_data.getIndices().size();
        int nb_col = col_data.getIndices().size();
 
        if(nb_row == 0) PetscFunctionReturn(0);
        if(nb_col == 0) PetscFunctionReturn(0);
 
 
        
        const FENumeredDofEntity *dof_ptr;
        ierr = getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(row_data.getIndices()[0],&dof_ptr); CHKERRQ(ierr);
        int rank = dof_ptr->getNbOfCoeffs();
        for(unsigned int gg = 0;gg<row_data.getN().size1();gg++) {
          double area;
          if(hoGeometry) {
            area = norm_2(getNormalsAtGaussPt(gg))*0.5;
          }   else {
            area = getArea();
          }
          double val = getGaussPts()(2,gg)*area;
          ierr = commonFunctions.shapeMat(rank, gg,  row_data, N_mat_row); CHKERRQ(ierr);
          ierr = commonFunctions.shapeMat(rank, gg,  col_data, N_mat_col); CHKERRQ(ierr);
 
//        cout<<"row_data.getN() "<<row_data.getN()<<endl;
//        cout<<"N_mat_row() "<<N_mat_row<<endl;
//
//        cout<<"col_data.getN() "<<col_data.getN()<<endl;
//        cout<<"N_mat_col "<<N_mat_col<<endl;
 
          if(gg==0){
            K =val* prod(trans(N_mat_row),N_mat_col);  //we don't need to define size K
          } else {
            K += val*prod(trans(N_mat_row),N_mat_col);
          }
        }
        CU=prod(K,       col_data.getFieldData());
        CTLam=prod(trans(K),row_data.getFieldData());
 
 
        Vec myF = F;
        if(F == PETSC_NULL) {
          myF = getFEMethod()->snes_f;
        }
        
        ierr = VecSetValues(myF,row_data.getIndices().size(),&row_data.getIndices()[0],&CU[0],ADD_VALUES); CHKERRQ(ierr);
        ierr = VecSetValues(myF,col_data.getIndices().size(),&col_data.getIndices()[0],&CTLam[0],ADD_VALUES); CHKERRQ(ierr);
 
      } catch (const std::exception& ex) {
        ostringstream ss;
        ss << "throw in method: " << ex.what() << endl;
        SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
      }
 
      PetscFunctionReturn(0);
    }
  };
  
  
  //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>::iterator sit = setOfRVEBC.begin();
    for(;sit!=setOfRVEBC.end();sit++) {
      
      //Lhs (calculate LHS)
      feRVEBCLhs.getOpPtrVector().push_back(
        new OpRVEBCsLhs(field_name,lagrang_field_name, aij, sit->second,ho_geometry)
      );
      
      
      //Rhs (calculate RHS)
      feRVEBCRhs.getOpPtrVector().push_back(
        new OpRVEBCsRhs(field_name,lagrang_field_name,f,given_strain, methodsOp, sit->second,ho_geometry)
      );
      
      
      //Rhs (F[6] used to calculate Homgenised stiffness matrix)
      feRVEBCRhsHomoC.getOpPtrVector().push_back(
         new OpRVEBCsRhsHomoC(field_name,lagrang_field_name,fvec,sit->second,ho_geometry)
       );
 
      //for internal forces (i.e. RHS)
      //calculate and assemble CU and CT Lam
      feRVEBCRhsResidual.getOpPtrVector().push_back(
        new OpRVEBCsRhsForceExternal_CU(field_name,lagrang_field_name,f,common_functions,sit->second,ho_geometry)
      );
      
      
     }
 
    PetscFunctionReturn(0);
  }
  
  //for linear problems
  PetscErrorCode setRVEBCsOperators(
      string field_name,
      string lagrang_field_name,
      string mesh_nodals_positions,
      Mat aij,vector<Vec> &fvec
      ) {
    PetscFunctionBegin;
    
    bool ho_geometry = false;
    if(mField.check_field(mesh_nodals_positions)) {
      ho_geometry = true;
    }
    map<int,RVEBC_Data>::iterator sit = setOfRVEBC.begin();
    for(;sit!=setOfRVEBC.end();sit++) {
      
      //Lhs (calculate LHS)
      feRVEBCLhs.getOpPtrVector().push_back(
      new OpRVEBCsLhs(field_name,lagrang_field_name, aij, sit->second,ho_geometry)
      );
 
      //Rhs (F[6] used to calculate Homgenised stiffness matrix)
      feRVEBCRhs.getOpPtrVector().push_back(
       new OpRVEBCsRhsHomoC(field_name,lagrang_field_name,fvec,sit->second,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 &data,
                    bool ho_geometry = false
                    ):
    OpDmatRhs(field_name,lagrang_field_name,data,ho_geometry),
    Stress_Homo(stress_homo) {
    }
    
    VectorDouble Stress_Homo_elem;
    
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSurcesCore::EntData &data) {
      PetscFunctionBegin;
      if(data.getIndices().size()==0) PetscFunctionReturn(0);
      if(dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt())==dAta.tRis.end()) PetscFunctionReturn(0);
      //        cout<<"OpRVEBCsRhs "<<endl;
      
      
      const FENumeredDofEntity *dof_ptr;
      ierr = getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(data.getIndices()[0],&dof_ptr); CHKERRQ(ierr);
      rank = dof_ptr->getNbOfCoeffs();
      nb_row_dofs = data.getIndices().size()/rank;
      
      ierr = calculateDmat(side,type,data); CHKERRQ(ierr);
      
      Stress_Homo_elem.resize(D_mat.size2());
      
      // data.getFieldData() is value of Lagrange multply
      // data.getFieldData() is reaction force (so multiply for -1 to get the force)
      noalias(Stress_Homo_elem) = prod(trans(D_mat), -data.getFieldData());
      
      const int Indices6[6] = {0, 1, 2, 3, 4, 5};
      const int Indices3[3] = {0, 1, 2};
      switch(rank) {
        case 3:
          ierr = VecSetValues(Stress_Homo,6,Indices6,&(Stress_Homo_elem.data())[0],ADD_VALUES); CHKERRQ(ierr);
          break;
        case 1:
          ierr = VecSetValues(Stress_Homo,3,Indices3,&(Stress_Homo_elem.data())[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>::iterator sit = setOfRVEBC.begin();
    for(;sit!=setOfRVEBC.end();sit++) {
      feRVEBCStress.getOpPtrVector().push_back(new OpRVEHomoStress(field_name, lagrang_field_name, Stress_Homo, sit->second, ho_geometry));
    }
    
    PetscFunctionReturn(0);
  }
};
 
#endif // __BCS_RVELAGRANGE_DISP_HPP