v0.14.0
NitschePeriodicMethod.hpp
Go to the documentation of this file.
1 /** \file NitschePeriodicMethod.hpp
2  * \ingroup nitsche_method
3  * \brief Implementation of Nitsche's method for periodic constrains.
4  */
5 
6 /*
7  * This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 #ifndef __NITSCHE_PERIODIC_METHOD_HPP__
22 #define __NITSCHE_PERIODIC_METHOD_HPP__
23 
24 /** \brief Periodic Nitsche method
25 * \ingroup nitsche_method
26 */
28 
29  struct CommonData {
30 
33  double eps;
34 
35  map<EntityHandle,vector<VectorDouble > > localCoordsMap;
36  map<EntityHandle,vector<int> > inTetFaceGaussPtsNumber;
37  map<EntityHandle,vector<int> > inTetTetGaussPtsNumber;
38  map<EntityHandle,double> dIstance;
39  vector<VectorDouble > dIsplacements;
40  vector<VectorDouble > coordsAtGaussPts;
41  vector<VectorDouble > hoCoordsAtGaussPts;
42  vector<MatrixDouble > stressJacobian;
43  vector<MatrixDouble > sTress;
44 
45  struct MultiIndexData {
46  int gG;
47  int sIde;
48  EntityType tYpe;
49  ublas::vector<int> iNdices;
50  ublas::vector<int> dofOrders;
54  int gg,int side,EntityType type
55  ):
56  gG(gg),
57  sIde(side),
58  tYpe(type) {
59  }
60  };
61 
62  typedef multi_index_container<
63  MultiIndexData,
64  indexed_by<
65  ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(MultiIndexData,int,MultiIndexData::gG) >,
66  ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(MultiIndexData,int,MultiIndexData::sIde) >,
67  ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(MultiIndexData,EntityType,MultiIndexData::tYpe) >,
68  ordered_unique<
69  composite_key<
70  MultiIndexData,
71  member<MultiIndexData,int,&MultiIndexData::gG>,
72  member<MultiIndexData,int,&MultiIndexData::sIde>,
73  member<MultiIndexData,EntityType,&MultiIndexData::tYpe>
74  > >
75  >
79 
80  };
81 
82 
84 
86  bool fieldDisp;
87 
88  OpGetFaceData(CommonData &common_data,bool field_disp = true):
90  commonData(common_data),
91  fieldDisp(field_disp) {
92  }
93 
94  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
95  PetscFunctionBegin;
96 
97  if(data.getIndices().size()==0) PetscFunctionReturn(0);
98  try {
99  EntityHandle this_face = getNumeredEntFiniteElementPtr()->getEnt();
100  int nb_gauss_pts_on_this_face = data.getN().size1();
101  for(int ffgg = 0;ffgg<nb_gauss_pts_on_this_face;ffgg++) {
102  int gg = commonData.inTetFaceGaussPtsNumber[this_face][ffgg];
103  CommonData::MultiIndexData gauss_pt_data(gg,side,type);
104  pair<CommonData::Container::iterator,bool> p;
105  p = commonData.facesContainer.insert(gauss_pt_data);
106  if(!p.second) {
107  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"data not inserted");
108  }
109  CommonData::MultiIndexData &p_data = const_cast<CommonData::MultiIndexData&>(*p.first);
110  VectorDouble &shape_fun = p_data.shapeFunctions;
111  int nb_shape_fun = data.getN().size2();
112  shape_fun.resize(nb_shape_fun);
113  cblas_dcopy(nb_shape_fun,&data.getN()(ffgg,0),1,&shape_fun[0],1);
114  p_data.iNdices = data.getIndices();
115  p_data.dofOrders.resize(data.getFieldDofs().size(),false);
116  for(unsigned int dd = 0;dd<data.getFieldDofs().size();dd++) {
117  p_data.dofOrders[dd] = data.getFieldDofs()[dd]->getDofOrder();
118  }
119  int nb_dofs = data.getFieldData().size();
120  if(type == MBVERTEX) {
121  commonData.dIsplacements[gg].resize(3,false);
122  commonData.coordsAtGaussPts[gg].resize(3,false);
123  commonData.hoCoordsAtGaussPts[gg].resize(3,false);
124  for(int rr = 0;rr<3;rr++) {
125  (commonData.dIsplacements[gg])[rr] = cblas_ddot(
126  nb_dofs/3,&data.getFieldData()[rr],3,&data.getN()(ffgg,0),1
127  );
128  if(!fieldDisp) {
129  (commonData.dIsplacements[gg])[rr] -= getCoordsAtGaussPts()(ffgg,rr);
130  }
131  commonData.coordsAtGaussPts[gg][rr] = getCoordsAtGaussPts()(ffgg,rr);
132  commonData.hoCoordsAtGaussPts[gg][rr] = getCoordsAtGaussPts()(ffgg,rr);
133  }
134  } else {
135  for(int rr = 0;rr<3;rr++) {
136  commonData.dIsplacements[gg][rr] += cblas_ddot(
137  nb_dofs/3,&data.getFieldData()[rr],3,&data.getN()(ffgg,0),1
138  );
139  }
140  }
141  }
142  } catch (const std::exception& ex) {
143  ostringstream ss;
144  ss << "throw in method: " << ex.what() << endl;
145  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
146  }
147 
148  //cerr << "done\n";
149  PetscFunctionReturn(0);
150  }
151 
152  };
153 
155 
157  PeriodicFace(MoFEM::Interface &m_field,CommonData &common_data):
159  commonData(common_data) {}
160 
161  int getRule(int order) { return -1; }
162 
163  PetscErrorCode setGaussPts(int order) {
164  PetscFunctionBegin;
165 
166  try {
167 
168  EntityHandle face = numeredEntFiniteElementPtr->getEnt();
169  gaussPts.resize(3,commonData.localCoordsMap[face].size());
170  //cerr << commonData.localCoordsMap[face].size() << endl;
171  for(unsigned int ffgg = 0;ffgg!=commonData.localCoordsMap[face].size();ffgg++) {
172  //cerr << commonData.localCoordsMap[face][ffgg] << endl;
173  gaussPts(0,ffgg) = commonData.localCoordsMap[face][ffgg][0];
174  gaussPts(1,ffgg) = commonData.localCoordsMap[face][ffgg][1];
175  gaussPts(2,ffgg) = 0; // not used, since filed on this face is not integrated
176  }
177 
178  } catch (const std::exception& ex) {
179  ostringstream ss;
180  ss << "throw in method: " << ex.what() << endl;
181  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
182  }
183  PetscFunctionReturn(0);
184  }
185 
186  };
187 
189 
192 
195  CommonData &periodic_common_data
196  ):
197  VolumeElementForcesAndSourcesCore::UserDataOperator("DISPLACEMENT",OPROW),
198  commonData(common_data),
199  periodicCommonData(periodic_common_data) {
200  }
201 
202  PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
203  PetscFunctionBegin;
204  if(data.getIndices().size()==0) PetscFunctionReturn(0);
205  try {
206  EntityHandle this_tet = getNumeredEntFiniteElementPtr()->getEnt();
207  int nb_gauss_pts = data.getN().size1();
208  for(int ffgg = 0;ffgg<nb_gauss_pts;ffgg++) {
209  int gg = periodicCommonData.inTetTetGaussPtsNumber[this_tet][ffgg];
210  CommonData::MultiIndexData gauss_pt_data(gg,side,type);
211  pair<CommonData::Container::iterator,bool> p;
212  p = periodicCommonData.volumesContainer.insert(gauss_pt_data);
213  if(!p.second) {
214  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"data not inserted");
215  }
216  CommonData::MultiIndexData &p_data = const_cast<CommonData::MultiIndexData&>(*p.first);
217  MatrixDouble &diff_shape_fun = p_data.diffShapeFunctions;
218  diff_shape_fun = data.getDiffN(ffgg);
219  p_data.iNdices = data.getIndices();
220  if(type == MBVERTEX) {
223  }
224  }
225  } catch (const std::exception& ex) {
226  ostringstream ss;
227  ss << "throw in method: " << ex.what() << endl;
228  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
229  }
230  PetscFunctionReturn(0);
231  }
232 
233  };
234 
235  struct PeriodicVolume: VolumeElementForcesAndSourcesCore {
236 
239 
241  VolumeElementForcesAndSourcesCore(m_field),
242  forFace(0),
243  commonData(common_data) {}
244  int getRule(int order) { return -1; }
245  PetscErrorCode setGaussPts(int order) {
246 
247  PetscFunctionBegin;
248  try {
249  gaussPts.resize(4,0,false);
250  EntityHandle tet = numeredEntFiniteElementPtr->getEnt();
251  int vffgg = 0;
252  for(int ff = 0;ff<4;ff++) {
253  EntityHandle face;
254  rval = mField.get_moab().side_element(tet,2,ff,face); CHKERRQ_MOAB(rval);
255  if(face!=forFace) continue;
256  if(commonData.localCoordsMap.find(face)!=commonData.localCoordsMap.end()) {
257  const double coords_tet[12] = { 0,0,0, 1,0,0, 0,1,0, 0,0,1 };
258  int nb_face_gauss_pts = commonData.localCoordsMap[face].size();
259  gaussPts.resize(4,nb_face_gauss_pts);
260  commonData.inTetTetGaussPtsNumber[tet].resize(nb_face_gauss_pts);
261  for(int ffgg = 0;ffgg!=nb_face_gauss_pts;ffgg++,vffgg++) {
263  double N[3];
264  VectorDouble &local_coords = commonData.localCoordsMap[face][ffgg];
265  ierr = ShapeMBTRI(N,&local_coords[0],&local_coords[1],1); CHKERRQ(ierr);
266  for(int dd = 0;dd<3;dd++) {
267  gaussPts(dd,vffgg) =
268  N[0]*coords_tet[3*dataH1.facesNodes(ff,0)+dd]+
269  N[1]*coords_tet[3*dataH1.facesNodes(ff,1)+dd]+
270  N[2]*coords_tet[3*dataH1.facesNodes(ff,2)+dd];
271  }
272  gaussPts(3,vffgg) = 0;
273  }
274  } else {
275  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"data inconstancy");
276  }
277  }
278  //cerr << gaussPts << endl;
279  } catch (const std::exception& ex) {
280  ostringstream ss;
281  ss << "throw in method: " << ex.what() << endl;
282  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
283  }
284  PetscFunctionReturn(0);
285  }
286 
287  };
288 
290 
294 
295  bool iNit;
296  AdaptiveKDTree tRee;
297 
299  MoFEM::Interface &m_field,
300  NitscheMethod::BlockData &block_data,
301  NitscheMethod::CommonData &common_data,
302  CommonData &periodic_common_data
303  ):
304  NitscheMethod::MyVolumeFE(m_field,block_data,common_data),
305  periodicCommonData(periodic_common_data),
306  periodicFace(m_field,periodic_common_data),
307  periodicVolume(m_field,periodic_common_data),
308  iNit(false),
309  tRee(&m_field.get_moab()) {
310  periodicFace.getOpPtrVector().push_back(
311  new OpGetFaceData(periodic_common_data)
312  );
313  }
314 
316  PetscErrorCode iNitialise() {
317  PetscFunctionBegin;
318 
319  treeRootSet = 0;
321  //cerr << blockData.fAces << endl;
322  PetscFunctionReturn(0);
323  }
324 
326  PetscFunctionBegin;
327 
328 
329 
330 
331  if(!iNit) {
332  ierr = iNitialise(); CHKERRQ(ierr);
333  iNit = true;
334  }
335 
336  //cerr << "part1\n";
337  try {
338  const double tol = 1e-18;
347  periodicCommonData.sTress.clear();
350  int gg = 0;
351  for(int ff = 0;ff<4;ff++) {
352  if(commonData.facesFePtr[ff]==NULL) continue;
353  for(unsigned int fgg = 0;fgg<commonData.faceGaussPts[ff].size2();fgg++,gg++) {
354  VectorDouble& ray = commonData.rAy[ff];
356  3,ublas::shallow_array_adaptor<double>(
357  3,&commonData.hoCoordsAtGaussPts[ff](fgg,0)
358  )
359  );
360  vector<EntityHandle> triangles_out;
361  vector<double> distance_out;
362  rval = tRee.ray_intersect_triangles(
363  treeRootSet,tol,&ray[0],&coords[0],triangles_out,distance_out
364  ); CHKERRQ_MOAB(rval);
365  double diffNTRI[6];
366  ierr = ShapeDiffMBTRI(diffNTRI); CHKERRQ(ierr);
367  vector<EntityHandle> triangles_out_on_other_side;
368  vector<double> distance_out_on_other_side;
369  {
370  vector<EntityHandle>::iterator hit = triangles_out.begin();
371  vector<double>::iterator dit = distance_out.begin();
372  for(;dit!=distance_out.end();dit++,hit++) {
373  const EntityHandle* conn;
374  int number_nodes;
375  rval = mField.get_moab().get_connectivity(*hit,conn,number_nodes,true); CHKERRQ_MOAB(rval);
376  double coords[9];
377  rval = mField.get_moab().get_coords(conn,number_nodes,coords); CHKERRQ_MOAB(rval);
378  double normal[3];
379  ierr = ShapeFaceNormalMBTRI(diffNTRI,coords,normal); CHKERRQ(ierr);
380  *dit *= cblas_ddot(3,normal,1,&ray[0],1)/cblas_dnrm2(3,normal,1);
381  if(*dit>tol) {
382  triangles_out_on_other_side.push_back(*hit);
383  distance_out_on_other_side.push_back(*dit);
384  }
385  }
386  }
387  vector<double>::iterator max_dit = max_element(
388  distance_out_on_other_side.begin(),distance_out_on_other_side.end()
389  );
390  EntityHandle other_face = triangles_out_on_other_side[
391  std::distance(distance_out_on_other_side.begin(),max_dit)
392  ];
393  VectorDouble other_coords = coords+ray*(*max_dit);
394  periodicCommonData.dIstance[other_face] = *max_dit;
395  VectorDouble local_coords;
396  local_coords.resize(2);
397  {
398  const EntityHandle* conn;
399  int number_nodes;
400  rval = mField.get_moab().get_connectivity(other_face,conn,number_nodes,false); CHKERRQ_MOAB(rval);
401 
402  {
403  double face_coords[18];
404  rval = mField.get_moab().get_coords(conn,number_nodes,face_coords); CHKERRQ_MOAB(rval);
405 
406  VectorDouble res,d_local_coords;
407  res.resize(3,false);
408  d_local_coords.resize(2,false);
409 
410  double N[number_nodes];
411  double diffNTRI[number_nodes*2];
412  double nrm2,nrm2_res;
413  const double eps = 1e-14;
414  const int max_it = 10;
415 
416  int it = 0;
417  local_coords.clear();
418  do {
419  if(number_nodes == 3) {
420  ierr = ShapeMBTRI(N,&local_coords[0],&local_coords[1],1); CHKERRQ(ierr);
421  } else if(number_nodes == 6) {
422  ierr = ShapeMBTRIQ(N,&local_coords[0],&local_coords[1],1); CHKERRQ(ierr);
423  } else {
424  SETERRQ(PETSC_COMM_SELF,MOFEM_NOT_IMPLEMENTED,"not implemented");
425  }
426  res.clear();
427  for(int dd = 0;dd<3;dd++) {
428  res[dd] += cblas_ddot(number_nodes,N,1,&face_coords[dd],3)-other_coords[dd];
429  }
430  nrm2_res = cblas_dnrm2(3,&res[0],1);
431  if(number_nodes == 3) {
432  ierr = ShapeDiffMBTRI(diffNTRI); CHKERRQ(ierr);
433  } else {
434  ierr = ShapeDiffMBTRIQ(diffNTRI,&local_coords[0],&local_coords[1],1); CHKERRQ(ierr);
435  }
436  MatrixDouble A;
437  A.resize(3,2);
438  for(int dd = 0;dd<3;dd++) {
439  A(dd,0) = cblas_ddot(number_nodes,&diffNTRI[0],2,&face_coords[dd],3);
440  A(dd,1) = cblas_ddot(number_nodes,&diffNTRI[1],2,&face_coords[dd],3);
441  }
442  MatrixDouble ATA;
443  ATA = prod(trans(A),A);
444  VectorDouble ATres;
445  ATres = prod(trans(A),res);
446  double det = (ATA(0,0)*ATA(1,1)-ATA(0,1)*ATA(1,0));
447  d_local_coords[0] = -(ATA(1,1)*ATres[0]-ATA(0,1)*ATres[1])/det;
448  d_local_coords[1] = -(ATA(0,0)*ATres[1]-ATA(1,0)*ATres[0])/det;
449  local_coords += d_local_coords;
450  nrm2 = cblas_dnrm2(2,&d_local_coords[0],1);
451  if((it++)>max_it) {
452  cerr << "Warrning: max_it reached in PeriodicNitscheConstrains:: ... ::doAdditionalJobWhenGuassPtsAreCalulated\n";
453  cerr << "nrm2 " << nrm2 << " nrm2_res " << nrm2_res << endl;
454  cerr << "res " << res << endl;
455  cerr << "other_coords " << other_coords << endl;
456  break;
457  }
458  } while(nrm2>eps || nrm2_res>eps);
459  if(local_coords[0]<-eps) {
460  SETERRQ1(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"local_coords[0]=%6.4e<eps",local_coords[0]);
461  }
462  if(local_coords[1]<-eps) {
463  SETERRQ1(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"local_coords[1]=%6.4e<eps",local_coords[1]);
464  }
465  if(local_coords[0]>1+eps) {
466  SETERRQ1(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"local_coords[0]=%6.4e>1+eps",local_coords[0]);
467  }
468  if(local_coords[1]>1+eps) {
469  SETERRQ1(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"local_coords[1]=%6.4e>1+eps",local_coords[1]);
470  }
471  if(local_coords[0]+local_coords[1]>1+eps) {
472  SETERRQ1(
473  PETSC_COMM_SELF,
475  "local_coords[0]+local_coords[1]>1",
476  local_coords[0]+local_coords[1]
477  );
478  }
479  }
480 
481  }
482  periodicCommonData.localCoordsMap[other_face].push_back(local_coords);
483  periodicCommonData.inTetFaceGaussPtsNumber[other_face].push_back(gg);
484  }
485  }
488  periodicCommonData.sTress.resize(gg);
491 
492  } catch (const std::exception& ex) {
493  ostringstream ss;
494  ss << "throw in method: " << ex.what() << endl;
495  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
496  }
497 
498  try {
499  for(
500  map<EntityHandle,vector<int> >::iterator mit = periodicCommonData.inTetFaceGaussPtsNumber.begin();
502  mit++
503  ) {
504  //cerr << mit->first << " " << mit->second.size() << endl;
505  NumeredEntFiniteElement_multiIndex::index<Composite_Name_And_Ent_mi_tag>::type::iterator it;
507  problemPtr->numeredFiniteElements)
508  .get<Composite_Name_And_Ent_mi_tag>()
509  .find(
510  boost::make_tuple(blockData.faceElemName, mit->first));
512  problemPtr->numeredFiniteElements)
513  .get<Composite_Name_And_Ent_mi_tag>()
514  .end()) {
515  SETERRQ1(
516  PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"No finite element found < %s >",
517  blockData.faceElemName.c_str()
518  );
519  }
520  //cerr << *it << endl;
521  boost::shared_ptr<const NumeredEntFiniteElement> face_ptr = *it;
522  periodicFace.copyBasicMethod(*this);
524  periodicFace.nInTheLoop = 0;
525  periodicFace.numeredEntFiniteElementPtr = face_ptr;
526  periodicFace.dataPtr = face_ptr->sPtr->data_dofs;
527  periodicFace.rowPtr = face_ptr->rows_dofs;
528  periodicFace.colPtr = face_ptr->cols_dofs;
529  ierr = periodicFace(); CHKERRQ(ierr);
530 
531  Range tets;
532  ierr = mField.get_adjacencies_equality(mit->first,3,tets); CHKERRQ(ierr);
533  if(tets.size()!=1) {
534  SETERRQ1(
535  PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"Expected one tet, but is %d",
536  tets.size()
537  );
538  }
540  problemPtr->numeredFiniteElements)
541  .get<Composite_Name_And_Ent_mi_tag>()
542  .find(boost::make_tuple(periodicCommonData.volumeElemName,
543  tets[0]));
545  problemPtr->numeredFiniteElements)
546  .get<Composite_Name_And_Ent_mi_tag>()
547  .end()) {
548  SETERRQ1(
549  PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"No finite element found < %s >",
551  );
552  }
553  boost::shared_ptr<const NumeredEntFiniteElement> tet_ptr = *it;
554  periodicVolume.copyBasicMethod(*this);
556  periodicVolume.nInTheLoop = 0;
557  periodicVolume.numeredEntFiniteElementPtr = tet_ptr;
558  periodicVolume.dataPtr = face_ptr->sPtr->data_dofs;
559  periodicVolume.rowPtr = face_ptr->rows_dofs;
560  periodicVolume.colPtr = face_ptr->cols_dofs;
561  periodicVolume.forFace = mit->first;
562  ierr = periodicVolume(); CHKERRQ(ierr);
563 
564  }
565  } catch (const std::exception& ex) {
566  ostringstream ss;
567  ss << "throw in method: " << ex.what() << endl;
568  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
569  }
570 
571  PetscFunctionReturn(0);
572  }
573 
574  };
575 
577 
579  vector<Vec> &F;
580 
582  string field_name,
583  NitscheMethod::BlockData &nitsche_block_data,
584  NitscheMethod::CommonData &nitsche_common_data,
587  CommonData &periodic_common_data,
588  vector<Vec> &f,
589  bool field_disp = true
590  ):
591  OpCommon(
592  field_name,
593  nitsche_block_data,
594  nitsche_common_data,
595  data,
596  common_data,
597  field_disp,
599  ),
600  periodicCommonData(periodic_common_data),
601  F(f) {
602  }
603 
605 
606 
610 
612  static PetscErrorCode calcualteDMatrix(VectorAdaptor coords,MatrixDouble &mat_d) {
613  PetscFunctionBegin;
614  mat_d.resize(3,6,false);
615  mat_d.clear();
616  mat_d(0,0) = coords[0]; mat_d(0,3) = coords[1]*0.5; mat_d(0,5) = coords[2]*0.5;
617  mat_d(1,1) = coords[1]; mat_d(1,3) = coords[0]*0.5; mat_d(1,4) = coords[2]*0.5;
618  mat_d(2,2) = coords[2]; mat_d(2,4) = coords[1]*0.5; mat_d(2,5) = coords[0]*0.5;
619  PetscFunctionReturn(0);
620  }
621 
622  PetscErrorCode doWork(
623  int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
624  ) {
625  PetscFunctionBegin;
626 
627  if(dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) == dAta.tEts.end()) {
628  PetscFunctionReturn(0);
629  }
630 
631  if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
632  int nb_dofs = row_data.getIndices().size();
633  const double gamma = nitscheBlockData.gamma;
634  const double phi = nitscheBlockData.phi;
635  double eps = periodicCommonData.eps;
636 
637  try {
638 
639  int gg = 0;
640  for(int ff = 0;ff<4;ff++) {
641  if(nitscheCommonData.facesFePtr[ff]==NULL) continue;
642  int nb_face_gauss_pts = nitscheCommonData.faceGaussPts[ff].size2();
643  ierr = getFaceRadius(ff); CHKERRQ(ierr);
644  for(int ii = 0;ii<6;ii++) {
645  nF0[ii].resize(nb_dofs,false);
646  nF0[ii].clear();
647  nF1[ii].resize(nb_dofs,false);
648  nF1[ii].clear();
649  }
650  for(int fgg = 0;fgg<nb_face_gauss_pts;fgg++,gg++) {
651  ierr = getGammaH(gamma,gg); CHKERRQ(ierr);
652  double val = getGaussPts()(3,gg);
653  ierr = getJac(row_data,gg,jAc_row); CHKERRQ(ierr);
654  ierr = getTractionVariance(gg,fgg,ff,jAc_row,tRac_v); CHKERRQ(ierr);
655  VectorAdaptor normal = VectorAdaptor(
656  3,ublas::shallow_array_adaptor<double>(
657  3,&nitscheCommonData.faceNormals[ff](fgg,0)
658  )
659  );
660  double area = cblas_dnrm2(3,&normal[0],1);
661 
662  voightStrain.resize(6,false);
663  dispOnThisSide.resize(3,false);
664  dispOnOtherSide.resize(3,false);
666  3,ublas::shallow_array_adaptor<double>(3,&getCoordsAtGaussPts()(gg,0))
667  );
669  VectorDouble other_x = VectorAdaptor(
670  3,ublas::shallow_array_adaptor<double>(3,&periodicCommonData.hoCoordsAtGaussPts[gg][0])
671  );
672  ierr = calcualteDMatrix(other_x,matD1);
673  const double err = 1e-10;
674  double delta = inner_prod(other_x,normal)+inner_prod(x,normal);
675  if(fabs(delta)>err) {
676  cerr << x << endl;
677  cerr << other_x << endl;
678  SETERRQ1(
679  PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"err = %6.4e",delta
680  );
681  }
682 
683  for(int ii = 0;ii<6;ii++) {
684  voightStrain.clear();
685  voightStrain[ii] = 1;
686  noalias(dispOnThisSide) = prod(matD0,voightStrain);
687  noalias(dispOnOtherSide) = prod(matD1,voightStrain);
688  for(int dd1 = 0;dd1<nb_dofs/3;dd1++) {
689  double n_val = row_data.getN(gg)[dd1];
690  for(int dd2 = 0;dd2<3;dd2++) {
691  nF0[ii][3*dd1+dd2] -= dispOnThisSide[dd2]*n_val*val*area;
692  nF0[ii][3*dd1+dd2] += (1-eps)*dispOnOtherSide[dd2]*n_val*val*area;
693  }
694  }
695  for(int dd = 0;dd<nb_dofs;dd++) {
696  double dot;
697  dot = cblas_ddot(
698  3,&dispOnThisSide[0],1,&tRac_v(0,dd),tRac_v.size2()
699  );
700  dot -= (1-eps)*cblas_ddot(
701  3,&dispOnOtherSide[0],1,&tRac_v(0,dd),tRac_v.size2()
702  );
703  nF1[ii][dd] += val*phi*dot;
704  }
705  }
706  }
707 
708  for(int ii = 0;ii<6;ii++) {
709  nF0[ii] *= -(1/gammaH);
710  ierr = VecSetValues(
711  F[ii],
712  nb_dofs,
713  &row_data.getIndices()[0],
714  &nF0[ii][0],
715  ADD_VALUES
716  ); CHKERRQ(ierr);
717  nF1[ii] *= -1;
718  ierr = VecSetValues(
719  F[ii],
720  nb_dofs,
721  &row_data.getIndices()[0],
722  &nF1[ii][0],
723  ADD_VALUES
724  ); CHKERRQ(ierr);
725  }
726 
727  }
728 
729  } catch (const std::exception& ex) {
730  ostringstream ss;
731  ss << "throw in method: " << ex.what() << endl;
732  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
733  }
734 
735  PetscFunctionReturn(0);
736  }
737 
738  };
739 
741 
743 
745  const string field_name,
746  NitscheMethod::BlockData &nitsche_block_data,
747  NitscheMethod::CommonData &nitsche_common_data,
750  CommonData &periodic_common_data,
751  bool field_disp = true
752  ):
753  OpCommon(
754  field_name,
755  nitsche_block_data,
756  nitsche_common_data,
757  data,
758  common_data,
759  field_disp,
761  ),
762  periodicCommonData(periodic_common_data) {
763  }
764 
767 
768  PetscErrorCode doWork(
769  int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
770  ) {
771  PetscFunctionBegin;
772 
773 
774  if(dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) == dAta.tEts.end()) {
775  PetscFunctionReturn(0);
776  }
777  if(row_data.getIndices().size()==0) PetscFunctionReturn(0);
778  int nb_dofs_row = row_data.getIndices().size();
779  const double gamma = nitscheBlockData.gamma;
780  const double phi = nitscheBlockData.phi;
781  double eps = periodicCommonData.eps;
782 
783  try {
784 
785  int gg = 0;
786  for(int ff = 0;ff<4;ff++) {
787  if(nitscheCommonData.facesFePtr[ff]==NULL) continue;
788  int nb_face_gauss_pts = nitscheCommonData.faceGaussPts[ff].size2();
789  vector<vector<MatrixDouble > > &P = nitscheCommonData.P;
790  P.resize(4);
791  P[ff].resize(nb_face_gauss_pts);
792  for(int fgg = 0;fgg<nb_face_gauss_pts;fgg++) {
793  P[ff][fgg].resize(3,3);
794  P[ff][fgg].clear();
795  for(int dd = 0;dd<3;dd++) {
796  P[ff][fgg](dd,dd) = 1;
797  }
798  }
799  ierr = getFaceRadius(ff); CHKERRQ(ierr);
800  for(int fgg = 0;fgg<nb_face_gauss_pts;fgg++,gg++) {
801  ierr = getGammaH(gamma,gg); CHKERRQ(ierr);
802  double val = getGaussPts()(3,gg);
803  ierr = getJac(row_data,gg,jAc_row); CHKERRQ(ierr);
804  ierr = getTractionVariance(gg,fgg,ff,jAc_row,tRac_v); CHKERRQ(ierr);
805  VectorAdaptor normal = VectorAdaptor(
806  3,ublas::shallow_array_adaptor<double>(3,&nitscheCommonData.faceNormals[ff](fgg,0))
807  );
808  double area = cblas_dnrm2(3,&normal[0],1);
809  const double alpha = 0.5; // mortar parameter
810 
811  CommonData::Container::nth_index<0>::type::iterator it,hi_it;
812  it = periodicCommonData.facesContainer.get<0>().lower_bound(gg);
813  hi_it = periodicCommonData.facesContainer.get<0>().upper_bound(gg);
814  for(;it!=hi_it;it++) {
815 
816  const ublas::vector<int> &indices = it->iNdices;
817  const VectorDouble &shape = it->shapeFunctions;
818 
819  int nb_dofs_col = indices.size();
820  if(indices.size()==0) {
821  continue;
822  }
823  kMatrix0.resize(nb_dofs_row,nb_dofs_col,false);
824  kMatrix0.clear();
825  for(int dd1 = 0;dd1<nb_dofs_row/3;dd1++) {
826  double n_row = row_data.getN()(gg,dd1);
827  for(int dd2 = 0;dd2<nb_dofs_col/3;dd2++) {
828  double n_col = shape[dd2];
829  for(int dd3 = 0;dd3<3;dd3++) {
830  kMatrix0(3*dd1+dd3,3*dd2+dd3) -= (1-eps)*val*n_row*n_col*area;
831  }
832  }
833  }
834  kMatrix1.resize(nb_dofs_row,nb_dofs_col,false);
835  kMatrix1.clear();
836  for(int dd1 = 0;dd1<nb_dofs_row;dd1++) {
837  for(int dd2 = 0;dd2<nb_dofs_col/3;dd2++) {
838  double n_col = shape[dd2];
839  for(int dd3 = 0;dd3<3;dd3++) {
840  kMatrix1(dd1,3*dd2+dd3) += (1-eps)*phi*val*tRac_v(dd3,dd1)*n_col;
841  }
842  }
843  }
844 
845  kMatrix0 /= gammaH;
846  ierr = MatSetValues(
847  getFEMethod()->snes_B,
848  nb_dofs_row,
849  &row_data.getIndices()[0],
850  nb_dofs_col,
851  &indices[0],
852  &kMatrix0(0,0),
853  ADD_VALUES
854  ); CHKERRQ(ierr);
855 
856  ierr = MatSetValues(
857  getFEMethod()->snes_B,
858  nb_dofs_row,
859  &row_data.getIndices()[0],
860  nb_dofs_col,
861  &indices[0],
862  &kMatrix1(0,0),
863  ADD_VALUES
864  ); CHKERRQ(ierr);
865 
866  kMatrix1.resize(nb_dofs_col,nb_dofs_row);
867  kMatrix1.clear();
868  for(int dd1 = 0;dd1<nb_dofs_col/3;dd1++) {
869  double n_col = shape[dd1];
870  for(int dd2 = 0;dd2<3;dd2++) {
871  for(int dd3 = 0;dd3<nb_dofs_row;dd3++) {
872  kMatrix1(3*dd1+dd2,dd3) += (1-eps)*val*n_col*tRac_v(dd2,dd3);
873  }
874  }
875  }
876 
877  kMatrix1 *= alpha;
878  ierr = MatSetValues(
879  getFEMethod()->snes_B,
880  nb_dofs_col,
881  &indices[0],
882  nb_dofs_row,
883  &row_data.getIndices()[0],
884  &kMatrix1(0,0),
885  ADD_VALUES
886  ); CHKERRQ(ierr);
887 
888  }
889 
890  MatrixDouble jac;
891  MatrixDouble trac;
892 
893  it = periodicCommonData.volumesContainer.get<0>().lower_bound(gg);
894  hi_it = periodicCommonData.volumesContainer.get<0>().upper_bound(gg);
895  for(;it!=hi_it;it++) {
896 
897  int nb = it->iNdices.size();
898  jac.resize(9,nb,false);
899  jac.clear();
900  const MatrixDouble diff_n = it->diffShapeFunctions;
901  const MatrixDouble jac_stress = periodicCommonData.stressJacobian[gg];
902  for(int dd = 0;dd<nb/3;dd++) {
903  for(int rr = 0;rr<3;rr++) {
904  for(int ii = 0;ii<9;ii++) {
905  for(int jj = 0;jj<3;jj++) {
906  jac(ii,3*dd+rr) += jac_stress(ii,3*rr+jj)*diff_n(dd,jj);
907  }
908  }
909  }
910  }
911  trac.resize(3,jac.size2());
912  trac.clear();
913  for(unsigned int dd2 = 0;dd2<jac.size2();dd2++) {
914  for(int nn = 0;nn<3;nn++) {
915  trac(nn,dd2) = -cblas_ddot(3,&jac(3*nn,dd2),jac.size2(),&normal[0],1);
916  }
917  }
918 
919  kMatrixOff.resize(nb_dofs_row,nb);
920  kMatrixOff.clear();
921  for(int dd1 = 0;dd1<nb_dofs_row/3;dd1++) {
922  double n_row = row_data.getN()(gg,dd1);
923  for(int dd2 = 0;dd2<3;dd2++) {
924  for(int dd3 = 0;dd3<nb;dd3++) {
925  kMatrixOff(3*dd1+dd2,dd3) += (1-eps)*val*n_row*trac(dd2,dd3);
926  }
927  }
928  }
929 
930  kMatrixOff *= (1-alpha);
931  ierr = MatSetValues(
932  getFEMethod()->snes_B,
933  nb_dofs_row,
934  &row_data.getIndices()[0],
935  nb,
936  &it->iNdices[0],
937  &kMatrixOff(0,0),
938  ADD_VALUES
939  ); CHKERRQ(ierr);
940 
941  }
942 
943  }
944  }
945 
946  } catch (const std::exception& ex) {
947  ostringstream ss;
948  ss << "throw in method: " << ex.what() << endl;
949  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
950  }
951 
952  PetscFunctionReturn(0);
953  }
954 
955  };
956 
958 
960 
962  const string field_name,
965  Vec stress_homo_ghost_vector
966  ):
968  field_name,data,common_data
969  ),
970  stressHomoGhostVector(stress_homo_ghost_vector) {
971  }
972 
974 
975  PetscErrorCode doWork(
976  int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
977  ) {
978  PetscFunctionBegin;
979 
980 
981 
982  if(row_type != MBVERTEX) PetscFunctionReturn(0);
983  if(dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) == dAta.tEts.end()) {
984  PetscFunctionReturn(0);
985  }
986 
987  try {
988 
989  homoStress.resize(6,false);
990  homoStress.clear();
991 
992  for(unsigned int gg = 0;gg<row_data.getN().size1();gg++) {
993  const MatrixDouble& stress = commonData.sTress[gg];
994  double val = getGaussPts()(3,gg)*getVolume();
995  if(getHoGaussPtsDetJac().size()>0) {
996  val *= getHoGaussPtsDetJac()[gg]; ///< higher order geometry
997  }
998  homoStress[0] += stress(0,0)*val;
999  homoStress[1] += stress(1,1)*val;
1000  homoStress[2] += stress(2,2)*val;
1001  homoStress[3] += (stress(0,1)+stress(1,0))*0.5*val;
1002  homoStress[4] += (stress(1,2)+stress(2,1))*0.5*val;
1003  homoStress[5] += (stress(0,2)+stress(2,0))*0.5*val;
1004  }
1005 
1006  const int indices_6[6] = {0, 1, 2, 3, 4, 5};
1007  ierr = VecSetValues(
1008  stressHomoGhostVector,6,indices_6,&homoStress[0],ADD_VALUES
1009  ); CHKERRQ(ierr);
1010 
1011  } catch (const std::exception& ex) {
1012  ostringstream ss;
1013  ss << "throw in method: " << ex.what() << endl;
1014  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1015  }
1016 
1017  PetscFunctionReturn(0);
1018  }
1019 
1020  };
1021 
1023 
1025 
1027  const string field_name,
1028  NitscheMethod::BlockData &nitsche_block_data,
1029  NitscheMethod::CommonData &nitsche_common_data,
1032  Vec stress_homo_ghost_vector,
1033  bool field_disp = 0
1034  ):
1036  field_name,
1037  nitsche_block_data,
1038  nitsche_common_data,
1039  data,
1040  common_data,
1041  field_disp,
1043  ),
1044  stressHomoGhostVector(stress_homo_ghost_vector) {
1045  }
1046 
1050 
1053 
1054  PetscErrorCode doWork(
1055  int row_side,EntityType row_type,DataForcesAndSourcesCore::EntData &row_data
1056  ) {
1057  PetscFunctionBegin;
1058 
1059 
1060 
1061  if(row_type != MBVERTEX) PetscFunctionReturn(0);
1062  if(dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) == dAta.tEts.end()) {
1063  PetscFunctionReturn(0);
1064  }
1065 
1066  try {
1067 
1068  homoStress.resize(6,false);
1069  homoStress.clear();
1070 
1071  int gg = 0;
1072  for(int ff = 0;ff<4;ff++) {
1073  if(nitscheCommonData.facesFePtr[ff]==NULL) continue;
1074  int nb_face_gauss_pts = nitscheCommonData.faceGaussPts[ff].size2();
1075  for(int fgg = 0;fgg<nb_face_gauss_pts;fgg++,gg++) {
1076  double val = getGaussPts()(3,gg);
1077  const MatrixDouble& stress = commonData.sTress[gg];
1078  VectorAdaptor normal = VectorAdaptor(
1079  3,ublas::shallow_array_adaptor<double>(
1080  3,&nitscheCommonData.faceNormals[ff](fgg,0)
1081  )
1082  );
1083 
1085  3,ublas::shallow_array_adaptor<double>(
1086  3,&getCoordsAtGaussPts()(gg,0)
1087  )
1088  );
1089 
1090  tRaction.resize(3,false);
1091  noalias(tRaction) = prod(normal,stress);
1092 
1093  homoStress[0] += x[0]*tRaction[0]*val;
1094  homoStress[1] += x[1]*tRaction[1]*val;
1095  homoStress[2] += x[2]*tRaction[2]*val;
1096  homoStress[3] += (x[0]*tRaction[1]+x[1]*tRaction[0])*0.5*val;
1097  homoStress[4] += (x[1]*tRaction[2]+x[2]*tRaction[1])*0.5*val;
1098  homoStress[5] += (x[0]*tRaction[2]+x[2]*tRaction[0])*0.5*val;
1099  }
1100  }
1101 
1102  const int indices_6[6] = {0, 1, 2, 3, 4, 5};
1103  ierr = VecSetValues(
1104  stressHomoGhostVector,6,indices_6,&homoStress[0],ADD_VALUES
1105  ); CHKERRQ(ierr);
1106 
1107  if(gg != (int)getGaussPts().size2()) {
1108  SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"wrong number of gauss pts");
1109  }
1110 
1111 
1112  } catch (const std::exception& ex) {
1113  ostringstream ss;
1114  ss << "throw in method: " << ex.what() << endl;
1115  SETERRQ(PETSC_COMM_SELF,1,ss.str().c_str());
1116  }
1117 
1118  PetscFunctionReturn(0);
1119  }
1120 
1121  };
1122 
1123 
1124 };
1125 
1126 #endif //__NITSCHE_PERIODIC_METHOD_HPP__
PeriodicNitscheConstrains
Periodic Nitsche method.
Definition: NitschePeriodicMethod.hpp:27
NitscheMethod::OpCommon::jAc_row
MatrixDouble jAc_row
Definition: NitscheMethod.hpp:449
PeriodicNitscheConstrains::CommonData::dIsplacements
vector< VectorDouble > dIsplacements
Definition: NitschePeriodicMethod.hpp:39
PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral
Definition: NitschePeriodicMethod.hpp:1022
NitscheMethod::OpBasicCommon::getFaceRadius
virtual MoFEMErrorCode getFaceRadius(int ff)
Definition: NitscheMethod.hpp:412
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
NitscheMethod::OpCommon::dAta
NonlinearElasticElement::BlockData & dAta
Definition: NitscheMethod.hpp:435
EntityHandle
PeriodicNitscheConstrains::OpRhsPeriodicNormal::dispOnOtherSide
VectorDouble dispOnOtherSide
Definition: NitschePeriodicMethod.hpp:609
NonlinearElasticElement::OpRhsPiolaKirchhoff::dAta
BlockData & dAta
Definition: NonLinearElasticElement.hpp:523
PeriodicNitscheConstrains::OpRhsPeriodicNormal::matD0
MatrixDouble matD0
Definition: NitschePeriodicMethod.hpp:611
NonlinearElasticElement::CommonData::jacStress
std::vector< MatrixDouble > jacStress
this is simply material tangent operator
Definition: NonLinearElasticElement.hpp:113
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
PeriodicNitscheConstrains::CommonData::eps
double eps
Definition: NitschePeriodicMethod.hpp:33
PeriodicNitscheConstrains::PeriodicVolume::PeriodicVolume
PeriodicVolume(MoFEM::Interface &m_field, CommonData &common_data)
Definition: NitschePeriodicMethod.hpp:240
PeriodicNitscheConstrains::CommonData::MultiIndexData::diffShapeFunctions
MatrixDouble diffShapeFunctions
Definition: NitschePeriodicMethod.hpp:52
NitscheMethod::CommonData
Common data shared between finite element operators.
Definition: NitscheMethod.hpp:91
NitscheMethod::BlockData
Block data for Nitsche method.
Definition: NitscheMethod.hpp:82
PeriodicNitscheConstrains::CommonData::inTetTetGaussPtsNumber
map< EntityHandle, vector< int > > inTetTetGaussPtsNumber
Definition: NitschePeriodicMethod.hpp:37
NitscheMethod::OpCommon::tRac_v
MatrixDouble tRac_v
Definition: NitscheMethod.hpp:451
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral::matD
MatrixDouble matD
Definition: NitschePeriodicMethod.hpp:1047
PeriodicNitscheConstrains::CommonData::MultiIndexData::iNdices
ublas::vector< int > iNdices
Definition: NitschePeriodicMethod.hpp:49
NitscheMethod::CommonData::P
std::vector< std::vector< MatrixDouble > > P
projection matrix
Definition: NitscheMethod.hpp:110
NonlinearElasticElement::OpRhsPiolaKirchhoff::commonData
CommonData & commonData
Definition: NonLinearElasticElement.hpp:524
PeriodicNitscheConstrains::OpRhsPeriodicNormal::dispOnThisSide
VectorDouble dispOnThisSide
Definition: NitschePeriodicMethod.hpp:608
NitscheMethod::MyVolumeFE::commonData
CommonData & commonData
Definition: NitscheMethod.hpp:243
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:30
PeriodicNitscheConstrains::OpGetVolumeData::periodicCommonData
CommonData & periodicCommonData
Definition: NitschePeriodicMethod.hpp:191
PeriodicNitscheConstrains::PeriodicFace::getRule
int getRule(int order)
Definition: NitschePeriodicMethod.hpp:161
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
PeriodicNitscheConstrains::PeriodicFace
Definition: NitschePeriodicMethod.hpp:154
PeriodicNitscheConstrains::CommonData::MultiIndexData
Definition: NitschePeriodicMethod.hpp:45
PeriodicNitscheConstrains::OpGetFaceData::commonData
CommonData & commonData
Definition: NitschePeriodicMethod.hpp:85
NitscheMethod::MyVolumeFE::blockData
BlockData & blockData
Definition: NitscheMethod.hpp:242
PeriodicNitscheConstrains::CommonData::facesContainer
Container facesContainer
Definition: NitschePeriodicMethod.hpp:77
PeriodicNitscheConstrains::OpLhsPeriodicNormal
Definition: NitschePeriodicMethod.hpp:740
PeriodicNitscheConstrains::PeriodicFace::PeriodicFace
PeriodicFace(MoFEM::Interface &m_field, CommonData &common_data)
Definition: NitschePeriodicMethod.hpp:157
PeriodicNitscheConstrains::OpGetVolumeData::doWork
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: NitschePeriodicMethod.hpp:202
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getVolume
double getVolume() const
element volume (linear geometry)
Definition: VolumeElementForcesAndSourcesCore.hpp:161
PeriodicNitscheConstrains::CommonData::coordsAtGaussPts
vector< VectorDouble > coordsAtGaussPts
Definition: NitschePeriodicMethod.hpp:40
NonlinearElasticElement::OpRhsPiolaKirchhoff::OpRhsPiolaKirchhoff
OpRhsPiolaKirchhoff(const std::string field_name, BlockData &data, CommonData &common_data)
Definition: NonLinearElasticElement.cpp:595
PeriodicNitscheConstrains::OpLhsPeriodicNormal::kMatrix1
MatrixDouble kMatrix1
Definition: NitschePeriodicMethod.hpp:765
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
NitscheMethod::BlockData::faceElemName
string faceElemName
name of element face
Definition: NitscheMethod.hpp:85
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::VolumeElementForcesAndSourcesCore::coords
VectorDouble coords
Definition: VolumeElementForcesAndSourcesCore.hpp:89
PeriodicNitscheConstrains::OpCalculateHomogenisedStressVolumeIntegral::stressHomoGhostVector
Vec stressHomoGhostVector
Definition: NitschePeriodicMethod.hpp:959
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1236
PeriodicNitscheConstrains::MyNitscheVolume::periodicCommonData
CommonData & periodicCommonData
Definition: NitschePeriodicMethod.hpp:291
PeriodicNitscheConstrains::PeriodicVolume::forFace
EntityHandle forFace
Definition: NitschePeriodicMethod.hpp:237
PeriodicNitscheConstrains::OpRhsPeriodicNormal::periodicCommonData
CommonData & periodicCommonData
Definition: NitschePeriodicMethod.hpp:578
PeriodicNitscheConstrains::PeriodicVolume
Definition: NitschePeriodicMethod.hpp:235
ShapeDiffMBTRI
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
PeriodicNitscheConstrains::OpGetVolumeData::OpGetVolumeData
OpGetVolumeData(NonlinearElasticElement::CommonData &common_data, CommonData &periodic_common_data)
Definition: NitschePeriodicMethod.hpp:193
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
PeriodicNitscheConstrains::OpGetFaceData
Definition: NitschePeriodicMethod.hpp:83
PeriodicNitscheConstrains::OpGetFaceData::OpGetFaceData
OpGetFaceData(CommonData &common_data, bool field_disp=true)
Definition: NitschePeriodicMethod.hpp:88
PeriodicNitscheConstrains::OpCalculateHomogenisedStressVolumeIntegral::OpCalculateHomogenisedStressVolumeIntegral
OpCalculateHomogenisedStressVolumeIntegral(const string field_name, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, Vec stress_homo_ghost_vector)
Definition: NitschePeriodicMethod.hpp:961
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
PeriodicNitscheConstrains::OpRhsPeriodicNormal::voightStrain
VectorDouble voightStrain
Definition: NitschePeriodicMethod.hpp:607
PeriodicNitscheConstrains::OpRhsPeriodicNormal::OpRhsPeriodicNormal
OpRhsPeriodicNormal(string field_name, NitscheMethod::BlockData &nitsche_block_data, NitscheMethod::CommonData &nitsche_common_data, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, CommonData &periodic_common_data, vector< Vec > &f, bool field_disp=true)
Definition: NitschePeriodicMethod.hpp:581
PeriodicNitscheConstrains::OpLhsPeriodicNormal::doWork
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: NitschePeriodicMethod.hpp:768
delta
static constexpr double delta
Definition: prism_elements_from_surface.cpp:18
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
PeriodicNitscheConstrains::OpRhsPeriodicNormal::matD1
MatrixDouble matD1
Definition: NitschePeriodicMethod.hpp:611
PeriodicNitscheConstrains::MyNitscheVolume::tRee
AdaptiveKDTree tRee
Definition: NitschePeriodicMethod.hpp:296
PeriodicNitscheConstrains::PeriodicFace::commonData
CommonData & commonData
Definition: NitschePeriodicMethod.hpp:156
PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral::stressHomoGhostVector
Vec stressHomoGhostVector
Definition: NitschePeriodicMethod.hpp:1024
PeriodicNitscheConstrains::OpRhsPeriodicNormal::F
vector< Vec > & F
Definition: NitschePeriodicMethod.hpp:579
convert.type
type
Definition: convert.py:64
PeriodicNitscheConstrains::CommonData::MultiIndexData::MultiIndexData
MultiIndexData(int gg, int side, EntityType type)
Definition: NitschePeriodicMethod.hpp:53
PeriodicNitscheConstrains::CommonData::sTress
vector< MatrixDouble > sTress
Definition: NitschePeriodicMethod.hpp:43
PeriodicNitscheConstrains::OpLhsPeriodicNormal::kMatrix0
MatrixDouble kMatrix0
Definition: NitschePeriodicMethod.hpp:765
NitscheMethod::OpCommon::commonData
NonlinearElasticElement::CommonData & commonData
Definition: NitscheMethod.hpp:436
PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral::homoStress
VectorDouble homoStress
Definition: NitschePeriodicMethod.hpp:1049
MoFEM::PetscData::A
Mat A
Definition: LoopMethods.hpp:49
MoFEM::VolumeElementForcesAndSourcesCore::conn
const EntityHandle * conn
Definition: VolumeElementForcesAndSourcesCore.hpp:101
ShapeMBTRIQ
PetscErrorCode ShapeMBTRIQ(double *N, const double *X, const double *Y, const int G_DIM)
Definition: fem_tools.c:780
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:197
ShapeFaceNormalMBTRI
PetscErrorCode ShapeFaceNormalMBTRI(double *diffN, const double *coords, double *normal)
Definition: fem_tools.c:229
PeriodicNitscheConstrains::OpRhsPeriodicNormal::doWork
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: NitschePeriodicMethod.hpp:622
PeriodicNitscheConstrains::CommonData::skinFaces
Range skinFaces
Definition: NitschePeriodicMethod.hpp:32
NitscheMethod::OpBasicCommon::nitscheBlockData
BlockData & nitscheBlockData
Definition: NitscheMethod.hpp:397
NitscheMethod::MyVolumeFE
Definition of volume element.
Definition: NitscheMethod.hpp:240
PeriodicNitscheConstrains::CommonData::Container
multi_index_container< MultiIndexData, indexed_by< ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(MultiIndexData, int, MultiIndexData::gG) >, ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(MultiIndexData, int, MultiIndexData::sIde) >, ordered_non_unique< BOOST_MULTI_INDEX_MEMBER(MultiIndexData, EntityType, MultiIndexData::tYpe) >, ordered_unique< composite_key< MultiIndexData, member< MultiIndexData, int,&MultiIndexData::gG >, member< MultiIndexData, int,&MultiIndexData::sIde >, member< MultiIndexData, EntityType,&MultiIndexData::tYpe > > > > > Container
Definition: NitschePeriodicMethod.hpp:76
PeriodicNitscheConstrains::OpLhsPeriodicNormal::OpLhsPeriodicNormal
OpLhsPeriodicNormal(const string field_name, NitscheMethod::BlockData &nitsche_block_data, NitscheMethod::CommonData &nitsche_common_data, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, CommonData &periodic_common_data, bool field_disp=true)
Definition: NitschePeriodicMethod.hpp:744
PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral::tRaction
VectorDouble tRaction
Definition: NitschePeriodicMethod.hpp:1048
MoFEM::ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Definition: ForcesAndSourcesCore.hpp:1000
returnNumeredEntFiniteElement_multiIndex
NumeredEntFiniteElement_multiIndex & returnNumeredEntFiniteElement_multiIndex(T &fes)
Definition: NitscheMethod.hpp:33
PeriodicNitscheConstrains::OpRhsPeriodicNormal::calcualteDMatrix
static PetscErrorCode calcualteDMatrix(VectorAdaptor coords, MatrixDouble &mat_d)
Definition: NitschePeriodicMethod.hpp:612
PeriodicNitscheConstrains::OpGetVolumeData::commonData
NonlinearElasticElement::CommonData & commonData
Definition: NitschePeriodicMethod.hpp:190
NitscheMethod::BlockData::gamma
double gamma
Penalty term, see .
Definition: NitscheMethod.hpp:83
PeriodicNitscheConstrains::CommonData::MultiIndexData::gG
int gG
Definition: NitschePeriodicMethod.hpp:46
MoFEM::ForcesAndSourcesCore::UserDataOperator::getCoordsAtGaussPts
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
Definition: ForcesAndSourcesCore.hpp:1265
PeriodicNitscheConstrains::MyNitscheVolume::doAdditionalJobWhenGuassPtsAreCalulated
PetscErrorCode doAdditionalJobWhenGuassPtsAreCalulated()
Definition: NitschePeriodicMethod.hpp:325
PeriodicNitscheConstrains::OpRhsPeriodicNormal
Definition: NitschePeriodicMethod.hpp:576
NonlinearElasticElement::BlockData
data for calculation heat conductivity and heat capacity elements
Definition: HookeElement.hpp:32
PeriodicNitscheConstrains::OpCalculateHomogenisedStressVolumeIntegral::doWork
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: NitschePeriodicMethod.hpp:975
PeriodicNitscheConstrains::OpLhsPeriodicNormal::periodicCommonData
CommonData & periodicCommonData
Definition: NitschePeriodicMethod.hpp:742
MoFEM::Types::VectorAdaptor
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
PeriodicNitscheConstrains::CommonData::volumeElemName
string volumeElemName
Definition: NitschePeriodicMethod.hpp:31
PeriodicNitscheConstrains::PeriodicVolume::getRule
int getRule(int order)
Definition: NitschePeriodicMethod.hpp:244
NitscheMethod::MyVolumeFE::MyVolumeFE
MyVolumeFE(MoFEM::Interface &m_field, BlockData &block_data, CommonData &common_data)
Definition: NitscheMethod.hpp:252
PeriodicNitscheConstrains::MyNitscheVolume::iNit
bool iNit
Definition: NitschePeriodicMethod.hpp:295
FaceElementForcesAndSourcesCore
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEMethod
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Definition: ForcesAndSourcesCore.hpp:1042
PeriodicNitscheConstrains::PeriodicVolume::commonData
CommonData & commonData
Definition: NitschePeriodicMethod.hpp:238
PeriodicNitscheConstrains::PeriodicFace::setGaussPts
PetscErrorCode setGaussPts(int order)
Definition: NitschePeriodicMethod.hpp:163
PeriodicNitscheConstrains::OpRhsPeriodicNormal::nF0
VectorDouble nF0[6]
Definition: NitschePeriodicMethod.hpp:604
PeriodicNitscheConstrains::MyNitscheVolume
Definition: NitschePeriodicMethod.hpp:289
N
const int N
Definition: speed_test.cpp:3
NonlinearElasticElement::CommonData::sTress
std::vector< MatrixDouble3by3 > sTress
Definition: NonLinearElasticElement.hpp:111
NitscheMethod::CommonData::facesFePtr
std::vector< boost::shared_ptr< const NumeredEntFiniteElement > > facesFePtr
Definition: NitscheMethod.hpp:94
Range
face_coords
static const double face_coords[4][9]
Definition: forces_and_sources_h1_continuity_check.cpp:13
NitscheMethod::OpBasicCommon::nitscheCommonData
CommonData & nitscheCommonData
Definition: NitscheMethod.hpp:398
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PeriodicNitscheConstrains::CommonData::MultiIndexData::shapeFunctions
VectorDouble shapeFunctions
Definition: NitschePeriodicMethod.hpp:51
NitscheMethod::CommonData::hoCoordsAtGaussPts
std::vector< MatrixDouble > hoCoordsAtGaussPts
Definition: NitscheMethod.hpp:100
PeriodicNitscheConstrains::OpLhsPeriodicNormal::kMatrixOff
MatrixDouble kMatrixOff
Definition: NitschePeriodicMethod.hpp:766
NitscheMethod::OpCommon::OpCommon
OpCommon(const std::string field_name, BlockData &nitsche_block_data, CommonData &nitsche_common_data, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, bool field_disp, const char type)
Definition: NitscheMethod.hpp:438
PeriodicNitscheConstrains::CommonData::stressJacobian
vector< MatrixDouble > stressJacobian
Definition: NitschePeriodicMethod.hpp:42
NonlinearElasticElement::OpRhsPiolaKirchhoff
Definition: NonLinearElasticElement.hpp:520
PeriodicNitscheConstrains::CommonData::localCoordsMap
map< EntityHandle, vector< VectorDouble > > localCoordsMap
Definition: NitschePeriodicMethod.hpp:35
PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral::OpCalculateHomogenisedStressSurfaceIntegral
OpCalculateHomogenisedStressSurfaceIntegral(const string field_name, NitscheMethod::BlockData &nitsche_block_data, NitscheMethod::CommonData &nitsche_common_data, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, Vec stress_homo_ghost_vector, bool field_disp=0)
Definition: NitschePeriodicMethod.hpp:1026
PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral::doWork
PetscErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: NitschePeriodicMethod.hpp:1054
PeriodicNitscheConstrains::OpGetFaceData::fieldDisp
bool fieldDisp
Definition: NitschePeriodicMethod.hpp:86
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral::resDisp
VectorDouble resDisp
Definition: NitschePeriodicMethod.hpp:1051
NitscheMethod::CommonData::rAy
std::vector< VectorDouble > rAy
Definition: NitscheMethod.hpp:101
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
NitscheMethod::OpBasicCommon::gammaH
double gammaH
Definition: NitscheMethod.hpp:422
NonlinearElasticElement::BlockData::tEts
Range tEts
constrains elements in block set
Definition: HookeElement.hpp:36
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::BasicMethod::problemPtr
const Problem * problemPtr
raw pointer to problem
Definition: LoopMethods.hpp:250
PeriodicNitscheConstrains::OpRhsPeriodicNormal::nF1
VectorDouble nF1[6]
Definition: NitschePeriodicMethod.hpp:604
PeriodicNitscheConstrains::PeriodicVolume::setGaussPts
PetscErrorCode setGaussPts(int order)
Definition: NitschePeriodicMethod.hpp:245
PeriodicNitscheConstrains::OpCalculateHomogenisedStressVolumeIntegral
Definition: NitschePeriodicMethod.hpp:957
NitscheMethod::CommonData::faceNormals
std::vector< MatrixDouble > faceNormals
Definition: NitscheMethod.hpp:96
PeriodicNitscheConstrains::MyNitscheVolume::MyNitscheVolume
MyNitscheVolume(MoFEM::Interface &m_field, NitscheMethod::BlockData &block_data, NitscheMethod::CommonData &common_data, CommonData &periodic_common_data)
Definition: NitschePeriodicMethod.hpp:298
PeriodicNitscheConstrains::CommonData::inTetFaceGaussPtsNumber
map< EntityHandle, vector< int > > inTetFaceGaussPtsNumber
Definition: NitschePeriodicMethod.hpp:36
NitscheMethod::CommonData::faceGaussPts
std::vector< MatrixDouble > faceGaussPts
Definition: NitscheMethod.hpp:97
PeriodicNitscheConstrains::MyNitscheVolume::periodicFace
PeriodicFace periodicFace
Definition: NitschePeriodicMethod.hpp:292
PeriodicNitscheConstrains::CommonData::hoCoordsAtGaussPts
vector< VectorDouble > hoCoordsAtGaussPts
Definition: NitschePeriodicMethod.hpp:41
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
PeriodicNitscheConstrains::OpGetVolumeData
Definition: NitschePeriodicMethod.hpp:188
PeriodicNitscheConstrains::MyNitscheVolume::periodicVolume
PeriodicVolume periodicVolume
Definition: NitschePeriodicMethod.hpp:293
PeriodicNitscheConstrains::MyNitscheVolume::iNitialise
PetscErrorCode iNitialise()
Definition: NitschePeriodicMethod.hpp:316
PeriodicNitscheConstrains::CommonData
Definition: NitschePeriodicMethod.hpp:29
ShapeDiffMBTRIQ
PetscErrorCode ShapeDiffMBTRIQ(double *diffN, const double *X, const double *Y, const int G_DIM)
Definition: fem_tools.c:795
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getJac
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
Definition: VolumeElementForcesAndSourcesCore.hpp:170
PeriodicNitscheConstrains::MyNitscheVolume::treeRootSet
EntityHandle treeRootSet
Definition: NitschePeriodicMethod.hpp:315
PeriodicNitscheConstrains::CommonData::MultiIndexData::tYpe
EntityType tYpe
Definition: NitschePeriodicMethod.hpp:48
PeriodicNitscheConstrains::CommonData::dIstance
map< EntityHandle, double > dIstance
Definition: NitschePeriodicMethod.hpp:38
NitscheMethod
Basic implementation of Nitsche's method.
Definition: NitscheMethod.hpp:61
PeriodicNitscheConstrains::OpCalculateHomogenisedStressVolumeIntegral::homoStress
VectorDouble homoStress
Definition: NitschePeriodicMethod.hpp:973
PeriodicNitscheConstrains::CommonData::volumesContainer
Container volumesContainer
Definition: NitschePeriodicMethod.hpp:78
NitscheMethod::OpCommon::getTractionVariance
MoFEMErrorCode getTractionVariance(int gg, int fgg, int ff, MatrixDouble &jac, MatrixDouble &trac)
Definition: NitscheMethod.hpp:486
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
PeriodicNitscheConstrains::CommonData::MultiIndexData::dofOrders
ublas::vector< int > dofOrders
Definition: NitschePeriodicMethod.hpp:50
NitscheMethod::OpCommon
Calculate jacobian and variation of tractions.
Definition: NitscheMethod.hpp:433
PeriodicNitscheConstrains::OpCalculateHomogenisedStressSurfaceIntegral::dispOnOtherSide
VectorDouble dispOnOtherSide
Definition: NitschePeriodicMethod.hpp:1052
NitscheMethod::OpBasicCommon::getGammaH
virtual MoFEMErrorCode getGammaH(double gamma, int gg)
Definition: NitscheMethod.hpp:423
PeriodicNitscheConstrains::OpGetFaceData::doWork
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: NitschePeriodicMethod.hpp:94
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
PeriodicNitscheConstrains::CommonData::MultiIndexData::sIde
int sIde
Definition: NitschePeriodicMethod.hpp:47
tol
double tol
Definition: mesh_smoothing.cpp:26
NitscheMethod::BlockData::phi
double phi
Nitsche method parameter, see .
Definition: NitscheMethod.hpp:84
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
ShapeMBTRI
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:182
NonlinearElasticElement::CommonData
common data used by volume elements
Definition: NonLinearElasticElement.hpp:105