v0.13.1
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
15namespace 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;
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);
37
38private:
39 boost::shared_ptr<CommonData> commonDataPtr;
41};
42
44 OpConstrainBoundaryRhs(const std::string field_name,
45 boost::shared_ptr<CommonData> common_data_ptr);
47
48private:
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
59private:
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
70private:
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);
80
81private:
83 boost::shared_ptr<CommonData> commonDataPtr;
84 bool uPdate;
85};
86
87
88template <typename T>
89inline 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
94template <typename T1, typename T2>
95inline double normal_traction(Tensor1<T1, 3> &t_traction,
96 Tensor1<T2, 3> &t_normal) {
97 return -t_traction(i) * t_normal(i);
98}
99
100inline 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
109inline double w(const double g, const double t) {
110 return g - (*cache).cn_cont * t;
111}
112
113inline double constraint(double g0, double g, double &&t) {
114 return (w(g, t) + std::abs(w(g, t))) / 2 + g0;
115};
116
117inline double diff_constraints_dtraction(double g, double &&t) {
118 return -(*cache).cn_cont * (1 + sign(w(g, t))) / 2;
119}
120
121inline 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)
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)
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) {
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),
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
591private:
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)
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
698 boost::shared_ptr<CommonData> common_data_ptr)
700 commonDataPtr(common_data_ptr) {
701 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
702 doEntities[MBTRI] = doEntities[MBQUAD] = true;
703 }
704
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
738private:
739 boost::shared_ptr<CommonData> commonDataPtr;
740};
741
743
744 OpGetContactArea(const std::string field_name,
745 boost::shared_ptr<CommonData> common_data_ptr)
747 commonDataPtr(common_data_ptr) {
748 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
749 doEntities[MBTRI] = doEntities[MBQUAD] = true;
750 }
751
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
800private:
801 boost::shared_ptr<CommonData> commonDataPtr;
802};
803
804}; // namespace OpContactTools
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:105
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
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:477
const double r
rate factor
double cn
Definition: plastic.cpp:102
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
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 & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
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)
Operator for linear form, usually to calculate values on right hand side.
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)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
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)
Operator for linear form, usually to calculate values on right hand side.
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)