v0.13.2
Loading...
Searching...
No Matches
ContactOps.hpp
Go to the documentation of this file.
1
2
3/**
4 * \file ContactOps.hpp
5 * \example ContactOps.hpp
6 */
7
8#ifndef __CONTACTOPS_HPP__
9#define __CONTACTOPS_HPP__
10
11namespace ContactOps {
12
13//! [Common data]
14struct CommonData : public boost::enable_shared_from_this<CommonData> {
15 // MatrixDouble contactStress;
16 MatrixDouble contactTraction;
17 MatrixDouble contactDisp;
18
19 static SmartPetscObj<Vec>
20 totalTraction; // User have to release and create vector when appropiate.
21
22 static auto createTotalTraction(MoFEM::Interface &m_field) {
23 constexpr int ghosts[] = {0, 1, 2};
25 createGhostVector(m_field.get_comm(),
26
27 (m_field.get_comm_rank() == 0) ? 3 : 0, 3,
28
29 (m_field.get_comm_rank() == 0) ? 0 : 3, ghosts);
30 return totalTraction;
31 }
32
35 const double *t_ptr;
36 CHK_THROW_MESSAGE(VecGetArrayRead(CommonData::totalTraction, &t_ptr),
37 "get array");
38 FTensor::Tensor1<double, 3> t{t_ptr[0], t_ptr[1], t_ptr[2]};
39 CHK_THROW_MESSAGE(VecRestoreArrayRead(CommonData::totalTraction, &t_ptr),
40 "restore array");
41 return t;
42 } else {
43 return FTensor::Tensor1<double, 3>{0., 0., 0.};
44 }
45 }
46
47 inline auto contactTractionPtr() {
48 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
50 }
51
52 inline auto contactDispPtr() {
53 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactDisp);
54 }
55};
56
57SmartPetscObj<Vec> CommonData::totalTraction;
58
59//! [Common data]
60
61//! [Surface distance function from python]
62#ifdef PYTHON_SFD
63struct SDFPython {
64 SDFPython() = default;
65 virtual ~SDFPython() = default;
66
67 MoFEMErrorCode sdfInit(const std::string py_file) {
69 try {
70
71 // create main module
72 auto main_module = bp::import("__main__");
73 mainNamespace = main_module.attr("__dict__");
74 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
75 // create a reference to python function
76 sdfFun = mainNamespace["sdf"];
77 sdfGradFun = mainNamespace["grad_sdf"];
78 sdfHessFun = mainNamespace["hess_sdf"];
79
80 } catch (bp::error_already_set const &) {
81 // print all other errors to stderr
82 PyErr_Print();
84 }
86 };
87
88 template <typename T>
89 inline std::vector<T>
90 py_list_to_std_vector(const boost::python::object &iterable) {
91 return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
92 boost::python::stl_input_iterator<T>());
93 }
94
95 MoFEMErrorCode evalSdf(
96
97 double t, double x, double y, double z, double tx, double ty, double tz,
98 double &sdf
99
100 ) {
102 try {
103
104 // call python function
105 sdf = bp::extract<double>(sdfFun(t, x, y, z, tx, ty, tz));
106
107 } catch (bp::error_already_set const &) {
108 // print all other errors to stderr
109 PyErr_Print();
111 }
113 }
114
115 MoFEMErrorCode evalGradSdf(
116
117 double t, double x, double y, double z, double tx, double ty, double tz,
118 std::vector<double> &grad_sdf
119
120 ) {
122 try {
123
124 // call python function
125 grad_sdf =
126 py_list_to_std_vector<double>(sdfGradFun(t, x, y, z, tx, ty, tz));
127
128 } catch (bp::error_already_set const &) {
129 // print all other errors to stderr
130 PyErr_Print();
132 }
133
134 if (grad_sdf.size() != 3) {
135 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Expected size 3");
136 }
137
139 }
140
141 MoFEMErrorCode evalHessSdf(
142
143 double t, double x, double y, double z, double tx, double ty, double tz,
144 std::vector<double> &hess_sdf
145
146 ) {
148 try {
149
150 // call python function
151 hess_sdf =
152 py_list_to_std_vector<double>(sdfHessFun(t, x, y, z, tx, ty, tz));
153
154 } catch (bp::error_already_set const &) {
155 // print all other errors to stderr
156 PyErr_Print();
158 }
159
160 if (hess_sdf.size() != 6) {
161 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Expected size 6");
162 }
163
165 }
166
167private:
168 bp::object mainNamespace;
169 bp::object sdfFun;
170 bp::object sdfGradFun;
171 bp::object sdfHessFun;
172};
173
174static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
175
176#endif
177//! [Surface distance function from python]
178
179using SurfaceDistanceFunction = boost::function<double(
180 double t, double x, double y, double z, double tx, double ty, double tz)>;
181
182using GradSurfaceDistanceFunction = boost::function<FTensor::Tensor1<double, 3>(
183 double t, double x, double y, double z, double tx, double ty, double tz)>;
184
186 boost::function<FTensor::Tensor2_symmetric<double, 3>(
187 double t, double x, double y, double z, double tx, double ty,
188 double tz)>;
189
190inline double surface_distance_function(double t, double x, double y, double z,
191 double tx, double ty, double tz) {
192
193#ifdef PYTHON_SFD
194 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
195 double sdf;
196 CHK_MOAB_THROW(sdf_ptr->evalSdf(t, x, y, z, tx, ty, tz, sdf),
197 "Failed python call");
198 return sdf;
199 }
200#endif
201 return y + 0.5;
202}
203
205grad_surface_distance_function(double t, double x, double y, double z,
206 double tx, double ty, double tz) {
207#ifdef PYTHON_SFD
208 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
209 std::vector<double> grad_sdf;
210 CHK_MOAB_THROW(sdf_ptr->evalGradSdf(t, x, y, z, tx, ty, tz, grad_sdf),
211 "Failed python call");
212 return FTensor::Tensor1<double, 3>{grad_sdf[0], grad_sdf[1], grad_sdf[2]};
213 }
214#endif
215 return FTensor::Tensor1<double, 3>{0., 1., 0.};
216}
217
219hess_surface_distance_function(double t, double x, double y, double z,
220 double tx, double ty, double tz) {
221#ifdef PYTHON_SFD
222 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
223 std::vector<double> hess_sdf;
224 CHK_MOAB_THROW(sdf_ptr->evalHessSdf(t, x, y, z, tx, ty, tz, hess_sdf),
225 "Failed python call");
226 return FTensor::Tensor2_symmetric<double, 3>{hess_sdf[0], hess_sdf[1],
227 hess_sdf[2], hess_sdf[3],
228 hess_sdf[4], hess_sdf[5]};
229 }
230#endif
231 return FTensor::Tensor2_symmetric<double, 3>{0., 0., 0., 0., 0., 0.};
232}
233
234template <int DIM, IntegrationType I, typename BoundaryEleOp>
236
237template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
239
240template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
242
243template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
245
246template <int DIM, typename BoundaryEleOp>
248 : public BoundaryEleOp {
250 boost::shared_ptr<CommonData> common_data_ptr);
251 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
252
253private:
254 boost::shared_ptr<CommonData> commonDataPtr;
255};
256
257template <int DIM, typename AssemblyBoundaryEleOp>
259 : public AssemblyBoundaryEleOp {
260 OpConstrainBoundaryRhsImpl(const std::string field_name,
261 boost::shared_ptr<CommonData> common_data_ptr);
262 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
263
265 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
267
268private:
269 boost::shared_ptr<CommonData> commonDataPtr;
270};
271
272template <int DIM, typename AssemblyBoundaryEleOp>
274 : public AssemblyBoundaryEleOp {
275 OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
276 const std::string col_field_name,
277 boost::shared_ptr<CommonData> common_data_ptr);
278 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
279 EntitiesFieldData::EntData &col_data);
280
282 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
284 HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
286
287 boost::shared_ptr<CommonData> commonDataPtr;
288};
289
290template <int DIM, typename AssemblyBoundaryEleOp>
292 : public AssemblyBoundaryEleOp {
294 const std::string row_field_name, const std::string col_field_name,
295 boost::shared_ptr<CommonData> common_data_ptr);
296 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
297 EntitiesFieldData::EntData &col_data);
298
300 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
302
303private:
304 boost::shared_ptr<CommonData> commonDataPtr;
305};
306
307template <typename BoundaryEleOp> struct ContactIntegrators {
308
309 template <int DIM, IntegrationType I>
312
313 template <AssemblyType A> struct Assembly {
314
316 typename FormsIntegrators<BoundaryEleOp>::template Assembly<A>::OpBase;
317
318 template <int DIM, IntegrationType I>
321
322 template <int DIM, IntegrationType I>
325
326 template <int DIM, IntegrationType I>
329 };
330};
331
332inline double sign(double x) {
333 if (x == 0)
334 return 0;
335 else if (x > 0)
336 return 1;
337 else
338 return -1;
339};
340
341inline double w(const double sdf, const double tn) {
342 return sdf - cn_contact * tn;
343}
344
345/**
346 * @brief constrain function
347 *
348 * return 1 if negative sdn or positive tn
349 *
350 * @param sdf signed distance
351 * @param tn traction
352 * @return double
353 */
354inline double constrain(double sdf, double tn) {
355 const auto s = sign(w(sdf, tn));
356 return (1 - s) / 2;
357}
358
359template <int DIM, typename BoundaryEleOp>
362 boost::shared_ptr<CommonData> common_data_ptr)
364 commonDataPtr(common_data_ptr) {}
365
366template <int DIM, typename BoundaryEleOp>
367MoFEMErrorCode
369 int side, EntityType type, EntData &data) {
371
373 FTensor::Tensor1<double, 3> t_sum_t{0., 0., 0.};
374
376 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
377
378 const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
379 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
380 const double alpha = t_w * BoundaryEleOp::getMeasure();
381 t_sum_t(i) += alpha * t_traction(i);
382 ++t_w;
383 ++t_traction;
384 }
385
386 constexpr int ind[] = {0, 1, 2};
387 CHKERR VecSetValues(commonDataPtr->totalTraction, 3, ind, &t_sum_t(0),
388 ADD_VALUES);
389
391}
392
393template <int DIM, typename AssemblyBoundaryEleOp>
396 boost::shared_ptr<CommonData> common_data_ptr)
398 commonDataPtr(common_data_ptr) {}
399
400template <int DIM, typename AssemblyBoundaryEleOp>
401MoFEMErrorCode
403 EntitiesFieldData::EntData &data) {
405
410
411 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
412
413 auto &nf = AssemblyBoundaryEleOp::locF;
414
415 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
416 auto t_total_traction = CommonData::getFTensor1TotalTraction();
417
418 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
419 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
420 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
421 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
422
423 size_t nb_base_functions = data.getN().size2() / 3;
424 auto t_base = data.getFTensor1N<3>();
425 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
426
428 t_normal(i) =
429 t_normal_at_pts(i) / std::sqrt(t_normal_at_pts(i) * t_normal_at_pts(i));
430
431 auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
432 const double alpha = t_w * AssemblyBoundaryEleOp::getMeasure();
433
434 FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
435 t_spatial_coords(i) = t_coords(i) + t_disp(i);
436
437 auto ts_a = AssemblyBoundaryEleOp::getTStime();
438
439 auto sdf = surfaceDistanceFunction(
440 ts_a, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
441 t_total_traction(0), t_total_traction(1), t_total_traction(2));
442
443 auto t_grad_sdf = gradSurfaceDistanceFunction(
444 ts_a, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
445 t_total_traction(0), t_total_traction(1), t_total_traction(2));
446
447 auto tn = -t_traction(i) * t_grad_sdf(i);
448 auto c = constrain(sdf, tn);
449
451 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
453 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
454
456 t_rhs(i) =
457
458 t_cQ(i, j) * (t_disp(j) - cn_contact * t_traction(j))
459
460 +
461
462 t_cP(i, j) * t_disp(j) +
463 c * (sdf * t_grad_sdf(i)); // add gap0 displacements
464
465 size_t bb = 0;
466 for (; bb != AssemblyBoundaryEleOp::nbRows / DIM; ++bb) {
467 const double beta = alpha * (t_base(i) * t_normal(i));
468 t_nf(i) -= beta * t_rhs(i);
469
470 ++t_nf;
471 ++t_base;
472 }
473 for (; bb < nb_base_functions; ++bb)
474 ++t_base;
475
476 ++t_disp;
477 ++t_traction;
478 ++t_coords;
479 ++t_w;
480 ++t_normal_at_pts;
481 }
482
484}
485
486template <int DIM, typename AssemblyBoundaryEleOp>
488 OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
489 const std::string col_field_name,
490 boost::shared_ptr<CommonData> common_data_ptr)
491 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
492 DomainEleOp::OPROWCOL),
493 commonDataPtr(common_data_ptr) {
494 AssemblyBoundaryEleOp::sYmm = false;
495}
496
497template <int DIM, typename AssemblyBoundaryEleOp>
498MoFEMErrorCode
500 EntitiesFieldData::EntData &row_data,
501 EntitiesFieldData::EntData &col_data) {
503
507
508 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
509 auto &locMat = AssemblyBoundaryEleOp::locMat;
510
511 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
512 auto t_total_traction = CommonData::getFTensor1TotalTraction();
513
514 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
515 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
516 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
517
518 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
519 auto t_row_base = row_data.getFTensor1N<3>();
520 size_t nb_face_functions = row_data.getN().size2() / 3;
521
522 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
523
524 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
525
527 t_normal(i) =
528 t_normal_at_pts(i) / std::sqrt(t_normal_at_pts(i) * t_normal_at_pts(i));
529
530 const double alpha = t_w * AssemblyBoundaryEleOp::getMeasure();
531
532 FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
533 t_spatial_coords(i) = t_coords(i) + t_disp(i);
534
535 auto ts_time = AssemblyBoundaryEleOp::getTStime();
536
537 auto sdf = surfaceDistanceFunction(
538 ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
539 t_total_traction(0), t_total_traction(1), t_total_traction(2));
540 auto t_grad_sdf = gradSurfaceDistanceFunction(
541 ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
542 t_total_traction(0), t_total_traction(1), t_total_traction(2));
543
544 auto tn = -t_traction(i) * t_grad_sdf(i);
545 auto c = constrain(sdf, tn);
546
548 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
550 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
551
553 t_res_dU(i, j) = kronecker_delta(i, j) + t_cP(i, j);
554
555 if (c > 0) {
556 auto t_hess_sdf = hessSurfaceDistanceFunction(
557 ts_time, t_spatial_coords(0), t_spatial_coords(1),
558 t_spatial_coords(2), t_total_traction(0), t_total_traction(1),
559 t_total_traction(2));
560 t_res_dU(i, j) +=
561 (c * cn_contact) *
562 (t_hess_sdf(i, j) * (t_grad_sdf(k) * t_traction(k)) +
563 t_grad_sdf(i) * t_hess_sdf(k, j) * t_traction(k)) +
564 c * sdf * t_hess_sdf(i, j);
565 }
566
567 size_t rr = 0;
568 for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
569
570 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
571
572 const double row_base = t_row_base(i) * t_normal(i);
573
574 auto t_col_base = col_data.getFTensor0N(gg, 0);
575 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
576 const double beta = alpha * row_base * t_col_base;
577
578 t_mat(i, j) -= beta * t_res_dU(i, j);
579
580 ++t_col_base;
581 ++t_mat;
582 }
583
584 ++t_row_base;
585 }
586 for (; rr < nb_face_functions; ++rr)
587 ++t_row_base;
588
589 ++t_disp;
590 ++t_traction;
591 ++t_coords;
592 ++t_w;
593 ++t_normal_at_pts;
594 }
595
597}
598
599template <int DIM, typename AssemblyBoundaryEleOp>
602 const std::string row_field_name, const std::string col_field_name,
603 boost::shared_ptr<CommonData> common_data_ptr)
604 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
605 DomainEleOp::OPROWCOL),
606 commonDataPtr(common_data_ptr) {
607 AssemblyBoundaryEleOp::sYmm = false;
608}
609
610template <int DIM, typename AssemblyBoundaryEleOp>
611MoFEMErrorCode
613 iNtegrate(EntitiesFieldData::EntData &row_data,
614 EntitiesFieldData::EntData &col_data) {
616
620
621 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
622 auto &locMat = AssemblyBoundaryEleOp::locMat;
623
624 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
625 auto t_total_traction = CommonData::getFTensor1TotalTraction();
626
627 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
628 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
629 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
630
631 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
632 auto t_row_base = row_data.getFTensor1N<3>();
633 size_t nb_face_functions = row_data.getN().size2() / 3;
634
635 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
636
638 t_normal(i) =
639 t_normal_at_pts(i) / std::sqrt(t_normal_at_pts(i) * t_normal_at_pts(i));
640
641 const double alpha = t_w * AssemblyBoundaryEleOp::getMeasure();
642
643 FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
644 t_spatial_coords(i) = t_coords(i) + t_disp(i);
645
646 auto ts_time = AssemblyBoundaryEleOp::getTStime();
647
648 auto sdf = surfaceDistanceFunction(
649 ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
650 t_total_traction(0), t_total_traction(1), t_total_traction(2));
651 auto t_grad_sdf = gradSurfaceDistanceFunction(
652 ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
653 t_total_traction(0), t_total_traction(1), t_total_traction(2));
654
655 auto tn = -t_traction(i) * t_grad_sdf(i);
656 auto c = constrain(sdf, tn);
657
659 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
661 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
662
664 t_res_dt(i, j) = -cn_contact * t_cQ(i, j);
665
666 size_t rr = 0;
667 for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
668
669 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
670 const double row_base = t_row_base(i) * t_normal(i);
671
672 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
673 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
674 const double col_base = t_col_base(i) * t_normal(i);
675 const double beta = alpha * row_base * col_base;
676
677 t_mat(i, j) -= beta * t_res_dt(i, j);
678
679 ++t_col_base;
680 ++t_mat;
681 }
682
683 ++t_row_base;
684 }
685 for (; rr < nb_face_functions; ++rr)
686 ++t_row_base;
687
688 ++t_disp;
689 ++t_traction;
690 ++t_coords;
691 ++t_w;
692 ++t_normal_at_pts;
693 }
694
696}
697
698template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
699MoFEMErrorCode opFactoryDomainRhs(
700 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
701 std::string sigma, std::string u) {
703 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
704 A>::template LinearForm<I>;
705 using OpMixDivURhs = typename B::template OpMixDivTimesU<3, DIM, DIM>;
706 using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
708 typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
710 typename B::template OpGradTimesTensor<1, DIM, DIM>;
711
712 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
713 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
714 auto div_stress_ptr = boost::make_shared<MatrixDouble>();
715 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
716
717 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
718 u, common_data_ptr->contactDispPtr()));
719 pip.push_back(
720 new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
721 pip.push_back(
722 new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
723
724 pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(u, mat_grad_ptr));
725
726 pip.push_back(
727 new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(),
728 [](double, double, double) constexpr { return 1; }));
729 pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr));
730 pip.push_back(new OpMixUTimesDivLambdaRhs(u, div_stress_ptr));
731 pip.push_back(new OpMixUTimesLambdaRhs(u, contact_stress_ptr));
732
734}
735
736template <typename OpMixLhs> struct OpMixLhsSide : public OpMixLhs {
737 using OpMixLhs::OpMixLhs;
738 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
739 EntityType col_type,
740 EntitiesFieldData::EntData &row_data,
741 EntitiesFieldData::EntData &col_data) {
743 auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
744 auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
745 // Only assemble side which correspond to edge entity on boundary
746 if (side_fe_entity == side_fe_data) {
747 CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
748 col_data);
749 }
751 }
752};
753
754template <int DIM, AssemblyType A, IntegrationType I, typename DomainEle>
756 MoFEM::Interface &m_field,
757 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
758 std::string fe_domain_name, std::string sigma, std::string u,
759 std::string geom, ForcesAndSourcesCore::RuleHookFun rule) {
761 auto op_loop_side = new OpLoopSide<DomainEle>(
762 m_field, fe_domain_name, DIM, Sev::noisy,
763 boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
764 pip.push_back(op_loop_side);
765
766 CHKERR AddHOOps<DIM, DIM, DIM>::add(op_loop_side->getOpPtrVector(),
767 {H1, HDIV}, geom);
768
769 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
770 A>::template BiLinearForm<I>;
771
772 using OpMixDivULhs = typename B::template OpMixDivTimesVec<DIM>;
773 using OpLambdaGraULhs = typename B::template OpMixTensorTimesGrad<DIM>;
774 using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
775 using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
776
777 auto unity = []() { return 1; };
778 op_loop_side->getOpPtrVector().push_back(
779 new OpMixDivULhsSide(sigma, u, unity, true));
780 op_loop_side->getOpPtrVector().push_back(
781 new OpLambdaGraULhsSide(sigma, u, unity, true));
782
783 op_loop_side->getSideFEPtr()->getRuleHook = rule;
785}
786
787template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
788MoFEMErrorCode opFactoryBoundaryLhs(
789 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
790 std::string sigma, std::string u) {
792
794
795 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
796
797 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
798 u, common_data_ptr->contactDispPtr()));
799 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
800 sigma, common_data_ptr->contactTractionPtr()));
801 pip.push_back(
802 new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
803 DIM, GAUSS>(sigma, u, common_data_ptr));
804 pip.push_back(new typename C::template Assembly<A>::
805 template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
806 sigma, sigma, common_data_ptr));
807
809}
810
811template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
812MoFEMErrorCode opFactoryBoundaryRhs(
813 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
814 std::string sigma, std::string u) {
816
818
819 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
820
821 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
822 u, common_data_ptr->contactDispPtr()));
823 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
824 sigma, common_data_ptr->contactTractionPtr()));
825 pip.push_back(
826 new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
827 DIM, GAUSS>(sigma, common_data_ptr));
828
830}
831
832template <int DIM, IntegrationType I, typename BoundaryEleOp>
834 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
835 std::string sigma) {
837
839
840 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
841 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
842 sigma, common_data_ptr->contactTractionPtr()));
843 pip.push_back(new typename C::template OpAssembleTotalContactTraction<DIM, I>(
844 common_data_ptr));
845
847}
848
849}; // namespace ContactOps
850
851#endif // __CONTACTOPS_HPP__
Kronecker Delta class.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
@ NOSPACE
Definition: definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#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
constexpr auto t_kd
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
[Only used for dynamics]
MoFEMErrorCode opFactoryBoundaryLhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u)
Definition: ContactOps.hpp:788
double cn_contact
Definition: contact.cpp:109
boost::function< double(double t, double x, double y, double z, double tx, double ty, double tz)> SurfaceDistanceFunction
[Common data]
Definition: ContactOps.hpp:180
double surface_distance_function(double t, double x, double y, double z, double tx, double ty, double tz)
Definition: ContactOps.hpp:190
boost::function< FTensor::Tensor1< double, 3 >(double t, double x, double y, double z, double tx, double ty, double tz)> GradSurfaceDistanceFunction
Definition: ContactOps.hpp:183
boost::function< FTensor::Tensor2_symmetric< double, 3 >(double t, double x, double y, double z, double tx, double ty, double tz)> HessSurfaceDistanceFunction
Definition: ContactOps.hpp:188
MoFEMErrorCode opFactoryDomainRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u)
Definition: ContactOps.hpp:699
double w(const double sdf, const double tn)
Definition: ContactOps.hpp:341
FTensor::Tensor2_symmetric< double, 3 > hess_surface_distance_function(double t, double x, double y, double z, double tx, double ty, double tz)
Definition: ContactOps.hpp:219
FTensor::Tensor1< double, 3 > grad_surface_distance_function(double t, double x, double y, double z, double tx, double ty, double tz)
Definition: ContactOps.hpp:205
MoFEMErrorCode opFactoryBoundaryToDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string fe_domain_name, std::string sigma, std::string u, std::string geom, ForcesAndSourcesCore::RuleHookFun rule)
Definition: ContactOps.hpp:755
EntitiesFieldData::EntData EntData
double constrain(double sdf, double tn)
constrain function
Definition: ContactOps.hpp:354
MoFEMErrorCode opFactoryCalculateTraction(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma)
Definition: ContactOps.hpp:833
MoFEMErrorCode opFactoryBoundaryRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u)
Definition: ContactOps.hpp:812
double sign(double x)
Definition: ContactOps.hpp:332
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
Definition: sdf.py:1
def grad_sdf(t, x, y, z, tx, ty, tz)
Definition: sdf.py:15
def hess_sdf(t, x, y, z, tx, ty, tz)
Definition: sdf.py:19
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr AssemblyType A
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
MatrixDouble contactDisp
Definition: ContactOps.hpp:17
static auto getFTensor1TotalTraction()
Definition: ContactOps.hpp:33
MatrixDouble contactTraction
Definition: ContactOps.hpp:16
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:20
static auto createTotalTraction(MoFEM::Interface &m_field)
Definition: ContactOps.hpp:22
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:738
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
auto getFTensor0IntegrationWeight()
Get integration weights.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element