v0.13.0
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);
65  ierr = mField.modify_finite_element_add_field_col(element_name,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);
68  ierr = mField.modify_finite_element_add_field_row(element_name,field_name); CHKERRQ(ierr);
69  //data
70  ierr = mField.modify_finite_element_add_field_data(element_name,lagrang_field_name); CHKERRQ(ierr);
71  ierr = mField.modify_finite_element_add_field_data(element_name,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,
170  RVEBC_Data_Periodic &data,
171  CommonDataPeriodic &common_data_periodic,
172  CommonFunctionsPeriodic &common_functions_periodic,
173  bool ho_geometry = false
174  ):
176  field_name,UserDataOperator::OPCOL
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 
189  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
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,
258  RVEBC_Data_Periodic &data,
259  CommonDataPeriodic &common_data_periodic,
260  CommonFunctionsPeriodic &common_functions_periodic,
261  bool ho_geometry = false
262  ):
264  field_name,UserDataOperator::OPROW
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 
277  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
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,
349  RVEBC_Data_Periodic &data,
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,
503  RVEBC_Data_Periodic &data,
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 
520  PetscErrorCode calculateDmat(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
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,
605  RVEBC_Data_Periodic &data,
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 
619  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
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,
681  RVEBC_Data_Periodic &data,
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 
694  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
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,
749  RVEBC_Data_Periodic &data,
750  CommonDataPeriodic &common_data_periodic,
751  CommonFunctionsPeriodic &common_functions_periodic,
752  bool ho_geometry = false
753  ):
755  field_name,UserDataOperator::OPCOL
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 
768  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
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,
835  RVEBC_Data_Periodic &data,
836  Vec f,
837  CommonDataPeriodic &common_data_periodic,
838  CommonFunctionsPeriodic &common_functions_periodic,
839  bool ho_geometry = false
840  ):
842  field_name,UserDataOperator::OPROW
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 
857  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
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,
916  RVEBC_Data_Periodic &data,
917  CommonDataPeriodic &common_data_periodic,
918  CommonFunctionsPeriodic &common_functions_periodic,
919  bool ho_geometry = false
920  ):
922  field_name,UserDataOperator::OPROW
923  ),
924  dAta(data),
925  hoGeometry(ho_geometry),
926  commonDataPeriodic(common_data_periodic),
927  commonFunctionsPeriodic(common_functions_periodic) {
928  }
929 
930 
933 
934 
935 
936  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
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,
991  RVEBC_Data_Periodic &data,
992  Vec f,
993  CommonDataPeriodic &common_data_periodic,
994  CommonFunctionsPeriodic &common_functions_periodic,
995  bool ho_geometry = false
996  ):
998  field_name,UserDataOperator::OPCOL
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 
1013  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
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(
1110  field_name,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
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(
1140  field_name,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
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(
1162  field_name,sit->second,f,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
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(
1193  field_name,sit->second,commonDataPeriodic,commonFunctionsPeriodic,ho_geometry
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(
1216  new OpRVEBCsPeriodicRhs(
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
1227  struct OpRVEHomoStress: public OpDmatRhs {
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 
1249  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
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:111
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
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
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
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.