v0.13.0
NitscheMethod.hpp
Go to the documentation of this file.
1 /** \file NitscheMethod.hpp
2  * \ingroup nitsche_method
3  * \brief Basic implementation of Nitsche method
4  *
5  * \note This is quite an old implementation, kept here for compatibility with
6  * older users modules. Currently, MoFEM has a generic implementation for
7  * integration over the skeleton, and that one should be used, instead one
8  * presented below.
9  *
10  */
11 
12 /*
13  * This file is part of MoFEM.
14  * MoFEM is free software: you can redistribute it and/or modify it under
15  * the terms of the GNU Lesser General Public License as published by the
16  * Free Software Foundation, either version 3 of the License, or (at your
17  * option) any later version.
18  *
19  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
20  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
22  * License for more details.
23  *
24  * You should have received a copy of the GNU Lesser General Public
25  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
26 
27 #ifndef __NITCHE_BOUNDARY_CONDITIONS_HPP__
28 #define __NITCHE_BOUNDARY_CONDITIONS_HPP__
29 
30 // This is fix to have back compatibility
31 template <typename T>
34  return fes;
35 }
36 
37 template <>
41  return fes;
42 }
43 
44 template <>
47  boost::shared_ptr<NumeredEntFiniteElement_multiIndex>>(
48  boost::shared_ptr<NumeredEntFiniteElement_multiIndex> &fes) {
49  return *fes;
50 }
51 
52 /** \brief Basic implementation of Nitsche's method
53  * \ingroup nitsche_method
54 
55  This implementation have been used to enforce periodic constrains, see \ref
56  nitsche_periodic. For theoretical foundations of the method explanations see
57  \cite nitsche_method_hal and \cite juntunen2009nitsche. For scratch book
58  derivations see <a href="nitsche_bc_for_ch.pdf" target="_blank"><b>link</b></a>.
59 
60 */
61 struct NitscheMethod {
62 
64  int addToRule;
67  int getRule(int order) { return order + addToRule; }
68  /*int getRule(int order) { return -1; }
69  MoFEMErrorCode setGaussPts(int order) {
70  MoFEMFunctionBeginHot;
71  int rule = order+addToRule+1;
72  int nb_gauss_pts = triangle_ncc_order_num(rule);
73  gaussPts.resize(3,nb_gauss_pts,false);
74  triangle_ncc_rule(rule,nb_gauss_pts,&gaussPts(0,0),&gaussPts(2,0));
75  MoFEMFunctionReturnHot(0);
76  }*/
77  };
78 
79  /** \brief Block data for Nitsche method
80  * \ingroup nitsche_method
81  */
82  struct BlockData {
83  double gamma; ///< Penalty term, see \cite nitsche_method_hal
84  double phi; ///< Nitsche method parameter, see \cite nitsche_method_hal
85  string faceElemName; ///< name of element face
86  Range fAces; ///< faces on which constrain is applied
87  };
88 
89  /** \brief Common data shared between finite element operators
90  */
91  struct CommonData {
93  std::vector<EntityHandle> fAces;
94  std::vector<boost::shared_ptr<const NumeredEntFiniteElement>> facesFePtr;
95  std::vector<VectorDouble> cOords;
96  std::vector<MatrixDouble> faceNormals;
97  std::vector<MatrixDouble> faceGaussPts;
98  std::vector<std::vector<int>> inTetFaceGaussPtsNumber;
99  std::vector<MatrixDouble> coordsAtGaussPts;
100  std::vector<MatrixDouble> hoCoordsAtGaussPts;
101  std::vector<VectorDouble> rAy;
103 
104  /** \brief projection matrix
105 
106  First index is face index, then integration point and
107  projection matrix 3x3
108 
109  */
110  std::vector<std::vector<MatrixDouble>> P;
111 
112  /** \brief derivative of projection matrix in respect DoFs
113  This is EntityType, face, gauss point at face.
114  Matrix has 3 rows (components of displacements)
115  Matrix has columns equal to number of DoFs on entity
116  */
117  std::map<EntityType, std::vector<std::vector<MatrixDouble>>> dP;
118 
120 
121  struct MultiIndexData {
122  int gG;
123  int sIde;
124  EntityType tYpe;
125  ublas::vector<int> iNdices;
126  ublas::vector<int> dofOrders;
129  MultiIndexData(int gg, int side, EntityType type)
130  : gG(gg), sIde(side), tYpe(type) {}
131  };
132 
133  typedef multi_index_container<
134  MultiIndexData,
135  indexed_by<
136  ordered_non_unique<BOOST_MULTI_INDEX_MEMBER(MultiIndexData, int,
138  ordered_non_unique<BOOST_MULTI_INDEX_MEMBER(MultiIndexData, int,
140  ordered_non_unique<BOOST_MULTI_INDEX_MEMBER(
141  MultiIndexData, EntityType, MultiIndexData::tYpe)>,
142  ordered_unique<composite_key<
143  MultiIndexData,
144  member<MultiIndexData, int, &MultiIndexData::gG>,
145  member<MultiIndexData, int, &MultiIndexData::sIde>,
146  member<MultiIndexData, EntityType, &MultiIndexData::tYpe>>>>>
149  };
150 
151  /** \brief Get integration pts data on face
152  */
155 
157 
158  OpGetFaceData(CommonData &common_data)
160  "DISPLACEMENT", OPROW),
161  commonData(common_data) {}
162 
163  MoFEMErrorCode doWork(int side, EntityType type,
166 
167  int faceInRespectToTet = getFEMethod()->nInTheLoop;
168  int nb_face_gauss_pts = getGaussPts().size2();
169 
170  try {
171  if (type == MBVERTEX) {
172  commonData.faceGaussPts.resize(4);
173  commonData.faceGaussPts[faceInRespectToTet] = getGaussPts();
174  for (int fgg = 0; fgg < nb_face_gauss_pts; fgg++) {
175  int gg = commonData.nbTetGaussPts;
176  commonData.inTetFaceGaussPtsNumber[faceInRespectToTet].push_back(
177  gg);
179  }
180  commonData.faceNormals.resize(4);
181  commonData.faceNormals[faceInRespectToTet] =
182  0.5 * getNormalsAtGaussPtss();
183  commonData.hoCoordsAtGaussPts.resize(4);
184  commonData.hoCoordsAtGaussPts[faceInRespectToTet] =
185  getCoordsAtGaussPts();
186  commonData.cOords.resize(4);
187  commonData.cOords[faceInRespectToTet] = getCoords();
188  commonData.coordsAtGaussPts.resize(4);
189  commonData.coordsAtGaussPts[faceInRespectToTet] =
190  getCoordsAtGaussPts();
191  commonData.rAy.resize(4);
192  commonData.rAy[faceInRespectToTet] = -getNormal();
193  commonData.rAy[faceInRespectToTet] /= norm_2(getNormal());
194  }
195  } catch (const std::exception &ex) {
196  std::ostringstream ss;
197  ss << "throw in method: " << ex.what() << std::endl;
198  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
199  }
200 
201  try {
202  for (int fgg = 0; fgg < nb_face_gauss_pts; fgg++) {
203  int gg = commonData.inTetFaceGaussPtsNumber[faceInRespectToTet][fgg];
204  // std::cerr << fgg << " " << gg << " " << side << " " << type <<
205  // std::endl;
206  CommonData::MultiIndexData gauss_pt_data(gg, side, type);
207  std::pair<CommonData::Container::iterator, bool> p;
208  p = commonData.facesContainer.insert(gauss_pt_data);
209  if (!p.second) {
210  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
211  "data not inserted");
212  }
214  const_cast<CommonData::MultiIndexData &>(*p.first);
215  VectorDouble &shape_fun = p_data.shapeFunctions;
216  int nb_shape_fun = data.getN().size2();
217  shape_fun.resize(nb_shape_fun);
218  // std::cerr << "nb_shape_fun " << nb_shape_fun << std::endl;
219  if (nb_shape_fun) {
220  cblas_dcopy(nb_shape_fun, &data.getN()(fgg, 0), 1, &shape_fun[0],
221  1);
222  }
223  p_data.iNdices = data.getIndices();
224  p_data.dofOrders.resize(data.getFieldDofs().size(), false);
225  for (unsigned int dd = 0; dd < data.getFieldDofs().size(); dd++) {
226  p_data.dofOrders[dd] = data.getFieldDofs()[dd]->getDofOrder();
227  }
228  // std::cerr << shape_fun << std::endl;
229  }
230  } catch (const std::exception &ex) {
231  std::ostringstream ss;
232  ss << "throw in method: " << ex.what() << std::endl;
233  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
234  }
236  }
237  };
238 
239  /// \brief Definition of volume element
241 
246 
250  }
251 
252  MyVolumeFE(MoFEM::Interface &m_field, BlockData &block_data,
253  CommonData &common_data)
255  blockData(block_data), commonData(common_data), faceFE(m_field),
256  addToRule(0) {
257  faceFE.getOpPtrVector().push_back(new OpGetFaceData(commonData));
258  }
259 
260  int getRule(int order) { return -1; };
261 
263 
265 
266  gaussPts.resize(4, 0, false);
267 
268  try {
270  commonData.fAces.resize(4);
272  for (int ff = 0; ff < 4; ff++) {
273  EntityHandle face;
274  rval = mField.get_moab().side_element(tet, 2, ff, face);
275  CHKERRG(rval);
276  if (blockData.fAces.find(face) != blockData.fAces.end()) {
277  commonData.fAces[ff] = face;
279  } else {
280  commonData.fAces[ff] = 0;
281  }
282  }
283  } catch (const std::exception &ex) {
284  std::ostringstream ss;
285  ss << "throw in method: " << ex.what() << std::endl;
286  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
287  }
288 
289  try {
290  commonData.facesFePtr.resize(4);
291  for (int ff = 0; ff < 4; ff++) {
292  if (commonData.fAces[ff] != 0) {
294  Composite_Name_And_Ent_mi_tag>::type::iterator it,
295  hi_it;
296 
298  problemPtr->numeredFiniteElements)
299  .get<Composite_Name_And_Ent_mi_tag>()
300  .lower_bound(boost::make_tuple(blockData.faceElemName,
301  commonData.fAces[ff]));
303  problemPtr->numeredFiniteElements)
304  .get<Composite_Name_And_Ent_mi_tag>()
305  .upper_bound(boost::make_tuple(blockData.faceElemName,
306  commonData.fAces[ff]));
308  problemPtr->numeredFiniteElements)
309  .get<Composite_Name_And_Ent_mi_tag>()
310  .end()) {
311  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
312  "No finite element found < %s >",
313  blockData.faceElemName.c_str());
314  }
315  commonData.facesFePtr[ff] = *it;
316  } else {
317  commonData.facesFePtr[ff].reset();
318  }
319  }
321  commonData.facesContainer.clear();
324  for (int ff = 0; ff < 4; ff++) {
325  if (commonData.facesFePtr[ff]) {
326  boost::shared_ptr<const NumeredEntFiniteElement> faceFEPtr =
328  faceFE.copyBasicMethod(*this);
329  faceFE.feName = blockData.faceElemName;
330  faceFE.nInTheLoop = ff;
331  faceFE.numeredEntFiniteElementPtr = faceFEPtr;
332  faceFE.dataPtr = faceFEPtr->sPtr->data_dofs;
333  faceFE.rowPtr = faceFEPtr->rows_dofs;
334  faceFE.colPtr = faceFEPtr->cols_dofs;
336  ierr = faceFE();
337  CHKERRG(ierr);
338  }
339  }
340  } catch (const std::exception &ex) {
341  std::ostringstream ss;
342  ss << "throw in method: " << ex.what() << std::endl;
343  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
344  }
345 
346  try {
347  int nb_gauss_pts = 0;
348  for (int ff = 0; ff < 4; ff++) {
349  if (!commonData.facesFePtr[ff])
350  continue;
351  nb_gauss_pts += commonData.faceGaussPts[ff].size2();
352  }
353  gaussPts.resize(4, nb_gauss_pts, false);
354  const double coords_tet[12] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
355  int gg = 0;
356  for (int ff = 0; ff < 4; ff++) {
357  if (!commonData.facesFePtr[ff])
358  continue;
359  int nb_gauss_face_pts = commonData.faceGaussPts[ff].size2();
360  // std::cerr << "nb_gauss_face_pts " << nb_gauss_face_pts <<
361  // std::endl;
362  for (int fgg = 0; fgg < nb_gauss_face_pts; fgg++, gg++) {
363  // std::cerr << ff << " gg " << gg << " fgg " << fgg << std::endl;
364  CommonData::Container::nth_index<3>::type::iterator sit;
365  sit = commonData.facesContainer.get<3>().find(
366  boost::make_tuple(gg, 0, MBVERTEX));
367  const VectorDouble &shape_fun = sit->shapeFunctions;
368  // std::cerr << shape_fun << std::endl;
369 
370  for (int dd = 0; dd < 3; dd++) {
371  gaussPts(dd, gg) =
372  shape_fun[0] * coords_tet[3 * dataH1.facesNodes(ff, 0) + dd] +
373  shape_fun[1] * coords_tet[3 * dataH1.facesNodes(ff, 1) + dd] +
374  shape_fun[2] * coords_tet[3 * dataH1.facesNodes(ff, 2) + dd];
375  }
376  gaussPts(3, gg) = commonData.faceGaussPts[ff](2, fgg);
377  }
378  }
379  } catch (const std::exception &ex) {
380  std::ostringstream ss;
381  ss << "throw in method: " << ex.what() << std::endl;
382  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
383  }
384 
386  CHKERRG(ierr);
387 
389  }
390  };
391 
392  /** \brief Basic operated shared between all Nitsche operators
393  */
396 
399  bool fieldDisp;
400 
401  OpBasicCommon(const std::string field_name, BlockData &nitsche_block_data,
402  CommonData &nitsche_common_data, bool field_disp,
403  const char type)
405  type),
406  nitscheBlockData(nitsche_block_data),
407  nitscheCommonData(nitsche_common_data), fieldDisp(field_disp) {
408  sYmm = false;
409  }
410 
411  double faceRadius;
412  virtual MoFEMErrorCode getFaceRadius(int ff) {
414  VectorDouble &coords = nitscheCommonData.cOords[ff];
415  double center[3];
416  tricircumcenter3d_tp(&coords[0], &coords[3], &coords[6], center, NULL,
417  NULL);
418  cblas_daxpy(3, -1, &coords[0], 1, center, 1);
419  faceRadius = cblas_dnrm2(3, center, 1);
421  }
422  double gammaH;
423  virtual MoFEMErrorCode getGammaH(double gamma, int gg) {
425  gammaH = gamma;
426  // gammaH*= faceRadius;
428  }
429  };
430 
431  /** \brief Calculate jacobian and variation of tractions
432  */
433  struct OpCommon : public OpBasicCommon {
434 
437 
438  OpCommon(const std::string field_name, BlockData &nitsche_block_data,
439  CommonData &nitsche_common_data,
441  NonlinearElasticElement::CommonData &common_data, bool field_disp,
442  const char type)
443  : OpBasicCommon(field_name, nitsche_block_data, nitsche_common_data,
444  field_disp, type),
445  dAta(data), commonData(common_data) {}
446 
453 
454  virtual MoFEMErrorCode calculateP(int gg, int fgg, int ff) {
457  }
458 
460  MatrixDouble &jac) {
462  try {
463  int nb = data.getFieldData().size();
464  jac.resize(9, nb, false);
465  jac.clear();
466  const MatrixAdaptor diffN = data.getDiffN(gg, nb / 3);
467  MatrixDouble &jac_stress = commonData.jacStress[gg];
468  for (int dd = 0; dd < nb / 3; dd++) {
469  for (int rr = 0; rr < 3; rr++) {
470  for (int ii = 0; ii < 9; ii++) {
471  for (int jj = 0; jj < 3; jj++) {
472  jac(ii, 3 * dd + rr) +=
473  jac_stress(ii, 3 * rr + jj) * diffN(dd, jj);
474  }
475  }
476  }
477  }
478  } catch (const std::exception &ex) {
479  std::ostringstream ss;
480  ss << "throw in method: " << ex.what() << std::endl;
481  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
482  }
484  }
485 
486  MoFEMErrorCode getTractionVariance(int gg, int fgg, int ff,
487  MatrixDouble &jac, MatrixDouble &trac) {
489  try {
491  3, ublas::shallow_array_adaptor<double>(
492  3, &nitscheCommonData.faceNormals[ff](fgg, 0)));
493  trac.resize(3, jac.size2());
494  trac.clear();
495  for (unsigned int dd2 = 0; dd2 < jac.size2(); dd2++) {
496  for (int nn = 0; nn < 3; nn++) {
497  trac(nn, dd2) =
498  cblas_ddot(3, &jac(3 * nn, dd2), jac.size2(), &normal[0], 1);
499  }
500  }
501  } catch (const std::exception &ex) {
502  std::ostringstream ss;
503  ss << "throw in method: " << ex.what() << std::endl;
504  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
505  }
507  }
508  };
509 
510  /** \brief Calculate Nitsche method terms on left hand side
511  * \ingroup nitsche_method
512  */
513  struct OpLhsNormal : public OpCommon {
514 
515  OpLhsNormal(const std::string field_name, BlockData &nitsche_block_data,
516  CommonData &nitsche_common_data,
519  bool field_disp)
520  : OpCommon(field_name, nitsche_block_data, nitsche_common_data, data,
521  common_data, field_disp, UserDataOperator::OPROWCOL) {}
522 
524  std::vector<MatrixDouble> kMatrixFace, kMatrixFace0, kMatrixFace1;
525 
526  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
527  EntityType col_type,
531 
532  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
533  dAta.tEts.end()) {
535  }
536  if (row_data.getIndices().size() == 0)
538  if (col_data.getIndices().size() == 0)
540  int nb_dofs_row = row_data.getIndices().size();
541  int nb_dofs_col = col_data.getIndices().size();
542  double gamma = nitscheBlockData.gamma;
543  double phi = nitscheBlockData.phi;
544 
545  kMatrix0.resize(nb_dofs_row, nb_dofs_col, false);
546  kMatrix1.resize(nb_dofs_row, nb_dofs_col, false);
547  kMatrix.resize(nb_dofs_row, nb_dofs_col, false);
548  kMatrix.clear();
549 
550  try {
551 
552  int gg = 0;
553  for (int ff = 0; ff < 4; ff++) {
554  if (!nitscheCommonData.facesFePtr[ff])
555  continue;
556  int nb_face_gauss_pts = nitscheCommonData.faceGaussPts[ff].size2();
557  ierr = getFaceRadius(ff);
558  CHKERRG(ierr);
559  kMatrix0.clear();
560  kMatrix1.clear();
561  for (int fgg = 0; fgg < nb_face_gauss_pts; fgg++, gg++) {
562  ierr = getGammaH(gamma, gg);
563  CHKERRG(ierr);
564  double val = getGaussPts()(3, gg);
565  ierr = getJac(row_data, gg, jAc_row);
566  CHKERRG(ierr);
567  ierr = getTractionVariance(gg, fgg, ff, jAc_row, tRac_v);
568  CHKERRG(ierr);
569  ierr = getJac(col_data, gg, jAc_col);
570  CHKERRG(ierr);
571  ierr = getTractionVariance(gg, fgg, ff, jAc_col, tRac_u);
572  CHKERRG(ierr);
574  3, ublas::shallow_array_adaptor<double>(
575  3, &nitscheCommonData.faceNormals[ff](fgg, 0)));
576  double area = cblas_dnrm2(3, &normal[0], 1);
577 
578  ierr = calculateP(gg, fgg, ff);
579  CHKERRG(ierr);
580  MatrixDouble &P = nitscheCommonData.P[ff][fgg];
581 
582  // P
583  for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
584  double n_row = row_data.getN()(gg, dd1);
585  for (int dd2 = 0; dd2 < nb_dofs_col / 3; dd2++) {
586  double n_col = col_data.getN()(gg, dd2);
587  for (int dd3 = 0; dd3 < 3; dd3++) {
588  for (int dd4 = 0; dd4 < 3; dd4++) {
589  if (!P(dd3, dd4))
590  continue;
591  kMatrix0(3 * dd1 + dd3, 3 * dd2 + dd4) +=
592  val * n_row * P(dd3, dd4) * n_col * area;
593  }
594  }
595  }
596  }
597  for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
598  double n_row = row_data.getN()(gg, dd1);
599  for (int dd2 = 0; dd2 < nb_dofs_col; dd2++) {
600  for (int dd3 = 0; dd3 < 3; dd3++) {
601  double t = cblas_ddot(3, &P(dd3, 0), 1, &tRac_u(0, dd2),
602  tRac_u.size2());
603  kMatrix1(3 * dd1 + dd3, dd2) -= val * n_row * t;
604  }
605  }
606  }
607  for (int dd1 = 0; dd1 < nb_dofs_row; dd1++) {
608  for (int dd2 = 0; dd2 < nb_dofs_col / 3; dd2++) {
609  double n_col = col_data.getN()(gg, dd2);
610  for (int dd3 = 0; dd3 < 3; dd3++) {
611  double t = cblas_ddot(3, &P(0, dd3), 3, &tRac_v(0, dd1),
612  tRac_v.size2());
613  kMatrix1(dd1, 3 * dd2 + dd3) -= phi * val * t * n_col;
614  }
615  }
616  }
617  // dP
618  if (nitscheCommonData.dP[col_type].empty() != 0) {
619  if (nitscheCommonData.dP[col_type].size() == 4) {
620  if (nitscheCommonData.dP[col_type][ff].size() ==
621  (unsigned int)nb_face_gauss_pts) {
622  MatrixDouble &dP = nitscheCommonData.dP[col_type][ff][fgg];
623  if (dP.size1() == 3 &&
624  dP.size2() == (unsigned int)nb_dofs_col) {
625  VectorDouble &u =
626  (commonData
628  for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
629  double n_row = row_data.getN()(gg, dd1);
630  for (int dd2 = 0; dd2 < nb_dofs_col / 3; dd2++) {
631  for (int dd3 = 0; dd3 < 3; dd3++) {
632  for (int dd4 = 0; dd4 < 3; dd4++) {
633  kMatrix0(3 * dd1 + dd3, 3 * dd2 + dd4) +=
634  val * n_row * u[dd4] * dP(dd4, 3 * dd2 + dd4);
635  }
636  }
637  }
638  }
639  VectorDouble t =
640  prod(trans(commonData.sTress[gg]),
641  normal); // this is Piola-Kirchhoff I stress
642  for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
643  double n_row = row_data.getN()(gg, dd1);
644  for (int dd2 = 0; dd2 < nb_dofs_col / 3; dd2++) {
645  for (int dd3 = 0; dd3 < 3; dd3++) {
646  for (int dd4 = 0; dd4 < 3; dd4++) {
647  kMatrix1(3 * dd1 + dd3, 3 * dd2 + dd4) +=
648  val * n_row * t[dd4] * dP(dd4, 3 * dd2 + dd4);
649  }
650  }
651  }
652  }
653  }
654  }
655  }
656  }
657  }
658 
659  kMatrix0 /= gammaH;
660  noalias(kMatrix) += kMatrix0 + kMatrix1;
661  }
662 
663  if (gg != (int)getGaussPts().size2()) {
664  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
665  "wrong number of gauss pts");
666  }
667 
668  ierr = MatSetValues(
669  getFEMethod()->snes_B, nb_dofs_row, &row_data.getIndices()[0],
670  nb_dofs_col, &col_data.getIndices()[0], &kMatrix(0, 0), ADD_VALUES);
671  CHKERRG(ierr);
672 
673  } catch (const std::exception &ex) {
674  std::ostringstream ss;
675  ss << "throw in method: " << ex.what() << std::endl;
676  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
677  }
678 
680  }
681  };
682 
683  /** \brief Calculate Nitsche method terms on right hand side
684  * \ingroup nitsche_method
685  */
686  struct OpRhsNormal : public OpCommon {
687 
688  OpRhsNormal(const std::string field_name, BlockData &nitsche_block_data,
689  CommonData &nitsche_common_data,
692  bool field_disp)
693  : OpCommon(field_name, nitsche_block_data, nitsche_common_data, data,
694  common_data, field_disp, UserDataOperator::OPROW) {}
695 
697 
698  MoFEMErrorCode doWork(int row_side, EntityType row_type,
701 
702  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
703  dAta.tEts.end()) {
705  }
706  if (row_data.getIndices().size() == 0)
708  int nb_dofs_row = row_data.getIndices().size();
709  double gamma = nitscheBlockData.gamma;
710  double phi = nitscheBlockData.phi;
711 
712  try {
713 
714  nF.resize(nb_dofs_row, false);
715  nF.clear();
716 
717  int gg = 0;
718  for (int ff = 0; ff < 4; ff++) {
719  if (!nitscheCommonData.facesFePtr[ff])
720  continue;
721  int nb_face_gauss_pts = nitscheCommonData.faceGaussPts[ff].size2();
722  ierr = getFaceRadius(ff);
723  CHKERRG(ierr);
724  for (int fgg = 0; fgg < nb_face_gauss_pts; fgg++, gg++) {
725 
726  ierr = getGammaH(gamma, gg);
727  CHKERRG(ierr);
728  double val = getGaussPts()(3, gg);
729  ierr = getJac(row_data, gg, jAc_row);
730  CHKERRG(ierr);
731  ierr = getTractionVariance(gg, fgg, ff, jAc_row, tRac_v);
732  CHKERRG(ierr);
734  3, ublas::shallow_array_adaptor<double>(
735  3, &nitscheCommonData.faceNormals[ff](fgg, 0)));
736  double area = cblas_dnrm2(3, &normal[0], 1);
737 
738  ierr = calculateP(gg, fgg, ff);
739  CHKERRG(ierr);
740  MatrixDouble &P = nitscheCommonData.P[ff][fgg];
741 
742  VectorDouble &u =
744  VectorDouble t = prod(trans(commonData.sTress[gg]),
745  normal); // this is Piola-Kirchhoff I stress
746 
747  for (int dd1 = 0; dd1 < nb_dofs_row / 3; dd1++) {
748  double n_row = row_data.getN()(gg, dd1);
749  for (int dd2 = 0; dd2 < 3; dd2++) {
750  nF[3 * dd1 + dd2] +=
751  gammaH *
752  (val * area * n_row *
753  (P(dd2, 0) * u[0] + P(dd2, 1) * u[1] + P(dd2, 2) * u[2]));
754  nF[3 * dd1 + dd2] +=
755  val * n_row *
756  (P(dd2, 0) * t[0] + P(dd2, 1) * t[1] + P(dd2, 2) * t[2]);
757  nF[3 * dd1 + dd2] +=
758  phi * tRac_v(dd2, 3 * dd1 + dd2) *
759  (P(dd2, 0) * u[0] + P(dd2, 1) * u[1] + P(dd2, 2) * u[2]);
760  }
761  }
762  }
763  }
764 
765  ierr = VecSetValues(getFEMethod()->snes_f, nb_dofs_row,
766  &row_data.getIndices()[0], &nF[0], ADD_VALUES);
767  CHKERRG(ierr);
768 
769  } catch (const std::exception &ex) {
770  std::ostringstream ss;
771  ss << "throw in method: " << ex.what() << std::endl;
772  SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
773  }
774 
776  }
777  };
778 };
779 
780 #endif // __NITCHE_BOUNDARY_CONDITIONS_HPP__
781 
782 /***************************************************************************/ /**
783 * \defgroup nitsche_method Nitsche Method
784 * \ingroup user_modules
785 ******************************************************************************/
static Index< 'p', 3 > p
ForcesAndSourcesCore::UserDataOperator UserDataOperator
NumeredEntFiniteElement_multiIndex & returnNumeredEntFiniteElement_multiIndex(T &fes)
NumeredEntFiniteElement_multiIndex & returnNumeredEntFiniteElement_multiIndex< NumeredEntFiniteElement_multiIndex >(NumeredEntFiniteElement_multiIndex &fes)
EntitiesFieldData::EntData EntData
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
multi_index_container< boost::shared_ptr< NumeredEntFiniteElement >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_EntFiniteElement, UId, &NumeredEntFiniteElement::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > >, ordered_non_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_RefEntity, EntityHandle, &NumeredEntFiniteElement::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< NumeredEntFiniteElement, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef >, const_mem_fun< NumeredEntFiniteElement::interface_type_RefEntity, EntityHandle, &NumeredEntFiniteElement::getEnt > > >, ordered_non_unique< tag< Composite_Name_And_Part_mi_tag >, composite_key< NumeredEntFiniteElement, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > > > > > NumeredEntFiniteElement_multiIndex
MultiIndex for entities for NumeredEntFiniteElement.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
const double T
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:77
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
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:126
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:143
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
constexpr double t
plate stiffness
Definition: plate.cpp:76
static double phi
const Problem * problemPtr
raw pointer to problem
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
MatrixInt facesNodes
nodes on finite element faces
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
MatrixDouble gaussPts
Matrix of integration points.
FTensor::Tensor2< double *, 3, 3 > & getJac()
get element Jacobian
Block data for Nitsche method.
double gamma
Penalty term, see .
string faceElemName
name of element face
double phi
Nitsche method parameter, see .
Range fAces
faces on which constrain is applied
MultiIndexData(int gg, int side, EntityType type)
Common data shared between finite element operators.
std::vector< std::vector< int > > inTetFaceGaussPtsNumber
std::vector< boost::shared_ptr< const NumeredEntFiniteElement > > facesFePtr
std::vector< MatrixDouble > faceNormals
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
std::vector< VectorDouble > rAy
std::vector< MatrixDouble > coordsAtGaussPts
std::map< EntityType, std::vector< std::vector< MatrixDouble > > > dP
derivative of projection matrix in respect DoFs This is EntityType, face, gauss point at face....
std::vector< std::vector< MatrixDouble > > P
projection matrix
std::vector< VectorDouble > cOords
std::vector< MatrixDouble > hoCoordsAtGaussPts
std::vector< MatrixDouble > faceGaussPts
std::vector< EntityHandle > fAces
int getRule(int order)
MyFace(MoFEM::Interface &m_field)
Definition of volume element.
MyVolumeFE(MoFEM::Interface &m_field, BlockData &block_data, CommonData &common_data)
virtual MoFEMErrorCode doAdditionalJobWhenGuassPtsAreCalulated()
MoFEMErrorCode setGaussPts(int order)
Basic operated shared between all Nitsche operators.
virtual MoFEMErrorCode getFaceRadius(int ff)
OpBasicCommon(const std::string field_name, BlockData &nitsche_block_data, CommonData &nitsche_common_data, bool field_disp, const char type)
virtual MoFEMErrorCode getGammaH(double gamma, int gg)
Calculate jacobian and variation of tractions.
NonlinearElasticElement::CommonData & commonData
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)
MoFEMErrorCode getJac(DataForcesAndSourcesCore::EntData &data, int gg, MatrixDouble &jac)
virtual MoFEMErrorCode calculateP(int gg, int fgg, int ff)
NonlinearElasticElement::BlockData & dAta
MoFEMErrorCode getTractionVariance(int gg, int fgg, int ff, MatrixDouble &jac, MatrixDouble &trac)
Get integration pts data on face.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpGetFaceData(CommonData &common_data)
Calculate Nitsche method terms on left hand side.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
std::vector< MatrixDouble > kMatrixFace
std::vector< MatrixDouble > kMatrixFace0
OpLhsNormal(const std::string field_name, BlockData &nitsche_block_data, CommonData &nitsche_common_data, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, bool field_disp)
std::vector< MatrixDouble > kMatrixFace1
Calculate Nitsche method terms on right hand side.
OpRhsNormal(const std::string field_name, BlockData &nitsche_block_data, CommonData &nitsche_common_data, NonlinearElasticElement::BlockData &data, NonlinearElasticElement::CommonData &common_data, bool field_disp)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Basic implementation of Nitsche's method.
data for calculation heat conductivity and heat capacity elements
Range tEts
constrains elements in block set
common data used by volume elements
std::vector< MatrixDouble > jacStress
this is simply material tangent operator
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::vector< MatrixDouble3by3 > sTress
void tricircumcenter3d_tp(a, b, c, circumcenter, double *xi, double *eta)