v0.13.0
ContactOperators.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 namespace OpContactTools {
16 
17 
18 //! [Common data]
20  boost::shared_ptr<MatrixDouble> contactStressPtr;
21  boost::shared_ptr<MatrixDouble> contactStressDivergencePtr;
22  boost::shared_ptr<MatrixDouble> contactTractionPtr;
23  boost::shared_ptr<MatrixDouble> contactNormalPtr;
24  boost::shared_ptr<VectorDouble> contactGapPtr;
25  boost::shared_ptr<MatrixDouble> contactGapDiffPtr;
26  boost::shared_ptr<MatrixDouble> contactNormalDiffPtr;
27  // boost::shared_ptr<MatrixDouble> mDispPtr;
28  Tag rollerTag;
29 };
30 //! [Common data]
31 
33  OpPostProcVertex(const std::string field_name,
34  boost::shared_ptr<CommonData> common_data_ptr,
35  moab::Interface *moab_inst);
36  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
37 
38 private:
39  boost::shared_ptr<CommonData> commonDataPtr;
41 };
42 
44  OpConstrainBoundaryRhs(const std::string field_name,
45  boost::shared_ptr<CommonData> common_data_ptr);
47 
48 private:
49  boost::shared_ptr<CommonData> commonDataPtr;
50 };
51 
53  OpConstrainBoundaryLhs_dU(const std::string row_field_name,
54  const std::string col_field_name,
55  boost::shared_ptr<CommonData> common_data_ptr);
56  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
57 
58 
59 private:
60  boost::shared_ptr<CommonData> commonDataPtr;
61 };
62 
65  const std::string row_field_name, const std::string col_field_name,
66  boost::shared_ptr<CommonData> common_data_ptr);
67  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
68 
69 
70 private:
71  boost::shared_ptr<CommonData> commonDataPtr;
72 };
73 
76  const std::string field_name,
77  boost::shared_ptr<CommonData> common_data_ptr,
78  bool update);
79  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
80 
81 private:
83  boost::shared_ptr<CommonData> commonDataPtr;
84  bool uPdate;
85 };
86 
87 
88 template <typename T>
89 inline double gap0(double g, FTensor::Tensor1<T, 3> &t_disp,
90  FTensor::Tensor1<T, 3> &t_normal) {
91  return -g + t_disp(i) * t_normal(i);
92 }
93 
94 template <typename T1, typename T2>
95 inline double normal_traction(Tensor1<T1, 3> &t_traction,
96  Tensor1<T2, 3> &t_normal) {
97  return -t_traction(i) * t_normal(i);
98 }
99 
100 inline double sign(double x) {
101  if (x == 0)
102  return 0;
103  else if (x > 0)
104  return 1;
105  else
106  return -1;
107 };
108 
109 inline double w(const double g, const double t) {
110  return g - (*cache).cn_cont * t;
111 }
112 
113 inline double constraint(double g0, double g, double &&t) {
114  return (w(g, t) + std::abs(w(g, t))) / 2 + g0;
115 };
116 
117 inline double diff_constraints_dtraction(double g, double &&t) {
118  return -(*cache).cn_cont * (1 + sign(w(g, t))) / 2;
119 }
120 
121 inline double diff_constraints_dgap(double g, double &&t) {
122  return (sign(w(g, t)) - 1) / 2;
123 }
124 
126  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
127  : BoundaryEleOpAssembly(field_name, field_name, BoundaryEleOp::OPROW),
128  commonDataPtr(common_data_ptr) {}
129 
132 
133  const size_t nb_gauss_pts = getGaussPts().size2();
134  const size_t nb_dofs = data.getIndices().size();
135 
136  if (nb_dofs) {
137 
138  auto t_normal = getFTensor1Normal();
139  const double ls = sqrt(t_normal(i) * t_normal(i));
140  t_normal(i) /= ls;
141 
142  auto t_w = getFTensor0IntegrationWeight();
143  auto t_disp = getFTensor1FromMat<3>(*(commonDataPtr->mDispPtr));
144  auto t_traction =
145  getFTensor1FromMat<3>(*(commonDataPtr->contactTractionPtr));
146  auto t_coords = getFTensor1CoordsAtGaussPts();
147 
148  size_t nb_base_functions = data.getN().size2() / 3;
149  auto t_base = data.getFTensor1N<3>();
150 
151  auto t_contact_normal =
152  getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
153  auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
154 
155  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
156 
157  auto t_nf = getFTensor1FromArray<3, 3>(locF);
158 
159  const double alpha = t_w * getMeasure();
160  if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
161  VectorDouble n = getNormalsAtGaussPts(gg);
162  auto t_n = getFTensor1FromPtr<3>(&*n.data().begin());
163  t_normal(i) = t_n(i) / sqrt(t_n(j) * t_n(j));
164  }
165  Tensor2<double, 3, 3> t_P;
166  t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
167 
168  Tensor2<double, 3, 3> t_Q;
169  t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
170 
171  const double gap_0 = gap0(t_gap, t_disp, t_contact_normal);
172  Tensor1<double, 3> t_rhs_constrains;
173  t_rhs_constrains(i) =
174  t_contact_normal(i) *
175  constraint(gap_0, t_gap,
176  normal_traction(t_traction, t_contact_normal));
177 
178  Tensor1<double, 3> t_rhs_tangent_disp, t_rhs_tangent_traction;
179  t_rhs_tangent_disp(i) = t_Q(i, j) * t_disp(j);
180  t_rhs_tangent_traction(i) = (*cache).cn_cont * t_Q(i, j) * t_traction(j);
181 
182  size_t bb = 0;
183  for (; bb != nb_dofs / 3; ++bb) {
184  const double beta = alpha * (t_base(i) * t_normal(i));
185 
186  t_nf(i) -= beta * t_rhs_constrains(i);
187  t_nf(i) -= beta * t_rhs_tangent_disp(i);
188  t_nf(i) += beta * t_rhs_tangent_traction(i);
189 
190  ++t_nf;
191  ++t_base;
192  }
193  for (; bb < nb_base_functions; ++bb)
194  ++t_base;
195 
196  ++t_contact_normal;
197  ++t_gap;
198  ++t_disp;
199  ++t_traction;
200  ++t_coords;
201  ++t_w;
202  }
203 
204  }
205 
207 }
208 
210  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr,
211  moab::Interface *moab_inst)
212  : BoundaryEleOp(field_name, field_name, BoundaryEleOp::OPROW),
213  commonDataPtr(common_data_ptr), moab(moab_inst) {
214  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
215  doEntities[MBTRI] = doEntities[MBTRI] = doEntities[MBQUAD] = true;
216 }
217 
219  EntData &data) {
221 
222  const size_t nb_gauss_pts = getGaussPts().size2();
223  const size_t nb_dofs = data.getIndices().size();
224 
225  bool empty = std::all_of(data.getIndices().begin(), data.getIndices().end(),
226  [](int i) { return i == -1; });
227  if (empty)
229 
230  EntityHandle ent = getFEEntityHandle();
231  auto get_tag = [&](const std::string name, size_t size) {
232  Tag th;
233  double def[3] = {0, 0, 0};
234  CHKERR moab->tag_get_handle(name.c_str(), size, MB_TYPE_DOUBLE, th,
235  MB_TAG_CREAT | MB_TAG_SPARSE, def);
236  return th;
237  };
238 
239  auto th_gap = get_tag("GAP", 1);
240  auto th_cons = get_tag("CONSTRAINT", 1);
241  auto th_roller = get_tag("ROLLER_ID", 1);
242  auto th_traction = get_tag("TRACTION", 3);
243  auto th_normal = get_tag("NORMAL", 3);
244  auto th_normal_tri = get_tag("NORMAL_TRI", 3);
245  auto th_disp = get_tag("DISPLACEMENT", 3);
246 
247  auto t_contact_normal =
248  getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
249  auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
250  auto t_traction = getFTensor1FromMat<3>(*(commonDataPtr->contactTractionPtr));
251  auto t_disp = getFTensor1FromMat<3>(*(commonDataPtr->mDispPtr));
252  auto t_coords = getFTensor1CoordsAtGaussPts();
253 
254  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
255 
256  double coords[] = {0, 0, 0};
257  EntityHandle new_vertex;
258  for (int dd = 0; dd != 3; dd++) {
259  coords[dd] = getCoordsAtGaussPts()(gg, dd);
260  }
261  CHKERR moab->create_vertex(&coords[0], new_vertex);
262 
263  Tensor1<double, 3> trac(t_traction(0), t_traction(1), t_traction(2));
264  Tensor1<double, 3> disp(t_disp(0), t_disp(1), t_disp(2));
265  Tensor1<double, 3> norm(t_contact_normal(0), t_contact_normal(1),
266  t_contact_normal(2));
267  trac(i) /= (*cache).young_modulus_inv;
268  double fin_gap = t_gap;
269  const double gap_0 = gap0(t_gap, t_disp, t_contact_normal);
270  double constra =
271  constraint(gap_0, t_gap, normal_traction(t_traction, t_contact_normal));
272  CHKERR moab->tag_set_data(th_gap, &new_vertex, 1, &fin_gap);
273  CHKERR moab->tag_set_data(th_cons, &new_vertex, 1, &constra);
274 
275  CHKERR moab->tag_set_data(th_traction, &new_vertex, 1, &trac(0));
276  CHKERR moab->tag_set_data(th_normal, &new_vertex, 1, &norm(0));
277  CHKERR moab->tag_set_data(th_disp, &new_vertex, 1, &disp(0));
278  auto el_normal = getNormal();
279  if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
280  VectorDouble n = getNormalsAtGaussPts(gg);
281  el_normal = n;
282  }
283  CHKERR moab->tag_set_data(th_normal_tri, &new_vertex, 1, &el_normal(0));
284 
285  double d_tag_data = (*cache).closest_roller_id;
286  CHKERR moab->tag_set_data(th_roller, &new_vertex, 1, &d_tag_data);
287 
288  ++t_contact_normal;
289  ++t_gap;
290  ++t_disp;
291  ++t_coords;
292  ++t_traction;
293  }
294 
296 }
297 
299  const std::string row_field_name, const std::string col_field_name,
300  boost::shared_ptr<CommonData> common_data_ptr)
301  : BoundaryEleOpAssembly(row_field_name, col_field_name, OPROWCOL),
302  commonDataPtr(common_data_ptr) {
303  sYmm = false;
304 }
306  EntData &col_data) {
308 
309  const size_t nb_gauss_pts = getGaussPts().size2();
310  const size_t row_nb_dofs = row_data.getIndices().size();
311  const size_t col_nb_dofs = col_data.getIndices().size();
312 
313  if (row_nb_dofs && col_nb_dofs) {
314 
315 
316  auto t_normal = getFTensor1Normal();
317  const double ls = sqrt(t_normal(i) * t_normal(i));
318  t_normal(i) /= ls;
319 
320  auto t_disp = getFTensor1FromMat<3>(*(commonDataPtr->mDispPtr));
321  auto t_traction =
322  getFTensor1FromMat<3>(*(commonDataPtr->contactTractionPtr));
323  auto t_coords = getFTensor1CoordsAtGaussPts();
324 
325  auto t_w = getFTensor0IntegrationWeight();
326  auto t_row_base = row_data.getFTensor1N<3>();
327  size_t nb_face_functions = row_data.getN().size2() / 3;
328 
329  auto t_contact_normal =
330  getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
331  auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
332  auto t_gap_diff =
333  getFTensor1FromMat<3>(*(commonDataPtr->contactGapDiffPtr));
334  auto t_diff_normal =
335  getFTensor2FromMat<3, 3>(*(commonDataPtr->contactNormalDiffPtr));
337  auto &cn = (*cache).cn_cont;
338 
339  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
340 
341  const double alpha = t_w * getMeasure();
342 
343  if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
344  VectorDouble n = getNormalsAtGaussPts(gg);
345  auto t_n = getFTensor1FromPtr<3>(&*n.data().begin());
346  t_normal(i) = t_n(i) / sqrt(t_n(j) * t_n(j));
347  }
348 
349  Tensor3<double, 3, 3, 3> t_diff_contact_tangent_tensor;
350  t_diff_contact_tangent_tensor(i, j, k) =
351  -t_diff_normal(i, k) * t_contact_normal(j) -
352  t_contact_normal(i) * t_diff_normal(j, k);
353 
354  Tensor2<double, 3, 3> t_P;
355  t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
356 
357  Tensor2<double, 3, 3> t_Q;
358  t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
359  const double gap_0 = gap0(t_gap, t_disp, t_contact_normal);
360 
361  t1(i) = t_contact_normal(i) - 0.5 * t_gap_diff(i) +
362  0.5 * cn * t_traction(j) * t_diff_normal(i, j) +
363  t_disp(j) * t_diff_normal(i, j) +
364  0.5 * sign(t_gap + cn * t_traction(k) * t_contact_normal(k)) *
365  (t_gap_diff(i) + cn * t_traction(j) * t_diff_normal(i, j));
366 
367  auto constrain = constraint(
368  gap_0, t_gap, normal_traction(t_traction, t_contact_normal));
369  size_t rr = 0;
370  for (; rr != row_nb_dofs / 3; ++rr) {
371 
372  auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
373  const double row_base = t_row_base(i) * t_normal(i);
374 
375  auto t_col_base = col_data.getFTensor0N(gg, 0);
376 
377  for (size_t cc = 0; cc != col_nb_dofs / 3; ++cc) {
378  const double beta = alpha * row_base * t_col_base;
379 
380  // linearisation of t_rhs_constrains(i)
381  t_mat(i, j) -= beta * t_contact_normal(i) * t1(j);
382  t_mat(i, j) -= (beta * constrain) * t_diff_normal(i, j);
383 
384  // // linearisation of t_rhs_tangent_disp(i)
385  t_mat(i, j) -= beta * t_Q(i, j);
386  t_mat(i, k) -=
387  beta * t_diff_contact_tangent_tensor(i, j, k) * t_disp(j);
388 
389  // // linearisation of t_rhs_tangent_traction(i);
390  t_mat(i, k) += beta * t_diff_contact_tangent_tensor(i, j, k) *
391  (*cache).cn_cont * t_traction(j);
392 
393  ++t_col_base;
394  ++t_mat;
395  }
396  ++t_row_base;
397  }
398  for (; rr < nb_face_functions; ++rr)
399  ++t_row_base;
400 
401  ++t_contact_normal;
402  ++t_gap;
403  ++t_gap_diff;
404  ++t_diff_normal;
405  ++t_disp;
406  ++t_traction;
407  ++t_coords;
408  ++t_w;
409  }
410 
411  }
412 
414 }
415 
417  const std::string row_field_name, const std::string col_field_name,
418  boost::shared_ptr<CommonData> common_data_ptr)
419  : BoundaryEleOpAssembly(row_field_name, col_field_name, OPROWCOL),
420  commonDataPtr(common_data_ptr) {
421  sYmm = false;
422 }
424  EntData &row_data, EntData &col_data) {
426 
427  const size_t nb_gauss_pts = getGaussPts().size2();
428  const size_t row_nb_dofs = row_data.getIndices().size();
429  const size_t col_nb_dofs = col_data.getIndices().size();
430 
431  if (row_nb_dofs && col_nb_dofs) {
432 
433  auto t_normal = getFTensor1Normal();
434  const double ls = sqrt(t_normal(i) * t_normal(i));
435  t_normal(i) /= ls;
436 
437  auto t_disp = getFTensor1FromMat<3>(*(commonDataPtr->mDispPtr));
438  auto t_traction =
439  getFTensor1FromMat<3>(*(commonDataPtr->contactTractionPtr));
440  auto t_coords = getFTensor1CoordsAtGaussPts();
441 
442  auto t_w = getFTensor0IntegrationWeight();
443  auto t_row_base = row_data.getFTensor1N<3>();
444  size_t nb_face_functions = row_data.getN().size2() / 3;
445 
446  auto t_contact_normal =
447  getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
448  auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
449 
450  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
451 
452  const double alpha = t_w * getMeasure();
453 
454  if (getNormalsAtGaussPts().size1() == nb_gauss_pts) {
455  VectorDouble n = getNormalsAtGaussPts(gg);
456  auto t_n = getFTensor1FromPtr<3>(&*n.data().begin());
457  t_normal(i) = t_n(i) / sqrt(t_n(j) * t_n(j));
458  }
459 
460  // auto t_contact_normal = normal(t_coords, t_disp);
461  Tensor2<double, 3, 3> t_P;
462  t_P(i, j) = t_contact_normal(i) * t_contact_normal(j);
463 
464  Tensor2<double, 3, 3> t_Q;
465  t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
466 
467  const double diff_traction = diff_constraints_dtraction(
468  t_gap, normal_traction(t_traction, t_contact_normal));
469 
470  size_t rr = 0;
471  for (; rr != row_nb_dofs / 3; ++rr) {
472  auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
473  const double row_base = t_row_base(i) * t_normal(i);
474 
475  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
476 
477  for (size_t cc = 0; cc != col_nb_dofs / 3; ++cc) {
478  const double col_base = t_col_base(i) * t_normal(i);
479  const double beta = alpha * row_base * col_base;
480 
481  t_mat(i, j) += (beta * diff_traction) * t_P(i, j);
482  t_mat(i, j) += beta * (*cache).cn_cont * t_Q(i, j);
483 
484  ++t_col_base;
485  ++t_mat;
486  }
487 
488  ++t_row_base;
489  }
490  for (; rr < nb_face_functions; ++rr)
491  ++t_row_base;
492 
493  ++t_contact_normal;
494  ++t_gap;
495  ++t_disp;
496  ++t_traction;
497  ++t_coords;
498  ++t_w;
499  }
500 
501  }
502 
504 }
505 
507  MoFEM::Interface &m_field, const std::string field_name,
508  boost::shared_ptr<CommonData> common_data_ptr, bool update)
509  : mField(m_field),
510  BoundaryEleOp(field_name, field_name, DomainEleOp::OPROW),
511  commonDataPtr(common_data_ptr), uPdate(update) {
512  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
513  doEntities[MBTRI] = doEntities[MBQUAD] = true;
514 }
515 
517  EntData &data) {
519 
520  const size_t nb_gauss_pts = getGaussPts().size2();
521  const size_t nb_dofs = data.getIndices().size();
522 
523  EntityHandle ent = getFEEntityHandle();
524  auto t_coords = getFTensor1CoordsAtGaussPts();
525  auto t_disp = getFTensor1FromMat<3>(*(commonDataPtr->mDispPtr));
526 
527  commonDataPtr->contactNormalPtr->resize(3, nb_gauss_pts, false);
528  commonDataPtr->contactNormalPtr->clear();
529  commonDataPtr->contactNormalDiffPtr->resize(9, nb_gauss_pts, false);
530  commonDataPtr->contactNormalDiffPtr->clear();
531  commonDataPtr->contactGapDiffPtr->resize(3, nb_gauss_pts, false);
532  commonDataPtr->contactGapDiffPtr->clear();
533  commonDataPtr->contactGapPtr->resize(nb_gauss_pts, false);
534  commonDataPtr->contactGapPtr->clear();
535 
536  auto t_contact_normal =
537  getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
538  auto t_diff_normal =
539  getFTensor2FromMat<3, 3>(*(commonDataPtr->contactNormalDiffPtr));
540  auto t_gap_diff = getFTensor1FromMat<3>(*(commonDataPtr->contactGapDiffPtr));
541  auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
542 
543  if (uPdate) {
544  // FIXME: takes first gauss pt
545  int roll_id = (*cache).getClosestRollerID(t_coords, t_disp);
546  CHKERR mField.get_moab().tag_set_data(commonDataPtr->rollerTag, &ent, 1,
547  &roll_id);
548  }
549  int tag_data;
550  CHKERR mField.get_moab().tag_get_data(commonDataPtr->rollerTag, &ent, 1,
551  &tag_data);
552  (*cache).closest_roller = &(*cache).rollerDataVec.at(tag_data);
553  (*cache).closest_roller_id = tag_data;
554 
555  const bool is_lhs = (int)getTSCtx() == 3;
556  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
557 
558  auto t_contact = (*cache).closest_roller->getNormal(t_coords, t_disp);
559  t_contact_normal(i) = t_contact(i);
560  t_gap = (*cache).closest_roller->getGap(t_coords);
561 
562  if (is_lhs) {
563  auto t_dn_du = ((*cache).closest_roller)
564  ->getDiffNormal(t_coords, t_disp, t_contact_normal);
565  t_diff_normal(i, j) = t_dn_du(i, j);
566  auto t_d_gap =
567  (*cache).closest_roller->getdGap(t_coords, t_contact_normal);
568  t_gap_diff(i) = t_d_gap(i);
569  }
570 
571  ++t_contact_normal;
572  ++t_gap;
573  ++t_gap_diff;
574  ++t_coords;
575  ++t_disp;
576  ++t_diff_normal;
577  }
578 
579  // } //dofs
580 
582 }
583 
585  OpPostProcContact(const std::string field_name,
586  moab::Interface &post_proc_mesh,
587  std::vector<EntityHandle> &map_gauss_pts,
588  boost::shared_ptr<CommonData> common_data_ptr);
589  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
590 
591 private:
593  std::vector<EntityHandle> &mapGaussPts;
594  boost::shared_ptr<CommonData> commonDataPtr;
595 };
596 //! [Operators definitions]
597 
599  const std::string field_name, moab::Interface &post_proc_mesh,
600  std::vector<EntityHandle> &map_gauss_pts,
601  boost::shared_ptr<CommonData> common_data_ptr)
602  : DomainEleOp(field_name, field_name, DomainEleOp::OPROW),
603  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
604  commonDataPtr(common_data_ptr) {
605  // Operator is only executed for vertices
606  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
607  doEntities[MBVERTEX] = true;
608 }
609 
610 //! [Postprocessing]
612  EntData &data) {
614 
615  auto get_tag_mat = [&](const std::string name) {
616  std::array<double, 9> def;
617  std::fill(def.begin(), def.end(), 0);
618  Tag th;
619  CHKERR postProcMesh.tag_get_handle(name.c_str(), 9, MB_TYPE_DOUBLE, th,
620  MB_TAG_CREAT | MB_TAG_SPARSE,
621  def.data());
622  return th;
623  };
624 
625  auto get_tag_vec = [&](const std::string name) {
626  std::array<double, 3> def;
627  std::fill(def.begin(), def.end(), 0);
628  Tag th;
629  CHKERR postProcMesh.tag_get_handle(name.c_str(), 3, MB_TYPE_DOUBLE, th,
630  MB_TAG_CREAT | MB_TAG_SPARSE,
631  def.data());
632  return th;
633  };
634 
635  MatrixDouble3by3 mat(3, 3);
636 
637  auto set_matrix = [&](auto &t) -> MatrixDouble3by3 & {
638  mat.clear();
639  for (size_t r = 0; r != 3; ++r)
640  for (size_t c = 0; c != 3; ++c)
641  mat(r, c) = t(r, c);
642  return mat;
643  };
644 
645  auto set_vector = [&](auto &t) -> MatrixDouble3by3 & {
646  mat.clear();
647  for (size_t r = 0; r != 3; ++r)
648  mat(0, r) = t(r);
649  return mat;
650  };
651 
652  auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
653  return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
654  &*mat.data().begin());
655  };
656 
657  auto th_stress = get_tag_mat("SIGMA");
658  auto th_div = get_tag_vec("DIV_SIGMA");
659 
660  size_t nb_gauss_pts = getGaussPts().size2();
661  auto t_stress = getFTensor2FromMat<3, 3>(*(commonDataPtr->contactStressPtr));
662  auto t_div =
663  getFTensor1FromMat<3>(*(commonDataPtr->contactStressDivergencePtr));
664 
665  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
666 
667  t_stress(i, j) /= (*cache).young_modulus_inv;
668  t_div(i) /= (*cache).young_modulus_inv;
669 
670  // FIXME: add option for more postprocessing
671  // CHKERR set_tag(th_stress, gg, set_matrix(t_stress));
672  CHKERR set_tag(th_div, gg, set_vector(t_div));
673  ++t_stress;
674  ++t_div;
675  }
676 
677  // Tag th_time;
678  // const EntityHandle root_meshset = postProcMesh.get_root_set();
679  // double def_t_val = 0;
680  // rval = postProcMesh.tag_get_handle("_TS_TIME_", 1, MB_TYPE_DOUBLE, th_time,
681  // MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH,
682  // &def_t_val);
683  // if (rval == MB_ALREADY_ALLOCATED) {
684  // // MOAB_THROW(rval);
685  // } else {
686  // def_t_val = getFEMethod()->ts_t;
687  // CHKERR postProcMesh.tag_set_data(th_time, &root_meshset, 1, &def_t_val);
688  // MOAB_THROW(rval);
689  // }
690 
692 }
693 //! [Postprocessing]
694 
696 
697  OpGetGaussPtsContactState(const std::string field_name,
698  boost::shared_ptr<CommonData> common_data_ptr)
699  : BoundaryEleOp(field_name, field_name, BoundaryEleOp::OPROW),
700  commonDataPtr(common_data_ptr) {
701  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
702  doEntities[MBTRI] = doEntities[MBQUAD] = true;
703  }
704 
705  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
707  const size_t nb_dofs = data.getIndices().size();
708  if (nb_dofs) {
709  bool empty =
710  std::all_of(data.getIndices().begin(), data.getIndices().end(),
711  [](int i) { return i == -1; });
712  if (empty)
714  size_t nb_gauss_pts = getGaussPts().size2();
715  auto t_traction =
716  getFTensor1FromMat<3>(*(commonDataPtr->contactTractionPtr));
717 
718  auto &data_vec = commonDataPtr->stateArrayCont;
719  auto t_contact_normal =
720  getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
721  auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
722  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
723 
724  const double trac = normal_traction(t_traction, t_contact_normal);
725 
726  if (t_gap <= (*cache).cn_cont * trac)
727  data_vec[0] += 1;
728  data_vec[1] += 1;
729 
730  ++t_contact_normal;
731  ++t_gap;
732  ++t_traction;
733  }
734  }
736  }
737 
738 private:
739  boost::shared_ptr<CommonData> commonDataPtr;
740 };
741 
743 
744  OpGetContactArea(const std::string field_name,
745  boost::shared_ptr<CommonData> common_data_ptr)
746  : BoundaryEleOp(field_name, field_name, BoundaryEleOp::OPROW),
747  commonDataPtr(common_data_ptr) {
748  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
749  doEntities[MBTRI] = doEntities[MBQUAD] = true;
750  }
751 
752  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
754  const size_t nb_dofs = data.getIndices().size();
755  if (nb_dofs) {
756  auto ent = getFEEntityHandle();
757  bool empty =
758  std::all_of(data.getIndices().begin(), data.getIndices().end(),
759  [](int i) { return i == -1; });
760  if (empty)
762 
763  size_t nb_gauss_pts = getGaussPts().size2();
764  auto t_w = getFTensor0IntegrationWeight();
765  auto t_traction =
766  getFTensor1FromMat<3>(*(commonDataPtr->contactTractionPtr));
767 
768  auto &data_vec = commonDataPtr->stateArrayArea;
769  auto t_contact_normal =
770  getFTensor1FromMat<3>(*(commonDataPtr->contactNormalPtr));
771  auto t_gap = getFTensor0FromVec(*(commonDataPtr->contactGapPtr));
772 
773  int id;
774  CHKERR getPtrFE()->mField.get_moab().tag_get_data(
775  commonDataPtr->rollerTag, &ent, 1, &id);
776  auto &l_reactions = commonDataPtr->bodyReactions;
777 
778  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
779  const double alpha = t_w * getMeasure();
780  const double trac = normal_traction(t_traction, t_contact_normal);
781  const auto sign_w = sign(w(t_gap, trac));
782 
783  for (int dd = 0; dd != 3; ++dd)
784  l_reactions[3 * id + dd] +=
785  t_traction(dd) * alpha / (*cache).young_modulus_inv;
786 
787  if (sign_w < 0)
788  data_vec[0] += alpha;
789  data_vec[1] += alpha;
790 
791  ++t_contact_normal;
792  ++t_gap;
793  ++t_traction;
794  ++t_w;
795  }
796  }
798  }
799 
800 private:
801  boost::shared_ptr<CommonData> commonDataPtr;
802 };
803 
804 }; // namespace OpContactTools
static Index< 'n', 3 > n
constexpr double alpha
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:116
FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > getFTensor1FromPtr< 3 >(double *ptr)
Definition: Templates.hpp:732
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
double constraint(double g0, double g, double &&t)
double w(const double g, const double t)
double gap0(double g, FTensor::Tensor1< T, 3 > &t_disp, FTensor::Tensor1< T, 3 > &t_normal)
double diff_constraints_dgap(double g, double &&t)
double normal_traction(Tensor1< T1, 3 > &t_traction, Tensor1< T2, 3 > &t_normal)
double diff_constraints_dtraction(double g, double &&t)
double sign(double x)
double constrain(long double dot_tau, long double f, long double sigma_y)
Definition: PlasticOps.hpp:491
const double r
rate factor
double cn
Definition: plastic.cpp:111
constexpr double t
plate stiffness
Definition: plate.cpp:76
constexpr double g
virtual moab::Interface & get_moab()=0
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
boost::shared_ptr< MatrixDouble > contactTractionPtr
boost::shared_ptr< MatrixDouble > contactGapDiffPtr
boost::shared_ptr< MatrixDouble > contactStressDivergencePtr
boost::shared_ptr< MatrixDouble > contactNormalDiffPtr
boost::shared_ptr< VectorDouble > contactGapPtr
boost::shared_ptr< MatrixDouble > contactStressPtr
boost::shared_ptr< MatrixDouble > contactNormalPtr
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
OpConstrainBoundaryLhs_dTraction(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
OpConstrainBoundaryLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode iNtegrate(EntData &data)
boost::shared_ptr< CommonData > commonDataPtr
OpConstrainBoundaryRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpGetClosestRigidBodyData(MoFEM::Interface &m_field, const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, bool update)
boost::shared_ptr< CommonData > commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
OpGetContactArea(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpGetGaussPtsContactState(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
moab::Interface & postProcMesh
OpPostProcContact(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
[Operators definitions]
std::vector< EntityHandle > & mapGaussPts
boost::shared_ptr< CommonData > commonDataPtr
[Common data]
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
moab::Interface * moab
boost::shared_ptr< CommonData > commonDataPtr
OpPostProcVertex(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, moab::Interface *moab_inst)