v0.9.2
ContactOps.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 //! [Common data]
19  boost::shared_ptr<MatrixDouble> contactStressPtr;
20  boost::shared_ptr<MatrixDouble> contactStressDivergencePtr;
21  boost::shared_ptr<MatrixDouble> contactTractionPtr;
22  boost::shared_ptr<MatrixDouble> contactDispPtr;
23 };
24 //! [Common data]
25 
30 
32  OpInternalDomainContactRhs(const std::string field_name,
33  boost::shared_ptr<CommonData> common_data_ptr);
34  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
35 
36 private:
37  boost::shared_ptr<CommonData> commonDataPtr;
38 };
39 
41  OpConstrainBoundaryRhs(const std::string field_name,
42  boost::shared_ptr<CommonData> common_data_ptr);
43  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
44 
45 private:
46  boost::shared_ptr<CommonData> commonDataPtr;
47 };
48 
50  OpConstrainBoundaryLhs_dU(const std::string row_field_name,
51  const std::string col_field_name,
52  boost::shared_ptr<CommonData> common_data_ptr);
53  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
54  EntityType col_type, EntData &row_data,
55  EntData &col_data);
56 
57 private:
58  boost::shared_ptr<CommonData> commonDataPtr;
60 };
61 
64  const std::string row_field_name, const std::string col_field_name,
65  boost::shared_ptr<CommonData> common_data_ptr);
66  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
67  EntityType col_type, EntData &row_data,
68  EntData &col_data);
69 
70 private:
71  boost::shared_ptr<CommonData> commonDataPtr;
73 };
74 
75 struct OpSpringRhs : public BoundaryEleOp {
76  OpSpringRhs(const std::string field_name,
77  boost::shared_ptr<CommonData> common_data_ptr);
78  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
79 
80 private:
81  boost::shared_ptr<CommonData> commonDataPtr;
82 };
83 
84 struct OpSpringLhs : public BoundaryEleOp {
85  OpSpringLhs(const std::string row_field_name,
86  const std::string col_field_name);
87  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
88  EntityType col_type, EntData &row_data,
89  EntData &col_data);
90 
91 private:
93 };
94 
96  OpConstrainBoundaryTraction(const std::string field_name,
97  boost::shared_ptr<CommonData> common_data_ptr);
98  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
99 
100 private:
101  boost::shared_ptr<CommonData> commonDataPtr;
102 };
103 
105  OpConstrainDomainRhs(const std::string field_name,
106  boost::shared_ptr<CommonData> common_data_ptr);
107  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
108 
109 private:
110  boost::shared_ptr<CommonData> commonDataPtr;
111 };
112 
114  OpConstrainDomainLhs_dU(const std::string row_field_name,
115  const std::string col_field_name,
116  boost::shared_ptr<CommonData> common_data_ptr);
117  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
118  EntityType col_type, EntData &row_data,
119  EntData &col_data);
120 
121 private:
122  boost::shared_ptr<CommonData> commonDataPtr;
125 };
126 
127 template <typename T1, typename T2>
129  FTensor::Tensor1<T2, 2> &t_disp) {
130  return FTensor::Tensor1<double, 2>{0., 1.};
131 }
132 
133 template <typename T>
134 inline double gap0(FTensor::Tensor1<T, 3> &t_coords,
135  FTensor::Tensor1<double, 2> &t_normal) {
136  return (-1 - t_coords(1)) * t_normal(1);
137 }
138 
139 template <typename T>
140 inline double gap(FTensor::Tensor1<T, 2> &t_disp,
141  FTensor::Tensor1<double, 2> &t_normal) {
142  return t_disp(i) * t_normal(i);
143 }
144 
145 template <typename T>
146 inline double normal_traction(FTensor::Tensor1<T, 2> &t_traction,
147  FTensor::Tensor1<double, 2> &t_normal) {
148  return t_traction(i) * t_normal(i);
149 }
150 
151 inline double sign(double x) {
152  if (x == 0)
153  return 0;
154  else if (x > 0)
155  return 1;
156  else
157  return -1;
158 };
159 
160 inline double w(const double g, const double t) {
161  return g - cn * t;
162 }
163 
164 inline double constrian(double &&g0, double &&g, double &&t) {
165  return (w(g - g0, t) + std::abs(w(g - g0, t))) / 2 + g0;
166 };
167 
168 inline double diff_constrains_dtraction(double &&g0, double &&g, double &&t) {
169  return -cn * (1 + sign(w(g - g0, t))) / 2;
170 }
171 
172 inline double diff_constrains_dgap(double &&g0, double &&g, double &&t) {
173  return (1 + sign(w(g - g0, t))) / 2;
174 }
175 
177  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
178  : DomainEleOp(field_name, DomainEleOp::OPROW),
179  commonDataPtr(common_data_ptr) {}
180 
182  EntData &data) {
184  const size_t nb_gauss_pts = getGaussPts().size2();
185  const size_t nb_dofs = data.getIndices().size();
186 
187  if (nb_dofs) {
188 
189  std::array<double, MAX_DOFS_ON_ENTITY> nf;
190  std::fill(&nf[0], &nf[nb_dofs], 0);
191 
192  const size_t nb_base_functions = data.getN().size2();
193  auto t_w = getFTensor0IntegrationWeight();
194  auto t_base = data.getFTensor0N();
195  auto t_diff_base = data.getFTensor1DiffN<2>();
196  auto t_stress =
197  getFTensor2FromMat<2, 2>(*(commonDataPtr->contactStressPtr));
198  auto t_div_stress =
199  getFTensor1FromMat<2>(*(commonDataPtr->contactStressDivergencePtr));
200 
201  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
202 
203  const double alpha = getMeasure() * t_w;
204  FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2> t_nf{&nf[0], &nf[1]};
205 
206  size_t bb = 0;
207  for (; bb != nb_dofs / 2; ++bb) {
208 
209  t_nf(i) += alpha * t_base * t_div_stress(i);
210  t_nf(i) += alpha * t_diff_base(j) * t_stress(i, j);
211 
212  ++t_nf;
213  ++t_base;
214  ++t_diff_base;
215  }
216  for (; bb < nb_base_functions; ++bb) {
217  ++t_base;
218  ++t_diff_base;
219  }
220 
221  ++t_div_stress;
222  ++t_stress;
223  ++t_w;
224  }
225 
226  CHKERR VecSetValues(getSNESf(), data, nf.data(), ADD_VALUES);
227  }
228 
230 }
231 
233  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
234  : BoundaryEleOp(field_name, DomainEleOp::OPROW),
235  commonDataPtr(common_data_ptr) {}
236 
238  EntData &data) {
240 
241  const size_t nb_gauss_pts = getGaussPts().size2();
242  const size_t nb_dofs = data.getIndices().size();
243 
244  if (nb_dofs) {
245 
246  std::array<double, MAX_DOFS_ON_ENTITY> nf;
247  std::fill(&nf[0], &nf[nb_dofs], 0);
248 
249  FTensor::Tensor1<double, 2> t_direction{getDirection()[0],
250  getDirection()[1]};
251  FTensor::Tensor1<double, 2> t_normal{-t_direction(1), t_direction(0)};
252  const double l = sqrt(t_normal(i) * t_normal(i));
253  t_normal(i) /= l;
254  t_direction(i) /= l;
255 
256  auto t_w = getFTensor0IntegrationWeight();
257  auto t_disp = getFTensor1FromMat<2>(*(commonDataPtr->contactDispPtr));
258  auto t_traction =
259  getFTensor1FromMat<2>(*(commonDataPtr->contactTractionPtr));
260  auto t_coords = getFTensor1CoordsAtGaussPts();
261 
262  size_t nb_base_functions = data.getN().size2() / 3;
263  auto t_base = data.getFTensor1N<3>();
264  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
265 
266  FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2> t_nf{&nf[0], &nf[1]};
267 
268  const double alpha = t_w * l;
269 
270  auto t_contact_normal = normal(t_coords, t_disp);
271  FTensor::Tensor2<double, 2, 2> t_contact_normal_tensor;
272  t_contact_normal_tensor(i, j) = t_contact_normal(i) * t_contact_normal(j);
273 
274  FTensor::Tensor2<double, 2, 2> t_contact_tangent_tensor;
275  t_contact_tangent_tensor(i, j) = t_contact_normal_tensor(i, j);
276  t_contact_tangent_tensor(0, 0) -= 1;
277  t_contact_tangent_tensor(1, 1) -= 1;
278 
279  FTensor::Tensor1<double, 2> t_rhs_constrains;
280 
281  t_rhs_constrains(i) =
282  t_contact_normal(i) *
283  constrian(gap0(t_coords, t_contact_normal),
284  gap(t_disp, t_contact_normal),
285  normal_traction(t_traction, t_contact_normal));
286 
287  FTensor::Tensor1<double, 2> t_rhs_tangent_disp, t_rhs_tangent_traction;
288  t_rhs_tangent_disp(i) = t_contact_tangent_tensor(i, j) * t_disp(j);
289  t_rhs_tangent_traction(i) =
290  cn * t_contact_tangent_tensor(i, j) * t_traction(j);
291 
292  size_t bb = 0;
293  for (; bb != nb_dofs / 2; ++bb) {
294  const double beta = alpha * (t_base(i) * t_normal(i));
295 
296  t_nf(i) += beta * t_rhs_constrains(i);
297  t_nf(i) += beta * t_rhs_tangent_disp(i);
298  t_nf(i) += beta * t_rhs_tangent_traction(i);
299 
300  ++t_nf;
301  ++t_base;
302  }
303  for (; bb < nb_base_functions; ++bb)
304  ++t_base;
305 
306  ++t_disp;
307  ++t_traction;
308  ++t_coords;
309  ++t_w;
310  }
311 
312  CHKERR VecSetValues(getSNESf(), data, nf.data(), ADD_VALUES);
313  }
314 
316 }
317 
319  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
320  : BoundaryEleOp(field_name, DomainEleOp::OPROW),
321  commonDataPtr(common_data_ptr) {
322  std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
323  doEntities[MBEDGE] = true;
324 }
325 
327  EntData &data) {
328 
330 
331  const size_t nb_gauss_pts = getGaussPts().size2();
332 
333  if (side == 0 && type == MBEDGE) {
334  commonDataPtr->contactTractionPtr->resize(2, nb_gauss_pts);
335  commonDataPtr->contactTractionPtr->clear();
336  }
337 
338  const size_t nb_dofs = data.getIndices().size();
339  if (nb_dofs) {
340 
341  auto t_direction = getFTensor1Direction();
342  FTensor::Tensor1<double, 2> t_normal{-t_direction(1), t_direction(0)};
343  const double l = sqrt(t_normal(i) * t_normal(i));
344  t_normal(i) /= l;
345 
346  auto t_traction =
347  getFTensor1FromMat<2>(*(commonDataPtr->contactTractionPtr));
348 
349  size_t nb_base_functions = data.getN().size2() / 3;
350  auto t_base = data.getFTensor1N<3>();
351 
352  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
353 
354  auto t_field_data = data.getFTensor1FieldData<2>();
355 
356  size_t bb = 0;
357  for (; bb != nb_dofs / 2; ++bb) {
358  t_traction(j) += (t_base(i) * t_normal(i)) * t_field_data(j);
359  ++t_field_data;
360  ++t_base;
361  }
362 
363  for (; bb < nb_base_functions; ++bb)
364  ++t_base;
365 
366  ++t_traction;
367  }
368  }
369 
371 }
372 
374  const std::string row_field_name, const std::string col_field_name,
375  boost::shared_ptr<CommonData> common_data_ptr)
376  : BoundaryEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
377  commonDataPtr(common_data_ptr) {
378  sYmm = false;
379 }
381  EntityType row_type,
382  EntityType col_type,
383  EntData &row_data,
384  EntData &col_data) {
386 
387  const size_t nb_gauss_pts = getGaussPts().size2();
388  const size_t row_nb_dofs = row_data.getIndices().size();
389  const size_t col_nb_dofs = col_data.getIndices().size();
390 
391  if (row_nb_dofs && col_nb_dofs) {
392 
393  locMat.resize(row_nb_dofs, col_nb_dofs, false);
394  locMat.clear();
395 
396  FTensor::Tensor1<double, 2> t_direction{getDirection()[0],
397  getDirection()[1]};
398  FTensor::Tensor1<double, 2> t_normal{-t_direction(1), t_direction(0)};
399  const double l = sqrt(t_normal(i) * t_normal(i));
400  t_normal(i) /= l;
401  t_direction(i) /= l;
402 
403  auto t_disp = getFTensor1FromMat<2>(*(commonDataPtr->contactDispPtr));
404  auto t_traction =
405  getFTensor1FromMat<2>(*(commonDataPtr->contactTractionPtr));
406  auto t_coords = getFTensor1CoordsAtGaussPts();
407 
408  auto t_w = getFTensor0IntegrationWeight();
409  auto t_row_base = row_data.getFTensor1N<3>();
410  size_t nb_face_functions = row_data.getN().size2() / 3;
411 
412  locMat.resize(row_nb_dofs, col_nb_dofs, false);
413  locMat.clear();
414 
415  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
416 
417  const double alpha = t_w * l;
418 
419  auto t_contact_normal = normal(t_coords, t_disp);
420  FTensor::Tensor2<double, 2, 2> t_contact_normal_tensor;
421  t_contact_normal_tensor(i, j) = t_contact_normal(i) * t_contact_normal(j);
422 
423  FTensor::Tensor2<double, 2, 2> t_contact_tangent_tensor;
424  t_contact_tangent_tensor(i, j) = t_contact_normal_tensor(i, j);
425  t_contact_tangent_tensor(0, 0) -= 1;
426  t_contact_tangent_tensor(1, 1) -= 1;
427 
428  auto diff_constrain = diff_constrains_dgap(
429  gap0(t_coords, t_contact_normal), gap(t_disp, t_contact_normal),
430  normal_traction(t_traction, t_contact_normal));
431 
432  size_t rr = 0;
433  for (; rr != row_nb_dofs / 2; ++rr) {
434 
436  &locMat(2 * rr + 0, 0), &locMat(2 * rr + 0, 1),
437  &locMat(2 * rr + 1, 0), &locMat(2 * rr + 1, 1));
438  const double row_base = t_row_base(i) * t_normal(i);
439 
440  auto t_col_base = col_data.getFTensor0N(gg, 0);
441  for (size_t cc = 0; cc != col_nb_dofs / 2; ++cc) {
442  const double beta = alpha * row_base * t_col_base;
443 
444  t_mat(i, j) += (beta * diff_constrain) * t_contact_normal_tensor(i, j);
445  t_mat(i, j) += beta * t_contact_tangent_tensor(i, j);
446 
447  ++t_col_base;
448  ++t_mat;
449  }
450 
451  ++t_row_base;
452  }
453  for (; rr < nb_face_functions; ++rr)
454  ++t_row_base;
455 
456  ++t_disp;
457  ++t_traction;
458  ++t_coords;
459  ++t_w;
460  }
461 
462  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
463  ADD_VALUES);
464  }
465 
467 }
468 
470  const std::string row_field_name, const std::string col_field_name,
471  boost::shared_ptr<CommonData> common_data_ptr)
472  : BoundaryEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
473  commonDataPtr(common_data_ptr) {
474  sYmm = false;
475 }
477  int row_side, int col_side, EntityType row_type, EntityType col_type,
478  EntData &row_data, EntData &col_data) {
480 
481  const size_t nb_gauss_pts = getGaussPts().size2();
482  const size_t row_nb_dofs = row_data.getIndices().size();
483  const size_t col_nb_dofs = col_data.getIndices().size();
484 
485  if (row_nb_dofs && col_nb_dofs) {
486 
487  locMat.resize(row_nb_dofs, col_nb_dofs, false);
488  locMat.clear();
489 
490  FTensor::Tensor1<double, 2> t_direction{getDirection()[0],
491  getDirection()[1]};
492  FTensor::Tensor1<double, 2> t_normal{-t_direction(1), t_direction(0)};
493  const double l = sqrt(t_normal(i) * t_normal(i));
494  t_normal(i) /= l;
495  t_direction(i) /= l;
496 
497  auto t_disp = getFTensor1FromMat<2>(*(commonDataPtr->contactDispPtr));
498  auto t_traction =
499  getFTensor1FromMat<2>(*(commonDataPtr->contactTractionPtr));
500  auto t_coords = getFTensor1CoordsAtGaussPts();
501 
502  auto t_w = getFTensor0IntegrationWeight();
503  auto t_row_base = row_data.getFTensor1N<3>();
504  size_t nb_face_functions = row_data.getN().size2() / 3;
505 
506  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
507 
508  const double alpha = t_w * l;
509 
510  auto t_contact_normal = normal(t_coords, t_disp);
511  FTensor::Tensor2<double, 2, 2> t_contact_normal_tensor;
512  t_contact_normal_tensor(i, j) = t_contact_normal(i) * t_contact_normal(j);
513 
514  FTensor::Tensor2<double, 2, 2> t_contact_tangent_tensor;
515  t_contact_tangent_tensor(i, j) = t_contact_normal_tensor(i, j);
516  t_contact_tangent_tensor(0, 0) -= 1;
517  t_contact_tangent_tensor(1, 1) -= 1;
518 
519  const double diff_traction = diff_constrains_dtraction(
520  gap0(t_coords, t_contact_normal), gap(t_disp, t_contact_normal),
521  normal_traction(t_traction, t_contact_normal));
522 
523  size_t rr = 0;
524  for (; rr != row_nb_dofs / 2; ++rr) {
526  &locMat(2 * rr + 0, 0), &locMat(2 * rr + 0, 1),
527  &locMat(2 * rr + 1, 0), &locMat(2 * rr + 1, 1));
528 
529  const double row_base = t_row_base(i) * t_normal(i);
530 
531  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
532  for (size_t cc = 0; cc != col_nb_dofs / 2; ++cc) {
533  const double col_base = t_col_base(i) * t_normal(i);
534  const double beta = alpha * row_base * col_base;
535 
536  t_mat(i, j) += (beta * diff_traction) * t_contact_normal_tensor(i, j);
537  t_mat(i, j) += beta * cn * t_contact_tangent_tensor(i, j);
538 
539  ++t_col_base;
540  ++t_mat;
541  }
542 
543  ++t_row_base;
544  }
545  for (; rr < nb_face_functions; ++rr)
546  ++t_row_base;
547 
548  ++t_disp;
549  ++t_traction;
550  ++t_coords;
551  ++t_w;
552  }
553 
554  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
555  ADD_VALUES);
556  }
557 
559 }
560 
562  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
563  : DomainEleOp(field_name, DomainEleOp::OPROW),
564  commonDataPtr(common_data_ptr) {}
565 
567  EntData &data) {
569  const size_t nb_gauss_pts = getGaussPts().size2();
570  const size_t nb_dofs = data.getIndices().size();
571 
572  if (nb_dofs) {
573 
574  std::array<double, MAX_DOFS_ON_ENTITY> nf;
575  std::fill(&nf[0], &nf[nb_dofs], 0);
576 
577  const size_t nb_base_functions = data.getN().size2() / 3;
578  auto t_w = getFTensor0IntegrationWeight();
579  auto t_base = data.getFTensor1N<3>();
580  auto t_diff_base = data.getFTensor2DiffN<3, 2>();
581  auto t_stress =
582  getFTensor2FromMat<2, 2>(*(commonDataPtr->contactStressPtr));
583  auto t_disp = getFTensor1FromMat<2>(*(commonDataPtr->contactDispPtr));
584  auto t_grad = getFTensor2FromMat<2, 2>(*(commonDataPtr->mGradPtr));
585 
586  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
587 
588  const double alpha = getMeasure() * t_w;
589  FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2> t_nf{&nf[0], &nf[1]};
590 
591  size_t bb = 0;
592  for (; bb != nb_dofs / 2; ++bb) {
593  const double t_div_base = t_diff_base(0, 0) + t_diff_base(1, 1);
594 
595  t_nf(i) += alpha * (t_base(j) * t_grad(i, j));
596  t_nf(i) += alpha * t_div_base * t_disp(i);
597 
598  ++t_nf;
599  ++t_base;
600  ++t_diff_base;
601  }
602  for (; bb < nb_base_functions; ++bb) {
603  ++t_base;
604  ++t_diff_base;
605  }
606 
607  ++t_stress;
608  ++t_disp;
609  ++t_grad;
610  ++t_w;
611  }
612 
613  CHKERR VecSetValues(getSNESf(), data, nf.data(), ADD_VALUES);
614  }
615 
617 }
618 
620  const std::string row_field_name, const std::string col_field_name,
621  boost::shared_ptr<CommonData> common_data_ptr)
622  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
623  commonDataPtr(common_data_ptr) {
624  sYmm = false;
625 }
626 
628  EntityType row_type,
629  EntityType col_type,
630  EntData &row_data,
631  EntData &col_data) {
633 
634  const size_t nb_gauss_pts = getGaussPts().size2();
635  const size_t row_nb_dofs = row_data.getIndices().size();
636  const size_t col_nb_dofs = col_data.getIndices().size();
637 
638  if (row_nb_dofs && col_nb_dofs) {
639 
640  auto t_w = getFTensor0IntegrationWeight();
641 
642  locMat.resize(row_nb_dofs, col_nb_dofs, false);
643  locMat.clear();
644  transLocMat.resize(col_nb_dofs, row_nb_dofs, false);
645 
646  size_t nb_base_functions = row_data.getN().size2() / 3;
647  auto t_row_base = row_data.getFTensor1N<3>();
648  auto t_row_diff_base = row_data.getFTensor2DiffN<3, 2>();
649 
650  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
651 
652  const double alpha = getMeasure() * t_w;
653 
654  size_t rr = 0;
655  for (; rr != row_nb_dofs / 2; ++rr) {
656 
658  &locMat(2 * rr + 0, 0), &locMat(2 * rr + 1, 1)};
660  &locMat(2 * rr + 0, 0), &locMat(2 * rr + 0, 1),
661  &locMat(2 * rr + 1, 0), &locMat(2 * rr + 1, 1)};
662 
663  const double t_row_div_base =
664  t_row_diff_base(0, 0) + t_row_diff_base(1, 1);
665 
666  auto t_col_base = col_data.getFTensor0N(gg, 0);
667  auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
668 
669  for (size_t cc = 0; cc != col_nb_dofs / 2; ++cc) {
670 
671  t_mat_diag(i) += alpha * t_row_base(j) * t_col_diff_base(j);
672  t_mat_diag(i) += alpha * t_row_div_base * t_col_base;
673 
674  ++t_col_base;
675  ++t_col_diff_base;
676  ++t_mat_diag;
677  ++t_mat;
678  }
679 
680  ++t_row_diff_base;
681  ++t_row_base;
682  }
683  for (; rr < nb_base_functions; ++rr) {
684  ++t_row_diff_base;
685  ++t_row_base;
686  }
687 
688  ++t_w;
689  }
690 
691  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
692  ADD_VALUES);
693  noalias(transLocMat) = trans(locMat);
694  CHKERR MatSetValues(getSNESB(), col_data, row_data,
695  &*transLocMat.data().begin(), ADD_VALUES);
696  }
697 
699 }
700 
701 OpSpringRhs::OpSpringRhs(const std::string field_name,
702  boost::shared_ptr<CommonData> common_data_ptr)
703  : BoundaryEleOp(field_name, DomainEleOp::OPROW),
704  commonDataPtr(common_data_ptr) {}
705 
706 MoFEMErrorCode OpSpringRhs::doWork(int side, EntityType type, EntData &data) {
708 
709  const size_t nb_gauss_pts = getGaussPts().size2();
710  const size_t nb_dofs = data.getIndices().size();
711 
712  if (nb_dofs) {
713 
714  std::array<double, MAX_DOFS_ON_ENTITY> nf;
715  std::fill(&nf[0], &nf[nb_dofs], 0);
716 
717  auto t_w = getFTensor0IntegrationWeight();
718  auto t_disp = getFTensor1FromMat<2>(*(commonDataPtr->contactDispPtr));
719 
720  size_t nb_base_functions = data.getN().size2();
721  auto t_base = data.getFTensor0N();
722  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
723 
724  FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2> t_nf{&nf[0], &nf[1]};
725  const double alpha = t_w * getMeasure();
726 
727  size_t bb = 0;
728  for (; bb != nb_dofs / 2; ++bb) {
729 
730  const double beta = alpha * t_base;
731  t_nf(i) += beta * spring_stiffness * t_disp(i);
732 
733  ++t_nf;
734  ++t_base;
735  }
736  for (; bb < nb_base_functions; ++bb)
737  ++t_base;
738 
739  ++t_disp;
740  ++t_w;
741  }
742 
743  CHKERR VecSetValues(getSNESf(), data, nf.data(), ADD_VALUES);
744  }
745 
747 }
748 
749 OpSpringLhs::OpSpringLhs(const std::string row_field_name,
750  const std::string col_field_name)
751  : BoundaryEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL) {
752  sYmm = false;
753 }
754 
755 MoFEMErrorCode OpSpringLhs::doWork(int row_side, int col_side,
756  EntityType row_type, EntityType col_type,
757  EntData &row_data, EntData &col_data) {
758 
760 
761  const size_t nb_gauss_pts = getGaussPts().size2();
762  const size_t row_nb_dofs = row_data.getIndices().size();
763  const size_t col_nb_dofs = col_data.getIndices().size();
764 
765  if (row_nb_dofs && col_nb_dofs) {
766 
767  locMat.resize(row_nb_dofs, col_nb_dofs, false);
768  locMat.clear();
769 
770  auto t_w = getFTensor0IntegrationWeight();
771  auto t_row_base = row_data.getFTensor0N();
772  size_t nb_face_functions = row_data.getN().size2();
773 
774  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
775 
776  const double alpha = t_w * getMeasure();
777 
778  size_t rr = 0;
779  for (; rr != row_nb_dofs / 2; ++rr) {
781  &locMat(2 * rr + 0, 0), &locMat(2 * rr + 1, 1));
782 
783  auto t_col_base = col_data.getFTensor0N(gg, 0);
784  for (size_t cc = 0; cc != col_nb_dofs / 2; ++cc) {
785  t_mat(i) += alpha * spring_stiffness * t_row_base * t_col_base;
786  ++t_col_base;
787  ++t_mat;
788  }
789 
790  ++t_row_base;
791  }
792  for (; rr < nb_face_functions; ++rr)
793  ++t_row_base;
794 
795  ++t_w;
796  }
797 
798  CHKERR MatSetValues(getSNESB(), row_data, col_data, &*locMat.data().begin(),
799  ADD_VALUES);
800  }
801 
803 }
804 
805 struct Monitor : public FEMethod {
806 
807  Monitor(SmartPetscObj<DM> &dm,
808  boost::shared_ptr<PostProcFaceOnRefinedMeshFor2D> &post_proc_fe,
809  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
810  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter)
811  : dM(dm), postProcFe(post_proc_fe), uXScatter(ux_scatter),
812  uYScatter(uy_scatter){};
813 
814  MoFEMErrorCode preProcess() { return 0; }
815  MoFEMErrorCode operator()() { return 0; }
816 
819 
820  auto make_vtk = [&]() {
823  CHKERR postProcFe->writeFile(
824  "out_contact_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
826  };
827 
828  auto print_max_min = [&](auto &tuple, const std::string msg) {
830  CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
831  INSERT_VALUES, SCATTER_FORWARD);
832  CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
833  INSERT_VALUES, SCATTER_FORWARD);
834  double max, min;
835  CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
836  CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
837  PetscPrintf(PETSC_COMM_WORLD, "%s time %3.4e min %3.4e max %3.4e\n",
838  msg.c_str(), ts_t, min, max);
840  };
841 
842  CHKERR make_vtk();
843  CHKERR print_max_min(uXScatter, "Ux");
844  CHKERR print_max_min(uYScatter, "Uy");
845 
847  }
848 
849 private:
850  SmartPetscObj<DM> dM;
851  boost::shared_ptr<PostProcFaceOnRefinedMeshFor2D> postProcFe;
852  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
853  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
854 };
855 
857  OpPostProcContact(const std::string field_name,
858  moab::Interface &post_proc_mesh,
859  std::vector<EntityHandle> &map_gauss_pts,
860  boost::shared_ptr<CommonData> common_data_ptr);
861  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
862 
863 private:
865  std::vector<EntityHandle> &mapGaussPts;
866  boost::shared_ptr<CommonData> commonDataPtr;
867 };
868 //! [Operators definitions]
869 
871  const std::string field_name, moab::Interface &post_proc_mesh,
872  std::vector<EntityHandle> &map_gauss_pts,
873  boost::shared_ptr<CommonData> common_data_ptr)
874  : DomainEleOp(field_name, DomainEleOp::OPROW), postProcMesh(post_proc_mesh),
875  mapGaussPts(map_gauss_pts), commonDataPtr(common_data_ptr) {
876  // Opetor is only executed for vertices
877  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
878 }
879 
880 //! [Postprocessing]
882  EntData &data) {
884 
885  auto get_tag_mat = [&](const std::string name) {
886  std::array<double, 9> def;
887  std::fill(def.begin(), def.end(), 0);
888  Tag th;
889  CHKERR postProcMesh.tag_get_handle(name.c_str(), 9, MB_TYPE_DOUBLE, th,
890  MB_TAG_CREAT | MB_TAG_SPARSE,
891  def.data());
892  return th;
893  };
894 
895  auto get_tag_vec = [&](const std::string name) {
896  std::array<double, 3> def;
897  std::fill(def.begin(), def.end(), 0);
898  Tag th;
899  CHKERR postProcMesh.tag_get_handle(name.c_str(), 3, MB_TYPE_DOUBLE, th,
900  MB_TAG_CREAT | MB_TAG_SPARSE,
901  def.data());
902  return th;
903  };
904 
905  MatrixDouble3by3 mat(3, 3);
906 
907  auto set_matrix = [&](auto &t) -> MatrixDouble3by3 & {
908  mat.clear();
909  for (size_t r = 0; r != 2; ++r)
910  for (size_t c = 0; c != 2; ++c)
911  mat(r, c) = t(r, c);
912  return mat;
913  };
914 
915  auto set_vector = [&](auto &t) -> MatrixDouble3by3 & {
916  mat.clear();
917  for (size_t r = 0; r != 2; ++r)
918  mat(0, r) = t(r);
919  return mat;
920  };
921 
922  auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
923  return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
924  &*mat.data().begin());
925  };
926 
927  auto th_stress = get_tag_mat("SIGMA");
928  auto th_div = get_tag_vec("DIV_SIGMA");
929 
930  size_t nb_gauss_pts = getGaussPts().size2();
931  auto t_stress = getFTensor2FromMat<2, 2>(*(commonDataPtr->contactStressPtr));
932  auto t_div =
933  getFTensor1FromMat<2>(*(commonDataPtr->contactStressDivergencePtr));
934 
935  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
936 
937  CHKERR set_tag(th_stress, gg, set_matrix(t_stress));
938  CHKERR set_tag(th_div, gg, set_vector(t_div));
939  ++t_stress;
940  ++t_div;
941  }
942 
944 }
945 //! [Postprocessing]
946 
947 }; // namespace OpContactTools
MoFEMErrorCode operator()()
Definition: ContactOps.hpp:815
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:58
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
constexpr double spring_stiffness
OpConstrainDomainLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:619
double gap0(FTensor::Tensor1< T, 3 > &t_coords, FTensor::Tensor1< double, 2 > &t_normal)
Definition: ContactOps.hpp:134
OpConstrainBoundaryRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:232
moab::Interface & postProcMesh
Definition: ContactOps.hpp:864
boost::shared_ptr< MatrixDouble > contactStressDivergencePtr
Definition: ContactOps.hpp:20
double normal_traction(FTensor::Tensor1< T, 2 > &t_traction, FTensor::Tensor1< double, 2 > &t_normal)
Definition: ContactOps.hpp:146
Definition: ContactOps.hpp:856
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: ContactOps.hpp:326
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:101
boost::shared_ptr< MatrixDouble > contactStressPtr
Definition: ContactOps.hpp:19
double diff_constrains_dtraction(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:168
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
FTensor::Tensor1< double, 2 > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, 2 > &t_disp)
Definition: ContactOps.hpp:128
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: ContactOps.hpp:476
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: ContactOps.hpp:755
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
MoFEMErrorCode preProcess()
Definition: ContactOps.hpp:814
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:71
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: ContactOps.hpp:380
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: ContactOps.hpp:566
OpConstrainDomainRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:561
const VectorInt & getIndices() const
Get global indices of dofs on entity.
OpSpringLhs(const std::string row_field_name, const std::string col_field_name)
Definition: ContactOps.hpp:749
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: ContactOps.hpp:852
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:122
bool sYmm
If true assume that matrix is symmetric structure.
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: ContactOps.hpp:181
boost::shared_ptr< MatrixDouble > contactTractionPtr
Definition: ContactOps.hpp:21
MoFEMErrorCode postProcess()
Definition: ContactOps.hpp:817
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
auto getFTensor0IntegrationWeight()
Get integration weights.
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
boost::shared_ptr< MatrixDouble > contactDispPtr
Definition: ContactOps.hpp:22
std::vector< EntityHandle > & mapGaussPts
Definition: ContactOps.hpp:865
constexpr double cn
double constrian(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:164
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:101
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:46
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
#define CHKERR
Inline error check.
Definition: definitions.h:602
double diff_constrains_dgap(double &&g0, double &&g, double &&t)
Definition: ContactOps.hpp:172
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: ContactOps.hpp:237
boost::shared_ptr< PostProcFaceOnRefinedMeshFor2D > postProcFe
Definition: ContactOps.hpp:851
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: ContactOps.hpp:706
OpInternalDomainContactRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:176
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
OpConstrainBoundaryLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:373
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
Data on single entity (This is passed as argument to DataOperator::doWork)
double sign(double x)
Definition: ContactOps.hpp:151
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:37
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]
Definition: ContactOps.hpp:870
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
SmartPetscObj< DM > dM
Definition: ContactOps.hpp:850
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: ContactOps.hpp:627
Monitor(SmartPetscObj< DM > &dm, boost::shared_ptr< PostProcFaceOnRefinedMeshFor2D > &post_proc_fe, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter >> ux_scatter, std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter >> uy_scatter)
Definition: ContactOps.hpp:807
OpConstrainBoundaryTraction(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:318
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:81
double gap(FTensor::Tensor1< T, 2 > &t_disp, FTensor::Tensor1< double, 2 > &t_normal)
Definition: ContactOps.hpp:140
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:866
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: ContactOps.hpp:853
OpSpringRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:701
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Index< 'l', 2 > l
Definition: ContactOps.hpp:29
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:110
double w(const double g, const double t)
Definition: ContactOps.hpp:160
OpConstrainBoundaryLhs_dTraction(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:469
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
Definition: ContactOps.hpp:881