BCs_RVELagrange_Trac.hpp

Go to the documentation of this file.

/* Copyright (C) 2014, Zahur Ullah (Zahur.Ullah AT glasgow.ac.uk)
 * --------------------------------------------------------------
 */
 
/* 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_TRAC_HPP
#define __BCS_RVELAGRANGE_TRAC_HPP
 
struct BCs_RVELagrange_Trac: public BCs_RVELagrange_Disp {
 
  BCs_RVELagrange_Trac(MoFEM::Interface &m_field): BCs_RVELagrange_Disp(m_field) {}
 
  struct CommonData {
    MatrixDouble D_mat;
  };
  CommonData commonData;
 
  struct CommonFunctions {
 
    PetscErrorCode shapeMat(
      int rank,unsigned int gg,DataForcesAndSourcesCore::EntData &col_data,MatrixDouble &n
    ) {
      PetscFunctionBegin;
 
      int shape_size = col_data.getN().size2();
 
      n.resize(rank,shape_size*rank);
      n.clear();
 
      int kk = 0;
      for(int ii=0; ii<shape_size; ii++){
        double val = col_data.getN()(gg,ii);
        for(int jj=0; jj<rank; jj++){
          n(jj,kk) = val;
          kk++;
        }
      }
 
      PetscFunctionReturn(0);
    }
 
 
    MatrixDouble H_mat_1Node;
 
    PetscErrorCode hMat(
      VectorDouble face_normal,
      int rank,
      MatrixDouble &n_mat,
      MatrixDouble &h_mat
    ) {
      PetscFunctionBegin;
 
      switch(rank) {
        //Mechanical problem
        case 3: {
          H_mat_1Node.resize(6,3);
          H_mat_1Node.clear(); //for one node
          if(face_normal[0]>0) {     //+X face of the RVE
            H_mat_1Node(0,0) = 1.0;
            H_mat_1Node(3,1) = 1.0;
            H_mat_1Node(5,2) = 1.0;
          } else if(face_normal[0]<0) {    //-X face of the RVE
            H_mat_1Node(0,0) = -1.0;
            H_mat_1Node(3,1) = -1.0;
            H_mat_1Node(5,2) = -1.0;
          } else if(face_normal[1]>0) {     //+Y face of the RVE
            H_mat_1Node(1,1) = 1.0;
            H_mat_1Node(3,0) = 1.0;
            H_mat_1Node(4,2) = 1.0;
          } else if(face_normal[1]<0) {    //-Y face of the RVE
            H_mat_1Node(1,1) = -1.0;
            H_mat_1Node(3,0) = -1.0;
            H_mat_1Node(4,2) = -1.0;
          } else if(face_normal[2]>0) {    //+Z face of the RVE
            H_mat_1Node(2,2) = 1.0;
            H_mat_1Node(4,1) = 1.0;
            H_mat_1Node(5,0) = 1.0;
          } else if(face_normal[2]<0) {    //-Z face of the RVE
            H_mat_1Node(2,2) = -1.0;
            H_mat_1Node(4,1) = -1.0;
            H_mat_1Node(5,0) = -1.0;
          } else {
            SETERRQ(PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"Can not be ?!");
          }
 
          int num_col = n_mat.size2();
          h_mat.resize(6,num_col);
 
          int cc1 = 0;
          for(int bb = 0; bb<num_col/3; bb++) {
            //blocks of 6x3
            for(int rr = 0; rr<6; rr++) {
              for(int cc = 0; cc<3; cc++) {
                h_mat(rr,(cc+cc1)) = H_mat_1Node(rr,cc);
              }
            }
            cc1+=3;
          }
 
        }
        break;
        case 1: //Moisture transport or thermal problem
        {
          H_mat_1Node.resize(3,1);
          H_mat_1Node.clear();
          if(face_normal[0]>0) {     //+X face of the RVE
            H_mat_1Node(0,0) = 1.0;
          }
          if(face_normal[0]<0) {    //-X face of the RVE
            H_mat_1Node(0,0) = -1.0;
          }
          if(face_normal[1]>0) {     //+Y face of the RVE
            H_mat_1Node(1,0) = 1.0;
          }
          if(face_normal[1]<0) {    //-Y face of the RVE
            H_mat_1Node(1,0) = -1.0;
          }
          if(face_normal[2]>0) {    //+Z face of the RVE
            H_mat_1Node(2,0) = 1.0;
          }
          if(face_normal[2]<0) {    //-Z face of the RVE
            H_mat_1Node(2,0) = -1.0;
          }
 
          int num_col = n_mat.size2();
          h_mat.resize(3,num_col);
 
          int cc1 = 0;
          for(int bb = 0; bb<num_col; bb++) {
            //blocks of 3x1
            for(int rr = 0; rr<3; rr++) {
              for(int cc = 0; cc<1; cc++) {
                h_mat(rr,(cc+cc1))=H_mat_1Node(rr,cc);
              }
            }
            cc1+=1;
          }
 
        }
        break;
        default:
        SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"not implemented");
      }
      PetscFunctionReturn(0);
    }
  };
 
  CommonFunctions commonFunctions;
 
  /// \biref operator to calculate the LHS for the RVE bounary conditions
  struct OpRVEBCsLhs:public FaceElementForcesAndSourcesCore::UserDataOperator  {
 
    RVEBC_Data &dAta;
    CommonFunctions &commonFunctions;
    bool hoGeometry;
    Mat Aij;
 
    OpRVEBCsLhs(const string field_name, const string lagrang_field_name,Mat aij,
      RVEBC_Data &data, CommonFunctions &common_functions, bool ho_geometry = false
    ):
    FaceElementForcesAndSourcesCore::UserDataOperator(lagrang_field_name, field_name, UserDataOperator::OPROWCOL),
    dAta(data),commonFunctions(common_functions),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, NTN;
    ublas::vector<MatrixDouble > N_mat_row;
    MatrixDouble N_mat_col,H_mat;
 
    PetscErrorCode doWork(
      int row_side,int col_side,
      EntityType row_type,EntityType col_type,
      DataForcesAndSourcesCore::EntData &row_data,
      DataForcesAndSourcesCore::EntData &col_data
    ) {
      PetscFunctionBegin;
 
      try {
        if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
        if(col_data.getIndices().size()==0) PetscFunctionReturn(0);
 
        int nb_row = row_data.getIndices().size();
        int nb_col = col_data.getIndices().size();
 
        auto weak_ptr_dof =
            getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(
                col_data.getIndices()[0]);
        const FENumeredDofEntity *dof_ptr;
        if (auto ptr = weak_ptr_dof.lock())
          dof_ptr = ptr.get();
        else
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Dof not found");
 
        int rank = dof_ptr->getNbOfCoeffs();
 
        K.resize(nb_row,nb_col);
        K.clear();
 
        N_mat_row.resize(col_data.getN().size1());
        for(unsigned int gg = 0;gg<col_data.getN().size1();gg++) {
          double area;
          if(hoGeometry) {
            area = norm_2(getNormalsAtGaussPts(gg))*0.5;
          }   else {
            area = getArea();
          }
          double val = getGaussPts()(2,gg)*area;
 
          if(col_type==MBVERTEX) {
            ierr = commonFunctions.shapeMat(rank,gg,col_data,N_mat_row[gg]); CHKERRQ(ierr);
          }
          ierr = commonFunctions.shapeMat(rank,gg,col_data,N_mat_col); CHKERRQ(ierr);
          ierr = commonFunctions.hMat(getNormalsAtGaussPts(gg),rank,N_mat_row(gg),H_mat); CHKERRQ(ierr);
 
          if(gg==0){
            NTN = prod(trans(N_mat_row(gg)),N_mat_col);  //we don't need to define size NTN
            K = val*prod(H_mat,NTN); //we don't need to defiend size of K
          } else {
            noalias(NTN) = prod(trans(N_mat_row(gg)),N_mat_col);
            K += val*prod(H_mat, NTN);
          }
        }
 
        // Matrix C
        int nb_rows=row_data.getIndices().size();
        int nb_cols=col_data.getIndices().size();
        ierr = MatSetValues(
          Aij,nb_rows,&row_data.getIndices()[0],nb_cols,&col_data.getIndices()[0],&K(0,0),ADD_VALUES
        ); CHKERRQ(ierr);
 
        // Matrix CT
        transK = trans(K);
        ierr = MatSetValues(
          Aij,nb_cols,&col_data.getIndices()[0],nb_rows,&row_data.getIndices()[0],&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,MOFEM_NOT_IMPLEMENTED,ss.str().c_str());
      }
      PetscFunctionReturn(0);
    }
 
  };
 
  /// \biref operator to calculate the RHS of the constrain for the RVE boundary conditions
  struct OpRVEBCsRhs_Cal:public FaceElementForcesAndSourcesCore::UserDataOperator {
 
    RVEBC_Data &dAta;
    bool hoGeometry;
    CommonData &commonData;
    CommonFunctions &commonFunctions;
 
    OpRVEBCsRhs_Cal(
      const string field_name, RVEBC_Data &data,CommonData &common_data,CommonFunctions &common_functions,bool ho_geometry = false
    ):
    FaceElementForcesAndSourcesCore::UserDataOperator(field_name, UserDataOperator::OPCOL),
    dAta(data),
    hoGeometry(ho_geometry),
    commonData(common_data),
    commonFunctions(common_functions) {
    }
 
    MatrixDouble X_mat, N_mat, NTX, H_mat;
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
 
      if(type!=MBVERTEX) PetscFunctionReturn(0);
      if(data.getIndices().size()==0) PetscFunctionReturn(0);
 
      auto weak_ptr_dof =
          getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(
              data.getIndices()[0]);
      const FENumeredDofEntity *dof_ptr;
      if (auto ptr = weak_ptr_dof.lock())
        dof_ptr = ptr.get();
      else
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Dof not found");
 
      int rank = dof_ptr->getNbOfCoeffs();
 
      double x,y,z;
      for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
        double area;
        if(hoGeometry) {
          area = norm_2(getNormalsAtGaussPts(gg))*0.5;
        } else {
          area = getArea();
        }
        double val = getGaussPts()(2,gg)*area;
 
        x = getCoordsAtGaussPts()(gg,0);
        y = getCoordsAtGaussPts()(gg,1);
        z = getCoordsAtGaussPts()(gg,2);
 
        switch(rank) {
          case 3: //mech problem
          X_mat.resize(rank,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(rank,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,MOFEM_NOT_IMPLEMENTED,"not implemented");
        }
        ierr = commonFunctions.shapeMat(rank, gg,  data, N_mat); CHKERRQ(ierr);
        ierr = commonFunctions.hMat(getNormalsAtGaussPts(gg), rank, N_mat, H_mat); CHKERRQ(ierr);
 
        if(gg==0) {
          NTX = val*prod(trans(N_mat),X_mat);
          commonData.D_mat = prod(H_mat, NTX);
        } else {
          noalias(NTX) = val*prod(trans(N_mat),X_mat);
          commonData.D_mat += prod(H_mat, NTX);
        }
      }
 
      PetscFunctionReturn(0);
    }
  };
 
  /// \biref operator to calculate the RHS of the constrain for the RVE boundary conditions
  struct OpRVEBCsRhs_Assemble: public FaceElementForcesAndSourcesCore::UserDataOperator {
 
    RVEBC_Data &dAta;
    CommonData &commonData;
    vector<Vec> &F;
 
    OpRVEBCsRhs_Assemble(const string lagrang_field_name,vector<Vec> &f,RVEBC_Data &data,CommonData &common_data):
    FaceElementForcesAndSourcesCore::UserDataOperator(lagrang_field_name,UserDataOperator::OPROW),
    dAta(data),commonData(common_data),F(f) {
    }
 
    
    VectorDouble applied_strain,f;
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
 
      if(data.getIndices().size()==0) PetscFunctionReturn(0);
 
      MatrixDouble &D_mat=commonData.D_mat;
      applied_strain.resize(D_mat.size2(),false);
 
      for(int ii = 0;ii<F.size();ii++) {
        applied_strain.clear();
        applied_strain(ii) = 1.0;
        f.resize(D_mat.size1(),false);
        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);
    }
  };
 
  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;
    }
    map<int,RVEBC_Data>::iterator sit = setOfRVEBC.begin();
    for(;sit!=setOfRVEBC.end();sit++) {
      //LHS
      feRVEBCLhs.getOpPtrVector().push_back(
        new OpRVEBCsLhs(field_name,lagrang_field_name,aij,sit->second,commonFunctions,ho_geometry)
      );
      //RHS
      feRVEBCRhs.getOpPtrVector().push_back(
        new OpRVEBCsRhs_Cal(field_name,sit->second,commonData,commonFunctions,ho_geometry)
      );
      feRVEBCRhs.getOpPtrVector().push_back(
        new OpRVEBCsRhs_Assemble(lagrang_field_name,f,sit->second,commonData)
      );
    }
    PetscFunctionReturn(0);
  }
 
  struct OpRVEHomoStress_Assemble:public FaceElementForcesAndSourcesCore::UserDataOperator {
 
    RVEBC_Data &dAta;
    CommonData &commonData;
    Vec stressHomo;
 
    OpRVEHomoStress_Assemble(
      const string lagrang_field_name,Vec stress_homo,RVEBC_Data &data,CommonData &common_data
    ):
    FaceElementForcesAndSourcesCore::UserDataOperator(lagrang_field_name,UserDataOperator::OPROW),
    dAta(data),commonData(common_data),stressHomo(stress_homo) {
    }
 
    
    VectorDouble stressVector;
 
    PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
      PetscFunctionBegin;
      if(data.getIndices().size()==0) PetscFunctionReturn(0);
      if(dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt())==dAta.tRis.end()) PetscFunctionReturn(0);
      stressVector.resize(commonData.D_mat.size2());
      //Lamda=data.getFieldData() is reaction force (so multiply for -1 to get the force)
      noalias(stressVector) = prod(trans(commonData.D_mat),-data.getFieldData());
      const int indices_6[6] = {0, 1, 2, 3, 4, 5};
      const int indices_3[3] = {0, 1, 2};
      switch(stressVector.size()) {
        case 6:
          ierr = VecSetValues(stressHomo,6,indices_6,&stressVector[0],ADD_VALUES); CHKERRQ(ierr);
          break;
        case 3:
          ierr = VecSetValues(stressHomo,3,indices_3,&stressVector[0],ADD_VALUES); CHKERRQ(ierr);
          break;
        default:
          SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"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 OpRVEBCsRhs_Cal(field_name,sit->second,commonData,commonFunctions)
      );
      feRVEBCStress.getOpPtrVector().push_back(
        new OpRVEHomoStress_Assemble(lagrang_field_name,stress_homo,sit->second,commonData)
      );
    }
    PetscFunctionReturn(0);
  }
 
};
 
#endif