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 ++t_coords;
269 }
270
271 return v_sdf;
272}
273
274inline MatrixDouble
275grad_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
276 MatrixDouble &m_spatial_coords,
277 MatrixDouble &m_normals_at_pts, int block_id) {
278#ifdef PYTHON_SDF
279 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
280
281 VectorDouble v_spatial_coords = m_spatial_coords.data();
282 VectorDouble v_normal_at_pts = m_normals_at_pts.data();
283
284 bp::list python_coords;
285 bp::list python_normals;
286
287 for (int idx = 0; idx < 3; ++idx) {
288 python_coords.append(
289 convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
290 python_normals.append(
291 convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
292 }
293
294 np::ndarray np_grad_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 3),
295 np::dtype::get_builtin<double>());
296 CHK_MOAB_THROW(sdf_ptr->evalGradSdf(
297 delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
298 bp::extract<np::ndarray>(python_coords[1]),
299 bp::extract<np::ndarray>(python_coords[2]),
300 bp::extract<np::ndarray>(python_normals[0]),
301 bp::extract<np::ndarray>(python_normals[1]),
302 bp::extract<np::ndarray>(python_normals[2]), block_id,
303 np_grad_sdf),
304 "Failed python call");
305
306 double *grad_ptr = reinterpret_cast<double *>(np_grad_sdf.get_data());
307
308 MatrixDouble m_grad_sdf;
309 m_grad_sdf.resize(3, nb_gauss_pts, false);
310 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
311 for (int idx = 0; idx < 3; ++idx)
312 m_grad_sdf(idx, gg) = *(grad_ptr + (3 * gg + idx));
313 }
314 return m_grad_sdf;
315 }
316#endif
317 MatrixDouble m_grad_sdf;
318 m_grad_sdf.resize(3, nb_gauss_pts, false);
319 FTensor::Index<'i', 3> i;
320 FTensor::Tensor1<double, 3> t_grad_sdf_set{0.0, 1.0, 0.0};
321 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
322
323 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
324 t_grad_sdf(i) = t_grad_sdf_set(i);
325 ++t_grad_sdf;
326 }
327
328 return m_grad_sdf;
329}
330
331inline MatrixDouble
332hess_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
333 MatrixDouble &m_spatial_coords,
334 MatrixDouble &m_normals_at_pts, int block_id) {
335#ifdef PYTHON_SDF
336 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
337
338 VectorDouble v_spatial_coords = m_spatial_coords.data();
339 VectorDouble v_normal_at_pts = m_normals_at_pts.data();
340
341 bp::list python_coords;
342 bp::list python_normals;
343
344 for (int idx = 0; idx < 3; ++idx) {
345 python_coords.append(
346 convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
347 python_normals.append(
348 convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
349 };
350
351 np::ndarray np_hess_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 6),
352 np::dtype::get_builtin<double>());
353 CHK_MOAB_THROW(sdf_ptr->evalHessSdf(
354 delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
355 bp::extract<np::ndarray>(python_coords[1]),
356 bp::extract<np::ndarray>(python_coords[2]),
357 bp::extract<np::ndarray>(python_normals[0]),
358 bp::extract<np::ndarray>(python_normals[1]),
359 bp::extract<np::ndarray>(python_normals[2]), block_id,
360 np_hess_sdf),
361 "Failed python call");
362
363 double *hess_ptr = reinterpret_cast<double *>(np_hess_sdf.get_data());
364
365 MatrixDouble m_hess_sdf;
366 m_hess_sdf.resize(6, nb_gauss_pts, false);
367 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
368 for (int idx = 0; idx < 6; ++idx)
369 m_hess_sdf(idx, gg) =
370 *(hess_ptr + (6 * gg + idx));
371 }
372 return m_hess_sdf;
373 }
374#endif
375 MatrixDouble m_hess_sdf;
376 m_hess_sdf.resize(6, nb_gauss_pts, false);
377 FTensor::Index<'i', 3> i;
378 FTensor::Index<'j', 3> j;
379 FTensor::Tensor2_symmetric<double, 3> t_hess_sdf_set{0., 0., 0., 0., 0., 0.};
380 auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
381
382 for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
383 t_hess_sdf(i, j) = t_hess_sdf_set(i, j);
384 ++t_hess_sdf;
385 }
386 return m_hess_sdf;
387}
388
389
390
391template <int DIM, IntegrationType I, typename BoundaryEleOp>
393
394template <int DIM, IntegrationType I, typename BoundaryEleOp>
396
397template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
399
400template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
402
403template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
405
406template <typename T1, typename T2, int DIM1, int DIM2>
409 size_t nb_gauss_pts) {
410 MatrixDouble m_spatial_coords(nb_gauss_pts, 3);
411 m_spatial_coords.clear();
412 auto t_spatial_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
413 FTensor::Index<'i', DIM2> i;
414 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
415 t_spatial_coords(i) = t_coords(i) + t_disp(i);
416 ++t_spatial_coords;
417 ++t_coords;
418 ++t_disp;
419 }
420 return m_spatial_coords;
421}
422
423template <typename T1, int DIM1>
425 size_t nb_gauss_pts) {
426 MatrixDouble m_normals_at_pts(3, nb_gauss_pts);
427 m_normals_at_pts.clear();
428 FTensor::Index<'i', DIM1> i;
429 auto t_set_normal = getFTensor1FromMat<3>(m_normals_at_pts);
430 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
431 t_set_normal(i) = t_normal_at_pts(i) / t_normal_at_pts.l2();
432 ++t_set_normal;
433 ++t_normal_at_pts;
434 }
435 return m_normals_at_pts;
436}
437
438template <int DIM, typename BoundaryEleOp>
440 : public BoundaryEleOp {
442 boost::shared_ptr<CommonData> common_data_ptr, double scale = 1,
443 bool is_axisymmetric = false);
444 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
445
446private:
447 boost::shared_ptr<CommonData> commonDataPtr;
448 const double scaleTraction;
450};
451
452template <int DIM, typename BoundaryEleOp>
453struct OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp> : public BoundaryEleOp {
454 OpEvaluateSDFImpl(boost::shared_ptr<CommonData> common_data_ptr);
455 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
456
457private:
458 boost::shared_ptr<CommonData> commonDataPtr;
459
461 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
463 HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
465};
466
467template <int DIM, typename AssemblyBoundaryEleOp>
469 : public AssemblyBoundaryEleOp {
470 OpConstrainBoundaryRhsImpl(const std::string field_name,
471 boost::shared_ptr<CommonData> common_data_ptr,
472 bool is_axisymmetric = false);
473 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
474
476 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
478
479private:
480 boost::shared_ptr<CommonData> commonDataPtr;
482};
483
484template <int DIM, typename AssemblyBoundaryEleOp>
486 : public AssemblyBoundaryEleOp {
487 OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
488 const std::string col_field_name,
489 boost::shared_ptr<CommonData> common_data_ptr,
490 bool is_axisymmetric = false);
491 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
492 EntitiesFieldData::EntData &col_data);
493
495 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
497 HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
499
500 boost::shared_ptr<CommonData> commonDataPtr;
502};
503
504template <int DIM, typename AssemblyBoundaryEleOp>
506 : public AssemblyBoundaryEleOp {
508 const std::string row_field_name, const std::string col_field_name,
509 boost::shared_ptr<CommonData> common_data_ptr,
510 bool is_axisymmetric = false);
511 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
512 EntitiesFieldData::EntData &col_data);
513
515 GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
517
518private:
519 boost::shared_ptr<CommonData> commonDataPtr;
521};
522
523template <typename BoundaryEleOp> struct ContactIntegrators {
524
525 template <int DIM, IntegrationType I>
528
529 template <int DIM, IntegrationType I>
531
532 template <AssemblyType A> struct Assembly {
533
535 typename FormsIntegrators<BoundaryEleOp>::template Assembly<A>::OpBase;
536
537 template <int DIM, IntegrationType I>
540
541 template <int DIM, IntegrationType I>
544
545 template <int DIM, IntegrationType I>
548 };
549};
550
551inline double sign(double x) {
552 constexpr auto eps = std::numeric_limits<float>::epsilon();
553 if (std::abs(x) < eps)
554 return 0;
555 else if (x > eps)
556 return 1;
557 else
558 return -1;
559};
560
561inline double w(const double sdf, const double tn) {
562 return sdf - cn_contact * tn;
563}
564
565/**
566 * @brief constrain function
567 *
568 * return 1 if negative sdf or positive tn
569 *
570 * @param sdf signed distance
571 * @param tn traction
572 * @return double
573 */
574inline double constrain(double sdf, double tn) {
575 const auto s = sign(w(sdf, tn));
576 return (1 - s) / 2;
577}
578
579template <int DIM, typename BoundaryEleOp>
582 boost::shared_ptr<CommonData> common_data_ptr, double scale,
583 bool is_axisymmetric)
585 commonDataPtr(common_data_ptr), scaleTraction(scale),
586 isAxisymmetric(is_axisymmetric) {}
587
588template <int DIM, typename BoundaryEleOp>
589MoFEMErrorCode
591 int side, EntityType type, EntData &data) {
593
594 FTensor::Index<'i', DIM> i;
595 FTensor::Tensor1<double, 3> t_sum_t{0., 0., 0.};
596
598 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
600
601 const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
602 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
603 double jacobian = 1.;
604 if (isAxisymmetric) {
605 jacobian = 2. * M_PI * t_coords(0);
606 }
607 const double alpha = t_w * jacobian * BoundaryEleOp::getMeasure();
608 t_sum_t(i) += alpha * t_traction(i);
609 ++t_w;
610 ++t_traction;
611 ++t_coords;
612 }
613
614 t_sum_t(i) *= scaleTraction;
615
616 constexpr int ind[] = {0, 1, 2};
617 CHKERR VecSetValues(commonDataPtr->totalTraction, 3, ind, &t_sum_t(0),
618 ADD_VALUES);
619
621}
622
623template <int DIM, typename BoundaryEleOp>
625 boost::shared_ptr<CommonData> common_data_ptr)
627 commonDataPtr(common_data_ptr) {}
628
629template <int DIM, typename BoundaryEleOp>
630MoFEMErrorCode
632 EntData &data) {
634
635 const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
636 auto &sdf_vec = commonDataPtr->sdfVals;
637 auto &grad_mat = commonDataPtr->gradsSdf;
638 auto &hess_mat = commonDataPtr->hessSdf;
639 auto &constraint_vec = commonDataPtr->constraintVals;
640 auto &contactTraction_mat = commonDataPtr->contactTraction;
641
642 sdf_vec.resize(nb_gauss_pts, false);
643 grad_mat.resize(DIM, nb_gauss_pts, false);
644 hess_mat.resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
645 constraint_vec.resize(nb_gauss_pts, false);
646
647 auto t_traction = getFTensor1FromMat<DIM>(contactTraction_mat);
648
649 auto t_sdf = getFTensor0FromVec(sdf_vec);
650 auto t_grad_sdf = getFTensor1FromMat<DIM>(grad_mat);
651 auto t_hess_sdf = getFTensor2SymmetricFromMat<DIM>(hess_mat);
652 auto t_constraint = getFTensor0FromVec(constraint_vec);
653
654 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
656 auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
657
658 FTensor::Index<'i', DIM> i;
659 FTensor::Index<'j', DIM> j;
660
661 auto ts_time = BoundaryEleOp::getTStime();
662 auto ts_time_step = BoundaryEleOp::getTStimeStep();
663
664 auto m_spatial_coords = get_spatial_coords(
666 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
667 auto m_normals_at_pts = get_normalize_normals(
669
670 // placeholder to pass boundary block id to python
671 int block_id = 0;
672
673 auto v_sdf =
674 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
675 m_spatial_coords, m_normals_at_pts, block_id);
676
677 auto m_grad_sdf =
678 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
679 m_spatial_coords, m_normals_at_pts, block_id);
680
681 auto m_hess_sdf =
682 hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
683 m_spatial_coords, m_normals_at_pts, block_id);
684
685 auto t_sdf_v = getFTensor0FromVec(v_sdf);
686 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
687 auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
688
689 auto next = [&]() {
690 ++t_sdf;
691 ++t_sdf_v;
692 ++t_grad_sdf;
693 ++t_grad_sdf_v;
694 ++t_hess_sdf;
695 ++t_hess_sdf_v;
696 ++t_disp;
697 ++t_traction;
698 ++t_constraint;
699 };
700
701 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
702
703 auto tn = -t_traction(i) * t_grad_sdf_v(i);
704 auto c = constrain(t_sdf_v, tn);
705
706 t_sdf = t_sdf_v;
707 t_grad_sdf(i) = t_grad_sdf_v(i);
708 t_hess_sdf(i, j) = t_hess_sdf_v(i, j);
709 t_constraint = c;
710
711 next();
712 }
713
715}
716
717template <int DIM, typename AssemblyBoundaryEleOp>
720 boost::shared_ptr<CommonData> common_data_ptr,
721 bool is_axisymmetric)
723 AssemblyBoundaryEleOp::OPROW),
724 commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {}
725
726template <int DIM, typename AssemblyBoundaryEleOp>
727MoFEMErrorCode
729 EntitiesFieldData::EntData &data) {
731
732 FTensor::Index<'i', DIM> i;
733 FTensor::Index<'j', DIM> j;
734 FTensor::Index<'k', DIM> k;
735 FTensor::Index<'l', DIM> l;
736
737 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
738
740
741 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
742
743 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
744 auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
745 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
746 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
747
748 size_t nb_base_functions = data.getN().size2() / 3;
749 auto t_base = data.getFTensor1N<3>();
750
751 auto m_spatial_coords = get_spatial_coords(
753 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
754 auto m_normals_at_pts = get_normalize_normals(
756
757 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
758
759 auto ts_time = AssemblyBoundaryEleOp::getTStime();
760 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
761
762 // placeholder to pass boundary block id to python
763 int block_id = 0;
764
765 auto v_sdf =
766 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
767 m_spatial_coords, m_normals_at_pts, block_id);
768
769 auto m_grad_sdf =
770 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
771 m_spatial_coords, m_normals_at_pts, block_id);
772
773 auto t_sdf = getFTensor0FromVec(v_sdf);
774 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
775
776 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
777
778 auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
779
780 double jacobian = 1.;
781 if (isAxisymmetric) {
782 jacobian = 2. * M_PI * t_coords(0);
783 }
784 const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
785
786 auto tn = -t_traction(i) * t_grad_sdf(i);
787 auto c = constrain(t_sdf, tn);
788
790 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
792 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
793
795 t_rhs(i) =
796
797 t_cQ(i, j) * (t_disp(j) - cn_contact * t_traction(j))
798
799 +
800
801 t_cP(i, j) * t_disp(j) +
802 c * (t_sdf * t_grad_sdf(i)); // add gap0 displacements
803
804 size_t bb = 0;
805 for (; bb != AssemblyBoundaryEleOp::nbRows / DIM; ++bb) {
806 const double beta = alpha * (t_base(i) * t_normal(i));
807 t_nf(i) -= beta * t_rhs(i);
808
809 ++t_nf;
810 ++t_base;
811 }
812 for (; bb < nb_base_functions; ++bb)
813 ++t_base;
814
815 ++t_disp;
816 ++t_traction;
817 ++t_coords;
818 ++t_w;
819 ++t_normal;
820 ++t_sdf;
821 ++t_grad_sdf;
822 }
823
825}
826
827template <int DIM, typename AssemblyBoundaryEleOp>
829 OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
830 const std::string col_field_name,
831 boost::shared_ptr<CommonData> common_data_ptr,
832 bool is_axisymmetric)
833 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
834 AssemblyBoundaryEleOp::OPROWCOL),
835 commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
836 AssemblyBoundaryEleOp::sYmm = false;
837}
838
839template <int DIM, typename AssemblyBoundaryEleOp>
840MoFEMErrorCode
842 EntitiesFieldData::EntData &row_data,
843 EntitiesFieldData::EntData &col_data) {
845
846 FTensor::Index<'i', DIM> i;
847 FTensor::Index<'j', DIM> j;
848 FTensor::Index<'k', DIM> k;
849
850 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
851 auto &locMat = AssemblyBoundaryEleOp::locMat;
852
853 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
854 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
855 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
856
857 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
858 auto t_row_base = row_data.getFTensor1N<3>();
859 size_t nb_face_functions = row_data.getN().size2() / 3;
860
861 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
862
863 auto m_spatial_coords = get_spatial_coords(
865 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
866 auto m_normals_at_pts = get_normalize_normals(
868
869 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
870
871 auto ts_time = AssemblyBoundaryEleOp::getTStime();
872 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
873
874 // placeholder to pass boundary block id to python
875 int block_id = 0;
876
877 auto v_sdf =
878 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
879 m_spatial_coords, m_normals_at_pts, block_id);
880
881 auto m_grad_sdf =
882 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
883 m_spatial_coords, m_normals_at_pts, block_id);
884
885 auto m_hess_sdf =
886 hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
887 m_spatial_coords, m_normals_at_pts, block_id);
888
889 auto t_sdf = getFTensor0FromVec(v_sdf);
890 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
891 auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
892
893 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
894
895 double jacobian = 1.;
896 if (isAxisymmetric) {
897 jacobian = 2. * M_PI * t_coords(0);
898 }
899 const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
900
901 auto tn = -t_traction(i) * t_grad_sdf(i);
902 auto c = constrain(t_sdf, tn);
903
905 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
907 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
908
910 t_res_dU(i, j) = kronecker_delta(i, j) + t_cP(i, j);
911
912 if (c > 0) {
913 t_res_dU(i, j) +=
914 (c * cn_contact) *
915 (t_hess_sdf(i, j) * (t_grad_sdf(k) * t_traction(k)) +
916 t_grad_sdf(i) * t_hess_sdf(k, j) * t_traction(k)) +
917 c * t_sdf * t_hess_sdf(i, j);
918 }
919
920 size_t rr = 0;
921 for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
922
923 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
924
925 const double row_base = t_row_base(i) * t_normal(i);
926
927 auto t_col_base = col_data.getFTensor0N(gg, 0);
928 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
929 const double beta = alpha * row_base * t_col_base;
930
931 t_mat(i, j) -= beta * t_res_dU(i, j);
932
933 ++t_col_base;
934 ++t_mat;
935 }
936
937 ++t_row_base;
938 }
939 for (; rr < nb_face_functions; ++rr)
940 ++t_row_base;
941
942 ++t_traction;
943 ++t_coords;
944 ++t_w;
945 ++t_normal;
946 ++t_sdf;
947 ++t_grad_sdf;
948 ++t_hess_sdf;
949 }
950
952}
953
954template <int DIM, typename AssemblyBoundaryEleOp>
957 const std::string row_field_name, const std::string col_field_name,
958 boost::shared_ptr<CommonData> common_data_ptr, bool is_axisymmetric)
959 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
960 AssemblyBoundaryEleOp::OPROWCOL),
961 commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
962 AssemblyBoundaryEleOp::sYmm = false;
963}
964
965template <int DIM, typename AssemblyBoundaryEleOp>
966MoFEMErrorCode
968 iNtegrate(EntitiesFieldData::EntData &row_data,
969 EntitiesFieldData::EntData &col_data) {
971
972 FTensor::Index<'i', DIM> i;
973 FTensor::Index<'j', DIM> j;
974 FTensor::Index<'k', DIM> k;
975
976 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
977 auto &locMat = AssemblyBoundaryEleOp::locMat;
978
979 auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
980 auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
981 auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
982
983 auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
984 auto t_row_base = row_data.getFTensor1N<3>();
985 size_t nb_face_functions = row_data.getN().size2() / 3;
986
987 auto m_spatial_coords = get_spatial_coords(
989 getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
990 auto m_normals_at_pts = get_normalize_normals(
992
993 auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
994
995 auto ts_time = AssemblyBoundaryEleOp::getTStime();
996 auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
997
998 // placeholder to pass boundary block id to python
999 int block_id = 0;
1000
1001 auto v_sdf =
1002 surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1003 m_spatial_coords, m_normals_at_pts, block_id);
1004
1005 auto m_grad_sdf =
1006 gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1007 m_spatial_coords, m_normals_at_pts, block_id);
1008
1009 auto t_sdf = getFTensor0FromVec(v_sdf);
1010 auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_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_dt(i, j) = -cn_contact * t_cQ(i, j);
1030
1031 size_t rr = 0;
1032 for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
1033
1034 auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
1035 const double row_base = t_row_base(i) * t_normal(i);
1036
1037 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
1038 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
1039 const double col_base = t_col_base(i) * t_normal(i);
1040 const double beta = alpha * row_base * col_base;
1041
1042 t_mat(i, j) -= beta * t_res_dt(i, j);
1043
1044 ++t_col_base;
1045 ++t_mat;
1046 }
1047
1048 ++t_row_base;
1049 }
1050 for (; rr < nb_face_functions; ++rr)
1051 ++t_row_base;
1052
1053 ++t_traction;
1054 ++t_coords;
1055 ++t_w;
1056 ++t_normal;
1057 ++t_sdf;
1058 ++t_grad_sdf;
1059 }
1060
1062}
1063
1064template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
1065MoFEMErrorCode opFactoryDomainRhs(
1066 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1067 std::string sigma, std::string u, bool is_axisymmetric = false) {
1069
1070 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
1071 A>::template LinearForm<I>;
1072 using OpMixDivURhs = typename B::template OpMixDivTimesU<3, DIM, DIM>;
1073 using OpMixDivUCylRhs =
1074 typename B::template OpMixDivTimesU<3, DIM, DIM, CYLINDRICAL>;
1075
1076 using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
1078 typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
1079 using OpMixUTimesLambdaRhs =
1080 typename B::template OpGradTimesTensor<1, DIM, DIM>;
1081
1082 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1083 auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
1084 auto div_stress_ptr = boost::make_shared<MatrixDouble>();
1085 auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
1086
1087 auto jacobian = [is_axisymmetric](const double r, const double,
1088 const double) {
1089 if (is_axisymmetric)
1090 return 2. * M_PI * r;
1091 else
1092 return 1.;
1093 };
1094
1095 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1096 u, common_data_ptr->contactDispPtr()));
1097 pip.push_back(
1098 new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
1099
1100 if (!is_axisymmetric) {
1101 pip.push_back(
1102 new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
1103 } else {
1104 pip.push_back(new OpCalculateHVecTensorDivergence<DIM, DIM, CYLINDRICAL>(
1105 sigma, div_stress_ptr));
1106 }
1107
1108 pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(u, mat_grad_ptr));
1109
1110 if (!is_axisymmetric) {
1111 pip.push_back(
1112 new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(), jacobian));
1113 } else {
1114 pip.push_back(new OpMixDivUCylRhs(sigma, common_data_ptr->contactDispPtr(),
1115 jacobian));
1116 }
1117
1118 pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr, jacobian));
1119 pip.push_back(new OpMixUTimesDivLambdaRhs(u, div_stress_ptr, jacobian));
1120 pip.push_back(new OpMixUTimesLambdaRhs(u, contact_stress_ptr, jacobian));
1121
1123}
1124
1125template <typename OpMixLhs> struct OpMixLhsSide : public OpMixLhs {
1126 using OpMixLhs::OpMixLhs;
1127 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
1128 EntityType col_type,
1129 EntitiesFieldData::EntData &row_data,
1130 EntitiesFieldData::EntData &col_data) {
1132 auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
1133 auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
1134 // Only assemble side which correspond to edge entity on boundary
1135 if (side_fe_entity == side_fe_data) {
1136 CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
1137 col_data);
1138 }
1140 }
1141};
1142
1143template <int DIM, AssemblyType A, IntegrationType I, typename DomainEle>
1145 MoFEM::Interface &m_field,
1146 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1147 std::string fe_domain_name, std::string sigma, std::string u,
1148 std::string geom, ForcesAndSourcesCore::RuleHookFun rule,
1149 bool is_axisymmetric = false) {
1151
1152 using DomainEleOp = typename DomainEle::UserDataOperator;
1153
1154 auto op_loop_side = new OpLoopSide<DomainEle>(
1155 m_field, fe_domain_name, DIM, Sev::noisy,
1156 boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
1157 pip.push_back(op_loop_side);
1158
1159 CHKERR AddHOOps<DIM, DIM, DIM>::add(op_loop_side->getOpPtrVector(),
1160 {H1, HDIV}, geom);
1161
1162 using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
1163 A>::template BiLinearForm<I>;
1164
1165 using OpMixDivULhs = typename B::template OpMixDivTimesVec<DIM>;
1166 using OpMixDivUCylLhs =
1167 typename B::template OpMixDivTimesVec<DIM, CYLINDRICAL>;
1168 using OpLambdaGraULhs = typename B::template OpMixTensorTimesGrad<DIM>;
1169
1170 using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
1171 using OpMixDivUCylLhsSide = OpMixLhsSide<OpMixDivUCylLhs>;
1172 using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
1173
1174 auto unity = []() { return 1; };
1175 auto jacobian = [is_axisymmetric](const double r, const double,
1176 const double) {
1177 if (is_axisymmetric)
1178 return 2. * M_PI * r;
1179 else
1180 return 1.;
1181 };
1182
1183 if (!is_axisymmetric) {
1184 op_loop_side->getOpPtrVector().push_back(
1185 new OpMixDivULhsSide(sigma, u, unity, jacobian, true));
1186 } else {
1187 op_loop_side->getOpPtrVector().push_back(
1188 new OpMixDivUCylLhsSide(sigma, u, unity, jacobian, true));
1189 }
1190 op_loop_side->getOpPtrVector().push_back(
1191 new OpLambdaGraULhsSide(sigma, u, unity, jacobian, true));
1192
1193 op_loop_side->getSideFEPtr()->getRuleHook = rule;
1195}
1196
1197template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
1199 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1200 std::string sigma, std::string u, bool is_axisymmetric = false) {
1202
1204
1205 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1206
1207 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1208 u, common_data_ptr->contactDispPtr()));
1209 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1210 sigma, common_data_ptr->contactTractionPtr()));
1211 pip.push_back(
1212 new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
1213 DIM, GAUSS>(sigma, u, common_data_ptr, is_axisymmetric));
1214 pip.push_back(new typename C::template Assembly<A>::
1215 template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
1216 sigma, sigma, common_data_ptr, is_axisymmetric));
1217
1219}
1220
1221template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
1223 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1224 std::string sigma, std::string u, bool is_axisymmetric = false) {
1226
1228
1229 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1230
1231 pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1232 u, common_data_ptr->contactDispPtr()));
1233 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1234 sigma, common_data_ptr->contactTractionPtr()));
1235 pip.push_back(
1236 new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
1237 DIM, GAUSS>(sigma, common_data_ptr, is_axisymmetric));
1238
1240}
1241
1242template <int DIM, IntegrationType I, typename BoundaryEleOp>
1244 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1245 std::string sigma, bool is_axisymmetric = false) {
1247
1249
1250 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1251 pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1252 sigma, common_data_ptr->contactTractionPtr()));
1253 pip.push_back(new typename C::template OpAssembleTotalContactTraction<DIM, I>(
1254 common_data_ptr, 1. / scale, is_axisymmetric));
1255
1257}
1258
1259}; // namespace ContactOps
1260
1261#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.
@ 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.
constexpr auto t_kd
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
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
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)
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]
EntitiesFieldData::EntData EntData
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:59
constexpr auto field_name
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.
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