BCs_RVELagrange_Trac_Rigid_Trans.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_TRAC_RIGID_TRANS_HPP
#define __BCS_RVELAGRANGE_TRAC_RIGID_TRANS_HPP
 
struct BCs_RVELagrange_Trac_Rigid_Trans: public BCs_RVELagrange_Disp {
 
  BCs_RVELagrange_Trac_Rigid_Trans(MoFEM::Interface &m_field):
  BCs_RVELagrange_Disp(m_field) {}
 
  /// \biref operator to calculate the LHS for the RVE bounary conditions
  struct OpRVEBCsLhs:public FaceElementForcesAndSourcesCore::UserDataOperator  {
 
    RVEBC_Data &dAta;
    Mat Aij;
    OpRVEBCsLhs(const string field_name, const string lagrang_field_name,Mat aij,RVEBC_Data &data):
    FaceElementForcesAndSourcesCore::UserDataOperator(
      lagrang_field_name,field_name,UserDataOperator::OPROWCOL
    ),
    dAta(data),Aij(aij) {
      //This will make sure to loop over all intities (e.g. for order=2 it will make doWork to loop 16 time)
      sYmm = false;
    }
 
    MatrixDouble Mat_face;
    MatrixDouble Mat_face_Tran;
 
    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);
        if(col_type==MBVERTEX) {
 
          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();
 
          Mat_face.resize(rank,3*rank);
          Mat_face.clear();
          Mat_face_Tran.resize(3*rank,rank);
          Mat_face_Tran.clear();
 
          switch(rank) {
            case 3:
            for(int nn=0; nn<3; nn++){
              Mat_face(0,3*nn+0) = 1.0;
              Mat_face(1,3*nn+1) = 1.0;
              Mat_face(2,3*nn+2) = 1.0;
            }
            break;
            case 1:
            Mat_face(0,0) = 1.0;
            Mat_face(0,1) = 1.0;
            Mat_face(0,2) = 1.0;
            break;
            default:
            SETERRQ(PETSC_COMM_SELF,1,"not implemented");
          }
 
          // Matrix C1
          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],&Mat_face(0,0),ADD_VALUES
          ); CHKERRQ(ierr);
 
          // Matrix C1T
          noalias(Mat_face_Tran) = trans(Mat_face);
          ierr = MatSetValues(
            Aij,nb_cols,&col_data.getIndices()[0],nb_rows,&row_data.getIndices()[0],&Mat_face_Tran(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);
    }
  };
 
  PetscErrorCode setRVEBCsRigidBodyTranOperators(
    string field_name,string lagrang_field_name,Mat aij,map<int,RVEBC_Data> setOfRVEBC
  ) {
    PetscFunctionBegin;
    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));
    }
    PetscFunctionReturn(0);
  }
 
};
 
#endif