v0.13.1
BCs_RVELagrange_Periodic.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
16#ifndef __BCS_RVELAGRANGE_PERIIODIC_HPP
17#define __BCS_RVELAGRANGE_PERIIODIC_HPP
18
20
21 boost::ptr_vector<MethodForForceScaling> methodsOp; //This is to apply strain vector in steps (for given steps) for nonlinear problem
22
24 Range pRisms; // All boundary surfaces
25 };
26 map<int,RVEBC_Data_Periodic> setOfRVEBCPrisms; ///< maps side set id with appropriate FluxData
27
30 int getRule(int order) { return order; };
31 };
32
33 MyPrismFE feRVEBCRhs; //To calculate the Rhs or RVE BCs
34 MyPrismFE feRVEBCLhs; //To calculate the Lhs or RVE BCs
35 MyPrismFE feRVEBCStress; //To calculate the Rhs or RVE BCs
36 MyPrismFE feRVEBCRhsResidual; //used to calculate the residual
37
40
43
45 BCs_RVELagrange_Disp(m_field),
46 feRVEBCRhs(m_field),
47 feRVEBCLhs(m_field),
48 feRVEBCRhsResidual(m_field),
49 feRVEBCStress(m_field) {
50 }
51
52 PetscErrorCode addLagrangiangElement(
53 const string element_name,
54 const string field_name,
55 const string lagrang_field_name,
56 const string mesh_nodals_positions,
57 Range &periodic_prisms
58 ) {
59 PetscFunctionBegin;
60
61 ierr = mField.add_finite_element(element_name,MF_ZERO); CHKERRQ(ierr);
62 //============================================================================================================
63 //C row as Lagrange_mul_disp and col as DISPLACEMENT
64 ierr = mField.modify_finite_element_add_field_row(element_name,lagrang_field_name); CHKERRQ(ierr);
66 //CT col as Lagrange_mul_disp and row as DISPLACEMENT
67 ierr = mField.modify_finite_element_add_field_col(element_name,lagrang_field_name); CHKERRQ(ierr);
69 //data
70 ierr = mField.modify_finite_element_add_field_data(element_name,lagrang_field_name); CHKERRQ(ierr);
72 //============================================================================================================
73
74 if(mField.check_field(mesh_nodals_positions)) { //for high order geometry
75 ierr = mField.modify_finite_element_add_field_data(element_name,mesh_nodals_positions); CHKERRQ(ierr);
76 }
77
78 //Adding Prisims to Element Lagrange_elem (to loop over these prisims)
79 setOfRVEBCPrisms[1].pRisms = periodic_prisms;
80 ierr = mField.add_ents_to_finite_element_by_type(periodic_prisms,MBPRISM,element_name); CHKERRQ(ierr);
81
82
83 Range tris;
84 rval = mField.get_moab().get_adjacencies(periodic_prisms,2,false,tris,moab::Interface::UNION); CHKERRQ_MOAB(rval);
85 setOfRVEBC[1].tRis = tris.subset_by_type(MBTRI);
86
87 PetscFunctionReturn(0);
88 }
89
91 PetscErrorCode shapeMat(
92 int rank, unsigned int gg, DataForcesAndSourcesCore::EntData &col_data, MatrixDouble &N_mat, int div=1
93 ) {
94 PetscFunctionBegin;
95
96 //we need it only for one face for periodic BC (col_data.getN() are shape functions inclding both faces, i.e. 6 shape functions for nodes)
97 int shape_size=col_data.getN().size2()/div;
98 // cout<<"shape_size = "<<shape_size<<endl;
99 // cout<<"rank = "<<rank<<endl;
100 N_mat.resize(rank,shape_size*rank); N_mat.clear();
101 int gg1=0; int kk=0;
102 for(int ii=0; ii<shape_size; ii++){ //number of shape functions
103 for(int jj=0; jj<rank; jj++){
104 // cout<<"ii "<<ii<<endl;
105 N_mat(jj,kk)=col_data.getN()(gg,gg1); kk++;
106 }
107 gg1++;
108 }
109 PetscFunctionReturn(0);
110 }
111
112
113
114 PetscErrorCode shapeMatNew(
115 int rank, unsigned int gg, MatrixDouble &RowN, MatrixDouble &N_mat, int div=1
116 ) {
117 PetscFunctionBegin;
118
119 //we need it only for one face for periodic BC (col_data.getN() are shape functions inclding both faces, i.e. 6 shape functions for nodes)
120 int shape_size=RowN.size2()/div;
121 // cout<<"shape_size = "<<shape_size<<endl;
122 // cout<<"rank = "<<rank<<endl;
123 N_mat.resize(rank,shape_size*rank); N_mat.clear();
124 int gg1=0; int kk=0;
125 for(int ii=0; ii<shape_size; ii++){ //number of shape functions
126 for(int jj=0; jj<rank; jj++){
127 // cout<<"ii "<<ii<<endl;
128 N_mat(jj,kk)=RowN(gg,gg1); kk++;
129 }
130 gg1++;
131 }
132 PetscFunctionReturn(0);
133 }
134
135 };
137
138
140 ublas::vector<MatrixDouble > D_mat;
141 ublas::vector<map<EntityType, map<int, map<EntityType, map<int, MatrixDouble > > > > > Cmat; //Cmat[ff][row_type][row_side][col_type][col_side]
142 ublas::vector<map<EntityType, map<int, ublas::vector<int> > > > RowInd; //RowInd[ff][row_sdie][row_side]
143 ublas::vector<map<EntityType, map<int, ublas::vector<int> > > > ColInd; //ColInd[ff][col_type][col_side]
144
145
146
147 ublas::vector<map<EntityType, map<int, MatrixDouble > > > RowN; //RowN[ff][col_type][col_side]
148 ublas::vector<map<EntityType, map<int, MatrixDouble > > > ColN; //ColN[ff][col_type][col_side]
149
150
151 //these are used in nonlineaer problem to caluclate residual vector
152 ublas::vector<ublas::vector<VectorDouble > > DispAtGaussPts; //DispAtGaussPts[ff][gg] -- different on both faces of prism
153 ublas::vector<VectorDouble > LagMulAtGaussPts; //LagMulAtGaussPts[gg] -- same on both faces of prism
154 };
156
157
158
159
160 //Operator to calculate column indices and colum shapue functions (N)
162
167
169 const string field_name,
171 CommonDataPeriodic &common_data_periodic,
172 CommonFunctionsPeriodic &common_functions_periodic,
173 bool ho_geometry = false
174 ):
177 ),
178 dAta(data),
179 hoGeometry(ho_geometry),
180 commonDataPeriodic(common_data_periodic),
181 commonFunctionsPeriodic(common_functions_periodic) {
182 }
183
184
185 ublas::vector<VectorDouble > field_data_nodes; //field_data_nodes[ff]
187
188
190 PetscFunctionBegin;
191
192 if(type==MBVERTEX){
193 commonDataPeriodic.ColInd.resize(2);
194 commonDataPeriodic.ColN.resize(2);
195 }
196 int nb;
197 //Global indices Columns
198 if(data.getIndices().size()==0) PetscFunctionReturn(0);
199
200 switch(type){
201 case MBVERTEX:
202 nb=data.getIndices().size()/2; //for one face for nodes only (total indieces are for both faces)
203 commonDataPeriodic.ColInd[0][type][side].resize(nb); //1st face of prism -- nodes
204 commonDataPeriodic.ColInd[1][type][side].resize(nb); //2nd face of prism -- nodes
205 //different indices for cols
206 for(int ii=0; ii<nb; ii++){
207 commonDataPeriodic.ColInd[0][type][side][ii]=data.getIndices()[ii];
208 commonDataPeriodic.ColInd[1][type][side][ii]=data.getIndices()[ii+nb];
209 }
210// cout<<"commonDataPeriodic.ColInd[0] "<<"type "<<type << " side "<< side << " = "<< commonDataPeriodic.ColInd[0][type][side]<<endl;
211// cout<<"commonDataPeriodic.ColInd[0] "<<"type "<<type << " side "<< side << " = "<< commonDataPeriodic.ColInd[1][type][side]<<endl;
212// cout<<"type "<<type << " side "<< side <<endl;
213 commonDataPeriodic.ColN[0][type][side]=data.getN();
214 commonDataPeriodic.ColN[1][type][side]=data.getN();
215 break;
216 case MBEDGE:
217 if(side<3){ //1st face eges 0, 1, 2 (cananical numbering)-- same indices for both faces
218 commonDataPeriodic.ColInd[0][type][side]=data.getIndices();
219 commonDataPeriodic.ColN[0][type][side]=data.getN();
220// cout<<"col_type "<<col_type << " col_side "<< col_side <<endl;
221// cout<<"col_data.getIndices "<<col_data.getIndices()<<endl;
222 }
223 else if (side >= 6){ //2nd face edges 6,7,8
224 commonDataPeriodic.ColInd[1][type][side-6]=data.getIndices();
225 commonDataPeriodic.ColN[1][type][side-6]=data.getN();
226// cout<<"col_type "<<col_type << " col_side "<< col_side <<endl;
227 }
228 break;
229 case MBTRI:
230 if(side==3){//1st face face 3
231 commonDataPeriodic.ColInd[0][type][side]=data.getIndices();
232 commonDataPeriodic.ColN[0][type][side]=data.getN();
233 }
234 else{ //2nd face face 4
235 commonDataPeriodic.ColInd[1][type][side-1]=data.getIndices();
236 commonDataPeriodic.ColN[1][type][side-1]=data.getN();
237 }
238 break;
239 default:
240 SETERRQ(PETSC_COMM_SELF,1,"data inconsitency");
241 }
242 PetscFunctionReturn(0);
243 }
244 };
245
246
247
248 //Operator to calculate row indices and row shapue functions (N)
250
255
257 const string field_name,
259 CommonDataPeriodic &common_data_periodic,
260 CommonFunctionsPeriodic &common_functions_periodic,
261 bool ho_geometry = false
262 ):
265 ),
266 dAta(data),
267 hoGeometry(ho_geometry),
268 commonDataPeriodic(common_data_periodic),
269 commonFunctionsPeriodic(common_functions_periodic) {
270 }
271
272
273 ublas::vector<VectorDouble > field_data_nodes; //field_data_nodes[ff]
275
276
278 PetscFunctionBegin;
279// cout<<"type "<<type << " side "<< side <<endl;
280 if(data.getIndices().size()==0) PetscFunctionReturn(0);
281 if(type == MBEDGE && side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
282 if(type == MBTRI && side == 4) PetscFunctionReturn(0); //ignore second triangel face (cananical number 4)
283
284
285 if(type==MBVERTEX){
286 commonDataPeriodic.RowInd.resize(2);
287 commonDataPeriodic.RowN.resize(2);
288 }
289
290 int nb;
291 //Global indices Rows
292 switch(type){
293 case MBVERTEX:
294 nb=data.getIndices().size()/2;
295 commonDataPeriodic.RowInd[0][type][side].resize(nb); //1st face of prism -- nodes
296 commonDataPeriodic.RowInd[1][type][side].resize(nb); //2nd face of prism -- nodes
297
298 //same indices for rows
299 for(int ii=0; ii<nb; ii++){
300 commonDataPeriodic.RowInd[0][type][side][ii]=data.getIndices()[ii];
301 commonDataPeriodic.RowInd[1][type][side][ii]=data.getIndices()[ii];
302 }
303// cout<<"type "<<type << " side "<< side <<endl;
304// cout<<"commonDataPeriodic.RowInd[0] "<<"type "<<type << " side "<< side << " = "<< commonDataPeriodic.RowInd[0][type][side]<<endl;
305 commonDataPeriodic.RowN[0][type][side]=data.getN();
306 commonDataPeriodic.RowN[1][type][side]=data.getN();
307 break;
308 case MBEDGE:
309 if(side<3){ //1st face eges 0, 1, 2 (cananical numbering)-- same indices for both faces
310 commonDataPeriodic.RowInd[0][type][side]=data.getIndices();
311 commonDataPeriodic.RowInd[1][type][side]=data.getIndices();
312 commonDataPeriodic.RowN[0][type][side]=data.getN();
313 commonDataPeriodic.RowN[1][type][side]=data.getN();
314// cout<<"commonDataPeriodic.RowInd[0] "<<"type "<<type << " side "<< side << " = "<< commonDataPeriodic.RowInd[0][type][side]<<endl;
315// cout<<"type "<<type << " side "<< side <<endl;
316 }
317 break;
318 case MBTRI:
319 if(side==3){//1st face face 3
320 commonDataPeriodic.RowInd[0][type][side]=data.getIndices();
321 commonDataPeriodic.RowInd[1][type][side]=data.getIndices();
322 commonDataPeriodic.RowN[0][type][side]=data.getN();
323 commonDataPeriodic.RowN[1][type][side]=data.getN();
324 }
325 break;
326 default:
327 SETERRQ(PETSC_COMM_SELF,1,"data inconsitency");
328 }
329
330// string aaa;
331// cin>>aaa;
332 PetscFunctionReturn(0);
333 }
334 };
335
336
337 /// \biref operator to calculate and assemble Cmat
339 Mat Aij;
344
346 const string field_name,
347 const string lagrang_field_name,
348 Mat aij,
350 CommonFunctionsPeriodic &common_functions_periodic,
351 CommonDataPeriodic &common_data_periodic,
352 bool ho_geometry = false
353 ):
354
356 lagrang_field_name, field_name, UserDataOperator::OPROWCOL
357 ),Aij(aij),
358 dAta(data),
359 commonFunctionsPeriodic(common_functions_periodic),
360 commonDataPeriodic(common_data_periodic),
361 hoGeometry(ho_geometry){
362 sYmm = false; //This will make sure to loop over all intities (e.g. for order=2 it will make doWork to loop 16 time)
363 }
364
366 ublas::vector<int> rrvec, ccvec;
368
369 PetscErrorCode doWork(
370 int row_side,int col_side,
371 EntityType row_type,EntityType col_type,
374 ) {
375 PetscFunctionBegin;
376
377 try {
378// cout<<"Hello from OpRVEBCsPeriodicLhsAssemCmat "<<endl;
379 if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
380 if(col_data.getIndices().size()==0) PetscFunctionReturn(0);
381 if(row_type == MBEDGE && row_side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
382 if(col_type == MBEDGE && col_side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
383 if(row_type == MBTRI && row_side == 4) PetscFunctionReturn(0); //ignore second triangel face (cananical number 4)
384 if(col_type == MBTRI && col_side == 4) PetscFunctionReturn(0); //ignore second triangel face (cananical number 4)
385
386
387 //Cmat calculation
388
389 int rank = col_data.getFieldDofs()[0]->getNbOfCoeffs();
390
391 MatrixDouble N_mat_row;
392 MatrixDouble N_mat_col;
393
394
395 //1st face
396 for(unsigned int gg = 0;gg<col_data.getN().size1();gg++) {
397 double area;
398 VectorDouble normal_f3 = getNormalsAtGaussPtsF3(gg);
399 area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
400 double val = getGaussPts()(2,gg)*area;
401
402 if(row_type == MBVERTEX) {
403 ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.RowN[0][row_type][row_side], N_mat_row, 2); CHKERRQ(ierr);
404// cerr << "N_mat_row "<< N_mat_row << endl;
405 }
406 else {
407 ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.RowN[0][row_type][row_side], N_mat_row); CHKERRQ(ierr);}
408
409 if(col_type == MBVERTEX) {
410 ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.ColN[0][col_type][col_side], N_mat_col, 2); CHKERRQ(ierr);
411// cerr << "N_mat_row "<< N_mat_row << endl;
412 }
413 else {
414 ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.ColN[0][col_type][col_side], N_mat_col); CHKERRQ(ierr);}
415 if(gg==0){
416// cerr << "row_type "<< row_type <<" row_side "<<row_side << " "<< endl;
417// cerr << "col_type "<< col_type <<" col_side "<< col_side<< " "<< endl;
418 NTN=prod(trans(N_mat_row),N_mat_col); //we don't need to define size of NTN
419 Kmat=-1.0*val*NTN; //we don't need to defien size of K
420 }else{
421 NTN=prod(trans(N_mat_row),N_mat_col);
422 Kmat+=-1.0*val*NTN;
423 }
424 }
425 rrvec=commonDataPeriodic.RowInd[0][row_type][row_side];
426 ccvec=commonDataPeriodic.ColInd[0][col_type][col_side];
427 ierr = MatSetValues(Aij,rrvec.size(),&rrvec[0],ccvec.size(),&ccvec[0],&Kmat(0,0),ADD_VALUES); CHKERRQ(ierr);
428 transKmat = trans(Kmat);
429 ierr = MatSetValues(Aij,ccvec.size(),&ccvec[0],rrvec.size(),&rrvec[0],&transKmat(0,0),ADD_VALUES); CHKERRQ(ierr);
430
431// cerr << "rrvec "<< rrvec << endl;
432// cerr << "ccvec "<< ccvec << endl;
433// cerr << "Kmat "<< Kmat << endl;
434
435
436 //2nd face
437 for(unsigned int gg = 0;gg<col_data.getN().size1();gg++) {
438 double area;
439 VectorDouble normal_f3 = getNormalsAtGaussPtsF3(gg);
440 area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
441 double val = getGaussPts()(2,gg)*area;
442
443 if(row_type == MBVERTEX) {
444 ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.RowN[1][row_type][row_side], N_mat_row, 2); CHKERRQ(ierr);
445// cerr << "N_mat_row "<< N_mat_row << endl;
446 }
447 else {
448 ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.RowN[1][row_type][row_side], N_mat_row); CHKERRQ(ierr);}
449
450 if(col_type == MBVERTEX) {
451 ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.ColN[1][col_type][col_side], N_mat_col, 2); CHKERRQ(ierr);
452// cerr << "N_mat_row "<< N_mat_row << endl;
453 }
454 else {
455 ierr = commonFunctionsPeriodic.shapeMatNew(rank, gg, commonDataPeriodic.ColN[1][col_type][col_side], N_mat_col); CHKERRQ(ierr);}
456 if(gg==0){
457// cerr << "row_type "<< row_type <<" row_side "<<row_side << " "<< endl;
458// cerr << "col_type "<< col_type <<" col_side "<< col_side<< " "<< endl;
459 NTN=prod(trans(N_mat_row),N_mat_col); //we don't need to define size of NTN
460 Kmat=val*NTN; //we don't need to defien size of K
461 }else{
462 NTN=prod(trans(N_mat_row),N_mat_col);
463 Kmat+=val*NTN;
464 }
465 }
466 rrvec=commonDataPeriodic.RowInd[1][row_type][row_side];
467 ccvec=commonDataPeriodic.ColInd[1][col_type][col_side];
468 ierr = MatSetValues(Aij,rrvec.size(),&rrvec[0],ccvec.size(),&ccvec[0],&Kmat(0,0),ADD_VALUES); CHKERRQ(ierr);
469 transKmat = trans(Kmat);
470 ierr = MatSetValues(Aij,ccvec.size(),&ccvec[0],rrvec.size(),&rrvec[0],&transKmat(0,0),ADD_VALUES); CHKERRQ(ierr);
471
472
473// cerr << "\n\n\nrrvec "<< rrvec << endl;
474// cerr << "ccvec "<< ccvec << endl;
475// cerr << "Kmat "<< Kmat << endl;
476
477
478 } catch (const std::exception& ex) {
479 ostringstream ss;
480 ss << "throw in method: " << ex.what() << endl;
481 SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
482 }
483// string wait;
484// cin>>wait;
485 PetscFunctionReturn(0);
486 }
487 };
488
489
490
491
492 /// \biref operator to calculate the RHS of the constrain for the RVE boundary conditions
494
499
501 const string field_name,
502 const string lagrang_field_name,
504 CommonDataPeriodic &common_data_periodic,
505 CommonFunctionsPeriodic &common_functions_periodic,
506 bool ho_geometry = false
507 ):
509 lagrang_field_name,UserDataOperator::OPROW
510 ),
511 dAta(data),
512 hoGeometry(ho_geometry),
513 commonDataPeriodic(common_data_periodic),
514 commonFunctionsPeriodic(common_functions_periodic) {
515 }
516
517 ublas::vector<MatrixDouble > X_mat; //Here X_mat are different for left and right face
519
521 PetscFunctionBegin;
522// if(data.getIndices().size()==0) PetscFunctionReturn(0);
523// if(type == MBEDGE && side >= 6) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
524// if(type == MBTRI && side == 4) PetscFunctionReturn(0); //ignore second triangel face (cananical number 4)
525
526 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
527
528 double x,y,z,x1,y1,z1;
529 for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
530 double area;
531 VectorDouble normal_f3 = getNormalsAtGaussPtsF3(gg);
532 area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
533 double val = getGaussPts()(2,gg)*area;
534
535 // Gauss point coordinates of the prism for both faces with HO geometry
536 x = getHOCoordsAtGaussPtsF3()(gg,0);
537 y = getHOCoordsAtGaussPtsF3()(gg,1);
538 z = getHOCoordsAtGaussPtsF3()(gg,2);
539
540 x1 = getHOCoordsAtGaussPtsF4()(gg,0);
541 y1 = getHOCoordsAtGaussPtsF4()(gg,1);
542 z1 = getHOCoordsAtGaussPtsF4()(gg,2);
543
544// cout<<"x "<<x << " y "<< y << " z "<< z <<endl;
545// cout<<"x1 "<<x1 << " y1 "<< y1 << " z1 "<< z1 <<endl;
546
547 commonDataPeriodic.D_mat.resize(2);
548 X_mat.resize(2);
549
550 switch(rank) {
551 case 3: //mech problem
552 X_mat[0].resize(rank,6,false); X_mat[0].clear();
553 X_mat[0](0,0)=2.0*x; X_mat[0](0,3)=y; X_mat[0](0,5)=z;
554 X_mat[0](1,1)=2.0*y; X_mat[0](1,3)=x; X_mat[0](1,4)=z;
555 X_mat[0](2,2)=2.0*z; X_mat[0](2,4)=y; X_mat[0](2,5)=x;
556 X_mat[0]=0.5*X_mat[0];
557
558 X_mat[1].resize(rank,6,false); X_mat[1].clear();
559 X_mat[1](0,0)=2.0*x1; X_mat[1](0,3)=y1; X_mat[1](0,5)=z1;
560 X_mat[1](1,1)=2.0*y1; X_mat[1](1,3)=x1; X_mat[1](1,4)=z1;
561 X_mat[1](2,2)=2.0*z1; X_mat[1](2,4)=y1; X_mat[1](2,5)=x1;
562 X_mat[1]=0.5*X_mat[1];
563 break;
564 case 1: //moisture transport or thermal problem
565 X_mat[0].resize(rank,3,false); X_mat[0].clear();
566 X_mat[0](0,0)=x; X_mat[0](0,1)=y; X_mat[0](0,2)=z;
567
568 X_mat[1].resize(rank,3,false); X_mat[1].clear();
569 X_mat[1](0,0)=x1; X_mat[1](0,1)=y1; X_mat[1](0,2)=z1;
570 break;
571 default:
572 SETERRQ(PETSC_COMM_SELF,1,"not implemented");
573 }
574
575 //
576 if(type == MBVERTEX) {
577 ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat, 2); CHKERRQ(ierr);
578 } else {
579 ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat); CHKERRQ(ierr);
580 }
581
582 if(gg==0) {
583 commonDataPeriodic.D_mat[0]=-1*val*prod(trans(N_mat),X_mat[0]);
584 commonDataPeriodic.D_mat[1]= val*prod(trans(N_mat),X_mat[1]);
585 } else {
586 commonDataPeriodic.D_mat[0]+=-1*val*prod(trans(N_mat),X_mat[0]);
587 commonDataPeriodic.D_mat[1]+= val*prod(trans(N_mat),X_mat[1]);
588 }
589 }
590
591// cout<<endl<<endl;
592
593 PetscFunctionReturn(0);
594 }
595 };
596
598
599 vector<Vec> &F;
600
602 const string field_name,
603 const string lagrang_field_name,
604 vector<Vec> &f,
606 CommonDataPeriodic &common_data_periodic,
607 CommonFunctionsPeriodic &common_functions_periodic,
608 bool ho_geometry = false
609 ):
610 OpDmatRhs(
611 field_name,lagrang_field_name,data,common_data_periodic,common_functions_periodic,ho_geometry
612 ),
613 F(f) {
614 }
615
618
620 PetscFunctionBegin;
621 if(data.getIndices().size()==0) PetscFunctionReturn(0);
622 if(type == MBEDGE && side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
623 if(type == MBTRI && side == 4) PetscFunctionReturn(0); //ignore second triangel face (cananical number 4)
624
625
626 ierr = calculateDmat(side,type,data); CHKERRQ(ierr);
627
628// if(type == MBVERTEX) {
629// cout<<"\n\nNew element "<<endl;
630// }
631// cout<<"type "<<type << " side "<< side <<endl;
632// cout<< "commonDataPeriodic.D_mat[0] "<<commonDataPeriodic.D_mat[0]<<endl;
633// cout<< "commonDataPeriodic.D_mat[1] "<<commonDataPeriodic.D_mat[1]<<endl<<endl;
634// string aaa;
635// cin >> aaa;
636
637
638 ublas::vector<int> rowvec;
639 if(type == MBVERTEX) {
640 int nb=data.getIndices().size()/2; //for one face
641 rowvec.resize(nb);
642 for(int ii=0; ii<nb; ii++){
643 rowvec[ii]=data.getIndices()[ii];
644 }
645 } else{
646 rowvec=data.getIndices();
647 }
648
649 appliedStrain.resize(commonDataPeriodic.D_mat[0].size2(),false);
650
651 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
652 int size = (rank == 3) ? 6 : 3;
653 for(int ii = 0;ii<size;ii++) {
654 appliedStrain.clear();
655 appliedStrain[ii] = 1.0;
657 ierr = VecSetValues(F[ii],rowvec.size(),&rowvec[0],&nF[0],ADD_VALUES); CHKERRQ(ierr);
659 ierr = VecSetValues(F[ii],rowvec.size(),&rowvec[0],&nF[0],ADD_VALUES); CHKERRQ(ierr);
660 }
661
662 PetscFunctionReturn(0);
663 }
664 };
665
666
667
668
670
673 boost::ptr_vector<MethodForForceScaling> &methodsOp;
674
676 const string field_name,
677 const string lagrang_field_name,
678 Vec f,
679 VectorDouble given_strain,
680 boost::ptr_vector<MethodForForceScaling> &methods_op,
682 CommonDataPeriodic &common_data_periodic,
683 CommonFunctionsPeriodic &common_functions_periodic,
684 bool ho_geometry = false
685 ):
686 OpDmatRhs(
687 field_name,lagrang_field_name,data,common_data_periodic,common_functions_periodic,ho_geometry
688 ),
689 F(f), givenStrain(given_strain),methodsOp(methods_op){
690 }
691
693
695 PetscFunctionBegin;
696 if(data.getIndices().size()==0) PetscFunctionReturn(0);
697 if(type == MBEDGE && side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
698 if(type == MBTRI && side == 4) PetscFunctionReturn(0); //ignore second triangel face (cananical number 4)
699
700
701 ierr = calculateDmat(side,type,data); CHKERRQ(ierr);
702
703 ublas::vector<int> rowvec;
704 if(type == MBVERTEX) {
705 int nb=data.getIndices().size()/2; //for one face
706 rowvec.resize(nb);
707 for(int ii=0; ii<nb; ii++){
708 rowvec[ii]=data.getIndices()[ii];
709 }
710 } else{
711 rowvec=data.getIndices();
712 }
713
714 //This is to scale givenStrain vector to apply the strian in steps
715 VectorDouble scaled_given_strain;
716 scaled_given_strain.resize(6);
717 scaled_given_strain=givenStrain;
718// cout <<"givenStrain Before = "<<scaled_given_strain<<endl;
719 ierr = MethodForForceScaling::applyScale(getFEMethod(),methodsOp,scaled_given_strain); CHKERRQ(ierr);
720// cout <<"givenStrain After = "<<scaled_given_strain<<endl;
721
722// SNES snes = getFEMethod()->snes;
723// int iter;
724// SNESGetIterationNumber(snes,&iter);
725// if(iter>0) PetscFunctionReturn(0);
726
727 nF = prod(commonDataPeriodic.D_mat[0], scaled_given_strain);
728 nF*=-1;
729 ierr = VecSetValues(F,rowvec.size(),&rowvec[0],&nF[0],ADD_VALUES); CHKERRQ(ierr);
730
731 nF = prod(commonDataPeriodic.D_mat[1], scaled_given_strain);
732 nF*=-1;
733 ierr = VecSetValues(F,rowvec.size(),&rowvec[0],&nF[0],ADD_VALUES); CHKERRQ(ierr);
734
735 PetscFunctionReturn(0);
736 }
737 };
738
739
741
746
748 const string field_name,
750 CommonDataPeriodic &common_data_periodic,
751 CommonFunctionsPeriodic &common_functions_periodic,
752 bool ho_geometry = false
753 ):
756 ),
757 dAta(data),
758 hoGeometry(ho_geometry),
759 commonDataPeriodic(common_data_periodic),
760 commonFunctionsPeriodic(common_functions_periodic) {
761 }
762
763
764 ublas::vector<VectorDouble > field_data_nodes; //field_data_nodes[ff]
766
767
769 PetscFunctionBegin;
770// cout<<" OpRVEBCsPeriodicCalDispAtGaussPts "<<endl;
771 if(data.getIndices().size()==0) PetscFunctionReturn(0);
772 int nb_gauss_pts = data.getN().size1();
773 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
774
775 //initialize
776 if(type == MBVERTEX) {
778 field_data_nodes.resize(2);
779
780 commonDataPeriodic.DispAtGaussPts[0].resize(nb_gauss_pts);
781 commonDataPeriodic.DispAtGaussPts[1].resize(nb_gauss_pts);
782 for(int gg = 0;gg<nb_gauss_pts;gg++) {
783 commonDataPeriodic.DispAtGaussPts[0][gg].resize(rank); commonDataPeriodic.DispAtGaussPts[0][gg].clear();
784 commonDataPeriodic.DispAtGaussPts[1][gg].resize(rank); commonDataPeriodic.DispAtGaussPts[1][gg].clear();
785 }
786
787 int nb=data.getFieldData().size()/2; //for one face for nodes only (total indieces are for both faces)
788 field_data_nodes[0].resize(nb); //1st face of prism -- nodes
789 field_data_nodes[1].resize(nb); //2nd face of prism -- nodes
790
791 for(int ii=0; ii<nb; ii++){
792 field_data_nodes[0][ii]=data.getFieldData()[ii];
793 field_data_nodes[1][ii]=data.getFieldData()[ii+nb];
794 }
795// cout<<"field_data_nodes[0] "<<field_data_nodes[0]<<endl;
796// cout<<"field_data_nodes[1] "<<field_data_nodes[1]<<endl;
797 }
798
799 for(unsigned int gg = 0; gg<nb_gauss_pts;gg++) {
800 if(type == MBVERTEX) {
801 ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat, 2); CHKERRQ(ierr);}
802 else {
803 ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat); CHKERRQ(ierr);}
804
805 if(type == MBVERTEX) {
808 else if(type == MBEDGE && side < 3){
809 commonDataPeriodic.DispAtGaussPts[0][gg]+=prod(N_mat, data.getFieldData());}
810 else if(type == MBEDGE && side >= 6){
811 commonDataPeriodic.DispAtGaussPts[1][gg]+=prod(N_mat, data.getFieldData());}
812 else if(type == MBTRI && side == 3){
813 commonDataPeriodic.DispAtGaussPts[0][gg]+=prod(N_mat, data.getFieldData());}
814 else if(type == MBTRI && side == 4){
815 commonDataPeriodic.DispAtGaussPts[1][gg]+=prod(N_mat, data.getFieldData());}
816 }
817
818 PetscFunctionReturn(0);
819 }
820 };
821
822
823
824
826
832
834 const string field_name,
836 Vec f,
837 CommonDataPeriodic &common_data_periodic,
838 CommonFunctionsPeriodic &common_functions_periodic,
839 bool ho_geometry = false
840 ):
843 ),
844 F(f),
845 dAta(data),
846 hoGeometry(ho_geometry),
847 commonDataPeriodic(common_data_periodic),
848 commonFunctionsPeriodic(common_functions_periodic) {
849 }
850
851
852 ublas::vector<int> row_ind;
855
856
858 PetscFunctionBegin;
859// cout<<" OpRVEBCsPeriodicCalCU "<<endl;
860 if(data.getIndices().size()==0) PetscFunctionReturn(0);
861 if(type == MBQUAD) PetscFunctionReturn(0);
862 if(type == MBEDGE && side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
863 if(type == MBTRI && side == 4) PetscFunctionReturn(0); //ignore second triangel face (cananical number 4)
864
865 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
866
867 if(type == MBVERTEX) {
868 int nb=data.getIndices().size()/2; //for one face for nodes only (total indieces are for both faces)
869 row_ind.resize(nb); //1st face of prism -- nodes
870 for(int ii=0; ii<nb; ii++){
871 row_ind[ii]=data.getIndices()[ii];
872 }
873 }
874
875 for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
876 double area;
877 VectorDouble normal_f3 = getNormalsAtGaussPtsF3(gg);
878 area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
879 double val = getGaussPts()(2,gg)*area;
880
881 if(type == MBVERTEX) {
882 ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat, 2); CHKERRQ(ierr);}
883 else {
884 ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat); CHKERRQ(ierr);}
885
886 if(gg==0){
887 fe1=-val*prod(trans(N_mat), commonDataPeriodic.DispAtGaussPts[0][gg]);
888 fe2= val*prod(trans(N_mat), commonDataPeriodic.DispAtGaussPts[1][gg]);}
889 else{
890 fe1+=-val*prod(trans(N_mat), commonDataPeriodic.DispAtGaussPts[0][gg]);
891 fe2+= val*prod(trans(N_mat), commonDataPeriodic.DispAtGaussPts[1][gg]);}
892 }
893
894 if(type == MBVERTEX) {
895 ierr = VecSetValues(F,row_ind.size(),&row_ind[0],&fe1[0],ADD_VALUES); CHKERRQ(ierr);
896 ierr = VecSetValues(F,row_ind.size(),&row_ind[0],&fe2[0],ADD_VALUES); CHKERRQ(ierr);}
897 else{
898 ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&fe1[0],ADD_VALUES); CHKERRQ(ierr);
899 ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&fe2[0],ADD_VALUES); CHKERRQ(ierr);}
900
901 PetscFunctionReturn(0);
902 }
903 };
904
905
906
908
913
915 const string field_name,
917 CommonDataPeriodic &common_data_periodic,
918 CommonFunctionsPeriodic &common_functions_periodic,
919 bool ho_geometry = false
920 ):
923 ),
924 dAta(data),
925 hoGeometry(ho_geometry),
926 commonDataPeriodic(common_data_periodic),
927 commonFunctionsPeriodic(common_functions_periodic) {
928 }
929
930
933
934
935
937 PetscFunctionBegin;
938
939
940 if(data.getIndices().size()==0) PetscFunctionReturn(0);
941 if(type == MBEDGE && side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
942 if(type == MBTRI && side == 4) PetscFunctionReturn(0); //ignore second triangel face (cananical number 4)
943
944 int nb_gauss_pts = data.getN().size1();
945 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
946
947 //initialize
948 if(type == MBVERTEX) {
949 commonDataPeriodic.LagMulAtGaussPts.resize(nb_gauss_pts);
950 for(int gg = 0;gg<nb_gauss_pts;gg++) {
952 }
953
954 int nb=data.getFieldData().size()/2; //for one face for nodes only (total indieces are for both faces)
955 field_data_nodes.resize(nb); //1st face of prism -- nodes
956
957 for(int ii=0; ii<nb; ii++){
958 field_data_nodes[ii]=data.getFieldData()[ii];
959 }
960 }
961
962
963 for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
964 if(type == MBVERTEX) {
965 ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat, 2); CHKERRQ(ierr);}
966 else {
967 ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat); CHKERRQ(ierr);}
968
969 if(type == MBVERTEX) {
971 else {
972 commonDataPeriodic.LagMulAtGaussPts[gg]+=prod(N_mat, data.getFieldData());}
973 }
974
975 PetscFunctionReturn(0);
976 }
977 };
978
979
980
982
988
990 const string field_name,
992 Vec f,
993 CommonDataPeriodic &common_data_periodic,
994 CommonFunctionsPeriodic &common_functions_periodic,
995 bool ho_geometry = false
996 ):
999 ),
1000 F(f),
1001 dAta(data),
1002 hoGeometry(ho_geometry),
1003 commonDataPeriodic(common_data_periodic),
1004 commonFunctionsPeriodic(common_functions_periodic) {
1005 }
1006
1007
1008 ublas::vector<ublas::vector<int> > col_ind; //indices for two faces
1010
1012
1014 PetscFunctionBegin;
1015// cout<<" OpRVEBCsPeriodicCalCTLam "<<endl;
1016 if(data.getIndices().size()==0) PetscFunctionReturn(0);
1017 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
1018 if(type == MBVERTEX) {
1019 col_ind.resize(2);
1020 int nb=data.getIndices().size()/2; //for one face for nodes only (total indieces are for both faces)
1021 col_ind[0].resize(nb); //1st face of prism -- nodes
1022 col_ind[1].resize(nb); //2nd face of prism -- nodes
1023 for(int ii=0; ii<nb; ii++){
1024 col_ind[0][ii]=data.getIndices()[ii];
1025 col_ind[1][ii]=data.getIndices()[ii+nb];
1026 }
1027 }
1028
1029
1030// cout<<"data.getN().size1() "<<data.getN().size1()<<endl;
1031 for(unsigned int gg = 0;gg<data.getN().size1();gg++) {
1032// cout<<"gg "<<gg<<endl;
1033 double area;
1034 VectorDouble normal_f3 = getNormalsAtGaussPtsF3(gg);
1035 area = cblas_dnrm2(3,&normal_f3[0],1)*0.5;
1036 double val = getGaussPts()(2,gg)*area;
1037
1038 if(type == MBVERTEX) {
1039 ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat, 2); CHKERRQ(ierr);}
1040 else {
1041 ierr = commonFunctionsPeriodic.shapeMat(rank, gg, data, N_mat); CHKERRQ(ierr);}
1042
1043
1044 if(gg==0){
1045 if(type == MBVERTEX) {
1046 fe1=-val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);
1047 fe2= val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
1048 else if(type == MBEDGE && side < 3){
1049 fe1=-val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
1050 else if(type == MBEDGE && side >= 6){
1051 fe2= val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
1052 else if(type == MBTRI && side == 3){
1053 fe1=-val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
1054 else if(type == MBTRI && side == 4){
1055 fe2= val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
1056 }else{
1057 if(type == MBVERTEX) {
1058 fe1+=-val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);
1059 fe2+= val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
1060 else if(type == MBEDGE && side < 3){
1061 fe1+=-val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
1062 else if(type == MBEDGE && side >= 6){
1063 fe2+= val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
1064 else if(type == MBTRI && side == 3){
1065 fe1+=-val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
1066 else if(type == MBTRI && side == 4){
1067 fe2+= val*prod(trans(N_mat), commonDataPeriodic.LagMulAtGaussPts[gg]);}
1068 }
1069
1070 }
1071// cout<<" Before assembley of F"<<endl;
1072
1073 if(type == MBVERTEX) {
1074 ierr = VecSetValues(F,col_ind[0].size(),&col_ind[0][0],&fe1[0],ADD_VALUES); CHKERRQ(ierr);
1075 ierr = VecSetValues(F,col_ind[1].size(),&col_ind[1][0],&fe2[0],ADD_VALUES); CHKERRQ(ierr);}
1076 else if(type == MBEDGE && side < 3){
1077 ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&fe1[0],ADD_VALUES); CHKERRQ(ierr);}
1078 else if(type == MBEDGE && side >= 6){
1079 ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&fe2[0],ADD_VALUES); CHKERRQ(ierr);}
1080 else if(type == MBTRI && side == 3){
1081 ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&fe1[0],ADD_VALUES); CHKERRQ(ierr);}
1082 else if(type == MBTRI && side == 4){
1083 ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&fe2[0],ADD_VALUES); CHKERRQ(ierr);}
1084 PetscFunctionReturn(0);
1085 }
1086 };
1087
1088
1089
1090 //for nonlinear problems
1092 string field_name,
1093 string lagrang_field_name,
1094 string mesh_nodals_positions,
1095 Mat aij,vector<Vec> &fvec,Vec f,VectorDouble given_strain
1096 ) {
1097 PetscFunctionBegin;
1098
1099 bool ho_geometry = false;
1100 if(mField.check_field(mesh_nodals_positions)) {
1101 ho_geometry = true;
1102 }
1103
1104 map<int,RVEBC_Data_Periodic>::iterator sit = setOfRVEBCPrisms.begin();
1105 for(;sit!=setOfRVEBCPrisms.end();sit++) {
1106
1107 //Col indices and Col N
1108 feRVEBCLhs.getOpPtrVector().push_back(
1111 )
1112 );
1113
1114
1115 //Row indices and Row N
1116 feRVEBCLhs.getOpPtrVector().push_back(
1118 lagrang_field_name,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
1119 )
1120 );
1121
1122 //Calculate and assemble Cmat
1123 feRVEBCLhs.getOpPtrVector().push_back(
1125 field_name,lagrang_field_name,aij,sit->second,commonFunctionsPeriodic, commonDataPeriodic,ho_geometry
1126 )
1127 );
1128
1129 //RHS
1130 feRVEBCRhs.getOpPtrVector().push_back(
1132 field_name,lagrang_field_name,f,given_strain, methodsOp, sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
1133 )
1134 );
1135
1136
1137 //Calculate displacement at gauss point
1138 feRVEBCRhsResidual.getOpPtrVector().push_back(
1141 )
1142 );
1143
1144 //Calculate and assemble CU = NT * u
1145 feRVEBCRhsResidual.getOpPtrVector().push_back(
1147 lagrang_field_name,sit->second,f,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
1148 )
1149 );
1150
1151
1152 //Calculate lagrange multipliers at gauss point
1153 feRVEBCRhsResidual.getOpPtrVector().push_back(
1155 lagrang_field_name,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
1156 )
1157 );
1158
1159 //Calculate and assemble CT Lam = NT * Lam
1160 feRVEBCRhsResidual.getOpPtrVector().push_back(
1163 )
1164 );
1165
1166// cout<<"Hi 1 from setRVEBCsOperators periodic "<<endl;
1167// string aaa;
1168// cin>>aaa;
1169 }
1170 PetscFunctionReturn(0);
1171 }
1172
1173
1174
1175 //for linear problems
1176 PetscErrorCode setRVEBCsOperators(
1177 string field_name,string lagrang_field_name,string mesh_nodals_positions,Mat aij,vector<Vec> &f
1178 ) {
1179 PetscFunctionBegin;
1180
1181 bool ho_geometry = false;
1182 if(mField.check_field(mesh_nodals_positions)) {
1183 ho_geometry = true;
1184 }
1185 // cout<<"Hi 1 from setRVEBCsOperators periodic "<<endl;
1186 map<int,RVEBC_Data_Periodic>::iterator sit = setOfRVEBCPrisms.begin();
1187 for(;sit!=setOfRVEBCPrisms.end();sit++) {
1188
1189 // LHS
1190 //Col indices and Col N
1191 feRVEBCLhs.getOpPtrVector().push_back(
1194 )
1195 );
1196
1197
1198 //Row indices and Row N
1199 feRVEBCLhs.getOpPtrVector().push_back(
1201 lagrang_field_name,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
1202 )
1203 );
1204
1205 //Calculate and assemble Cmat
1206 feRVEBCLhs.getOpPtrVector().push_back(
1208 field_name,lagrang_field_name,aij,sit->second,commonFunctionsPeriodic, commonDataPeriodic,ho_geometry
1209 )
1210 );
1211
1212
1213
1214 //RHS
1215 feRVEBCRhs.getOpPtrVector().push_back(
1217 field_name,lagrang_field_name,f,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
1218 )
1219 );
1220
1221 }
1222 PetscFunctionReturn(0);
1223 }
1224
1225
1226 /// \biref operator to calculate the RVE homogenised stress
1228
1230
1232 const string field_name,
1233 const string lagrang_field_name,
1234 Vec stress_homo,
1235 RVEBC_Data_Periodic &data,
1236 CommonDataPeriodic &common_data_periodic,
1237 CommonFunctionsPeriodic &common_functions_periodic,
1238 bool ho_geometry = false
1239 ):
1240 OpDmatRhs(
1241 field_name,lagrang_field_name,data,common_data_periodic,common_functions_periodic,ho_geometry
1242 ),
1243 Stress_Homo(stress_homo) {
1244 }
1245
1246 ublas::vector<VectorDouble > Stress_Homo_elem;
1248
1250 PetscFunctionBegin;
1251
1252 if(data.getIndices().size()==0) PetscFunctionReturn(0);
1253 if(type == MBEDGE && side >= 3) PetscFunctionReturn(0); //ignore second triangel edges (cananical number 6 to 8)
1254 if(type == MBTRI && side == 4) PetscFunctionReturn(0); //ignore second triangel face (cananical number 4)
1255
1256
1257 ierr = calculateDmat(side,type,data); CHKERRQ(ierr);
1258
1259// cerr << commonDataPeriodic.D_mat[0] << endl;
1260// cerr << data.getFieldData() << endl;
1261// cerr << rowFieldName << endl;
1262// cerr << data.getFieldDofs()[0]->getName() << endl;
1263// cerr << data.getFieldDofs()[0]->getNbOfCoeffs() << endl;
1264// cerr << data.getFieldDofs().size() << endl;
1265// cerr << data.getIndices() << endl;
1266// */
1267
1268
1269 Stress_Homo_elem.resize(2);
1270 if(type == MBVERTEX) {
1271 int nb=data.getFieldData().size()/2; //for one face for nodes only (total indieces are for both faces)
1272 field_data.resize(nb);
1273
1274 for(int ii=0; ii<nb; ii++){
1275 field_data[ii]=data.getFieldData()[ii];
1276 }
1277 Stress_Homo_elem[0] = prod(trans(commonDataPeriodic.D_mat[0]), -1*field_data); //Lamda=data.getFieldData() is reaction force (so multiply for -1 to get the force)
1278 Stress_Homo_elem[1] = prod(trans(commonDataPeriodic.D_mat[1]), -1*field_data); //Lamda=data.getFieldData() is reaction force (so multiply for -1 to get the force)
1279 }
1280 else{
1281 Stress_Homo_elem[0] = prod(trans(commonDataPeriodic.D_mat[0]), -1*data.getFieldData()); //Lamda=data.getFieldData() is reaction force (so multiply for -1 to get the force)
1282 Stress_Homo_elem[1] = prod(trans(commonDataPeriodic.D_mat[1]), -1*data.getFieldData()); //Lamda=data.getFieldData() is reaction force (so multiply for -1 to get the force)
1283 }
1284
1285
1286// cerr << "field_data "<<field_data << endl;
1287// cerr << "Stress_Homo_elem[0] "<<Stress_Homo_elem[0] << endl;
1288// cerr << "Stress_Homo_elem[1] "<<Stress_Homo_elem[1] << endl;
1289// cerr << "commonDataPeriodic.D_mat[0] "<<commonDataPeriodic.D_mat[0] << endl;
1290
1291 int Indices6[6]={0, 1, 2, 3, 4, 5};
1292 int Indices3[3]={0, 1, 2};
1293
1294 int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
1295 switch(rank) {
1296 case 3:
1297 ierr = VecSetValues(Stress_Homo,6,Indices6,&Stress_Homo_elem[0][0],ADD_VALUES); CHKERRQ(ierr);
1298 ierr = VecSetValues(Stress_Homo,6,Indices6,&Stress_Homo_elem[1][0],ADD_VALUES); CHKERRQ(ierr);
1299 break;
1300 case 1:
1301 ierr = VecSetValues(Stress_Homo,3,Indices3,&Stress_Homo_elem[0][0],ADD_VALUES); CHKERRQ(ierr);
1302 ierr = VecSetValues(Stress_Homo,3,Indices3,&Stress_Homo_elem[1][0],ADD_VALUES); CHKERRQ(ierr);
1303 break;
1304 default:
1305 SETERRQ(PETSC_COMM_SELF,1,"not implemented");
1306 }
1307
1308 PetscFunctionReturn(0);
1309 }
1310 };
1311
1313 string field_name,
1314 string lagrang_field_name,
1315 string mesh_nodals_positions,
1316 Vec stress_homo
1317 ) {
1318 PetscFunctionBegin;
1319 bool ho_geometry = false;
1320 if(mField.check_field(mesh_nodals_positions)) {
1321 ho_geometry = true;
1322 }
1323 map<int,RVEBC_Data_Periodic>::iterator sit = setOfRVEBCPrisms.begin();
1324 for(;sit!=setOfRVEBCPrisms.end();sit++) {
1325 feRVEBCStress.getOpPtrVector().push_back(
1326 new OpRVEHomoStress(
1327 field_name,lagrang_field_name,stress_homo,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
1328 )
1329 );
1330 }
1331 PetscFunctionReturn(0);
1332 }
1333
1334};
1335
1336#endif
ForcesAndSourcesCore::UserDataOperator UserDataOperator
EntitiesFieldData::EntData EntData
@ MF_ZERO
Definition: definitions.h:98
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
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
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
constexpr auto field_name
map< int, RVEBC_Data > setOfRVEBC
maps side set id with appropriate FluxData
MoFEM::Interface & mField
ublas::vector< map< EntityType, map< int, ublas::vector< int > > > > RowInd
ublas::vector< map< EntityType, map< int, MatrixDouble > > > ColN
ublas::vector< map< EntityType, map< int, MatrixDouble > > > RowN
ublas::vector< map< EntityType, map< int, map< EntityType, map< int, MatrixDouble > > > > > Cmat
ublas::vector< ublas::vector< VectorDouble > > DispAtGaussPts
ublas::vector< map< EntityType, map< int, ublas::vector< int > > > > ColInd
PetscErrorCode shapeMatNew(int rank, unsigned int gg, MatrixDouble &RowN, MatrixDouble &N_mat, int div=1)
PetscErrorCode shapeMat(int rank, unsigned int gg, DataForcesAndSourcesCore::EntData &col_data, MatrixDouble &N_mat, int div=1)
\biref operator to calculate the RHS of the constrain for the RVE boundary conditions
OpDmatRhs(const string field_name, const string lagrang_field_name, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
PetscErrorCode calculateDmat(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
\biref operator to calculate and assemble Cmat
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
OpRVEBCsPeriodicCalAssemCmat(const string field_name, const string lagrang_field_name, Mat aij, RVEBC_Data_Periodic &data, CommonFunctionsPeriodic &common_functions_periodic, CommonDataPeriodic &common_data_periodic, bool ho_geometry=false)
OpRVEBCsPeriodicCalCTLam(const string field_name, RVEBC_Data_Periodic &data, Vec f, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpRVEBCsPeriodicCalCU(const string field_name, RVEBC_Data_Periodic &data, Vec f, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpRVEBCsPeriodicCalDispAtGaussPts(const string field_name, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
OpRVEBCsPeriodicCalLagMulAtGaussPts(const string field_name, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpRVEBCsPeriodicColInd(const string field_name, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpRVEBCsPeriodicRhs_givenStrain(const string field_name, const string lagrang_field_name, Vec f, VectorDouble given_strain, boost::ptr_vector< MethodForForceScaling > &methods_op, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
OpRVEBCsPeriodicRhs(const string field_name, const string lagrang_field_name, vector< Vec > &f, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpRVEBCsPeriodicRowInd(const string field_name, RVEBC_Data_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
\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_Periodic &data, CommonDataPeriodic &common_data_periodic, CommonFunctionsPeriodic &common_functions_periodic, bool ho_geometry=false)
map< int, RVEBC_Data_Periodic > setOfRVEBCPrisms
maps side set id with appropriate FluxData
boost::ptr_vector< MethodForForceScaling > methodsOp
PetscErrorCode setRVEBCsHomoStressOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Vec stress_homo)
CommonFunctionsPeriodic commonFunctionsPeriodic
BCs_RVELagrange_Periodic(MoFEM::Interface &m_field)
PetscErrorCode addLagrangiangElement(const string element_name, const string field_name, const string lagrang_field_name, const string mesh_nodals_positions, Range &periodic_prisms)
PetscErrorCode setRVEBCsOperators(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &f)
PetscErrorCode setRVEBCsOperatorsNonlinear(string field_name, string lagrang_field_name, string mesh_nodals_positions, Mat aij, vector< Vec > &fvec, Vec f, VectorDouble given_strain)
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
virtual moab::Interface & get_moab()=0
Deprecated interface functions.