v0.14.0
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 VectorDouble sdfVals; ///< size is equal to number of gauss points on element
20 MatrixDouble gradsSdf; ///< nb of rows is equals to dimension, and nb of cols
21 ///< is equals to number of gauss points on element
22 MatrixDouble hessSdf; ///< nb of rows is equals to nb of element of symmetric
23 ///< matrix, and nb of cols is equals to number of gauss
24 ///< points on element
25 VectorDouble constraintVals;
26
27 static SmartPetscObj<Vec>
28 totalTraction; // User have to release and create vector when appropiate.
29
30 static auto createTotalTraction(MoFEM::Interface &m_field) {
31 constexpr int ghosts[] = {0, 1, 2};
33 createGhostVector(m_field.get_comm(),
34
35 (m_field.get_comm_rank() == 0) ? 3 : 0, 3,
36
37 (m_field.get_comm_rank() == 0) ? 0 : 3, ghosts);
38 return totalTraction;
39 }
40
43 const double *t_ptr;
44 CHK_THROW_MESSAGE(VecGetArrayRead(CommonData::totalTraction, &t_ptr),
45 "get array");
46 FTensor::Tensor1<double, 3> t{t_ptr[0], t_ptr[1], t_ptr[2]};
47 CHK_THROW_MESSAGE(VecRestoreArrayRead(CommonData::totalTraction, &t_ptr),
48 "restore array");
49 return t;
50 } else {
51 return FTensor::Tensor1<double, 3>{0., 0., 0.};
52 }
53 }
54
55 inline auto contactTractionPtr() {
56 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
58 }
59
60 inline auto contactDispPtr() {
61 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactDisp);
62 }
63
64 inline auto sdfPtr() {
65 return boost::shared_ptr<VectorDouble>(shared_from_this(), &sdfVals);
66 }
67
68 inline auto gradSdfPtr() {
69 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &gradsSdf);
70 }
71
72 inline auto hessSdfPtr() {
73 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &hessSdf);
74 }
75
76 inline auto constraintPtr() {
77 return boost::shared_ptr<VectorDouble>(shared_from_this(), &constraintVals);
78 }
79};
80
81SmartPetscObj<Vec> CommonData::totalTraction;
82
83//! [Common data]
84
85//! [Surface distance function from python]
86#ifdef PYTHON_SDF
87struct SDFPython {
88 SDFPython() = default;
89 virtual ~SDFPython() = default;
90
91 MoFEMErrorCode sdfInit(const std::string py_file) {
93 try {
94
95 // create main module
96 auto main_module = bp::import("__main__");
97 mainNamespace = main_module.attr("__dict__");
98 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
99 // create a reference to python function
100 sdfFun = mainNamespace["sdf"];
101 sdfGradFun = mainNamespace["grad_sdf"];
102 sdfHessFun = mainNamespace["hess_sdf"];
103
104 } catch (bp::error_already_set const &) {
105 // print all other errors to stderr
106 PyErr_Print();
108 }
110 };
111
112 template <typename T>
113 inline std::vector<T>
114 py_list_to_std_vector(const boost::python::object &iterable) {
115 return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
116 boost::python::stl_input_iterator<T>());
117 }
118
119 MoFEMErrorCode evalSdf(
120
121 double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
122 np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
123 np::ndarray &sdf
124
125 ) {
127 try {
128
129 // call python function
130 sdf = bp::extract<np::ndarray>(
131 sdfFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
132
133 } catch (bp::error_already_set const &) {
134 // print all other errors to stderr
135 PyErr_Print();
137 }
139 }
140
141 MoFEMErrorCode evalGradSdf(
142
143 double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
144 np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
145 np::ndarray &grad_sdf
146
147 ) {
149 try {
150
151 // call python function
152 grad_sdf = bp::extract<np::ndarray>(
153 sdfGradFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
154
155 } catch (bp::error_already_set const &) {
156 // print all other errors to stderr
157 PyErr_Print();
159 }
161 }
162
163 MoFEMErrorCode evalHessSdf(
164
165 double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
166 np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
167 np::ndarray &hess_sdf
168
169 ) {
171 try {
172
173 // call python function
174 hess_sdf = bp::extract<np::ndarray>(
175 sdfHessFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
176
177 } catch (bp::error_already_set const &) {
178 // print all other errors to stderr
179 PyErr_Print();
181 }
183 }
184
185private:
186 bp::object mainNamespace;
187 bp::object sdfFun;
188 bp::object sdfGradFun;
189 bp::object sdfHessFun;
190};
191
192static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
193
194inline np::ndarray convert_to_numpy(VectorDouble &data, int nb_gauss_pts,
195 int id) {
196 auto dtype = np::dtype::get_builtin<double>();
197 auto size = bp::make_tuple(nb_gauss_pts);
198 auto stride = bp::make_tuple(3 * sizeof(double));
199 return (np::from_data(&data[id], dtype, size, stride, bp::object()));
200};
201#endif
202//! [Surface distance function from python]
203
204using SurfaceDistanceFunction = boost::function<VectorDouble(
205 double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
206 MatrixDouble &normals_at_pts, int block_id)>;
207
208using GradSurfaceDistanceFunction = boost::function<MatrixDouble(
209 double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
210 MatrixDouble &normals_at_pts, int block_id)>;
211
212using HessSurfaceDistanceFunction = boost::function<MatrixDouble(
213 double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
214 MatrixDouble &normals_at_pts, int block_id)>;
215
216inline VectorDouble surface_distance_function(double delta_t, double t,
217 int nb_gauss_pts,
218 MatrixDouble &m_spatial_coords,
219 MatrixDouble &m_normals_at_pts,
220 int block_id) {
221
222#ifdef PYTHON_SDF
223 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
224
225 VectorDouble v_spatial_coords = m_spatial_coords.data();
226 VectorDouble v_normal_at_pts = m_normals_at_pts.data();
227
228 bp::list python_coords;
229 bp::list python_normals;
230
231 for (int idx = 0; idx < 3; ++idx) {
232 python_coords.append(
233 convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
234 python_normals.append(
235 convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
236 }
237
238 np::ndarray np_sdf = np::empty(bp::make_tuple(nb_gauss_pts),
239 np::dtype::get_builtin<double>());
240 CHK_MOAB_THROW(sdf_ptr->evalSdf(delta_t, t,
241 bp::extract<np::ndarray>(python_coords[0]),
242 bp::extract<np::ndarray>(python_coords[1]),
243 bp::extract<np::ndarray>(python_coords[2]),
244 bp::extract<np::ndarray>(python_normals[0]),
245 bp::extract<np::ndarray>(python_normals[1]),
246 bp::extract<np::ndarray>(python_normals[2]),
247 block_id, np_sdf),
248 "Failed python call");
249
250 double *sdf_val_ptr = reinterpret_cast<double *>(np_sdf.get_data());
251
252 VectorDouble v_sdf;
253 v_sdf.resize(nb_gauss_pts, false);
254
255 for (size_t gg = 0; gg < nb_gauss_pts; ++gg)
256 v_sdf[gg] = *(sdf_val_ptr + gg);
257
258 return v_sdf;
259 }
260#endif
261 VectorDouble v_sdf;
262 v_sdf.resize(nb_gauss_pts, false);
263 m_spatial_coords.resize(3, nb_gauss_pts);
264 auto t_coords = getFTensor1FromMat<3>(m_spatial_coords);
265
266 for (size_t gg = 0; gg < nb_gauss_pts; ++gg)
267 v_sdf[gg] = t_coords(1) + 0.5;
268
269 return v_sdf;
270}
271
272inline MatrixDouble
273grad_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
274 MatrixDouble &m_spatial_coords,
275 MatrixDouble &m_normals_at_pts, int block_id) {
276#ifdef PYTHON_SDF
277 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
278
279 VectorDouble v_spatial_coords = m_spatial_coords.data();
280 VectorDouble v_normal_at_pts = m_normals_at_pts.data();
281
282 bp::list python_coords;
283 bp::list python_normals;
284
285 for (int idx = 0; idx < 3; ++idx) {
286 python_coords.append(
287 convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
288 python_normals.append(
289 convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
290 }
291
292 np::ndarray np_grad_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 3),
293 np::dtype::get_builtin<double>());
294 CHK_MOAB_THROW(sdf_ptr->evalGradSdf(
295 delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
296 bp::extract<np::ndarray>(python_coords[1]),
297 bp::extract<np::ndarray>(python_coords[2]),
298 bp::extract<np::ndarray>(python_normals[0]),
299 bp::extract<np::ndarray>(python_normals[1]),
300 bp::extract<np::ndarray>(python_normals[2]), block_id,
301 np_grad_sdf),
302 "Failed python call");
303
304 double *grad_ptr = reinterpret_cast<double *>(np_grad_sdf.get_data());
305
306 MatrixDouble m_grad_sdf;
307 m_grad_sdf.resize(3, nb_gauss_pts, false);
308 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
309 for (int idx = 0; idx < 3; ++idx)
310 m_grad_sdf(idx, gg) = *(grad_ptr + (3 * gg + idx));
311 }
312 return m_grad_sdf;
313 }
314#endif
315 MatrixDouble m_grad_sdf;
316 m_grad_sdf.resize(3, nb_gauss_pts, false);
318 FTensor::Tensor1<double, 3> t_grad_sdf_set{0.0, 0.0, 1.0};
319 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
320
321 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
322 t_grad_sdf(i) = t_grad_sdf_set(i);
323 ++t_grad_sdf;
324 }
325
326 return m_grad_sdf;
327}
328
329inline MatrixDouble
330hess_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
331 MatrixDouble &m_spatial_coords,
332 MatrixDouble &m_normals_at_pts, int block_id) {
333#ifdef PYTHON_SDF
334 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
335
336 VectorDouble v_spatial_coords = m_spatial_coords.data();
337 VectorDouble v_normal_at_pts = m_normals_at_pts.data();
338
339 bp::list python_coords;
340 bp::list python_normals;
341
342 for (int idx = 0; idx < 3; ++idx) {
343 python_coords.append(
344 convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
345 python_normals.append(
346 convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
347 };
348
349 np::ndarray np_hess_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 6),
350 np::dtype::get_builtin<double>());
351 CHK_MOAB_THROW(sdf_ptr->evalHessSdf(
352 delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
353 bp::extract<np::ndarray>(python_coords[1]),
354 bp::extract<np::ndarray>(python_coords[2]),
355 bp::extract<np::ndarray>(python_normals[0]),
356 bp::extract<np::ndarray>(python_normals[1]),
357 bp::extract<np::ndarray>(python_normals[2]), block_id,
358 np_hess_sdf),
359 "Failed python call");
360
361 double *hess_ptr = reinterpret_cast<double *>(np_hess_sdf.get_data());
362
363 MatrixDouble m_hess_sdf;
364 m_hess_sdf.resize(6, nb_gauss_pts, false);
365 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
366 for (int idx = 0; idx < 6; ++idx)
367 m_hess_sdf(idx, gg) =
368 *(hess_ptr + (6 * gg + idx));
369 }
370 return m_hess_sdf;
371 }
372#endif
373 MatrixDouble m_hess_sdf;
374 m_hess_sdf.resize(6, nb_gauss_pts, false);
377 FTensor::Tensor2_symmetric<double, 3> t_hess_sdf_set{0., 0., 0., 0., 0., 0.};
378 auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
379
380 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
381 t_hess_sdf(i, j) = t_hess_sdf_set(i, j);
382 ++t_hess_sdf;
383 }
384 return m_hess_sdf;
385}
386
387
388
389template <int DIM, IntegrationType I, typename BoundaryEleOp>
391
392template <int DIM, IntegrationType I, typename BoundaryEleOp>
394
395template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
397
398template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
400
401template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
403
404template <typename T1, typename T2, int DIM1, int DIM2>
407 size_t nb_gauss_pts) {
408 MatrixDouble m_spatial_coords(nb_gauss_pts, 3);
409 m_spatial_coords.clear();
410 auto t_spatial_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
412 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
413 t_spatial_coords(i) = t_coords(i) + t_disp(i);
414 ++t_spatial_coords;
415 ++t_coords;
416 ++t_disp;
417 }
418 return m_spatial_coords;
419}
420
421template <typename T1, int DIM1>
423 size_t nb_gauss_pts) {
424 MatrixDouble m_normals_at_pts(3, nb_gauss_pts);
425 m_normals_at_pts.clear();
427 auto t_set_normal = getFTensor1FromMat<3>(m_normals_at_pts);
428 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
429 t_set_normal(i) = t_normal_at_pts(i) / t_normal_at_pts.l2();
430 ++t_set_normal;
431 ++t_normal_at_pts;
432 }
433 return m_normals_at_pts;
434}
435
436template <int DIM, typename BoundaryEleOp>
438 : public BoundaryEleOp {
440 boost::shared_ptr<CommonData> common_data_ptr, double scale = 1,
441 bool is_axisymmetric = false);
442 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
443
444private:
445 boost::shared_ptr<CommonData> commonDataPtr;
446 const double scaleTraction;
448};
449
450template <int DIM, typename BoundaryEleOp>
451struct OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp> : public BoundaryEleOp {
452 OpEvaluateSDFImpl(boost::shared_ptr<CommonData> common_data_ptr);
453 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
454
455private:
456 boost::shared_ptr<CommonData> commonDataPtr;
457
459 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
461 HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
463};
464
465template <int DIM, typename AssemblyBoundaryEleOp>
467 : public AssemblyBoundaryEleOp {
468 OpConstrainBoundaryRhsImpl(const std::string field_name,
469 boost::shared_ptr<CommonData> common_data_ptr,
470 bool is_axisymmetric = false);
471 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
472
474 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
476
477private:
478 boost::shared_ptr<CommonData> commonDataPtr;
480};
481
482template <int DIM, typename AssemblyBoundaryEleOp>
484 : public AssemblyBoundaryEleOp {
485 OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
486 const std::string col_field_name,
487 boost::shared_ptr<CommonData> common_data_ptr,
488 bool is_axisymmetric = false);
489 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
490 EntitiesFieldData::EntData &col_data);
491
493 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
495 HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
497
498 boost::shared_ptr<CommonData> commonDataPtr;
500};
501
502template <int DIM, typename AssemblyBoundaryEleOp>
504 : public AssemblyBoundaryEleOp {
506 const std::string row_field_name, const std::string col_field_name,
507 boost::shared_ptr<CommonData> common_data_ptr,
508 bool is_axisymmetric = false);
509 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
510 EntitiesFieldData::EntData &col_data);
511
513 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
515
516private:
517 boost::shared_ptr<CommonData> commonDataPtr;
519};
520
521template <typename BoundaryEleOp> struct ContactIntegrators {
522
523 template <int DIM, IntegrationType I>
526
527 template <int DIM, IntegrationType I>
529
530 template <AssemblyType A> struct Assembly {
531
533 typename FormsIntegrators<BoundaryEleOp>::template Assembly<A>::OpBase;
534
535 template <int DIM, IntegrationType I>
538
539 template <int DIM, IntegrationType I>
542
543 template <int DIM, IntegrationType I>
546 };
547};
548
549inline double sign(double x) {
550 constexpr auto eps = std::numeric_limits<float>::epsilon();
551 if (std::abs(x) < eps)
552 return 0;
553 else if (x > eps)
554 return 1;
555 else
556 return -1;
557};
558
559inline double w(const double sdf, const double tn) {
560 return sdf - cn_contact * tn;
561}
562
563/**
564 * @brief constrain function
565 *
566 * return 1 if negative sdf or positive tn
567 *
568 * @param sdf signed distance
569 * @param tn traction
570 * @return double
571 */
572inline double constrain(double sdf, double tn) {
573 const auto s = sign(w(sdf, tn));
574 return (1 - s) / 2;
575}
576
577template <int DIM, typename BoundaryEleOp>
580 boost::shared_ptr<CommonData> common_data_ptr, double scale,
581 bool is_axisymmetric)
583 commonDataPtr(common_data_ptr), scaleTraction(scale),
584 isAxisymmetric(is_axisymmetric) {}
585
586template <int DIM, typename BoundaryEleOp>
587MoFEMErrorCode
589 int side, EntityType type, EntData &data) {
591
593 FTensor::Tensor1<double, 3> t_sum_t{0., 0., 0.};
594
596 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
598
599 const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
600 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
601 double jacobian = 1.;
602 if (isAxisymmetric) {
603 jacobian = 2. * M_PI * t_coords(0);
604 }
605 const double alpha = t_w * jacobian * BoundaryEleOp::getMeasure();
606 t_sum_t(i) += alpha * t_traction(i);
607 ++t_w;
608 ++t_traction;
609 ++t_coords;
610 }
611
612 t_sum_t(i) *= scaleTraction;
613
614 constexpr int ind[] = {0, 1, 2};
615 CHKERR VecSetValues(commonDataPtr->totalTraction, 3, ind, &t_sum_t(0),
616 ADD_VALUES);
617
619}
620
621template <int DIM, typename BoundaryEleOp>
623 boost::shared_ptr<CommonData> common_data_ptr)
625 commonDataPtr(common_data_ptr) {}
626
627template <int DIM, typename BoundaryEleOp>
628MoFEMErrorCode
630 EntData &data) {
632
633 const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
634 auto &sdf_vec = commonDataPtr->sdfVals;
635 auto &grad_mat = commonDataPtr->gradsSdf;
636 auto &hess_mat = commonDataPtr->hessSdf;
637 auto &constraint_vec = commonDataPtr->constraintVals;
638 auto &contactTraction_mat = commonDataPtr->contactTraction;
639
640 sdf_vec.resize(nb_gauss_pts, false);
641 grad_mat.resize(DIM, nb_gauss_pts, false);
642 hess_mat.resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
643 constraint_vec.resize(nb_gauss_pts, false);
644
645 auto t_traction = getFTensor1FromMat<DIM>(contactTraction_mat);
646
647 auto t_sdf = getFTensor0FromVec(sdf_vec);
648 auto t_grad_sdf = getFTensor1FromMat<DIM>(grad_mat);
649 auto t_hess_sdf = getFTensor2SymmetricFromMat<DIM>(hess_mat);
650 auto t_constraint = getFTensor0FromVec(constraint_vec);
651
652 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
654 auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
655
658
659 auto ts_time = BoundaryEleOp::getTStime();
660 auto ts_time_step = BoundaryEleOp::getTStimeStep();
661
662 auto m_spatial_coords = get_spatial_coords(
664 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
665 auto m_normals_at_pts = get_normalize_normals(
667
668 // placeholder to pass boundary block id to python
669 int block_id = 0;
670
671 auto v_sdf =
672 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
673 m_spatial_coords, m_normals_at_pts, block_id);
674
675 auto m_grad_sdf =
676 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
677 m_spatial_coords, m_normals_at_pts, block_id);
678
679 auto m_hess_sdf =
680 hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
681 m_spatial_coords, m_normals_at_pts, block_id);
682
683 auto t_sdf_v = getFTensor0FromVec(v_sdf);
684 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
685 auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
686
687 auto next = [&]() {
688 ++t_sdf;
689 ++t_sdf_v;
690 ++t_grad_sdf;
691 ++t_grad_sdf_v;
692 ++t_hess_sdf;
693 ++t_hess_sdf_v;
694 ++t_disp;
695 ++t_traction;
696 ++t_constraint;
697 };
698
699 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
700
701 auto tn = -t_traction(i) * t_grad_sdf_v(i);
702 auto c = constrain(t_sdf_v, tn);
703
704 t_sdf = t_sdf_v;
705 t_grad_sdf(i) = t_grad_sdf_v(i);
706 t_hess_sdf(i, j) = t_hess_sdf_v(i, j);
707 t_constraint = c;
708
709 next();
710 }
711
713}
714
715template <int DIM, typename AssemblyBoundaryEleOp>
718 boost::shared_ptr<CommonData> common_data_ptr,
719 bool is_axisymmetric)
721 AssemblyBoundaryEleOp::OPROW),
722 commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {}
723
724template <int DIM, typename AssemblyBoundaryEleOp>
725MoFEMErrorCode
727 EntitiesFieldData::EntData &data) {
729
734
735 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
736
738
739 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
740
741 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
742 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
743 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
744 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
745
746 size_t nb_base_functions = data.getN().size2() / 3;
747 auto t_base = data.getFTensor1N<3>();
748
749 auto m_spatial_coords = get_spatial_coords(
751 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
752 auto m_normals_at_pts = get_normalize_normals(
754
755 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
756
757 auto ts_time = AssemblyBoundaryEleOp::getTStime();
758 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
759
760 // placeholder to pass boundary block id to python
761 int block_id = 0;
762
763 auto v_sdf =
764 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
765 m_spatial_coords, m_normals_at_pts, block_id);
766
767 auto m_grad_sdf =
768 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
769 m_spatial_coords, m_normals_at_pts, block_id);
770
771 auto t_sdf = getFTensor0FromVec(v_sdf);
772 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
773
774 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
775
776 auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
777
778 double jacobian = 1.;
779 if (isAxisymmetric) {
780 jacobian = 2. * M_PI * t_coords(0);
781 }
782 const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
783
784 auto tn = -t_traction(i) * t_grad_sdf(i);
785 auto c = constrain(t_sdf, tn);
786
788 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
790 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
791
793 t_rhs(i) =
794
795 t_cQ(i, j) * (t_disp(j) - cn_contact * t_traction(j))
796
797 +
798
799 t_cP(i, j) * t_disp(j) +
800 c * (t_sdf * t_grad_sdf(i)); // add gap0 displacements
801
802 size_t bb = 0;
803 for (; bb != AssemblyBoundaryEleOp::nbRows / DIM; ++bb) {
804 const double beta = alpha * (t_base(i) * t_normal(i));
805 t_nf(i) -= beta * t_rhs(i);
806
807 ++t_nf;
808 ++t_base;
809 }
810 for (; bb < nb_base_functions; ++bb)
811 ++t_base;
812
813 ++t_disp;
814 ++t_traction;
815 ++t_coords;
816 ++t_w;
817 ++t_normal;
818 ++t_sdf;
819 ++t_grad_sdf;
820 }
821
823}
824
825template <int DIM, typename AssemblyBoundaryEleOp>
827 OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
828 const std::string col_field_name,
829 boost::shared_ptr<CommonData> common_data_ptr,
830 bool is_axisymmetric)
831 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
832 AssemblyBoundaryEleOp::OPROWCOL),
833 commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
834 AssemblyBoundaryEleOp::sYmm = false;
835}
836
837template <int DIM, typename AssemblyBoundaryEleOp>
838MoFEMErrorCode
840 EntitiesFieldData::EntData &row_data,
841 EntitiesFieldData::EntData &col_data) {
843
847
848 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
849 auto &locMat = AssemblyBoundaryEleOp::locMat;
850
851 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
852 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
853 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
854
855 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
856 auto t_row_base = row_data.getFTensor1N<3>();
857 size_t nb_face_functions = row_data.getN().size2() / 3;
858
859 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
860
861 auto m_spatial_coords = get_spatial_coords(
863 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
864 auto m_normals_at_pts = get_normalize_normals(
866
867 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
868
869 auto ts_time = AssemblyBoundaryEleOp::getTStime();
870 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
871
872 // placeholder to pass boundary block id to python
873 int block_id = 0;
874
875 auto v_sdf =
876 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
877 m_spatial_coords, m_normals_at_pts, block_id);
878
879 auto m_grad_sdf =
880 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
881 m_spatial_coords, m_normals_at_pts, block_id);
882
883 auto m_hess_sdf =
884 hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
885 m_spatial_coords, m_normals_at_pts, block_id);
886
887 auto t_sdf = getFTensor0FromVec(v_sdf);
888 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
889 auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
890
891 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
892
893 double jacobian = 1.;
894 if (isAxisymmetric) {
895 jacobian = 2. * M_PI * t_coords(0);
896 }
897 const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
898
899 auto tn = -t_traction(i) * t_grad_sdf(i);
900 auto c = constrain(t_sdf, tn);
901
903 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
905 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
906
908 t_res_dU(i, j) = kronecker_delta(i, j) + t_cP(i, j);
909
910 if (c > 0) {
911 t_res_dU(i, j) +=
912 (c * cn_contact) *
913 (t_hess_sdf(i, j) * (t_grad_sdf(k) * t_traction(k)) +
914 t_grad_sdf(i) * t_hess_sdf(k, j) * t_traction(k)) +
915 c * t_sdf * t_hess_sdf(i, j);
916 }
917
918 size_t rr = 0;
919 for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
920
921 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
922
923 const double row_base = t_row_base(i) * t_normal(i);
924
925 auto t_col_base = col_data.getFTensor0N(gg, 0);
926 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
927 const double beta = alpha * row_base * t_col_base;
928
929 t_mat(i, j) -= beta * t_res_dU(i, j);
930
931 ++t_col_base;
932 ++t_mat;
933 }
934
935 ++t_row_base;
936 }
937 for (; rr < nb_face_functions; ++rr)
938 ++t_row_base;
939
940 ++t_traction;
941 ++t_coords;
942 ++t_w;
943 ++t_normal;
944 ++t_sdf;
945 ++t_grad_sdf;
946 ++t_hess_sdf;
947 }
948
950}
951
952template <int DIM, typename AssemblyBoundaryEleOp>
955 const std::string row_field_name, const std::string col_field_name,
956 boost::shared_ptr<CommonData> common_data_ptr, bool is_axisymmetric)
957 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
958 AssemblyBoundaryEleOp::OPROWCOL),
959 commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
960 AssemblyBoundaryEleOp::sYmm = false;
961}
962
963template <int DIM, typename AssemblyBoundaryEleOp>
964MoFEMErrorCode
966 iNtegrate(EntitiesFieldData::EntData &row_data,
967 EntitiesFieldData::EntData &col_data) {
969
973
974 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
975 auto &locMat = AssemblyBoundaryEleOp::locMat;
976
977 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
978 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
979 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
980
981 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
982 auto t_row_base = row_data.getFTensor1N<3>();
983 size_t nb_face_functions = row_data.getN().size2() / 3;
984
985 auto m_spatial_coords = get_spatial_coords(
987 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
988 auto m_normals_at_pts = get_normalize_normals(
990
991 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
992
993 auto ts_time = AssemblyBoundaryEleOp::getTStime();
994 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
995
996 // placeholder to pass boundary block id to python
997 int block_id = 0;
998
999 auto v_sdf =
1000 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1001 m_spatial_coords, m_normals_at_pts, block_id);
1002
1003 auto m_grad_sdf =
1004 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1005 m_spatial_coords, m_normals_at_pts, block_id);
1006
1007 auto t_sdf = getFTensor0FromVec(v_sdf);
1008 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
1009
1010 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1011
1012 double jacobian = 1.;
1013 if (isAxisymmetric) {
1014 jacobian = 2. * M_PI * t_coords(0);
1015 }
1016 const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
1017
1018 auto tn = -t_traction(i) * t_grad_sdf(i);
1019 auto c = constrain(t_sdf, tn);
1020
1022 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
1024 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
1025
1027 t_res_dt(i, j) = -cn_contact * t_cQ(i, j);
1028
1029 size_t rr = 0;
1030 for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
1031
1032 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
1033 const double row_base = t_row_base(i) * t_normal(i);
1034
1035 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
1036 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
1037 const double col_base = t_col_base(i) * t_normal(i);
1038 const double beta = alpha * row_base * col_base;
1039
1040 t_mat(i, j) -= beta * t_res_dt(i, j);
1041
1042 ++t_col_base;
1043 ++t_mat;
1044 }
1045
1046 ++t_row_base;
1047 }
1048 for (; rr < nb_face_functions; ++rr)
1049 ++t_row_base;
1050
1051 ++t_traction;
1052 ++t_coords;
1053 ++t_w;
1054 ++t_normal;
1055 ++t_sdf;
1056 ++t_grad_sdf;
1057 }
1058
1060}
1061
1062template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
1063MoFEMErrorCode opFactoryDomainRhs(
1064 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1065 std::string sigma, std::string u, bool is_axisymmetric = false) {
1067
1068 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
1069 A>::template LinearForm<I>;
1070 using OpMixDivURhs = typename B::template OpMixDivTimesU<3, DIM, DIM>;
1071 using OpMixDivUCylRhs =
1072 typename B::template OpMixDivTimesU<3, DIM, DIM, CYLINDRICAL>;
1073
1074 using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
1076 typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
1077 using OpMixUTimesLambdaRhs =
1078 typename B::template OpGradTimesTensor<1, DIM, DIM>;
1079
1080 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1081 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
1082 auto div_stress_ptr = boost::make_shared<MatrixDouble>();
1083 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
1084
1085 auto jacobian = [is_axisymmetric](const double r, const double,
1086 const double) {
1087 if (is_axisymmetric)
1088 return 2. * M_PI * r;
1089 else
1090 return 1.;
1091 };
1092
1093 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1094 u, common_data_ptr->contactDispPtr()));
1095 pip.push_back(
1096 new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
1097
1098 if (!is_axisymmetric) {
1099 pip.push_back(
1100 new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
1101 } else {
1102 pip.push_back(new OpCalculateHVecTensorDivergence<DIM, DIM, CYLINDRICAL>(
1103 sigma, div_stress_ptr));
1104 }
1105
1106 pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(u, mat_grad_ptr));
1107
1108 if (!is_axisymmetric) {
1109 pip.push_back(
1110 new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(), jacobian));
1111 } else {
1112 pip.push_back(new OpMixDivUCylRhs(sigma, common_data_ptr->contactDispPtr(),
1113 jacobian));
1114 }
1115
1116 pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr, jacobian));
1117 pip.push_back(new OpMixUTimesDivLambdaRhs(u, div_stress_ptr, jacobian));
1118 pip.push_back(new OpMixUTimesLambdaRhs(u, contact_stress_ptr, jacobian));
1119
1121}
1122
1123template <typename OpMixLhs> struct OpMixLhsSide : public OpMixLhs {
1124 using OpMixLhs::OpMixLhs;
1125 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
1126 EntityType col_type,
1127 EntitiesFieldData::EntData &row_data,
1128 EntitiesFieldData::EntData &col_data) {
1130 auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
1131 auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
1132 // Only assemble side which correspond to edge entity on boundary
1133 if (side_fe_entity == side_fe_data) {
1134 CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
1135 col_data);
1136 }
1138 }
1139};
1140
1141template <int DIM, AssemblyType A, IntegrationType I, typename DomainEle>
1143 MoFEM::Interface &m_field,
1144 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1145 std::string fe_domain_name, std::string sigma, std::string u,
1146 std::string geom, ForcesAndSourcesCore::RuleHookFun rule,
1147 bool is_axisymmetric = false) {
1149
1150 using DomainEleOp = typename DomainEle::UserDataOperator;
1151
1152 auto op_loop_side = new OpLoopSide<DomainEle>(
1153 m_field, fe_domain_name, DIM, Sev::noisy,
1154 boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
1155 pip.push_back(op_loop_side);
1156
1157 CHKERR AddHOOps<DIM, DIM, DIM>::add(op_loop_side->getOpPtrVector(),
1158 {H1, HDIV}, geom);
1159
1160 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
1161 A>::template BiLinearForm<I>;
1162
1163 using OpMixDivULhs = typename B::template OpMixDivTimesVec<DIM>;
1164 using OpMixDivUCylLhs =
1165 typename B::template OpMixDivTimesVec<DIM, CYLINDRICAL>;
1166 using OpLambdaGraULhs = typename B::template OpMixTensorTimesGrad<DIM>;
1167
1168 using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
1169 using OpMixDivUCylLhsSide = OpMixLhsSide<OpMixDivUCylLhs>;
1170 using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
1171
1172 auto unity = []() { return 1; };
1173 auto jacobian = [is_axisymmetric](const double r, const double,
1174 const double) {
1175 if (is_axisymmetric)
1176 return 2. * M_PI * r;
1177 else
1178 return 1.;
1179 };
1180
1181 if (!is_axisymmetric) {
1182 op_loop_side->getOpPtrVector().push_back(
1183 new OpMixDivULhsSide(sigma, u, unity, jacobian, true));
1184 } else {
1185 op_loop_side->getOpPtrVector().push_back(
1186 new OpMixDivUCylLhsSide(sigma, u, unity, jacobian, true));
1187 }
1188 op_loop_side->getOpPtrVector().push_back(
1189 new OpLambdaGraULhsSide(sigma, u, unity, jacobian, true));
1190
1191 op_loop_side->getSideFEPtr()->getRuleHook = rule;
1193}
1194
1195template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
1197 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1198 std::string sigma, std::string u, bool is_axisymmetric = false) {
1200
1202
1203 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1204
1205 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1206 u, common_data_ptr->contactDispPtr()));
1207 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1208 sigma, common_data_ptr->contactTractionPtr()));
1209 pip.push_back(
1210 new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
1211 DIM, GAUSS>(sigma, u, common_data_ptr, is_axisymmetric));
1212 pip.push_back(new typename C::template Assembly<A>::
1213 template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
1214 sigma, sigma, common_data_ptr, is_axisymmetric));
1215
1217}
1218
1219template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
1221 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1222 std::string sigma, std::string u, bool is_axisymmetric = false) {
1224
1226
1227 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1228
1229 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1230 u, common_data_ptr->contactDispPtr()));
1231 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1232 sigma, common_data_ptr->contactTractionPtr()));
1233 pip.push_back(
1234 new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
1235 DIM, GAUSS>(sigma, common_data_ptr, is_axisymmetric));
1236
1238}
1239
1240template <int DIM, IntegrationType I, typename BoundaryEleOp>
1242 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1243 std::string sigma, bool is_axisymmetric = false) {
1245
1247
1248 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1249 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1250 sigma, common_data_ptr->contactTractionPtr()));
1251 pip.push_back(new typename C::template OpAssembleTotalContactTraction<DIM, I>(
1252 common_data_ptr, 1. / scale, is_axisymmetric));
1253
1255}
1256
1257}; // namespace ContactOps
1258
1259#endif // __CONTACTOPS_HPP__
static const double eps
Kronecker Delta class.
PetscBool is_axisymmetric
Definition: contact.cpp:137
#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
#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
MoFEMErrorCode opFactoryDomainRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
auto get_normalize_normals(FTensor::Tensor1< T1, DIM1 > &&t_normal_at_pts, size_t nb_gauss_pts)
Definition: ContactOps.hpp:422
MatrixDouble grad_surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:273
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, bool is_axisymmetric=false)
MatrixDouble hess_surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:330
MoFEMErrorCode opFactoryBoundaryRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
double w(const double sdf, const double tn)
Definition: ContactOps.hpp:559
boost::function< MatrixDouble(double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> HessSurfaceDistanceFunction
Definition: ContactOps.hpp:214
EntitiesFieldData::EntData EntData
auto get_spatial_coords(FTensor::Tensor1< T1, DIM1 > &&t_coords, FTensor::Tensor1< T2, DIM2 > &&t_disp, size_t nb_gauss_pts)
Definition: ContactOps.hpp:405
double constrain(double sdf, double tn)
constrain function
Definition: ContactOps.hpp:572
VectorDouble surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:216
boost::function< VectorDouble(double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> SurfaceDistanceFunction
[Common data]
Definition: ContactOps.hpp:206
MoFEMErrorCode opFactoryBoundaryLhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
boost::function< MatrixDouble(double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> GradSurfaceDistanceFunction
Definition: ContactOps.hpp:210
MoFEMErrorCode opFactoryCalculateTraction(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, bool is_axisymmetric=false)
double sign(double x)
Definition: ContactOps.hpp:549
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
VectorDouble sdfVals
size is equal to number of gauss points on element
Definition: ContactOps.hpp:19
MatrixDouble contactTraction
Definition: ContactOps.hpp:16
static auto getFTensor1TotalTraction()
Definition: ContactOps.hpp:41
VectorDouble constraintVals
Definition: ContactOps.hpp:25
MatrixDouble contactDisp
Definition: ContactOps.hpp:17
MatrixDouble hessSdf
Definition: ContactOps.hpp:22
MatrixDouble gradsSdf
Definition: ContactOps.hpp:20
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:28
static auto createTotalTraction(MoFEM::Interface &m_field)
Definition: ContactOps.hpp:30
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column