v0.13.0
BCs_RVELagrange_Trac.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_HPP
20 #define __BCS_RVELAGRANGE_TRAC_HPP
21 
23 
25 
26  struct CommonData {
28  };
30 
31  struct CommonFunctions {
32 
33  PetscErrorCode shapeMat(
34  int rank,unsigned int gg,DataForcesAndSourcesCore::EntData &col_data,MatrixDouble &n
35  ) {
36  PetscFunctionBegin;
37 
38  int shape_size = col_data.getN().size2();
39 
40  n.resize(rank,shape_size*rank);
41  n.clear();
42 
43  int kk = 0;
44  for(int ii=0; ii<shape_size; ii++){
45  double val = col_data.getN()(gg,ii);
46  for(int jj=0; jj<rank; jj++){
47  n(jj,kk) = val;
48  kk++;
49  }
50  }
51 
52  PetscFunctionReturn(0);
53  }
54 
55 
57 
58  PetscErrorCode hMat(
59  VectorDouble face_normal,
60  int rank,
61  MatrixDouble &n_mat,
62  MatrixDouble &h_mat
63  ) {
64  PetscFunctionBegin;
65 
66  switch(rank) {
67  //Mechanical problem
68  case 3: {
69  H_mat_1Node.resize(6,3);
70  H_mat_1Node.clear(); //for one node
71  if(face_normal[0]>0) { //+X face of the RVE
72  H_mat_1Node(0,0) = 1.0;
73  H_mat_1Node(3,1) = 1.0;
74  H_mat_1Node(5,2) = 1.0;
75  } else if(face_normal[0]<0) { //-X face of the RVE
76  H_mat_1Node(0,0) = -1.0;
77  H_mat_1Node(3,1) = -1.0;
78  H_mat_1Node(5,2) = -1.0;
79  } else if(face_normal[1]>0) { //+Y face of the RVE
80  H_mat_1Node(1,1) = 1.0;
81  H_mat_1Node(3,0) = 1.0;
82  H_mat_1Node(4,2) = 1.0;
83  } else if(face_normal[1]<0) { //-Y face of the RVE
84  H_mat_1Node(1,1) = -1.0;
85  H_mat_1Node(3,0) = -1.0;
86  H_mat_1Node(4,2) = -1.0;
87  } else if(face_normal[2]>0) { //+Z face of the RVE
88  H_mat_1Node(2,2) = 1.0;
89  H_mat_1Node(4,1) = 1.0;
90  H_mat_1Node(5,0) = 1.0;
91  } else if(face_normal[2]<0) { //-Z face of the RVE
92  H_mat_1Node(2,2) = -1.0;
93  H_mat_1Node(4,1) = -1.0;
94  H_mat_1Node(5,0) = -1.0;
95  } else {
96  SETERRQ(PETSC_COMM_SELF,MOFEM_IMPOSIBLE_CASE,"Can not be ?!");
97  }
98 
99  int num_col = n_mat.size2();
100  h_mat.resize(6,num_col);
101 
102  int cc1 = 0;
103  for(int bb = 0; bb<num_col/3; bb++) {
104  //blocks of 6x3
105  for(int rr = 0; rr<6; rr++) {
106  for(int cc = 0; cc<3; cc++) {
107  h_mat(rr,(cc+cc1)) = H_mat_1Node(rr,cc);
108  }
109  }
110  cc1+=3;
111  }
112 
113  }
114  break;
115  case 1: //Moisture transport or thermal problem
116  {
117  H_mat_1Node.resize(3,1);
118  H_mat_1Node.clear();
119  if(face_normal[0]>0) { //+X face of the RVE
120  H_mat_1Node(0,0) = 1.0;
121  }
122  if(face_normal[0]<0) { //-X face of the RVE
123  H_mat_1Node(0,0) = -1.0;
124  }
125  if(face_normal[1]>0) { //+Y face of the RVE
126  H_mat_1Node(1,0) = 1.0;
127  }
128  if(face_normal[1]<0) { //-Y face of the RVE
129  H_mat_1Node(1,0) = -1.0;
130  }
131  if(face_normal[2]>0) { //+Z face of the RVE
132  H_mat_1Node(2,0) = 1.0;
133  }
134  if(face_normal[2]<0) { //-Z face of the RVE
135  H_mat_1Node(2,0) = -1.0;
136  }
137 
138  int num_col = n_mat.size2();
139  h_mat.resize(3,num_col);
140 
141  int cc1 = 0;
142  for(int bb = 0; bb<num_col; bb++) {
143  //blocks of 3x1
144  for(int rr = 0; rr<3; rr++) {
145  for(int cc = 0; cc<1; cc++) {
146  h_mat(rr,(cc+cc1))=H_mat_1Node(rr,cc);
147  }
148  }
149  cc1+=1;
150  }
151 
152  }
153  break;
154  default:
155  SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"not implemented");
156  }
157  PetscFunctionReturn(0);
158  }
159  };
160 
162 
163  /// \biref operator to calculate the LHS for the RVE bounary conditions
165 
169  Mat Aij;
170 
171  OpRVEBCsLhs(const string field_name, const string lagrang_field_name,Mat aij,
172  RVEBC_Data &data, CommonFunctions &common_functions, bool ho_geometry = false
173  ):
174  FaceElementForcesAndSourcesCore::UserDataOperator(lagrang_field_name, field_name, UserDataOperator::OPROWCOL),
175  dAta(data),commonFunctions(common_functions),hoGeometry(ho_geometry),Aij(aij) {
176  //This will make sure to loop over all entities
177  // (e.g. for order=2 it will make doWork to loop 16 time)
178  sYmm = false;
179  }
180 
182  ublas::vector<MatrixDouble > N_mat_row;
184 
185  PetscErrorCode doWork(
186  int row_side,int col_side,
187  EntityType row_type,EntityType col_type,
190  ) {
191  PetscFunctionBegin;
192 
193  try {
194  if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
195  if(col_data.getIndices().size()==0) PetscFunctionReturn(0);
196 
197  int nb_row = row_data.getIndices().size();
198  int nb_col = col_data.getIndices().size();
199 
200  auto weak_ptr_dof =
201  getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(
202  col_data.getIndices()[0]);
203  const FENumeredDofEntity *dof_ptr;
204  if (auto ptr = weak_ptr_dof.lock())
205  dof_ptr = ptr.get();
206  else
207  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Dof not found");
208 
209  int rank = dof_ptr->getNbOfCoeffs();
210 
211  K.resize(nb_row,nb_col);
212  K.clear();
213 
214  N_mat_row.resize(col_data.getN().size1());
215  for(unsigned int gg = 0;gg<col_data.getN().size1();gg++) {
216  double area;
217  if(hoGeometry) {
218  area = norm_2(getNormalsAtGaussPts(gg))*0.5;
219  } else {
220  area = getArea();
221  }
222  double val = getGaussPts()(2,gg)*area;
223 
224  if(col_type==MBVERTEX) {
225  ierr = commonFunctions.shapeMat(rank,gg,col_data,N_mat_row[gg]); CHKERRQ(ierr);
226  }
227  ierr = commonFunctions.shapeMat(rank,gg,col_data,N_mat_col); CHKERRQ(ierr);
228  ierr = commonFunctions.hMat(getNormalsAtGaussPts(gg),rank,N_mat_row(gg),H_mat); CHKERRQ(ierr);
229 
230  if(gg==0){
231  NTN = prod(trans(N_mat_row(gg)),N_mat_col); //we don't need to define size NTN
232  K = val*prod(H_mat,NTN); //we don't need to defiend size of K
233  } else {
234  noalias(NTN) = prod(trans(N_mat_row(gg)),N_mat_col);
235  K += val*prod(H_mat, NTN);
236  }
237  }
238 
239  // Matrix C
240  int nb_rows=row_data.getIndices().size();
241  int nb_cols=col_data.getIndices().size();
242  ierr = MatSetValues(
243  Aij,nb_rows,&row_data.getIndices()[0],nb_cols,&col_data.getIndices()[0],&K(0,0),ADD_VALUES
244  ); CHKERRQ(ierr);
245 
246  // Matrix CT
247  transK = trans(K);
248  ierr = MatSetValues(
249  Aij,nb_cols,&col_data.getIndices()[0],nb_rows,&row_data.getIndices()[0],&transK(0,0),ADD_VALUES
250  ); CHKERRQ(ierr);
251 
252  } catch (const std::exception& ex) {
253  ostringstream ss;
254  ss << "throw in method: " << ex.what() << endl;
255  SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,ss.str().c_str());
256  }
257  PetscFunctionReturn(0);
258  }
259 
260  };
261 
262  /// \biref operator to calculate the RHS of the constrain for the RVE boundary conditions
264 
269 
271  const string field_name, RVEBC_Data &data,CommonData &common_data,CommonFunctions &common_functions,bool ho_geometry = false
272  ):
274  dAta(data),
275  hoGeometry(ho_geometry),
276  commonData(common_data),
278  }
279 
281 
282  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
283  PetscFunctionBegin;
284 
285  if(type!=MBVERTEX) PetscFunctionReturn(0);
286  if(data.getIndices().size()==0) PetscFunctionReturn(0);
287 
288  auto weak_ptr_dof =
289  getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(
290  data.getIndices()[0]);
291  const FENumeredDofEntity *dof_ptr;
292  if (auto ptr = weak_ptr_dof.lock())
293  dof_ptr = ptr.get();
294  else
295  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Dof not found");
296 
297  int rank = dof_ptr->getNbOfCoeffs();
298 
299  double x,y,z;
300  for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
301  double area;
302  if(hoGeometry) {
303  area = norm_2(getNormalsAtGaussPts(gg))*0.5;
304  } else {
305  area = getArea();
306  }
307  double val = getGaussPts()(2,gg)*area;
308 
309  x = getCoordsAtGaussPts()(gg,0);
310  y = getCoordsAtGaussPts()(gg,1);
311  z = getCoordsAtGaussPts()(gg,2);
312 
313  switch(rank) {
314  case 3: //mech problem
315  X_mat.resize(rank,6,false);
316  X_mat.clear();
317  X_mat(0,0) = 2.0*x; X_mat(0,3) = y; X_mat(0,5) = z;
318  X_mat(1,1) = 2.0*y; X_mat(1,3) = x; X_mat(1,4) = z;
319  X_mat(2,2) = 2.0*z; X_mat(2,4) = y; X_mat(2,5) = x;
320  X_mat=0.5*X_mat;
321  break;
322  case 1: //moisture transport or thermal problem
323  X_mat.resize(rank,3,false);
324  X_mat.clear();
325  X_mat(0,0) = x;
326  X_mat(0,1) = y;
327  X_mat(0,2) = z;
328  break;
329  default:
330  SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"not implemented");
331  }
332  ierr = commonFunctions.shapeMat(rank, gg, data, N_mat); CHKERRQ(ierr);
333  ierr = commonFunctions.hMat(getNormalsAtGaussPts(gg), rank, N_mat, H_mat); CHKERRQ(ierr);
334 
335  if(gg==0) {
336  NTX = val*prod(trans(N_mat),X_mat);
337  commonData.D_mat = prod(H_mat, NTX);
338  } else {
339  noalias(NTX) = val*prod(trans(N_mat),X_mat);
340  commonData.D_mat += prod(H_mat, NTX);
341  }
342  }
343 
344  PetscFunctionReturn(0);
345  }
346  };
347 
348  /// \biref operator to calculate the RHS of the constrain for the RVE boundary conditions
350 
353  vector<Vec> &F;
354 
355  OpRVEBCsRhs_Assemble(const string lagrang_field_name,vector<Vec> &f,RVEBC_Data &data,CommonData &common_data):
357  dAta(data),commonData(common_data),F(f) {
358  }
359 
360 
362 
363  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
364  PetscFunctionBegin;
365 
366  if(data.getIndices().size()==0) PetscFunctionReturn(0);
367 
369  applied_strain.resize(D_mat.size2(),false);
370 
371  for(int ii = 0;ii<F.size();ii++) {
372  applied_strain.clear();
373  applied_strain(ii) = 1.0;
374  f.resize(D_mat.size1(),false);
375  noalias(f) = prod(D_mat,applied_strain);
376  ierr = VecSetValues(F[ii],data.getIndices().size(),&data.getIndices()[0],&f[0],ADD_VALUES); CHKERRQ(ierr);
377  }
378 
379  PetscFunctionReturn(0);
380  }
381  };
382 
383  PetscErrorCode setRVEBCsOperators(
384  string field_name,string lagrang_field_name,string mesh_nodals_positions,Mat aij,vector<Vec> &f
385  ) {
386  PetscFunctionBegin;
387  bool ho_geometry = false;
388  if(mField.check_field(mesh_nodals_positions)) {
389  ho_geometry = true;
390  }
391  map<int,RVEBC_Data>::iterator sit = setOfRVEBC.begin();
392  for(;sit!=setOfRVEBC.end();sit++) {
393  //LHS
394  feRVEBCLhs.getOpPtrVector().push_back(
395  new OpRVEBCsLhs(field_name,lagrang_field_name,aij,sit->second,commonFunctions,ho_geometry)
396  );
397  //RHS
398  feRVEBCRhs.getOpPtrVector().push_back(
399  new OpRVEBCsRhs_Cal(field_name,sit->second,commonData,commonFunctions,ho_geometry)
400  );
401  feRVEBCRhs.getOpPtrVector().push_back(
402  new OpRVEBCsRhs_Assemble(lagrang_field_name,f,sit->second,commonData)
403  );
404  }
405  PetscFunctionReturn(0);
406  }
407 
409 
413 
415  const string lagrang_field_name,Vec stress_homo,RVEBC_Data &data,CommonData &common_data
416  ):
418  dAta(data),commonData(common_data),stressHomo(stress_homo) {
419  }
420 
421 
423 
424  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
425  PetscFunctionBegin;
426  if(data.getIndices().size()==0) PetscFunctionReturn(0);
427  if(dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt())==dAta.tRis.end()) PetscFunctionReturn(0);
428  stressVector.resize(commonData.D_mat.size2());
429  //Lamda=data.getFieldData() is reaction force (so multiply for -1 to get the force)
430  noalias(stressVector) = prod(trans(commonData.D_mat),-data.getFieldData());
431  const int indices_6[6] = {0, 1, 2, 3, 4, 5};
432  const int indices_3[3] = {0, 1, 2};
433  switch(stressVector.size()) {
434  case 6:
435  ierr = VecSetValues(stressHomo,6,indices_6,&stressVector[0],ADD_VALUES); CHKERRQ(ierr);
436  break;
437  case 3:
438  ierr = VecSetValues(stressHomo,3,indices_3,&stressVector[0],ADD_VALUES); CHKERRQ(ierr);
439  break;
440  default:
441  SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"not implemented");
442  }
443  PetscFunctionReturn(0);
444  }
445 
446  };
447 
449  string field_name,string lagrang_field_name,string mesh_nodals_positions,Vec stress_homo
450  ) {
451  PetscFunctionBegin;
452  bool ho_geometry = false;
453  if(mField.check_field(mesh_nodals_positions)) {
454  ho_geometry = true;
455  }
456  map<int,RVEBC_Data>::iterator sit = setOfRVEBC.begin();
457  for(;sit!=setOfRVEBC.end();sit++) {
458  feRVEBCStress.getOpPtrVector().push_back(
459  new OpRVEBCsRhs_Cal(field_name,sit->second,commonData,commonFunctions)
460  );
461  feRVEBCStress.getOpPtrVector().push_back(
462  new OpRVEHomoStress_Assemble(lagrang_field_name,stress_homo,sit->second,commonData)
463  );
464  }
465  PetscFunctionReturn(0);
466  }
467 
468 };
469 
470 #endif
static Index< 'n', 3 > n
ForcesAndSourcesCore::UserDataOperator UserDataOperator
EntitiesFieldData::EntData EntData
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:48
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
virtual bool check_field(const std::string &name) const =0
check if field is in database
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
CommonFunctions common_functions
MoFEM::Interface & mField
PetscErrorCode hMat(VectorDouble face_normal, int rank, MatrixDouble &n_mat, MatrixDouble &h_mat)
PetscErrorCode shapeMat(int rank, unsigned int gg, DataForcesAndSourcesCore::EntData &col_data, MatrixDouble &n)
\biref operator to calculate the LHS for the RVE bounary conditions
ublas::vector< MatrixDouble > N_mat_row
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpRVEBCsLhs(const string field_name, const string lagrang_field_name, Mat aij, RVEBC_Data &data, CommonFunctions &common_functions, bool ho_geometry=false)
\biref operator to calculate the RHS of the constrain for the RVE boundary conditions
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpRVEBCsRhs_Assemble(const string lagrang_field_name, vector< Vec > &f, RVEBC_Data &data, CommonData &common_data)
\biref operator to calculate the RHS of the constrain for the RVE boundary conditions
OpRVEBCsRhs_Cal(const string field_name, RVEBC_Data &data, CommonData &common_data, CommonFunctions &common_functions, bool ho_geometry=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpRVEHomoStress_Assemble(const string lagrang_field_name, Vec stress_homo, RVEBC_Data &data, CommonData &common_data)
BCs_RVELagrange_Trac(MoFEM::Interface &m_field)
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &f)
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
Deprecated interface functions.