v0.13.0
BCs_RVELagrange_Trac_Rigid_Rot.hpp
Go to the documentation of this file.
1 /* Copyright (C) 2014, Zahur Ullah (Zahur.Ullah AT glasgow.ac.uk)
2  * --------------------------------------------------------------
3  */
4 
5 /* This file is part of MoFEM.
6  * MoFEM is free software: you can redistribute it and/or modify it under
7  * the terms of the GNU Lesser General Public License as published by the
8  * Free Software Foundation, either version 3 of the License, or (at your
9  * option) any later version.
10  *
11  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14  * License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
18 
19 #ifndef __BCS_RVELAGRANGE_TRAC_RIGID_ROT_HPP
20 #define __BCS_RVELAGRANGE_TRAC_RIGID_ROT_HPP
21 
23 
25 
26  /// \biref operator to calculate the LHS for the RVE boundary conditions
28 
30  Mat Aij;
32 
34  MoFEM::Interface &_mField, const string field_name, const string lagrang_field_name, Mat aij, RVEBC_Data &data
35  ):
37  lagrang_field_name, field_name, UserDataOperator::OPROWCOL),
38  mField(_mField),Aij(aij), dAta(data){
39  //This will make sure to loop over all entities
40  //(e.g. for order=2 it will make doWork to loop 16 time)
41  sYmm = false;
42  }
43 
46 
47  PetscErrorCode doWork(
48  int row_side,int col_side,
49  EntityType row_type,EntityType col_type,
52  ) {
53  PetscFunctionBegin;
54 
55  try {
56  if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
57  if(col_data.getIndices().size()==0) PetscFunctionReturn(0);
58  if(col_type==MBVERTEX){
59 
60 
61  Mat_face.resize(3,9);
62  Mat_face.clear();
63  Mat_face_Tran.resize(9,3);
64  Mat_face_Tran.clear();
65 
66  int num_nodes;
67  const EntityHandle* conn_face;
68  double coords_face[9];
69  EntityHandle face_tri = getNumeredEntFiniteElementPtr()->getEnt(); //handle of finite element
70 
71 
72  rval = mField.get_moab().get_connectivity(face_tri,conn_face,num_nodes,true); CHKERRQ_MOAB(rval);
73  rval = mField.get_moab().get_coords(conn_face,num_nodes,coords_face); CHKERRQ_MOAB(rval);
74 
75  for(int nn=0; nn<3; nn++){
76  Mat_face(0,3*nn+1)=-coords_face[3*nn+2]; Mat_face(0,3*nn+2)= coords_face[3*nn+1];
77  Mat_face(1,3*nn+0)= coords_face[3*nn+2]; Mat_face(1,3*nn+2)=-coords_face[3*nn+0];
78  Mat_face(2,3*nn+0)=-coords_face[3*nn+1]; Mat_face(2,3*nn+1)= coords_face[3*nn+0];
79  }
80 
81  // Matrix C1
82  int nb_rows=row_data.getIndices().size();
83  int nb_cols=col_data.getIndices().size();
85  Aij,nb_rows,&row_data.getIndices()[0],nb_cols,&col_data.getIndices()[0],&Mat_face(0,0),ADD_VALUES
86  ); CHKERRQ(ierr);
87 
88  // Matrix C1T
89  noalias(Mat_face_Tran) = trans(Mat_face);
91  Aij,nb_cols,&col_data.getIndices()[0],nb_rows,&row_data.getIndices()[0],&Mat_face_Tran(0,0),ADD_VALUES
92  ); CHKERRQ(ierr);
93 
94  }
95 
96  } catch (const std::exception& ex) {
97  ostringstream ss;
98  ss << "throw in method: " << ex.what() << endl;
99  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
100  }
101  PetscFunctionReturn(0);
102  }
103  };
104 
106  string field_name,string lagrang_field_name,Mat aij,map<int,RVEBC_Data> setOfRVEBC
107  ) {
108  PetscFunctionBegin;
109  map<int,RVEBC_Data>::iterator sit = setOfRVEBC.begin();
110  for(;sit!=setOfRVEBC.end();sit++) {
111  //LHS
112  feRVEBCLhs.getOpPtrVector().push_back(
113  new OpRVEBCsLhs(mField, field_name,lagrang_field_name,aij,sit->second)
114  );
115  }
116  PetscFunctionReturn(0);
117  }
118 
119 };
120 
121 #endif
ForcesAndSourcesCore::UserDataOperator UserDataOperator
EntitiesFieldData::EntData EntData
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
MoFEM::Interface & mField
\biref operator to calculate the LHS for the RVE boundary conditions
OpRVEBCsLhs(MoFEM::Interface &_mField, const string field_name, const string lagrang_field_name, Mat aij, RVEBC_Data &data)
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
PetscErrorCode setRVEBCsRigidBodyRotOperators(string field_name, string lagrang_field_name, Mat aij, map< int, RVEBC_Data > setOfRVEBC)
BCs_RVELagrange_Trac_Rigid_Rot(MoFEM::Interface &m_field)
virtual moab::Interface & get_moab()=0
Deprecated interface functions.