v0.13.0
BCs_RVELagrange_Disp.hpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2  * MoFEM is free software: you can redistribute it and/or modify it under
3  * the terms of the GNU Lesser General Public License as published by the
4  * Free Software Foundation, either version 3 of the License, or (at your
5  * option) any later version.
6  *
7  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10  * License for more details.
11  *
12  * You should have received a copy of the GNU Lesser General Public
13  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14 
15 #ifndef __BCS_RVELAGRANGE_DISP_HPP
16 #define __BCS_RVELAGRANGE_DISP_HPP
17 
19 
22  int getRule(int order) { return order; };
23  };
24 
25  boost::ptr_vector<MethodForForceScaling> methodsOp; //This is to apply strain vector in steps (for given steps) for nonlinear problem
26 
27 
29  MyTriFE feRVEBCRhs; //To calculate the Rhs or RVE BCs
30  MyTriFE feRVEBCLhs; //To calculate the Lhs or RVE BCs
31  MyTriFE feRVEBCStress; //To calculate the Rhs or RVE BCs
32  MyTriFE feRVEBCRhsResidual; //To calculate the Rhs or RVE BCs
33  MyTriFE feRVEBCRhsHomoC; //To be used for the calculation of homgenised stiffness matrix
34 
37 
40 
42 
43 
45  mField(m_field),
46  feRVEBCRhs(m_field),
47  feRVEBCLhs(m_field),
48  feRVEBCRhsResidual(m_field),
49  feRVEBCStress(m_field),
50  feRVEBCRhsHomoC(m_field){
51  }
52 
53  struct RVEBC_Data {
54  Range tRis; // All boundary surfaces
55  };
56  map<int,RVEBC_Data> setOfRVEBC; ///< maps side set id with appropriate FluxData
57 
58  PetscErrorCode addLagrangiangElement(
59  const string element_name,
60  const string field_name,
61  const string lagrang_field_name,
62  const string mesh_nodals_positions
63  ) {
64  PetscFunctionBegin;
65 
66 
67 
68 
69  ierr = mField.add_finite_element(element_name,MF_ZERO); CHKERRQ(ierr);
70 
71  //============================================================================================================
72  //C row as Lagrange_mul_disp and col as DISPLACEMENT
73  ierr = mField.modify_finite_element_add_field_row(element_name,lagrang_field_name); CHKERRQ(ierr);
74  ierr = mField.modify_finite_element_add_field_col(element_name,field_name); CHKERRQ(ierr);
75  //CT col as Lagrange_mul_disp and row as DISPLACEMENT
76  ierr = mField.modify_finite_element_add_field_col(element_name,lagrang_field_name); CHKERRQ(ierr);
77  ierr = mField.modify_finite_element_add_field_row(element_name,field_name); CHKERRQ(ierr);
78  //data
79  ierr = mField.modify_finite_element_add_field_data(element_name,lagrang_field_name); CHKERRQ(ierr);
80  ierr = mField.modify_finite_element_add_field_data(element_name,field_name); CHKERRQ(ierr);
81  //============================================================================================================
82 
83  if(mField.check_field(mesh_nodals_positions)) { //for high order geometry
84  ierr = mField.modify_finite_element_add_field_data(element_name,mesh_nodals_positions); CHKERRQ(ierr);
85  }
86 
87  //this is alternative method for setting boundary conditions, to bypass bu in cubit file reader.
88  //not elegant, but good enough
90  if(it->getName().compare(0,12,"AllBoundSurf") == 0 || it->getMeshsetId() == 103) {
91  // cout<<"Hi from addLagrangiangElement "<<endl;
92  rval = mField.get_moab().get_entities_by_type(it->meshset,MBTRI,setOfRVEBC[it->getMeshsetId()].tRis,true); CHKERRQ_MOAB(rval);
93  ierr = mField.add_ents_to_finite_element_by_type(setOfRVEBC[it->getMeshsetId()].tRis,MBTRI,element_name); CHKERRQ(ierr);
94  }
95  }
96 
97  PetscFunctionReturn(0);
98  }
99 
100  /// \biref operator to calculate the LHS for the RVE bounary conditions
102 
105  Mat Aij;
107  const string field_name, const string lagrang_field_name, Mat aij,
108  RVEBC_Data &data,bool ho_geometry = false
109  ):
111  lagrang_field_name, field_name, UserDataOperator::OPROWCOL
112  ),
113  dAta(data),hoGeometry(ho_geometry),Aij(aij) {
114  // This will make sure to loop over all entities
115  // (e.g. for order=2 it will make doWork to loop 16 time)
116  sYmm = false;
117  }
118 
120 
121  /** \brief calculate thermal convection term in the lhs of equations
122  *
123  * C = intS N^T N dS
124  */
125  PetscErrorCode doWork(
126  int row_side,int col_side,
127  EntityType row_type,EntityType col_type,
130  ) {
131  PetscFunctionBegin;
132 
133  try {
134  int nb_row = row_data.getIndices().size();
135  int nb_col = col_data.getIndices().size();
136 
137  if(nb_row == 0) PetscFunctionReturn(0);
138  if(nb_col == 0) PetscFunctionReturn(0);
139 
140  auto weak_ptr_dof =
141  getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(
142  row_data.getIndices()[0]);
143  const FENumeredDofEntity *dof_ptr;
144  if(auto ptr = weak_ptr_dof.lock())
145  dof_ptr = ptr.get();
146  else
147  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Dof not found");
148 
149  int rank = dof_ptr->getNbOfCoeffs();
150 
151  K.resize(nb_row/rank,nb_col/rank);
152  K.clear();
153 
154  for(unsigned int gg = 0;gg<row_data.getN().size1();gg++) {
155  double area;
156  if(hoGeometry) {
157  area = norm_2(getNormalsAtGaussPts(gg))*0.5;
158  } else {
159  area = getArea();
160  }
161  double val = getGaussPts()(2,gg)*area;
162  noalias(K) += val*outer_prod(
163  row_data.getN(gg,nb_row/rank),col_data.getN(gg,nb_col/rank)
164  );
165  }
166 
167  VectorInt row_indices,col_indices;
168  row_indices.resize(nb_row/rank);
169  col_indices.resize(nb_col/rank);
170 
171  for(int rr = 0;rr < rank; rr++) {
172  unsigned int nb_rows;
173  unsigned int nb_cols;
174  int *rows;
175  int *cols;
176  if(rank > 1) {
177 
178  ublas::noalias(row_indices) = ublas::vector_slice<VectorInt >
179  (row_data.getIndices(), ublas::slice(rr, rank, row_data.getIndices().size()/rank));
180 
181  ublas::noalias(col_indices) = ublas::vector_slice<VectorInt >
182  (col_data.getIndices(), ublas::slice(rr, rank, col_data.getIndices().size()/rank));
183 
184  nb_rows = row_indices.size();
185  nb_cols = col_indices.size();
186  rows = &*row_indices.data().begin();
187  cols = &*col_indices.data().begin();
188 
189  } else {
190 
191  nb_rows = row_data.getIndices().size();
192  nb_cols = col_data.getIndices().size();
193  rows = &*row_data.getIndices().data().begin();
194  cols = &*col_data.getIndices().data().begin();
195 
196  }
197 
198  // Matrix C
199  ierr = MatSetValues(Aij,nb_rows,rows,nb_cols,cols,&K(0,0),ADD_VALUES); CHKERRQ(ierr);
200 
201  // Matrix CT
202  transK.resize(nb_col/rank,nb_row/rank);
203  noalias(transK) = trans(K);
204  ierr = MatSetValues(Aij,nb_cols,cols,nb_rows,rows,&transK(0,0),ADD_VALUES); CHKERRQ(ierr);
205 
206  }
207 
208  } catch (const std::exception& ex) {
209  ostringstream ss;
210  ss << "throw in method: " << ex.what() << endl;
211  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
212  }
213 
214  PetscFunctionReturn(0);
215  }
216  };
217 
219 
222 
224  const string field_name, const string lagrang_field_name,
225  RVEBC_Data &data,bool ho_geometry = false
226  ):
228  lagrang_field_name,UserDataOperator::OPROW
229  ),
230  dAta(data),hoGeometry(ho_geometry) {
231  }
232 
235 
237 
238  PetscErrorCode calculateDmat(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
239  PetscFunctionBegin;
240 
241  double x,y,z;
242  for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
243 
244  double area;
245  if(hoGeometry) {
246  area = norm_2(getNormalsAtGaussPts(gg))*0.5;
247  } else {
248  area = getArea();
249  }
250  double val = getGaussPts()(2,gg)*area;
251 
252  x = getCoordsAtGaussPts()(gg,0);
253  y = getCoordsAtGaussPts()(gg,1);
254  z = getCoordsAtGaussPts()(gg,2);
255 
256  switch(rank) {
257  case 3: //mech problem
258  X_mat.resize(3,6,false);
259  X_mat.clear();
260  X_mat(0,0) = 2.0*x; X_mat(0,3)=y; X_mat(0,5)=z;
261  X_mat(1,1) = 2.0*y; X_mat(1,3)=x; X_mat(1,4)=z;
262  X_mat(2,2) = 2.0*z; X_mat(2,4)=y; X_mat(2,5)=x;
263  X_mat = 0.5*X_mat;
264  break;
265  case 1: //moisture transport or thermal problem
266  X_mat.resize(3,3,false);
267  X_mat.clear();
268  X_mat(0,0)= x; X_mat(0,1)= y; X_mat(0,2) = z;
269  break;
270  default:
271  SETERRQ(PETSC_COMM_SELF,1,"not implemented");
272  }
273 
274  int shape_size = data.getN().size2();
275  N_mat.resize(rank,shape_size*rank);
276  N_mat.clear();
277 
278  // FIXME: this is inefficient implementation
279  {
280  int kk=0;
281  for(int ii=0; ii<shape_size; ii++){
282  //number of shape functions
283  double val = data.getN()(gg,ii);
284  for(int jj=0; jj<rank; jj++) {
285  N_mat(jj,kk) = val;
286  kk++;
287  }
288  }
289  if(gg==0) {
290  D_mat = val*prod(trans(N_mat),X_mat);
291  } else{
292  D_mat += val*prod(trans(N_mat),X_mat);
293  }
294  }
295 
296  }
297 
298  PetscFunctionReturn(0);
299  }
300 
301  };
302 
303 
304  /// \biref operator to calculate the RHS of the constrain for the RVE boundary conditions
305  struct OpRVEBCsRhs:public OpDmatRhs {
308  boost::ptr_vector<MethodForForceScaling> &methodsOp;
309 
311  const string field_name, const string lagrang_field_name,
312  Vec f, VectorDouble given_strain, boost::ptr_vector<MethodForForceScaling> &methods_op, RVEBC_Data &data,bool ho_geometry = false
313  ):
314  OpDmatRhs(field_name,lagrang_field_name,data,ho_geometry),F(f),givenStrain(given_strain),methodsOp(methods_op) {
315  }
316 
317  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
318  PetscFunctionBegin;
319 
320  nb_row_dofs = data.getIndices().size();
321 
322  if(nb_row_dofs ==0) PetscFunctionReturn(0);
323  if(dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt())==dAta.tRis.end()) PetscFunctionReturn(0);
324 
325  auto weak_ptr_dof =
326  getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(
327  data.getIndices()[0]);
328  const FENumeredDofEntity *dof_ptr;
329  if (auto ptr = weak_ptr_dof.lock())
330  dof_ptr = ptr.get();
331  else
332  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Dof not found");
333 
334  rank = dof_ptr->getNbOfCoeffs();
335  nb_row_dofs /= rank;
336 
337  ierr = calculateDmat(side,type,data); CHKERRQ(ierr);
338 
339  //This is to scale givenStrain vector to apply the strian in steps
340  VectorDouble scaled_given_strain;
341  scaled_given_strain.resize(6);
342  scaled_given_strain=givenStrain;
343 // cout <<"givenStrain Before = "<<scaled_given_strain<<endl;
344  ierr = MethodForForceScaling::applyScale(getFEMethod(),methodsOp,scaled_given_strain); CHKERRQ(ierr);
345 // cout <<"givenStrain After = "<<scaled_given_strain<<endl;
346  f.resize(D_mat.size1(),false);
347  noalias(f) = prod(D_mat, scaled_given_strain);
348  f*=-1;
349 
350 // cout <<"f = "<<f<<endl;
351 // string wait;
352 // cin>>wait;
353 
354  Vec myF = F;
355  if(F == PETSC_NULL) {
356  myF = getFEMethod()->snes_f;
357  }
358 
359  ierr = VecSetValues(myF,data.getIndices().size(),&data.getIndices()[0],&f[0],ADD_VALUES); CHKERRQ(ierr);
360  PetscFunctionReturn(0);
361  }
362  };
363 
364  /// \biref operator to calculate the RHS for the calculation of the homgoensied stiffness matrix
365  struct OpRVEBCsRhsHomoC:public OpDmatRhs {
366 
367  vector<Vec> &F;
368 
370  const string field_name, const string lagrang_field_name,
371  vector<Vec> &f,RVEBC_Data &data,bool ho_geometry = false
372  ):
373  OpDmatRhs(field_name,lagrang_field_name,data,ho_geometry),F(f) {
374  }
375 
376  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
377  PetscFunctionBegin;
378 
379  nb_row_dofs = data.getIndices().size();
380 
381  if(nb_row_dofs ==0) PetscFunctionReturn(0);
382  if(dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt())==dAta.tRis.end()) PetscFunctionReturn(0);
383 
384  auto weak_ptr_dof =
385  getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(
386  data.getIndices()[0]);
387  const FENumeredDofEntity *dof_ptr;
388  if (auto ptr = weak_ptr_dof.lock())
389  dof_ptr = ptr.get();
390  else
391  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Dof not found");
392 
393  rank = dof_ptr->getNbOfCoeffs();
394  nb_row_dofs /= rank;
395 
396  ierr = calculateDmat(side,type,data); CHKERRQ(ierr);
397 
398  applied_strain.resize(D_mat.size2(),false);
399  f.resize(D_mat.size1(),false);
400 
401  int size = (rank == 3) ? 6 : 3;
402  for(int ii = 0;ii<size;ii++) {
403  applied_strain.clear();
404  applied_strain[ii] = 1.0;
405  noalias(f) = prod(D_mat, applied_strain);
406  ierr = VecSetValues(
407  F[ii],data.getIndices().size(),&data.getIndices()[0],&f[0],ADD_VALUES
408  ); CHKERRQ(ierr);
409  }
410 
411  PetscFunctionReturn(0);
412 
413  }
414  };
415 
417  PetscErrorCode shapeMat(int rank, unsigned int gg, DataForcesAndSourcesCore::EntData &col_data, MatrixDouble &N_mat) {
418  PetscFunctionBegin;
419  int shape_size=col_data.getN().size2();
420  // cout<<"shape_size = "<<shape_size<<endl;
421  // cout<<"rank = "<<rank<<endl;
422  N_mat.resize(rank,shape_size*rank); N_mat.clear();
423  int gg1=0; int kk=0;
424  for(int ii=0; ii<shape_size; ii++){ //number of shape funcitons
425  for(int jj=0; jj<rank; jj++){
426  // cout<<"ii "<<ii<<endl;
427  N_mat(jj,kk)=col_data.getN()(gg,gg1); kk++;
428  }
429  gg1++;
430  }
431  PetscFunctionReturn(0);
432  }
433  };
435 
436 
437  /// \biref operator to calculate and assemble CU
439 
444 
445 
447  const string field_name, const string lagrang_field_name, Vec f,
448  CommonFunctions &_common_functions, RVEBC_Data &data, bool ho_geometry = false
449  ):
451  lagrang_field_name, field_name, UserDataOperator::OPROWCOL
452  ),
453  dAta(data),hoGeometry(ho_geometry),F(f), commonFunctions(_common_functions){
454  // This will make sure to loop over all entities
455  // (e.g. for order=2 it will make doWork to loop 16 time)
456  sYmm = false;
457  }
458 
461 
462  PetscErrorCode doWork(
463  int row_side,int col_side,
464  EntityType row_type,EntityType col_type,
467  ) {
468  PetscFunctionBegin;
469 
470  try {
471  int nb_row = row_data.getIndices().size();
472  int nb_col = col_data.getIndices().size();
473 
474  if(nb_row == 0) PetscFunctionReturn(0);
475  if(nb_col == 0) PetscFunctionReturn(0);
476 
477  auto weak_ptr_dof =
478  getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(
479  row_data.getIndices()[0]);
480  const FENumeredDofEntity *dof_ptr;
481  if (auto ptr = weak_ptr_dof.lock())
482  dof_ptr = ptr.get();
483  else
484  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Dof not found");
485 
486  int rank = dof_ptr->getNbOfCoeffs();
487 
488  for(unsigned int gg = 0;gg<row_data.getN().size1();gg++) {
489  double area;
490  if(hoGeometry) {
491  area = norm_2(getNormalsAtGaussPts(gg))*0.5;
492  } else {
493  area = getArea();
494  }
495  double val = getGaussPts()(2,gg)*area;
496  ierr = commonFunctions.shapeMat(rank, gg, row_data, N_mat_row); CHKERRQ(ierr);
497  ierr = commonFunctions.shapeMat(rank, gg, col_data, N_mat_col); CHKERRQ(ierr);
498 
499 // cout<<"row_data.getN() "<<row_data.getN()<<endl;
500 // cout<<"N_mat_row() "<<N_mat_row<<endl;
501 //
502 // cout<<"col_data.getN() "<<col_data.getN()<<endl;
503 // cout<<"N_mat_col "<<N_mat_col<<endl;
504 
505  if(gg==0){
506  K =val* prod(trans(N_mat_row),N_mat_col); //we don't need to define size K
507  } else {
508  K += val*prod(trans(N_mat_row),N_mat_col);
509  }
510  }
511  CU=prod(K, col_data.getFieldData());
512  CTLam=prod(trans(K),row_data.getFieldData());
513 
514 
515  Vec myF = F;
516  if(F == PETSC_NULL) {
517  myF = getFEMethod()->snes_f;
518  }
519 
520  ierr = VecSetValues(myF,row_data.getIndices().size(),&row_data.getIndices()[0],&CU[0],ADD_VALUES); CHKERRQ(ierr);
521  ierr = VecSetValues(myF,col_data.getIndices().size(),&col_data.getIndices()[0],&CTLam[0],ADD_VALUES); CHKERRQ(ierr);
522 
523  } catch (const std::exception& ex) {
524  ostringstream ss;
525  ss << "throw in method: " << ex.what() << endl;
526  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
527  }
528 
529  PetscFunctionReturn(0);
530  }
531  };
532 
533 
534  //for nonlinear problems
536  string field_name,
537  string lagrang_field_name,
538  string mesh_nodals_positions,
539  Mat aij,vector<Vec> &fvec,Vec f,VectorDouble given_strain
540  ) {
541  PetscFunctionBegin;
542 
543  bool ho_geometry = false;
544  if(mField.check_field(mesh_nodals_positions)) {
545  ho_geometry = true;
546  }
547  map<int,RVEBC_Data>::iterator sit = setOfRVEBC.begin();
548  for(;sit!=setOfRVEBC.end();sit++) {
549 
550  //Lhs (calculate LHS)
551  feRVEBCLhs.getOpPtrVector().push_back(
552  new OpRVEBCsLhs(field_name,lagrang_field_name, aij, sit->second,ho_geometry)
553  );
554 
555 
556  //Rhs (calculate RHS)
557  feRVEBCRhs.getOpPtrVector().push_back(
558  new OpRVEBCsRhs(field_name,lagrang_field_name,f,given_strain, methodsOp, sit->second,ho_geometry)
559  );
560 
561 
562  //Rhs (F[6] used to calculate Homgenised stiffness matrix)
563  feRVEBCRhsHomoC.getOpPtrVector().push_back(
564  new OpRVEBCsRhsHomoC(field_name,lagrang_field_name,fvec,sit->second,ho_geometry)
565  );
566 
567  //for internal forces (i.e. RHS)
568  //calculate and assemble CU and CT Lam
569  feRVEBCRhsResidual.getOpPtrVector().push_back(
570  new OpRVEBCsRhsForceExternal_CU(field_name,lagrang_field_name,f,common_functions,sit->second,ho_geometry)
571  );
572 
573 
574  }
575 
576  PetscFunctionReturn(0);
577  }
578 
579  //for linear problems
580  PetscErrorCode setRVEBCsOperators(
581  string field_name,
582  string lagrang_field_name,
583  string mesh_nodals_positions,
584  Mat aij,vector<Vec> &fvec
585  ) {
586  PetscFunctionBegin;
587 
588  bool ho_geometry = false;
589  if(mField.check_field(mesh_nodals_positions)) {
590  ho_geometry = true;
591  }
592  map<int,RVEBC_Data>::iterator sit = setOfRVEBC.begin();
593  for(;sit!=setOfRVEBC.end();sit++) {
594 
595  //Lhs (calculate LHS)
596  feRVEBCLhs.getOpPtrVector().push_back(
597  new OpRVEBCsLhs(field_name,lagrang_field_name, aij, sit->second,ho_geometry)
598  );
599 
600  //Rhs (F[6] used to calculate Homgenised stiffness matrix)
601  feRVEBCRhs.getOpPtrVector().push_back(
602  new OpRVEBCsRhsHomoC(field_name,lagrang_field_name,fvec,sit->second,ho_geometry)
603  );
604  }
605  PetscFunctionReturn(0);
606  }
607 
608 
609  /// \biref operator to calculate the RVE homogenised stress
610  struct OpRVEHomoStress:public OpDmatRhs {
611 
613 
615  const string field_name,
616  const string lagrang_field_name,
617  Vec stress_homo,
618  RVEBC_Data &data,
619  bool ho_geometry = false
620  ):
621  OpDmatRhs(field_name,lagrang_field_name,data,ho_geometry),
622  Stress_Homo(stress_homo) {
623  }
624 
626 
627  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
628  PetscFunctionBegin;
629  if(data.getIndices().size()==0) PetscFunctionReturn(0);
630  if(dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt())==dAta.tRis.end()) PetscFunctionReturn(0);
631  // cout<<"OpRVEBCsRhs "<<endl;
632 
633 
634  auto weak_ptr_dof =
635  getNumeredEntFiniteElementPtr()->getRowDofsByPetscGlobalDofIdx(
636  data.getIndices()[0]);
637  const FENumeredDofEntity *dof_ptr;
638  if (auto ptr = weak_ptr_dof.lock())
639  dof_ptr = ptr.get();
640  else
641  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Dof not found");
642 
643  rank = dof_ptr->getNbOfCoeffs();
644  nb_row_dofs = data.getIndices().size()/rank;
645 
646  ierr = calculateDmat(side,type,data); CHKERRQ(ierr);
647 
648  Stress_Homo_elem.resize(D_mat.size2());
649 
650  // data.getFieldData() is value of Lagrange multply
651  // data.getFieldData() is reaction force (so multiply for -1 to get the force)
652  noalias(Stress_Homo_elem) = prod(trans(D_mat), -data.getFieldData());
653 
654  const int Indices6[6] = {0, 1, 2, 3, 4, 5};
655  const int Indices3[3] = {0, 1, 2};
656  switch(rank) {
657  case 3:
658  ierr = VecSetValues(Stress_Homo,6,Indices6,&(Stress_Homo_elem.data())[0],ADD_VALUES); CHKERRQ(ierr);
659  break;
660  case 1:
661  ierr = VecSetValues(Stress_Homo,3,Indices3,&(Stress_Homo_elem.data())[0],ADD_VALUES); CHKERRQ(ierr);
662  break;
663  default:
664  SETERRQ(PETSC_COMM_SELF,1,"not implemented");
665  }
666 
667  PetscFunctionReturn(0);
668  }
669  };
670 
672  string field_name,
673  string lagrang_field_name,
674  string mesh_nodals_positions,
675  Vec Stress_Homo
676  ) {
677  PetscFunctionBegin;
678 
679  bool ho_geometry = false;
680  if(mField.check_field(mesh_nodals_positions)) {
681  ho_geometry = true;
682  }
683 
684  map<int,RVEBC_Data>::iterator sit = setOfRVEBC.begin();
685  for(;sit!=setOfRVEBC.end();sit++) {
686  feRVEBCStress.getOpPtrVector().push_back(new OpRVEHomoStress(field_name, lagrang_field_name, Stress_Homo, sit->second, ho_geometry));
687  }
688 
689  PetscFunctionReturn(0);
690  }
691 };
692 
693 #endif // __BCS_RVELAGRANGE_DISP_HPP
ForcesAndSourcesCore::UserDataOperator UserDataOperator
EntitiesFieldData::EntData EntData
@ MF_ZERO
Definition: definitions.h:111
@ SIDESET
Definition: definitions.h:160
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
UBlasVector< int > VectorInt
Definition: Types.hpp:78
PetscErrorCode shapeMat(int rank, unsigned int gg, DataForcesAndSourcesCore::EntData &col_data, MatrixDouble &N_mat)
MyTriFE(MoFEM::Interface &_mField)
PetscErrorCode calculateDmat(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpDmatRhs(const string field_name, const string lagrang_field_name, RVEBC_Data &data, bool ho_geometry=false)
\biref operator to calculate the LHS for the RVE bounary conditions
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
calculate thermal convection term in the lhs of equations
OpRVEBCsLhs(const string field_name, const string lagrang_field_name, Mat aij, RVEBC_Data &data, bool ho_geometry=false)
\biref operator to calculate and assemble CU
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpRVEBCsRhsForceExternal_CU(const string field_name, const string lagrang_field_name, Vec f, CommonFunctions &_common_functions, RVEBC_Data &data, bool ho_geometry=false)
\biref operator to calculate the RHS for the calculation of the homgoensied stiffness matrix
OpRVEBCsRhsHomoC(const string field_name, const string lagrang_field_name, vector< Vec > &f, RVEBC_Data &data, bool ho_geometry=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
\biref operator to calculate the RHS of the constrain for the RVE boundary conditions
OpRVEBCsRhs(const string field_name, const string lagrang_field_name, Vec f, VectorDouble given_strain, boost::ptr_vector< MethodForForceScaling > &methods_op, RVEBC_Data &data, bool ho_geometry=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
boost::ptr_vector< MethodForForceScaling > & methodsOp
\biref operator to calculate the RVE homogenised stress
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpRVEHomoStress(const string field_name, const string lagrang_field_name, Vec stress_homo, RVEBC_Data &data, bool ho_geometry=false)
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
BCs_RVELagrange_Disp(MoFEM::Interface &m_field)
boost::ptr_vector< MethodForForceScaling > methodsOp
CommonFunctions common_functions
PetscErrorCode setRVEBCsOperatorsNonlinear(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &fvec, Vec f, VectorDouble given_strain)
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec Stress_Homo)
MoFEM::Interface & mField
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &fvec)
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions)
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
virtual moab::Interface & get_moab()=0
Deprecated interface functions.