BCs_RVELagrange_Trac_Rigid_Rot.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_RIGID_ROT_HPP
#define __BCS_RVELAGRANGE_TRAC_RIGID_ROT_HPP
 
struct BCs_RVELagrange_Trac_Rigid_Rot: public BCs_RVELagrange_Disp {
 
  BCs_RVELagrange_Trac_Rigid_Rot(MoFEM::Interface &m_field): BCs_RVELagrange_Disp(m_field){}
 
  /// \biref operator to calculate the LHS for the RVE boundary conditions
  struct OpRVEBCsLhs:public FaceElementForcesAndSourcesCore::UserDataOperator  {
 
    MoFEM::Interface &mField;
    Mat Aij;
    RVEBC_Data &dAta;
 
    OpRVEBCsLhs(
      MoFEM::Interface &_mField, const string field_name, const string lagrang_field_name, Mat aij, RVEBC_Data &data
    ):
    FaceElementForcesAndSourcesCore::UserDataOperator(
      lagrang_field_name, field_name, UserDataOperator::OPROWCOL),
      mField(_mField),Aij(aij), dAta(data){
      //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 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){
          
 
          Mat_face.resize(3,9);
          Mat_face.clear();
          Mat_face_Tran.resize(9,3);
          Mat_face_Tran.clear();
 
          int num_nodes;
          const EntityHandle* conn_face;
          double coords_face[9];
          EntityHandle face_tri = getNumeredEntFiniteElementPtr()->getEnt(); //handle of finite element
 
          
          rval = mField.get_moab().get_connectivity(face_tri,conn_face,num_nodes,true); CHKERRQ_MOAB(rval);
          rval = mField.get_moab().get_coords(conn_face,num_nodes,coords_face); CHKERRQ_MOAB(rval);
 
          for(int nn=0; nn<3; nn++){
            Mat_face(0,3*nn+1)=-coords_face[3*nn+2];    Mat_face(0,3*nn+2)= coords_face[3*nn+1];
            Mat_face(1,3*nn+0)= coords_face[3*nn+2];    Mat_face(1,3*nn+2)=-coords_face[3*nn+0];
            Mat_face(2,3*nn+0)=-coords_face[3*nn+1];    Mat_face(2,3*nn+1)= coords_face[3*nn+0];
          }
 
          // 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 setRVEBCsRigidBodyRotOperators(
    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(mField, field_name,lagrang_field_name,aij,sit->second)
      );
    }
    PetscFunctionReturn(0);
  }
 
};
 
#endif