v0.15.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 MatrixDouble contactDispGrad;
19
20 VectorDouble sdfVals; ///< size is equal to number of gauss points on element
21 MatrixDouble gradsSdf; ///< nb of rows is equals to dimension, and nb of cols
22 ///< is equals to number of gauss points on element
23 MatrixDouble hessSdf; ///< nb of rows is equals to nb of element of symmetric
24 ///< matrix, and nb of cols is equals to number of gauss
25 ///< points on element
26 VectorDouble constraintVals;
27
28 static SmartPetscObj<Vec>
29 totalTraction; // User have to release and create vector when appropiate.
30
31 static auto createTotalTraction(MoFEM::Interface &m_field) {
32 constexpr int ghosts[] = {0, 1, 2, 3, 4};
34 createGhostVector(m_field.get_comm(),
35
36 (m_field.get_comm_rank() == 0) ? 5 : 0, 5,
37
38 (m_field.get_comm_rank() == 0) ? 0 : 5, ghosts);
39 return totalTraction;
40 }
41
44 const double *t_ptr;
45 CHK_THROW_MESSAGE(VecGetArrayRead(CommonData::totalTraction, &t_ptr),
46 "get array");
47 FTensor::Tensor1<double, 5> t{t_ptr[0], t_ptr[1], t_ptr[2], t_ptr[3],
48 t_ptr[4]};
49 CHK_THROW_MESSAGE(VecRestoreArrayRead(CommonData::totalTraction, &t_ptr),
50 "restore array");
51 return t;
52 } else {
53 return FTensor::Tensor1<double, 5>{0., 0., 0., 0., 0.};
54 }
55 }
56
57 inline auto contactTractionPtr() {
58 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
60 }
61
62 inline auto contactDispPtr() {
63 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactDisp);
64 }
65
66 inline auto contactDispGradPtr() {
67 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
69 }
70
71 inline auto sdfPtr() {
72 return boost::shared_ptr<VectorDouble>(shared_from_this(), &sdfVals);
73 }
74
75 inline auto gradSdfPtr() {
76 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &gradsSdf);
77 }
78
79 inline auto hessSdfPtr() {
80 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &hessSdf);
81 }
82
83 inline auto constraintPtr() {
84 return boost::shared_ptr<VectorDouble>(shared_from_this(), &constraintVals);
85 }
86};
87
88SmartPetscObj<Vec> CommonData::totalTraction;
89
90//! [Common data]
91
92//! [Surface distance function from python]
93#ifdef ENABLE_PYTHON_BINDING
94struct SDFPython {
95 SDFPython() = default;
96 virtual ~SDFPython() = default;
97
98 MoFEMErrorCode sdfInit(const std::string py_file) {
100 try {
101
102 // create main module
103 auto main_module = bp::import("__main__");
104 mainNamespace = main_module.attr("__dict__");
105 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
106 // create a reference to python function
107 sdfFun = mainNamespace["sdf"];
108 sdfGradFun = mainNamespace["grad_sdf"];
109 sdfHessFun = mainNamespace["hess_sdf"];
110
111 } catch (bp::error_already_set const &) {
112 // print all other errors to stderr
113 PyErr_Print();
115 }
117 };
118
119 template <typename T>
120 inline std::vector<T>
121 py_list_to_std_vector(const boost::python::object &iterable) {
122 return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
123 boost::python::stl_input_iterator<T>());
124 }
125
126 MoFEMErrorCode evalSdf(
127
128 double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
129 np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
130 np::ndarray &sdf
131
132 ) {
134 try {
135
136 // call python function
137 sdf = bp::extract<np::ndarray>(
138 sdfFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
139
140 } catch (bp::error_already_set const &) {
141 // print all other errors to stderr
142 PyErr_Print();
144 }
146 }
147
148 MoFEMErrorCode evalGradSdf(
149
150 double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
151 np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
152 np::ndarray &grad_sdf
153
154 ) {
156 try {
157
158 // call python function
159 grad_sdf = bp::extract<np::ndarray>(
160 sdfGradFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
161
162 } catch (bp::error_already_set const &) {
163 // print all other errors to stderr
164 PyErr_Print();
166 }
168 }
169
170 MoFEMErrorCode evalHessSdf(
171
172 double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
173 np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
174 np::ndarray &hess_sdf
175
176 ) {
178 try {
179
180 // call python function
181 hess_sdf = bp::extract<np::ndarray>(
182 sdfHessFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
183
184 } catch (bp::error_already_set const &) {
185 // print all other errors to stderr
186 PyErr_Print();
188 }
190 }
191
192private:
193 bp::object mainNamespace;
194 bp::object sdfFun;
195 bp::object sdfGradFun;
196 bp::object sdfHessFun;
197};
198
199static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
200
201inline np::ndarray convert_to_numpy(VectorDouble &data, int nb_gauss_pts,
202 int id) {
203 auto dtype = np::dtype::get_builtin<double>();
204 auto size = bp::make_tuple(nb_gauss_pts);
205 auto stride = bp::make_tuple(3 * sizeof(double));
206 return (np::from_data(&data[id], dtype, size, stride, bp::object()));
207};
208#endif
209//! [Surface distance function from python]
210
211using SurfaceDistanceFunction = boost::function<VectorDouble(
212 double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
213 MatrixDouble &normals_at_pts, int block_id)>;
214
215using GradSurfaceDistanceFunction = boost::function<MatrixDouble(
216 double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
217 MatrixDouble &normals_at_pts, int block_id)>;
218
219using HessSurfaceDistanceFunction = boost::function<MatrixDouble(
220 double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
221 MatrixDouble &normals_at_pts, int block_id)>;
222
223inline VectorDouble surface_distance_function(double delta_t, double t,
224 int nb_gauss_pts,
225 MatrixDouble &m_spatial_coords,
226 MatrixDouble &m_normals_at_pts,
227 int block_id) {
228
229#ifdef ENABLE_PYTHON_BINDING
230 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
231
232 VectorDouble v_spatial_coords = m_spatial_coords.data();
233 VectorDouble v_normal_at_pts = m_normals_at_pts.data();
234
235 bp::list python_coords;
236 bp::list python_normals;
237
238 for (int idx = 0; idx < 3; ++idx) {
239 python_coords.append(
240 convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
241 python_normals.append(
242 convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
243 }
244
245 np::ndarray np_sdf = np::empty(bp::make_tuple(nb_gauss_pts),
246 np::dtype::get_builtin<double>());
247 CHK_MOAB_THROW(sdf_ptr->evalSdf(delta_t, t,
248 bp::extract<np::ndarray>(python_coords[0]),
249 bp::extract<np::ndarray>(python_coords[1]),
250 bp::extract<np::ndarray>(python_coords[2]),
251 bp::extract<np::ndarray>(python_normals[0]),
252 bp::extract<np::ndarray>(python_normals[1]),
253 bp::extract<np::ndarray>(python_normals[2]),
254 block_id, np_sdf),
255 "Failed python call");
256
257 double *sdf_val_ptr = reinterpret_cast<double *>(np_sdf.get_data());
258
259 VectorDouble v_sdf;
260 v_sdf.resize(nb_gauss_pts, false);
261
262 for (size_t gg = 0; gg < nb_gauss_pts; ++gg)
263 v_sdf[gg] = *(sdf_val_ptr + gg);
264
265 return v_sdf;
266 }
267#endif
268 VectorDouble v_sdf;
269 v_sdf.resize(nb_gauss_pts, false);
270 auto t_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
271
272 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
273 v_sdf[gg] = -t_coords(2) - 0.1;
274 ++t_coords;
275 }
276
277 return v_sdf;
278}
279
280inline MatrixDouble
281grad_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
282 MatrixDouble &m_spatial_coords,
283 MatrixDouble &m_normals_at_pts, int block_id) {
284#ifdef ENABLE_PYTHON_BINDING
285 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
286
287 VectorDouble v_spatial_coords = m_spatial_coords.data();
288 VectorDouble v_normal_at_pts = m_normals_at_pts.data();
289
290 bp::list python_coords;
291 bp::list python_normals;
292
293 for (int idx = 0; idx < 3; ++idx) {
294 python_coords.append(
295 convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
296 python_normals.append(
297 convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
298 }
299
300 np::ndarray np_grad_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 3),
301 np::dtype::get_builtin<double>());
302 CHK_MOAB_THROW(sdf_ptr->evalGradSdf(
303 delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
304 bp::extract<np::ndarray>(python_coords[1]),
305 bp::extract<np::ndarray>(python_coords[2]),
306 bp::extract<np::ndarray>(python_normals[0]),
307 bp::extract<np::ndarray>(python_normals[1]),
308 bp::extract<np::ndarray>(python_normals[2]), block_id,
309 np_grad_sdf),
310 "Failed python call");
311
312 double *grad_ptr = reinterpret_cast<double *>(np_grad_sdf.get_data());
313
314 MatrixDouble m_grad_sdf;
315 m_grad_sdf.resize(3, nb_gauss_pts, false);
316 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
317 for (int idx = 0; idx < 3; ++idx)
318 m_grad_sdf(idx, gg) = *(grad_ptr + (3 * gg + idx));
319 }
320 return m_grad_sdf;
321 }
322#endif
323 MatrixDouble m_grad_sdf;
324 m_grad_sdf.resize(3, nb_gauss_pts, false);
325 FTensor::Index<'i', 3> i;
326 FTensor::Tensor1<double, 3> t_grad_sdf_set{0.0, 0.0, -1.0};
327 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
328
329 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
330 t_grad_sdf(i) = t_grad_sdf_set(i);
331 ++t_grad_sdf;
332 }
333
334 return m_grad_sdf;
335}
336
337inline MatrixDouble
338hess_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
339 MatrixDouble &m_spatial_coords,
340 MatrixDouble &m_normals_at_pts, int block_id) {
341#ifdef ENABLE_PYTHON_BINDING
342 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
343
344 VectorDouble v_spatial_coords = m_spatial_coords.data();
345 VectorDouble v_normal_at_pts = m_normals_at_pts.data();
346
347 bp::list python_coords;
348 bp::list python_normals;
349
350 for (int idx = 0; idx < 3; ++idx) {
351 python_coords.append(
352 convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
353 python_normals.append(
354 convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
355 };
356
357 np::ndarray np_hess_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 6),
358 np::dtype::get_builtin<double>());
359 CHK_MOAB_THROW(sdf_ptr->evalHessSdf(
360 delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
361 bp::extract<np::ndarray>(python_coords[1]),
362 bp::extract<np::ndarray>(python_coords[2]),
363 bp::extract<np::ndarray>(python_normals[0]),
364 bp::extract<np::ndarray>(python_normals[1]),
365 bp::extract<np::ndarray>(python_normals[2]), block_id,
366 np_hess_sdf),
367 "Failed python call");
368
369 double *hess_ptr = reinterpret_cast<double *>(np_hess_sdf.get_data());
370
371 MatrixDouble m_hess_sdf;
372 m_hess_sdf.resize(6, nb_gauss_pts, false);
373 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
374 for (int idx = 0; idx < 6; ++idx)
375 m_hess_sdf(idx, gg) =
376 *(hess_ptr + (6 * gg + idx));
377 }
378 return m_hess_sdf;
379 }
380#endif
381 MatrixDouble m_hess_sdf;
382 m_hess_sdf.resize(6, nb_gauss_pts, false);
383 FTensor::Index<'i', 3> i;
384 FTensor::Index<'j', 3> j;
385 FTensor::Tensor2_symmetric<double, 3> t_hess_sdf_set{0., 0., 0., 0., 0., 0.};
386 auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
387
388 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
389 t_hess_sdf(i, j) = t_hess_sdf_set(i, j);
390 ++t_hess_sdf;
391 }
392 return m_hess_sdf;
393}
394
395template <int DIM, IntegrationType I, typename BoundaryEleOp>
397
398template <int DIM, IntegrationType I, typename BoundaryEleOp>
400
401template <int DIM, IntegrationType I, typename BoundaryEleOp>
403
404template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
406
407template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
409
410template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
412
413template <typename T1, typename T2, int DIM1, int DIM2>
416 size_t nb_gauss_pts) {
417 MatrixDouble m_spatial_coords(nb_gauss_pts, 3);
418 m_spatial_coords.clear();
419 auto t_spatial_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
420 FTensor::Index<'i', DIM2> i;
421 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
422 t_spatial_coords(i) = t_coords(i) + t_disp(i);
423 ++t_spatial_coords;
424 ++t_coords;
425 ++t_disp;
426 }
427 return m_spatial_coords;
428}
429
430template <typename T1, int DIM1>
432 size_t nb_gauss_pts) {
433 MatrixDouble m_normals_at_pts(3, nb_gauss_pts);
434 m_normals_at_pts.clear();
435 FTensor::Index<'i', DIM1> i;
436 auto t_set_normal = getFTensor1FromMat<3>(m_normals_at_pts);
437 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
438 t_set_normal(i) = t_normal_at_pts(i) / t_normal_at_pts.l2();
439 ++t_set_normal;
440 ++t_normal_at_pts;
441 }
442 return m_normals_at_pts;
443}
444
445template <int DIM, typename BoundaryEleOp>
447 : public BoundaryEleOp {
449 boost::shared_ptr<CommonData> common_data_ptr, double scale = 1,
450 bool is_axisymmetric = false);
451 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
452
453private:
454 boost::shared_ptr<CommonData> commonDataPtr;
455 const double scaleTraction;
457};
458
459template <int DIM, typename BoundaryEleOp>
461 : public BoundaryEleOp {
463 boost::shared_ptr<CommonData> common_data_ptr,
464 bool is_axisymmetric = false,
465 boost::shared_ptr<Range> contact_range_ptr = nullptr);
466 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
467
469 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
471
472private:
473 boost::shared_ptr<CommonData> commonDataPtr;
475 boost::shared_ptr<Range> contactRange;
476};
477
478template <int DIM, typename BoundaryEleOp>
479struct OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp> : public BoundaryEleOp {
480 OpEvaluateSDFImpl(boost::shared_ptr<CommonData> common_data_ptr);
481 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
482
483private:
484 boost::shared_ptr<CommonData> commonDataPtr;
485
487 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
489 HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
491};
492
493template <int DIM, typename AssemblyBoundaryEleOp>
495 : public AssemblyBoundaryEleOp {
496 OpConstrainBoundaryRhsImpl(const std::string field_name,
497 boost::shared_ptr<CommonData> common_data_ptr,
498 bool is_axisymmetric = false);
499 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
500
502 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
504
505private:
506 boost::shared_ptr<CommonData> commonDataPtr;
508};
509
510template <int DIM, typename AssemblyBoundaryEleOp>
512 : public AssemblyBoundaryEleOp {
513 OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
514 const std::string col_field_name,
515 boost::shared_ptr<CommonData> common_data_ptr,
516 bool is_axisymmetric = false);
517 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
518 EntitiesFieldData::EntData &col_data);
519
521 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
523 HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
525
526 boost::shared_ptr<CommonData> commonDataPtr;
528};
529
530template <int DIM, typename AssemblyBoundaryEleOp>
532 : public AssemblyBoundaryEleOp {
534 const std::string row_field_name, const std::string col_field_name,
535 boost::shared_ptr<CommonData> common_data_ptr,
536 bool is_axisymmetric = false);
537 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
538 EntitiesFieldData::EntData &col_data);
539
541 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
543
544private:
545 boost::shared_ptr<CommonData> commonDataPtr;
547};
548
549template <typename BoundaryEleOp> struct ContactIntegrators {
550 template <int DIM, IntegrationType I>
553
554 template <int DIM, IntegrationType I>
557
558 template <int DIM, IntegrationType I>
560
561 template <AssemblyType A> struct Assembly {
562
564 typename FormsIntegrators<BoundaryEleOp>::template Assembly<A>::OpBase;
565
566 template <int DIM, IntegrationType I>
569
570 template <int DIM, IntegrationType I>
573
574 template <int DIM, IntegrationType I>
577 };
578};
579
580inline double sign(double x) {
581 constexpr auto eps = std::numeric_limits<float>::epsilon();
582 if (std::abs(x) < eps)
583 return 0;
584 else if (x > eps)
585 return 1;
586 else
587 return -1;
588};
589
590inline double w(const double sdf, const double tn) {
591 return sdf - cn_contact * tn;
592}
593
594/**
595 * @brief constrain function
596 *
597 * return 1 if negative sdf or positive tn
598 *
599 * @param sdf signed distance
600 * @param tn traction
601 * @return double
602 */
603inline double constrain(double sdf, double tn) {
604 const auto s = sign(w(sdf, tn));
605 return (1 - s) / 2;
606}
607
608template <int DIM, typename BoundaryEleOp>
611 boost::shared_ptr<CommonData> common_data_ptr, double scale,
612 bool is_axisymmetric)
614 commonDataPtr(common_data_ptr), scaleTraction(scale),
615 isAxisymmetric(is_axisymmetric) {}
616
617template <int DIM, typename BoundaryEleOp>
618MoFEMErrorCode
620 int side, EntityType type, EntData &data) {
622
623 FTensor::Index<'i', DIM> i;
624 FTensor::Tensor1<double, 3> t_sum_t{0., 0., 0.};
625
626 auto t_w = BoundaryEleOp::getFTensor0IntegrationWeight();
627 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
628 auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
629
630 const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
631 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
632 double jacobian = 1.;
633 if (isAxisymmetric) {
634 jacobian = 2. * M_PI * t_coords(0);
635 }
636 const double alpha = t_w * jacobian * BoundaryEleOp::getMeasure();
637 t_sum_t(i) += alpha * t_traction(i);
638 ++t_w;
639 ++t_traction;
640 ++t_coords;
641 }
642
643 t_sum_t(i) *= scaleTraction;
644
645 constexpr int ind[] = {0, 1, 2};
646 CHKERR VecSetValues(commonDataPtr->totalTraction, 3, ind, &t_sum_t(0),
647 ADD_VALUES);
648
650}
651template <int DIM, typename BoundaryEleOp>
654 boost::shared_ptr<CommonData> common_data_ptr, bool is_axisymmetric,
655 boost::shared_ptr<Range> contact_range_ptr)
657 commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric),
658 contactRange(contact_range_ptr) {}
659
660template <int DIM, typename BoundaryEleOp>
661MoFEMErrorCode
663 int side, EntityType type, EntData &data) {
665
666 auto fe_type = BoundaryEleOp::getFEType();
667
668 const auto fe_ent = BoundaryEleOp::getFEEntityHandle();
669
670 if (contactRange->find(fe_ent) != contactRange->end()) {
671 FTensor::Index<'i', DIM> i;
672 FTensor::Index<'j', DIM> j;
673 FTensor::Tensor1<double, 2> t_sum_a{0., 0.};
674
675 auto t_w = BoundaryEleOp::getFTensor0IntegrationWeight();
676 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
677 auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
678
679 auto t_grad = getFTensor2FromMat<DIM, DIM>(commonDataPtr->contactDispGrad);
680 auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
681
682 const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
683 auto m_spatial_coords = get_spatial_coords(
684 BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
685 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
686 auto m_normals_at_pts = get_normalize_normals(
687 BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
688
689 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
690 auto ts_time = BoundaryEleOp::getTStime();
691 auto ts_time_step = BoundaryEleOp::getTStimeStep();
692 int block_id = 0;
693 auto v_sdf =
694 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
695 m_spatial_coords, m_normals_at_pts, block_id);
696 auto m_grad_sdf = gradSurfaceDistanceFunction(
697 ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
698 block_id);
699 auto t_sdf = getFTensor0FromVec(v_sdf);
700 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
701 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
702 double jacobian = 1.;
703 if (isAxisymmetric) {
704 jacobian = 2. * M_PI * t_coords(0); // Axisymmetric Jacobian
705 }
706 auto tn = -t_traction(i) * t_grad_sdf(i);
707 auto c = constrain(t_sdf, tn);
708 double alpha = t_w * jacobian;
709
712 FTensor::Tensor1<double, DIM> t_normal_current;
713
714 F(i, j) = t_grad(i, j) + kronecker_delta(i, j);
715 auto det = determinantTensor(F);
716 CHKERR invertTensor(F, det, invF);
717 t_normal_current(i) = det * (invF(j, i) * t_normal_at_pts(j));
718
719 alpha *= sqrt(t_normal_current(i) * t_normal_current(i));
720
721 if (fe_type == MBTRI) {
722 alpha /= 2;
723 }
724 if (c > 1e-12) {
725 t_sum_a(0) += alpha; // real area
726 }
727 t_sum_a(1) += alpha; // Potential area
728 ++t_w;
729 ++t_traction;
730 ++t_coords;
731 ++t_sdf;
732 ++t_grad_sdf;
733
734 ++t_grad;
735 ++t_normal_at_pts;
736 }
737 constexpr int ind[] = {3, 4};
738 CHKERR VecSetValues(commonDataPtr->totalTraction, 2, ind, &t_sum_a(0),
739 ADD_VALUES);
740 }
742}
743
744template <int DIM, typename BoundaryEleOp>
746 boost::shared_ptr<CommonData> common_data_ptr)
748 commonDataPtr(common_data_ptr) {}
749
750template <int DIM, typename BoundaryEleOp>
751MoFEMErrorCode
753 EntData &data) {
755
756 const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
757 auto &sdf_vec = commonDataPtr->sdfVals;
758 auto &grad_mat = commonDataPtr->gradsSdf;
759 auto &hess_mat = commonDataPtr->hessSdf;
760 auto &constraint_vec = commonDataPtr->constraintVals;
761 auto &contactTraction_mat = commonDataPtr->contactTraction;
762
763 sdf_vec.resize(nb_gauss_pts, false);
764 grad_mat.resize(DIM, nb_gauss_pts, false);
765 hess_mat.resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
766 constraint_vec.resize(nb_gauss_pts, false);
767
768 auto t_traction = getFTensor1FromMat<DIM>(contactTraction_mat);
769
770 auto t_sdf = getFTensor0FromVec(sdf_vec);
771 auto t_grad_sdf = getFTensor1FromMat<DIM>(grad_mat);
772 auto t_hess_sdf = getFTensor2SymmetricFromMat<DIM>(hess_mat);
773 auto t_constraint = getFTensor0FromVec(constraint_vec);
774
775 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
776 auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
777 auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
778
779 FTensor::Index<'i', DIM> i;
780 FTensor::Index<'j', DIM> j;
781
782 auto ts_time = BoundaryEleOp::getTStime();
783 auto ts_time_step = BoundaryEleOp::getTStimeStep();
784
785 auto m_spatial_coords = get_spatial_coords(
786 BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
787 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
788 auto m_normals_at_pts = get_normalize_normals(
789 BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
790
791 // placeholder to pass boundary block id to python
792 int block_id = 0;
793
794 auto v_sdf =
795 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
796 m_spatial_coords, m_normals_at_pts, block_id);
797
798 auto m_grad_sdf =
799 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
800 m_spatial_coords, m_normals_at_pts, block_id);
801
802 auto m_hess_sdf =
803 hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
804 m_spatial_coords, m_normals_at_pts, block_id);
805
806 auto t_sdf_v = getFTensor0FromVec(v_sdf);
807 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
808 auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
809
810 auto next = [&]() {
811 ++t_sdf;
812 ++t_sdf_v;
813 ++t_grad_sdf;
814 ++t_grad_sdf_v;
815 ++t_hess_sdf;
816 ++t_hess_sdf_v;
817 ++t_disp;
818 ++t_traction;
819 ++t_constraint;
820 };
821
822 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
823
824 auto tn = -t_traction(i) * t_grad_sdf_v(i);
825 auto c = constrain(t_sdf_v, tn);
826
827 t_sdf = t_sdf_v;
828 t_grad_sdf(i) = t_grad_sdf_v(i);
829 t_hess_sdf(i, j) = t_hess_sdf_v(i, j);
830 t_constraint = c;
831
832 next();
833 }
834
836}
837
838template <int DIM, typename AssemblyBoundaryEleOp>
841 boost::shared_ptr<CommonData> common_data_ptr,
842 bool is_axisymmetric)
844 AssemblyBoundaryEleOp::OPROW),
845 commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {}
846
847template <int DIM, typename AssemblyBoundaryEleOp>
848MoFEMErrorCode
850 EntitiesFieldData::EntData &data) {
852
853 FTensor::Index<'i', DIM> i;
854 FTensor::Index<'j', DIM> j;
855 FTensor::Index<'k', DIM> k;
856 FTensor::Index<'l', DIM> l;
857
858 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
859
860 auto &nf = AssemblyBoundaryEleOp::locF;
861
862 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
863
864 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
865 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
866 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
867 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
868
869 size_t nb_base_functions = data.getN().size2() / 3;
870 auto t_base = data.getFTensor1N<3>();
871
872 auto m_spatial_coords = get_spatial_coords(
873 BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
874 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
875 auto m_normals_at_pts = get_normalize_normals(
876 BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
877
878 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
879
880 auto ts_time = AssemblyBoundaryEleOp::getTStime();
881 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
882
883 // placeholder to pass boundary block id to python
884 int block_id = 0;
885
886 auto v_sdf =
887 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
888 m_spatial_coords, m_normals_at_pts, block_id);
889
890 auto m_grad_sdf =
891 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
892 m_spatial_coords, m_normals_at_pts, block_id);
893
894 auto t_sdf = getFTensor0FromVec(v_sdf);
895 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
896
897 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
898
899 auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
900
901 double jacobian = 1.;
902 if (isAxisymmetric) {
903 jacobian = 2. * M_PI * t_coords(0);
904 }
905 const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
906
907 auto tn = -t_traction(i) * t_grad_sdf(i);
908 auto c = constrain(t_sdf, tn);
909
911 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
913 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
914
916 t_rhs(i) =
917
918 t_cQ(i, j) * (t_disp(j) - cn_contact * t_traction(j))
919
920 +
921
922 t_cP(i, j) * t_disp(j) +
923 c * (t_sdf * t_grad_sdf(i)); // add gap0 displacements
924
925 size_t bb = 0;
926 for (; bb != AssemblyBoundaryEleOp::nbRows / DIM; ++bb) {
927 const double beta = alpha * (t_base(i) * t_normal(i));
928 t_nf(i) -= beta * t_rhs(i);
929
930 ++t_nf;
931 ++t_base;
932 }
933 for (; bb < nb_base_functions; ++bb)
934 ++t_base;
935
936 ++t_disp;
937 ++t_traction;
938 ++t_coords;
939 ++t_w;
940 ++t_normal;
941 ++t_sdf;
942 ++t_grad_sdf;
943 }
944
946}
947
948template <int DIM, typename AssemblyBoundaryEleOp>
950 OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
951 const std::string col_field_name,
952 boost::shared_ptr<CommonData> common_data_ptr,
953 bool is_axisymmetric)
954 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
955 AssemblyBoundaryEleOp::OPROWCOL),
956 commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
957 AssemblyBoundaryEleOp::sYmm = false;
958}
959
960template <int DIM, typename AssemblyBoundaryEleOp>
961MoFEMErrorCode
963 EntitiesFieldData::EntData &row_data,
964 EntitiesFieldData::EntData &col_data) {
966
967 FTensor::Index<'i', DIM> i;
968 FTensor::Index<'j', DIM> j;
969 FTensor::Index<'k', DIM> k;
970
971 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
972 auto &locMat = AssemblyBoundaryEleOp::locMat;
973
974 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
975 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
976 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
977
978 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
979 auto t_row_base = row_data.getFTensor1N<3>();
980 size_t nb_face_functions = row_data.getN().size2() / 3;
981
982 auto m_spatial_coords = get_spatial_coords(
983 BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
984 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
985 auto m_normals_at_pts = get_normalize_normals(
986 BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
987
988 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
989
990 auto ts_time = AssemblyBoundaryEleOp::getTStime();
991 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
992
993 // placeholder to pass boundary block id to python
994 int block_id = 0;
995
996 auto v_sdf =
997 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
998 m_spatial_coords, m_normals_at_pts, block_id);
999
1000 auto m_grad_sdf =
1001 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1002 m_spatial_coords, m_normals_at_pts, block_id);
1003
1004 auto m_hess_sdf =
1005 hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1006 m_spatial_coords, m_normals_at_pts, block_id);
1007
1008 auto t_sdf = getFTensor0FromVec(v_sdf);
1009 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
1010 auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
1011
1012 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1013
1014 double jacobian = 1.;
1015 if (isAxisymmetric) {
1016 jacobian = 2. * M_PI * t_coords(0);
1017 }
1018 const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
1019
1020 auto tn = -t_traction(i) * t_grad_sdf(i);
1021 auto c = constrain(t_sdf, tn);
1022
1024 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
1026 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
1027
1029 t_res_dU(i, j) = kronecker_delta(i, j) + t_cP(i, j);
1030
1031 if (c > 0) {
1032 t_res_dU(i, j) +=
1033 (c * cn_contact) *
1034 (t_hess_sdf(i, j) * (t_grad_sdf(k) * t_traction(k)) +
1035 t_grad_sdf(i) * t_hess_sdf(k, j) * t_traction(k)) +
1036 c * t_sdf * t_hess_sdf(i, j);
1037 }
1038
1039 size_t rr = 0;
1040 for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
1041
1042 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
1043
1044 const double row_base = t_row_base(i) * t_normal(i);
1045
1046 auto t_col_base = col_data.getFTensor0N(gg, 0);
1047 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
1048 const double beta = alpha * row_base * t_col_base;
1049
1050 t_mat(i, j) -= beta * t_res_dU(i, j);
1051
1052 ++t_col_base;
1053 ++t_mat;
1054 }
1055
1056 ++t_row_base;
1057 }
1058 for (; rr < nb_face_functions; ++rr)
1059 ++t_row_base;
1060
1061 ++t_traction;
1062 ++t_coords;
1063 ++t_w;
1064 ++t_normal;
1065 ++t_sdf;
1066 ++t_grad_sdf;
1067 ++t_hess_sdf;
1068 }
1069
1071}
1072
1073template <int DIM, typename AssemblyBoundaryEleOp>
1076 const std::string row_field_name, const std::string col_field_name,
1077 boost::shared_ptr<CommonData> common_data_ptr, bool is_axisymmetric)
1078 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
1079 AssemblyBoundaryEleOp::OPROWCOL),
1080 commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
1081 AssemblyBoundaryEleOp::sYmm = false;
1082}
1083
1084template <int DIM, typename AssemblyBoundaryEleOp>
1085MoFEMErrorCode
1087 iNtegrate(EntitiesFieldData::EntData &row_data,
1088 EntitiesFieldData::EntData &col_data) {
1090
1091 FTensor::Index<'i', DIM> i;
1092 FTensor::Index<'j', DIM> j;
1093 FTensor::Index<'k', DIM> k;
1094
1095 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
1096 auto &locMat = AssemblyBoundaryEleOp::locMat;
1097
1098 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
1099 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
1100 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
1101
1102 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
1103 auto t_row_base = row_data.getFTensor1N<3>();
1104 size_t nb_face_functions = row_data.getN().size2() / 3;
1105
1106 auto m_spatial_coords = get_spatial_coords(
1107 BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
1108 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
1109 auto m_normals_at_pts = get_normalize_normals(
1110 BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
1111
1112 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
1113
1114 auto ts_time = AssemblyBoundaryEleOp::getTStime();
1115 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
1116
1117 // placeholder to pass boundary block id to python
1118 int block_id = 0;
1119
1120 auto v_sdf =
1121 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1122 m_spatial_coords, m_normals_at_pts, block_id);
1123
1124 auto m_grad_sdf =
1125 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1126 m_spatial_coords, m_normals_at_pts, block_id);
1127
1128 auto t_sdf = getFTensor0FromVec(v_sdf);
1129 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
1130
1131 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1132
1133 double jacobian = 1.;
1134 if (isAxisymmetric) {
1135 jacobian = 2. * M_PI * t_coords(0);
1136 }
1137 const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
1138
1139 auto tn = -t_traction(i) * t_grad_sdf(i);
1140 auto c = constrain(t_sdf, tn);
1141
1143 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
1145 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
1146
1148 t_res_dt(i, j) = -cn_contact * t_cQ(i, j);
1149
1150 size_t rr = 0;
1151 for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
1152
1153 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
1154 const double row_base = t_row_base(i) * t_normal(i);
1155
1156 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
1157 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
1158 const double col_base = t_col_base(i) * t_normal(i);
1159 const double beta = alpha * row_base * col_base;
1160
1161 t_mat(i, j) -= beta * t_res_dt(i, j);
1162
1163 ++t_col_base;
1164 ++t_mat;
1165 }
1166
1167 ++t_row_base;
1168 }
1169 for (; rr < nb_face_functions; ++rr)
1170 ++t_row_base;
1171
1172 ++t_traction;
1173 ++t_coords;
1174 ++t_w;
1175 ++t_normal;
1176 ++t_sdf;
1177 ++t_grad_sdf;
1178 }
1179
1181}
1182
1183template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
1184MoFEMErrorCode opFactoryDomainRhs(
1185 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1186 std::string sigma, std::string u, bool is_axisymmetric = false) {
1188
1189 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
1190 A>::template LinearForm<I>;
1191 using OpMixDivURhs = typename B::template OpMixDivTimesU<3, DIM, DIM>;
1192 using OpMixDivUCylRhs =
1193 typename B::template OpMixDivTimesU<3, DIM, DIM, CYLINDRICAL>;
1194
1195 using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
1196 using OpMixUTimesDivLambdaRhs =
1197 typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
1198 using OpMixUTimesLambdaRhs =
1199 typename B::template OpGradTimesTensor<1, DIM, DIM>;
1200
1201 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1202 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
1203 auto div_stress_ptr = boost::make_shared<MatrixDouble>();
1204 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
1205
1206 auto jacobian = [is_axisymmetric](const double r, const double,
1207 const double) {
1208 if (is_axisymmetric)
1209 return 2. * M_PI * r;
1210 else
1211 return 1.;
1212 };
1213
1214 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1215 u, common_data_ptr->contactDispPtr()));
1216 pip.push_back(
1217 new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
1218
1219 if (!is_axisymmetric) {
1220 pip.push_back(
1221 new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
1222 } else {
1223 pip.push_back(new OpCalculateHVecTensorDivergence<DIM, DIM, CYLINDRICAL>(
1224 sigma, div_stress_ptr));
1225 }
1226
1227 pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(u, mat_grad_ptr));
1228
1229 if (!is_axisymmetric) {
1230 pip.push_back(
1231 new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(), jacobian));
1232 } else {
1233 pip.push_back(new OpMixDivUCylRhs(sigma, common_data_ptr->contactDispPtr(),
1234 jacobian));
1235 }
1236
1237 pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr, jacobian));
1238 pip.push_back(new OpMixUTimesDivLambdaRhs(u, div_stress_ptr, jacobian));
1239 pip.push_back(new OpMixUTimesLambdaRhs(u, contact_stress_ptr, jacobian));
1240
1242}
1243
1244template <typename OpMixLhs> struct OpMixLhsSide : public OpMixLhs {
1245 using OpMixLhs::OpMixLhs;
1246 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
1247 EntityType col_type,
1248 EntitiesFieldData::EntData &row_data,
1249 EntitiesFieldData::EntData &col_data) {
1251 auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
1252 auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
1253 // Only assemble side which correspond to edge entity on boundary
1254 if (side_fe_entity == side_fe_data) {
1255 CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
1256 col_data);
1257 }
1259 }
1260};
1261
1262template <int DIM, AssemblyType A, IntegrationType I, typename DomainEle>
1264 MoFEM::Interface &m_field,
1265 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1266 std::string fe_domain_name, std::string sigma, std::string u,
1267 std::string geom, ForcesAndSourcesCore::RuleHookFun rule,
1268 bool is_axisymmetric = false) {
1270
1271 using DomainEleOp = typename DomainEle::UserDataOperator;
1272
1273 auto op_loop_side = new OpLoopSide<DomainEle>(
1274 m_field, fe_domain_name, DIM, Sev::noisy,
1275 boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
1276 pip.push_back(op_loop_side);
1277
1278 CHKERR AddHOOps<DIM, DIM, DIM>::add(op_loop_side->getOpPtrVector(),
1279 {H1, HDIV}, geom);
1280
1281 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
1282 A>::template BiLinearForm<I>;
1283
1284 using OpMixDivULhs = typename B::template OpMixDivTimesVec<DIM>;
1285 using OpMixDivUCylLhs =
1286 typename B::template OpMixDivTimesVec<DIM, CYLINDRICAL>;
1287 using OpLambdaGraULhs = typename B::template OpMixTensorTimesGrad<DIM>;
1288
1289 using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
1290 using OpMixDivUCylLhsSide = OpMixLhsSide<OpMixDivUCylLhs>;
1291 using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
1292
1293 auto unity = []() { return 1; };
1294 auto jacobian = [is_axisymmetric](const double r, const double,
1295 const double) {
1296 if (is_axisymmetric)
1297 return 2. * M_PI * r;
1298 else
1299 return 1.;
1300 };
1301
1302 if (!is_axisymmetric) {
1303 op_loop_side->getOpPtrVector().push_back(
1304 new OpMixDivULhsSide(sigma, u, unity, jacobian, true));
1305 } else {
1306 op_loop_side->getOpPtrVector().push_back(
1307 new OpMixDivUCylLhsSide(sigma, u, unity, jacobian, true));
1308 }
1309 op_loop_side->getOpPtrVector().push_back(
1310 new OpLambdaGraULhsSide(sigma, u, unity, jacobian, true));
1311
1312 op_loop_side->getSideFEPtr()->getRuleHook = rule;
1314}
1315
1316template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
1318 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1319 std::string sigma, std::string u, bool is_axisymmetric = false) {
1321
1323
1324 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1325
1326 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1327 u, common_data_ptr->contactDispPtr()));
1328 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1329 sigma, common_data_ptr->contactTractionPtr()));
1330 pip.push_back(
1331 new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
1332 DIM, GAUSS>(sigma, u, common_data_ptr, is_axisymmetric));
1333 pip.push_back(new typename C::template Assembly<A>::
1334 template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
1335 sigma, sigma, common_data_ptr, is_axisymmetric));
1336
1338}
1339
1340template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
1342 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1343 std::string sigma, std::string u, bool is_axisymmetric = false) {
1345
1347
1348 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1349
1350 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1351 u, common_data_ptr->contactDispPtr()));
1352 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1353 sigma, common_data_ptr->contactTractionPtr()));
1354 pip.push_back(
1355 new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
1356 DIM, GAUSS>(sigma, common_data_ptr, is_axisymmetric));
1357
1359}
1360
1361template <int DIM, IntegrationType I, typename BoundaryEleOp>
1363 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1364 std::string sigma, bool is_axisymmetric = false) {
1366
1368
1369 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1370 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1371 sigma, common_data_ptr->contactTractionPtr()));
1372 pip.push_back(new typename C::template OpAssembleTotalContactTraction<DIM, I>(
1373 common_data_ptr, 1. / scale, is_axisymmetric));
1374
1376}
1377
1378template <int DIM, IntegrationType I, typename BoundaryEleOp>
1380 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1381 OpLoopSide<SideEle> *op_loop_side, std::string sigma, std::string u,
1382 bool is_axisymmetric = false,
1383 boost::shared_ptr<Range> contact_range_ptr = nullptr) {
1386
1387 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1388
1389 op_loop_side->getOpPtrVector().push_back(
1390 new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
1391 "U", common_data_ptr->contactDispGradPtr()));
1392
1393 if (contact_range_ptr) {
1394 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1395 u, common_data_ptr->contactDispPtr()));
1396 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1397 sigma, common_data_ptr->contactTractionPtr()));
1398 pip.push_back(op_loop_side);
1399 pip.push_back(new typename C::template OpAssembleTotalContactArea<DIM, I>(
1400 common_data_ptr, is_axisymmetric, contact_range_ptr));
1401 }
1403}
1404
1405}; // namespace ContactOps
1406
1407#endif // __CONTACTOPS_HPP__
static const double eps
EntitiesFieldData::EntData EntData
Data on entities.
PetscBool is_axisymmetric
Definition contact.cpp:93
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
@ F
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
constexpr int DIM2
Definition level_set.cpp:22
constexpr int DIM1
Definition level_set.cpp:21
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
MoFEMErrorCode opFactoryDomainRhs(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)> HessSurfaceDistanceFunction
auto get_normalize_normals(FTensor::Tensor1< T1, DIM1 > &&t_normal_at_pts, size_t nb_gauss_pts)
double cn_contact
Definition contact.cpp:99
MoFEMErrorCode opFactoryCalculateArea(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, OpLoopSide< SideEle > *op_loop_side, std::string sigma, std::string u, bool is_axisymmetric=false, boost::shared_ptr< Range > contact_range_ptr=nullptr)
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)
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)
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)
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]
auto get_spatial_coords(FTensor::Tensor1< T1, DIM1 > &&t_coords, FTensor::Tensor1< T2, DIM2 > &&t_disp, size_t nb_gauss_pts)
boost::function< MatrixDouble( double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> GradSurfaceDistanceFunction
double constrain(double sdf, double tn)
constrain function
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)
MoFEMErrorCode opFactoryBoundaryLhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
MoFEMErrorCode opFactoryCalculateTraction(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, bool is_axisymmetric=false)
double sign(double x)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition sdf.py:1
hess_sdf(t, x, y, z, tx, ty, tz)
Definition sdf.py:19
grad_sdf(t, x, y, z, tx, ty, tz)
Definition sdf.py:15
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
constexpr double t
plate stiffness
Definition plate.cpp:58
constexpr auto field_name
MatrixDouble contactDispGrad
MatrixDouble contactDisp
static auto getFTensor1TotalTraction()
VectorDouble constraintVals
MatrixDouble contactTraction
static SmartPetscObj< Vec > totalTraction
static auto createTotalTraction(MoFEM::Interface &m_field)
VectorDouble sdfVals
size is equal to number of gauss points on element
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.