v0.13.2
Loading...
Searching...
No Matches
ContactOps.hpp
Go to the documentation of this file.
1
2
3/**
4 * \file ContactOps.hpp
5 * \example ContactOps.hpp
6 */
7
8namespace ContactOps {
9
10//! [Common data]
11struct CommonData : public boost::enable_shared_from_this<CommonData> {
12 MatrixDouble mD;
13 MatrixDouble mGrad;
14
15 MatrixDouble contactStress;
17 MatrixDouble contactTraction;
18 MatrixDouble contactDisp;
19
20 inline auto mDPtr() {
21 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &mD);
22 }
23
24 inline auto mGradPtr() {
25 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &mGrad);
26 }
27
28 inline auto contactStressPtr() {
29 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactStress);
30 }
31
33 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
35 }
36
37 inline auto contactTractionPtr() {
38 return boost::shared_ptr<MatrixDouble>(shared_from_this(),
40 }
41
42 inline auto contactDispPtr() {
43 return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactDisp);
44 }
45};
46//! [Common data]
47
48//! [Surface distance function from python]
49#ifdef PYTHON_SFD
50struct SDFPython {
51 SDFPython() = default;
52 virtual ~SDFPython() = default;
53
54 MoFEMErrorCode sdfInit(const std::string py_file) {
56 try {
57
58 // create main module
59 auto main_module = bp::import("__main__");
60 mainNamespace = main_module.attr("__dict__");
61 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
62 // create a reference to python function
63 sdfFun = mainNamespace["sdf"];
64 sdfGradFun = mainNamespace["grad_sdf"];
65 sdfHessFun = mainNamespace["hess_sdf"];
66
67 } catch (bp::error_already_set const &) {
68 // print all other errors to stderr
69 PyErr_Print();
71 }
73 };
74
75 template <typename T>
76 inline std::vector<T>
77 py_list_to_std_vector(const boost::python::object &iterable) {
78 return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
79 boost::python::stl_input_iterator<T>());
80 }
81
82 MoFEMErrorCode evalSdf(
83
84 double t, double x, double y, double z, double tx, double ty, double tz,
85 double &sdf
86
87 ) {
89 try {
90
91 // call python function
92 sdf = bp::extract<double>(sdfFun(t, x, y, z, tx, ty, tz));
93
94 } catch (bp::error_already_set const &) {
95 // print all other errors to stderr
96 PyErr_Print();
98 }
100 }
101
102 MoFEMErrorCode evalGradSdf(
103
104 double t, double x, double y, double z, double tx, double ty, double tz,
105 std::vector<double> &grad_sdf
106
107 ) {
109 try {
110
111 // call python function
112 grad_sdf =
113 py_list_to_std_vector<double>(sdfGradFun(t, x, y, z, tx, ty, tz));
114
115 } catch (bp::error_already_set const &) {
116 // print all other errors to stderr
117 PyErr_Print();
119 }
120
121 if (grad_sdf.size() != 3) {
122 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Expected size 3");
123 }
124
126 }
127
128 MoFEMErrorCode evalHessSdf(
129
130 double t, double x, double y, double z, double tx, double ty, double tz,
131 std::vector<double> &hess_sdf
132
133 ) {
135 try {
136
137 // call python function
138 hess_sdf =
139 py_list_to_std_vector<double>(sdfHessFun(t, x, y, z, tx, ty, tz));
140
141 } catch (bp::error_already_set const &) {
142 // print all other errors to stderr
143 PyErr_Print();
145 }
146
147 if (hess_sdf.size() != 6) {
148 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Expected size 6");
149 }
150
152 }
153
154private:
155 bp::object mainNamespace;
156 bp::object sdfFun;
157 bp::object sdfGradFun;
158 bp::object sdfHessFun;
159};
160
161static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
162
163#endif
164//! [Surface distance function from python]
165
170
172 OpConstrainBoundaryRhs(const std::string field_name,
173 boost::shared_ptr<CommonData> common_data_ptr);
174 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
175
176private:
177 boost::shared_ptr<CommonData> commonDataPtr;
178};
179
181 OpConstrainBoundaryLhs_dU(const std::string row_field_name,
182 const std::string col_field_name,
183 boost::shared_ptr<CommonData> common_data_ptr);
184 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
185 EntitiesFieldData::EntData &col_data);
186
187private:
188 boost::shared_ptr<CommonData> commonDataPtr;
189};
190
193 const std::string row_field_name, const std::string col_field_name,
194 boost::shared_ptr<CommonData> common_data_ptr);
195 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
196 EntitiesFieldData::EntData &col_data);
197
198private:
199 boost::shared_ptr<CommonData> commonDataPtr;
200};
201
202template <typename T1, typename T2>
203inline double
206#ifdef PYTHON_SFD
207 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
208 double sdf;
209 CHK_MOAB_THROW(sdf_ptr->evalSdf(t, t_coords(0), t_coords(1), t_coords(2),
210 t_traction(0), t_traction(1),
211 (SPACE_DIM == 3) ? t_traction(2) : 0, sdf),
212 "Failed python call");
213 return sdf;
214 }
215#endif
216 return t_coords(1) + 0.5;
217};
218
219template <typename T1, typename T2>
223#ifdef PYTHON_SFD
224 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
225 std::vector<double> grad_sdf;
227 sdf_ptr->evalGradSdf(t, t_coords(0), t_coords(1), t_coords(2),
228 t_traction(0), t_traction(1),
229 (SPACE_DIM == 3) ? t_traction(2) : 0, grad_sdf),
230 "Failed python call");
231 return FTensor::Tensor1<double, 3>{grad_sdf[0], grad_sdf[1], grad_sdf[2]};
232 }
233#endif
234 return FTensor::Tensor1<double, 3>{0., 1., 0.};
235};
236
237template <typename T1, typename T2>
241#ifdef PYTHON_SFD
242 if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
243 std::vector<double> hess_sdf;
245 sdf_ptr->evalHessSdf(t, t_coords(0), t_coords(1), t_coords(2),
246 t_traction(0), t_traction(1),
247 (SPACE_DIM == 3) ? t_traction(2) : 0, hess_sdf),
248 "Failed python call");
249 return FTensor::Tensor2_symmetric<double, 3>{hess_sdf[0], hess_sdf[1],
250 hess_sdf[2], hess_sdf[3],
251 hess_sdf[4], hess_sdf[5]};
252 }
253#endif
254 return FTensor::Tensor2_symmetric<double, 3>{0., 0., 0., 0., 0., 0.};
255};
256
257inline double sign(double x) {
258 if (x == 0)
259 return 0;
260 else if (x > 0)
261 return 1;
262 else
263 return -1;
264};
265
266inline double w(const double sdf, const double tn) { return sdf - cn * tn; }
267
268inline double constrain(double sdf, double tn) {
269 const auto s = sign(w(sdf, tn));
270 return (1 - s) / 2;
271}
272
274 const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
276 commonDataPtr(common_data_ptr) {}
277
278MoFEMErrorCode
279OpConstrainBoundaryRhs::iNtegrate(EntitiesFieldData::EntData &data) {
281
282 const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
283
284 auto &nf = AssemblyBoundaryEleOp::locF;
285
286 auto t_normal = getFTensor1Normal();
287 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
288
289 auto t_w = getFTensor0IntegrationWeight();
290 auto t_disp = getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactDisp);
291 auto t_traction =
292 getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactTraction);
293 auto t_coords = getFTensor1CoordsAtGaussPts();
294
295 size_t nb_base_functions = data.getN().size2() / 3;
296 auto t_base = data.getFTensor1N<3>();
297 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
298
299 auto t_nf = getFTensor1FromPtr<SPACE_DIM>(&nf[0]);
300 const double alpha = t_w * getMeasure();
301
302 FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
303 t_spatial_coords(i) = t_coords(i) + t_disp(i);
304
305 auto sdf =
306 surface_distance_function(getTStime(), t_spatial_coords, t_traction);
307 auto t_grad_sdf = grad_surface_distance_function(
308 getTStime(), t_spatial_coords, t_traction);
309 auto tn = -t_traction(i) * t_grad_sdf(i);
310 auto c = constrain(sdf, tn);
311
313 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
315 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
316
318 t_rhs(i) =
319
320 t_cQ(i, j) * (t_disp(j) - cn * t_traction(j))
321
322 +
323
324 t_cP(i, j) * t_disp(j) +
325 c * (sdf * t_grad_sdf(i)); // add gap0 displacements
326
327 size_t bb = 0;
328 for (; bb != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++bb) {
329 const double beta = alpha * (t_base(i) * t_normal(i));
330 t_nf(i) -= beta * t_rhs(i);
331
332 ++t_nf;
333 ++t_base;
334 }
335 for (; bb < nb_base_functions; ++bb)
336 ++t_base;
337
338 ++t_disp;
339 ++t_traction;
340 ++t_coords;
341 ++t_w;
342 }
343
345}
346
348 const std::string row_field_name, const std::string col_field_name,
349 boost::shared_ptr<CommonData> common_data_ptr)
350 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
351 DomainEleOp::OPROWCOL),
352 commonDataPtr(common_data_ptr) {
353 sYmm = false;
354}
355
356MoFEMErrorCode
357OpConstrainBoundaryLhs_dU::iNtegrate(EntitiesFieldData::EntData &row_data,
358 EntitiesFieldData::EntData &col_data) {
360
361 const size_t nb_gauss_pts = getGaussPts().size2();
362 auto &locMat = AssemblyBoundaryEleOp::locMat;
363
364 auto t_normal = getFTensor1Normal();
365 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
366
367 auto t_disp = getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactDisp);
368 auto t_traction =
369 getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactTraction);
370 auto t_coords = getFTensor1CoordsAtGaussPts();
371
372 auto t_w = getFTensor0IntegrationWeight();
373 auto t_row_base = row_data.getFTensor1N<3>();
374 size_t nb_face_functions = row_data.getN().size2() / 3;
375
376 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
377
378 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
379
380 const double alpha = t_w * getMeasure();
381
382 FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
383 t_spatial_coords(i) = t_coords(i) + t_disp(i);
384
385 auto sdf =
386 surface_distance_function(getTStime(), t_spatial_coords, t_traction);
387 auto t_grad_sdf = grad_surface_distance_function(
388 getTStime(), t_spatial_coords, t_traction);
389
390 auto tn = -t_traction(i) * t_grad_sdf(i);
391 auto c = constrain(sdf, tn);
392
394 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
396 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
397
399 t_res_dU(i, j) =
400 kronecker_delta(i, j) + t_cP(i, j);
401
402 if(c>0) {
403 auto t_hess_sdf = hess_surface_distance_function(
404 getTStime(), t_spatial_coords, t_traction);
405 t_res_dU(i, j) +=
406 (c * cn) * (t_hess_sdf(i, j) * (t_grad_sdf(k) * t_traction(k)) +
407 t_grad_sdf(i) * t_hess_sdf(k, j) * t_traction(k)) +
408 c * sdf * t_hess_sdf(i, j);
409 }
410
411 size_t rr = 0;
412 for (; rr != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++rr) {
413
414 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
415 locMat, SPACE_DIM * rr);
416
417 const double row_base = t_row_base(i) * t_normal(i);
418
419 auto t_col_base = col_data.getFTensor0N(gg, 0);
420 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / SPACE_DIM;
421 ++cc) {
422 const double beta = alpha * row_base * t_col_base;
423
424 t_mat(i, j) -= beta * t_res_dU(i, j);
425
426 ++t_col_base;
427 ++t_mat;
428 }
429
430 ++t_row_base;
431 }
432 for (; rr < nb_face_functions; ++rr)
433 ++t_row_base;
434
435 ++t_disp;
436 ++t_traction;
437 ++t_coords;
438 ++t_w;
439 }
440
442}
443
445 const std::string row_field_name, const std::string col_field_name,
446 boost::shared_ptr<CommonData> common_data_ptr)
447 : AssemblyBoundaryEleOp(row_field_name, col_field_name,
448 DomainEleOp::OPROWCOL),
449 commonDataPtr(common_data_ptr) {
450 sYmm = false;
451}
452
454 EntitiesFieldData::EntData &row_data,
455 EntitiesFieldData::EntData &col_data) {
457
458 const size_t nb_gauss_pts = getGaussPts().size2();
459 auto &locMat = AssemblyBoundaryEleOp::locMat;
460
461 auto t_normal = getFTensor1Normal();
462 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
463
464 auto t_disp = getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactDisp);
465 auto t_traction =
466 getFTensor1FromMat<SPACE_DIM>(commonDataPtr->contactTraction);
467 auto t_coords = getFTensor1CoordsAtGaussPts();
468
469 auto t_w = getFTensor0IntegrationWeight();
470 auto t_row_base = row_data.getFTensor1N<3>();
471 size_t nb_face_functions = row_data.getN().size2() / 3;
472
473 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
474
475 const double alpha = t_w * getMeasure();
476
477 FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
478 t_spatial_coords(i) = t_coords(i) + t_disp(i);
479
480 auto sdf =
481 surface_distance_function(getTStime(), t_spatial_coords, t_traction);
482 auto t_grad_sdf = grad_surface_distance_function(
483 getTStime(), t_spatial_coords, t_traction);
484 auto tn = -t_traction(i) * t_grad_sdf(i);
485 auto c = constrain(sdf, tn);
486
488 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
490 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
491
493 t_res_dt(i, j) = -cn * t_cQ(i, j);
494
495 size_t rr = 0;
496 for (; rr != AssemblyBoundaryEleOp::nbRows / SPACE_DIM; ++rr) {
497
498 auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
499 locMat, SPACE_DIM * rr);
500 const double row_base = t_row_base(i) * t_normal(i);
501
502 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
503 for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / SPACE_DIM;
504 ++cc) {
505 const double col_base = t_col_base(i) * t_normal(i);
506 const double beta = alpha * row_base * col_base;
507
508 t_mat(i, j) -= beta * t_res_dt(i, j);
509
510 ++t_col_base;
511 ++t_mat;
512 }
513
514 ++t_row_base;
515 }
516 for (; rr < nb_face_functions; ++rr)
517 ++t_row_base;
518
519 ++t_disp;
520 ++t_traction;
521 ++t_coords;
522 ++t_w;
523 }
524
526}
527
528}; // namespace ContactOps
constexpr int SPACE_DIM
Kronecker Delta class.
double cn
Definition: contact.cpp:124
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
constexpr auto t_kd
const double c
speed of light (cm/ns)
FTensor::Index< 'k', SPACE_DIM > k
Definition: ContactOps.hpp:168
FTensor::Tensor1< double, 3 > grad_surface_distance_function(double t, FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_traction)
Definition: ContactOps.hpp:221
FTensor::Index< 'j', SPACE_DIM > j
Definition: ContactOps.hpp:167
FTensor::Index< 'i', SPACE_DIM > i
[Common data]
Definition: ContactOps.hpp:166
FTensor::Index< 'l', SPACE_DIM > l
Definition: ContactOps.hpp:169
double w(const double sdf, const double tn)
Definition: ContactOps.hpp:266
FTensor::Tensor2_symmetric< double, 3 > hess_surface_distance_function(double t, FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_traction)
Definition: ContactOps.hpp:239
double surface_distance_function(double t, FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_traction)
Definition: ContactOps.hpp:204
double constrain(double sdf, double tn)
Definition: ContactOps.hpp:268
double sign(double x)
Definition: ContactOps.hpp:257
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:17
def hess_sdf(t, x, y, z, tx, ty, tz)
Definition: sdf.py:22
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
MatrixDouble contactDisp
Definition: ContactOps.hpp:18
MatrixDouble contactTraction
Definition: ContactOps.hpp:17
auto contactStressDivergencePtr()
Definition: ContactOps.hpp:32
MatrixDouble contactStressDivergence
Definition: ContactOps.hpp:16
MatrixDouble contactStress
Definition: ContactOps.hpp:15
OpConstrainBoundaryLhs_dTraction(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:444
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:199
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:453
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:188
OpConstrainBoundaryLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:347
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:357
OpConstrainBoundaryRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
Definition: ContactOps.hpp:273
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:177
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Definition: ContactOps.hpp:279