v0.13.1
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);
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);
78 //data
79 ierr = mField.modify_finite_element_add_field_data(element_name,lagrang_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
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
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
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
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);
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
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
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:98
@ SIDESET
Definition: definitions.h:147
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
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
auto f
Definition: HenckyOps.hpp:5
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
UBlasVector< int > VectorInt
Definition: Types.hpp:67
constexpr auto field_name
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.