v0.14.0
Loading...
Searching...
No Matches
EshelbianContact.cpp
Go to the documentation of this file.
1/**
2 * @file EshelbianContact.cpp
3 * @brief
4 * @date 2023-05-13
5 *
6 * @copyright Copyright (c) 2023
7 *
8 */
9
10extern "C" {
11void tricircumcenter3d_tp(double a[3], double b[3], double c[3],
12 double circumcenter[3], double *xi, double *eta);
13}
14
15template <>
18 [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
19 const EntitiesFieldData::EntData &row_data,
20 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
21 return MatSetValues<AssemblyTypeSelector<A>>(
22 op_ptr->getKSPA(), row_data, col_data, m, ADD_VALUES);
23 };
24
25namespace EshelbianPlasticity {
26
28 const std::string row_field_name,
29 boost::shared_ptr<ContactOps::CommonData> common_data_ptr)
30 : ContactOps::AssemblyBoundaryEleOp(row_field_name, row_field_name,
32 commonDataPtr(common_data_ptr) {}
33
34MoFEMErrorCode
35OpConstrainBoundaryRhs::iNtegrate(EntitiesFieldData::EntData &data) {
37
42
43 const size_t nb_gauss_pts = getGaussPts().size2();
44
45 // #ifndef NDEBUGbase
46 if (commonDataPtr->contactDisp.size2() != nb_gauss_pts) {
47 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
48 "Wrong number of integration pts %d != %d",
49 commonDataPtr->contactDisp.size2(), nb_gauss_pts);
50 }
51 // #endif // !NDEBUG
52
53 auto &nf = locF;
54 locF.clear();
55
56 auto t_w = getFTensor0IntegrationWeight();
57 auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
58 auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
59 auto t_coords = getFTensor1CoordsAtGaussPts();
60 auto t_normal_at_pts = getFTensor1NormalsAtGaussPts();
61
62 auto next = [&]() {
63 ++t_w;
64 ++t_disp;
65 ++t_traction;
66 ++t_coords;
67 ++t_normal_at_pts;
68 };
69
70 size_t nb_base_functions = data.getN().size2() / 3;
71 auto t_base = data.getFTensor1N<3>();
72 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
73
75 t_normal(i) = t_normal_at_pts(i) / t_normal_at_pts.l2();
76
77 auto t_nf = getFTensor1FromPtr<3>(&nf[0]);
78 const double alpha = t_w * getMeasure();
79
80 FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
81 t_spatial_coords(i) = t_coords(i) + t_disp(i);
82
83 auto ts_t = getTStime();
84
86 ts_t, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
87 t_normal_at_pts(0), t_normal_at_pts(1), t_normal_at_pts(2));
88
89 auto t_grad_sdf = gradSurfaceDistanceFunction(
90 ts_t, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
91 t_normal_at_pts(0), t_normal_at_pts(1), t_normal_at_pts(2));
92
93 auto tn = t_traction(i) * t_grad_sdf(i);
94 auto c = ContactOps::constrain(sdf, tn);
95
97 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
99 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
100
102 t_disp_e(i) = t_disp(i) - ContactOps::cn_contact * t_traction(i);
103
105 t_rhs(i) = -c * sdf * t_grad_sdf(i)
106
107 +
108
109 t_cP(i, j) * t_disp(j)
110
111 +
112
113 t_cQ(i, j) * t_disp_e(j);
114
115 size_t bb = 0;
116 for (; bb != nbRows / 3; ++bb) {
117 const double row_base = t_base(i) * t_normal(i);
118 const double beta = alpha * row_base;
119 t_nf(i) -= beta * t_rhs(i);
120 ++t_nf;
121 ++t_base;
122 }
123 for (; bb < nb_base_functions; ++bb)
124 ++t_base;
125
126 next();
127 }
128
130}
131
133 const std::string row_field_name,
134 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
135 boost::shared_ptr<EntitiesFieldData::EntData> col_data_ptr)
136 : ContactOps::AssemblyBoundaryEleOp(row_field_name, row_field_name,
137 ContactOps::BoundaryEleOp::OPROW),
138 commonDataPtr(common_data_ptr), colDataPtr(col_data_ptr) {}
139
140MoFEMErrorCode
141OpConstrainBoundaryLhs_dU::iNtegrate(EntitiesFieldData::EntData &row_data) {
143
144 using namespace ContactOps;
145
149
150 auto nb_rows = row_data.getIndices().size();
151 auto nb_cols = colDataPtr->getIndices().size();
152
153 auto &locMat = AssemblyBoundaryEleOp::locMat;
154 locMat.resize(nb_rows, nb_cols, false);
155 locMat.clear();
156
157 if (nb_cols && nb_rows) {
158
159 const size_t nb_gauss_pts = getGaussPts().size2();
160
161 // #ifndef NDEBUG
162 if (row_data.getN().size1() != nb_gauss_pts) {
163 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
164 "Wrong number of integration pts %d != %d",
165 row_data.getN().size1(), nb_gauss_pts);
166 }
167 if (row_data.getN().size2() < nb_rows / 3) {
168 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
169 "Wrong number of base functions", row_data.getN().size2(),
170 nb_rows / 3);
171 }
172 if (colDataPtr->getN().size1() != nb_gauss_pts) {
173 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
174 "Wrong number of integration pts %d != %d",
175 colDataPtr->getN().size1(), nb_gauss_pts);
176 }
177 if (colDataPtr->getN().size2() < nb_cols / 3) {
178 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
179 "Wrong number of base functions", colDataPtr->getN().size2(),
180 nb_cols / 3);
181 }
182 // #endif
183
184 auto t_w = getFTensor0IntegrationWeight();
185 auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
186 auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
187 auto t_coords = getFTensor1CoordsAtGaussPts();
188 auto t_normal_at_pts = getFTensor1NormalsAtGaussPts();
189
190 auto next = [&]() {
191 ++t_w;
192 ++t_disp;
193 ++t_traction;
194 ++t_coords;
195 ++t_normal_at_pts;
196 };
197
198 auto t_row_base = row_data.getFTensor1N<3>();
199 auto nb_face_functions = row_data.getN().size2() / 3;
200
201 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
202
203 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
204
206 t_normal(i) = t_normal_at_pts(i) / t_normal_at_pts.l2();
207
208 const double alpha = t_w * getMeasure();
209
210 FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
211 t_spatial_coords(i) = t_coords(i) + t_disp(i);
212
213 auto ts_time = getTStime();
214
216 ts_time, t_spatial_coords(0), t_spatial_coords(1),
217 t_spatial_coords(2), t_normal_at_pts(0), t_normal_at_pts(1),
218 t_normal_at_pts(2));
219 auto t_grad_sdf = gradSurfaceDistanceFunction(
220 ts_time, t_spatial_coords(0), t_spatial_coords(1),
221 t_spatial_coords(2), t_normal_at_pts(0), t_normal_at_pts(1),
222 t_normal_at_pts(2));
223
224 auto t_hess_sdf = hessSurfaceDistanceFunction(
225 ts_time, t_spatial_coords(0), t_spatial_coords(1),
226 t_spatial_coords(2), t_normal_at_pts(0), t_normal_at_pts(1),
227 t_normal_at_pts(2));
228
229 auto tn = t_traction(i) * t_grad_sdf(i);
230 auto c = constrain(sdf, tn);
231
233 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
235 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
236
238 t_disp_e(i) = t_disp(i) - ContactOps::cn_contact * t_traction(i);
239
241 t_res_dU(i, j) =
242
243 -c * sdf * t_hess_sdf(i, j) +
244 c * (t_hess_sdf(i, j) * t_grad_sdf(k) * t_disp(k) +
245 t_grad_sdf(i) * t_hess_sdf(k, j) * t_disp(k))
246
247 +
248
249 t_cQ(i, j)
250
251 -
252
253 c * (t_hess_sdf(i, j) * t_grad_sdf(k) * t_disp_e(k) +
254 t_grad_sdf(i) * t_hess_sdf(k, j) * t_disp_e(k));
255
256 size_t rr = 0;
257 for (; rr != nb_rows / 3; ++rr) {
258
259 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
260 const double row_base = t_row_base(i) * t_normal(i);
261 auto t_col_base = colDataPtr->getFTensor0N(gg, 0);
262
263 for (size_t cc = 0; cc != nb_cols / 3; ++cc) {
264 const double beta = alpha * row_base * t_col_base;
265 t_mat(i, j) -= beta * t_res_dU(i, j);
266 ++t_col_base;
267 ++t_mat;
268 }
269
270 ++t_row_base;
271 }
272 for (; rr < nb_face_functions; ++rr)
273 ++t_row_base;
274
275 next();
276 }
277 }
278
280}
281
282MoFEMErrorCode
283OpConstrainBoundaryLhs_dU::aSsemble(EntitiesFieldData::EntData &row_data) {
285 auto &locMat = ContactOps::AssemblyBoundaryEleOp::locMat;
286 CHKERR matSetValuesHook(this, row_data, *colDataPtr, locMat);
288}
289
291 const std::string row_field_name,
292 boost::shared_ptr<ContactOps::CommonData> common_data_ptr)
293 : ContactOps::AssemblyBoundaryEleOp(row_field_name, row_field_name,
294 ContactOps::BoundaryEleOp::OPROWCOL),
295 commonDataPtr(common_data_ptr) {
296 sYmm = false;
297}
298
299MoFEMErrorCode
300OpConstrainBoundaryLhs_dP::iNtegrate(EntitiesFieldData::EntData &row_data,
301 EntitiesFieldData::EntData &col_data) {
303
307
308 const size_t nb_gauss_pts = getGaussPts().size2();
309
310 auto t_w = getFTensor0IntegrationWeight();
311 auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
312 auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
313 auto t_coords = getFTensor1CoordsAtGaussPts();
314 auto t_normal_at_pts = getFTensor1NormalsAtGaussPts();
315
316 auto next = [&]() {
317 ++t_w;
318 ++t_disp;
319 ++t_traction;
320 ++t_coords;
321 ++t_normal_at_pts;
322 };
323
324 auto t_row_base = row_data.getFTensor1N<3>();
325 size_t nb_face_functions = row_data.getN().size2() / 3;
326
327 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
328
330 t_normal(i) = t_normal_at_pts(i) / t_normal_at_pts.l2();
331
332 const double alpha = t_w * getMeasure();
333
334 FTensor::Tensor1<double, 3> t_spatial_coords{0., 0., 0.};
335 t_spatial_coords(i) = t_coords(i) + t_disp(i);
336
337 auto ts_time = getTStime();
338
339 auto sdf = surfaceDistanceFunction(ts_time, t_spatial_coords(0),
340 t_spatial_coords(1), t_spatial_coords(2),
341 0, 0, 0);
342 auto t_grad_sdf = gradSurfaceDistanceFunction(
343 ts_time, t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2),
344 t_normal_at_pts(0), t_normal_at_pts(1), t_normal_at_pts(2));
345
346 auto tn = t_traction(i) * t_grad_sdf(i);
347 auto c = ContactOps::constrain(sdf, tn);
348
350 t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
352 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
353
355 t_res_dt(i, j) = -ContactOps::cn_contact * t_cQ(i, j);
356
357 size_t rr = 0;
358 for (; rr != nbRows / 3; ++rr) {
359
360 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
361 const double row_base = t_row_base(i) * t_normal(i);
362
363 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
364 for (size_t cc = 0; cc != nbCols / 3; ++cc) {
365 const double col_base = t_col_base(i) * t_normal(i);
366 const double beta = alpha * row_base * col_base;
367
368 t_mat(i, j) -= beta * t_res_dt(i, j);
369
370 ++t_col_base;
371 ++t_mat;
372 }
373
374 ++t_row_base;
375 }
376 for (; rr < nb_face_functions; ++rr)
377 ++t_row_base;
378
379 next();
380 }
381
383}
384
386 const std::string field_name,
387 boost::shared_ptr<EntitiesFieldData::EntData> col_data_ptr)
388 : UDO(field_name, UDO::OPCOL), colDataPtr(col_data_ptr) {
389 std::fill(&UDO::doEntities[MBVERTEX], &UDO::doEntities[MBMAXTYPE], false);
390 for (EntityType t = CN::TypeDimensionMap[3].first;
391 t <= CN::TypeDimensionMap[3].second; ++t) {
392 UDO::doEntities[t] = true;
393 }
394}
395
396MoFEMErrorCode
398 EntitiesFieldData::EntData &data) {
400 colDataPtr->getFieldEntities() = data.getFieldEntities();
401 colDataPtr->getIndices() = data.getIndices();
402 colDataPtr->getFieldData() = data.getFieldData();
403 colDataPtr->getN() = data.getN();
405}
406
408 boost::shared_ptr<moab::Core> core_mesh_ptr,
409 int max_order, std::map<int, Range> &&body_map)
410 : Base(m_field, core_mesh_ptr, "contact"), maxOrder(max_order),
411 bodyMap(body_map) {
412
413 auto ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTRI>>();
414 ref_ele_ptr->hoNodes =
415 PETSC_FALSE; ///< So far only linear geometry is implememted for contact
416
417 CHK_THROW_MESSAGE(ref_ele_ptr->getOptions(optionsPrefix), "getOptions");
418 CHK_THROW_MESSAGE(ref_ele_ptr->generateReferenceElementMesh(),
419 "Error when generating reference element");
420
421 // MOFEM_LOG("EP", Sev::inform) << "Contact hoNodes " << ref_ele_ptr->hoNodes
422 // ? "true"
423 // : "false";
424 // MOFEM_LOG("EP", Sev::inform)
425 // << "Contact maxOrder " << ref_ele_ptr->defMaxLevel;
426
427 refElementsMap[MBTRI] = ref_ele_ptr;
428
429 int def_ele_id = -1;
430 CHKERR getPostProcMesh().tag_get_handle("ELE_ID", 1, MB_TYPE_INTEGER, thEleId,
431 MB_TAG_CREAT | MB_TAG_DENSE,
432 &def_ele_id);
433 CHKERR getPostProcMesh().tag_get_handle("BODY_ID", 1, MB_TYPE_INTEGER,
434 thBodyId, MB_TAG_CREAT | MB_TAG_DENSE,
435 &def_ele_id);
436
437 const auto nb_nodes = (refElementsMap.at(MBTRI)->hoNodes) ? 6 : 3;
438 const auto nb_bases = NBVOLUMETET_L2(maxOrder);
439
440 std::vector<int> def_ids(3 * nb_bases, 0);
441 CHKERR getPostProcMesh().tag_get_handle("IDS", 3 * nb_bases, MB_TYPE_INTEGER,
442 thIds, MB_TAG_CREAT | MB_TAG_DENSE,
443 &*def_ids.begin());
444 std::vector<double> def_coeffs(3 * nb_bases, 0);
445 CHKERR getPostProcMesh().tag_get_handle(
446 "COEFFS", 3 * nb_bases, MB_TYPE_DOUBLE, thCoeffs,
447 MB_TAG_CREAT | MB_TAG_DENSE, &*def_coeffs.begin());
448 std::vector<double> def_basses(nb_bases, 0);
449 CHKERR getPostProcMesh().tag_get_handle("BASES", nb_bases, MB_TYPE_DOUBLE,
450 thBases, MB_TAG_CREAT | MB_TAG_DENSE,
451 &*def_basses.begin());
452
453 std::vector<double> def_small_x_h1(3, 0);
454 CHKERR getPostProcMesh().tag_get_handle("U_H1", 3, MB_TYPE_DOUBLE, thSmallXH1,
455 MB_TAG_CREAT | MB_TAG_DENSE,
456 &*def_small_x_h1.begin());
457}
458
459MoFEMErrorCode ContactTree::preProcess() {
461 CHKERR Base::preProcess();
463}
464
465MoFEMErrorCode ContactTree::postProcess() {
467
468 CHKERR Base::postProcess();
469 shadowDataMap.clear();
470
471 PetscBarrier(nullptr);
472
473 ParallelComm *pcomm_post_proc_mesh =
474 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
475 if (pcomm_post_proc_mesh == nullptr)
476 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "PComm not allocated");
477
478 // MOFEM_LOG("EP", Sev::verbose) << "ContactTree brodact entities";
479
480 auto brodacts = [&](auto &brodacts_ents) {
482 Range recived_ents;
483 for (auto rr = 0; rr != mField.get_comm_size(); ++rr) {
484 pcomm_post_proc_mesh->broadcast_entities(
485 rr, rr == mField.get_comm_rank() ? brodacts_ents : recived_ents,
486 false, true);
487 }
489 };
490
491 Range brodacts_ents;
492 CHKERR getPostProcMesh().get_entities_by_dimension(0, 2, brodacts_ents, true);
493 CHKERR brodacts(brodacts_ents);
494
495 // MOFEM_LOG("EP", Sev::verbose) << "ContactTree build tree";
496
497 Range ents;
498 CHKERR getPostProcMesh().get_entities_by_dimension(0, 2, ents, true);
499 CHKERR buildTree(ents);
500
501 CHKERR writeFile("debug_tree.h5m");
502
504}
505
506MoFEMErrorCode ContactTree::buildTree(Range &ents) {
508 // MOFEM_LOG("SELF", Sev::inform) << ents << endl;
509
510 treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
511 new OrientedBoxTreeTool(&getPostProcMesh(), "ROOTSETSURF", true));
512 CHKERR treeSurfPtr->build(ents, rootSetSurf);
513
514 // CHKERR treeSurfPtr->stats(rootSetSurf, std::cerr);
515
517}
518
520 boost::shared_ptr<ContactTree> contact_tree_ptr,
521 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
522 boost::shared_ptr<MatrixDouble> u_h1_ptr,
523 boost::shared_ptr<EntitiesFieldData::EntData> col_data_ptr)
524 : UOP(NOSPACE, UOP::OPSPACE), contactTreePtr(contact_tree_ptr),
525 uH1Ptr(u_h1_ptr), commonDataPtr(common_data_ptr),
526 colDataPtr(col_data_ptr) {}
527
528MoFEMErrorCode OpMoveNode::doWork(int side, EntityType type,
529 EntitiesFieldData::EntData &data) {
531 auto &moab_post_proc_mesh = contactTreePtr->getPostProcMesh();
532 auto &post_proc_ents = contactTreePtr->getPostProcElements();
533
534 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
535 auto fe_id = id_from_handle(fe_ent);
536 auto &map_gauss_pts = contactTreePtr->getMapGaussPts();
537
538 auto get_body_id = [&](auto fe_ent) {
539 for (auto &m : contactTreePtr->bodyMap) {
540 if (m.second.find(fe_ent) != m.second.end()) {
541 return m.first;
542 }
543 }
544 return -1;
545 };
546
547 CHKERR moab_post_proc_mesh.tag_clear_data(contactTreePtr->thEleId,
548 post_proc_ents, &fe_id);
549 auto body_id = get_body_id(fe_ent);
550 CHKERR moab_post_proc_mesh.tag_clear_data(contactTreePtr->thBodyId,
551 post_proc_ents, &body_id);
552
553 const auto nb_bases = NBVOLUMETET_L2(contactTreePtr->maxOrder);
554 VectorInt ids(3 * nb_bases, 0);
555 std::copy(colDataPtr->getIndices().begin(), colDataPtr->getIndices().end(),
556 ids.begin());
557 CHKERR moab_post_proc_mesh.tag_clear_data(contactTreePtr->thIds,
558 post_proc_ents, &*ids.begin());
559 VectorDouble coeffs(3 * nb_bases, 0);
560 std::copy(colDataPtr->getFieldData().begin(),
561 colDataPtr->getFieldData().end(), coeffs.begin());
562 CHKERR moab_post_proc_mesh.tag_clear_data(contactTreePtr->thCoeffs,
563 post_proc_ents, &*coeffs.begin());
564
565 auto nb_gauss_pts = getGaussPts().size2();
566 auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
567 auto t_u_h1 = getFTensor1FromMat<3>(*uH1Ptr);
568 auto t_coords = getFTensor1CoordsAtGaussPts();
569
570 MatrixDouble x_l2(nb_gauss_pts, 3);
572 &x_l2(0, 0), &x_l2(0, 1), &x_l2(0, 2)};
573 MatrixDouble x_h1(nb_gauss_pts, 3);
575 &x_h1(0, 0), &x_h1(0, 1), &x_h1(0, 2)};
576
578 auto nb_cols = colDataPtr->getIndices().size();
579 auto t_col_base = colDataPtr->getFTensor0N();
580
581 VectorDouble bases(nb_bases, 0);
582 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
583 t_x_l2(i) = t_coords(i) + t_disp(i);
584 t_x_h1(i) = t_coords(i) + t_u_h1(i);
585
586 for (size_t cc = 0; cc != nb_cols / 3; ++cc) {
587 bases[cc] = t_col_base;
588 ++t_col_base;
589 }
590 CHKERR moab_post_proc_mesh.tag_set_data(
591 contactTreePtr->thBases, &map_gauss_pts[gg], 1, &*bases.begin());
592
593 ++t_coords;
594 ++t_disp;
595 ++t_u_h1;
596 ++t_x_l2;
597 ++t_x_h1;
598 }
599
600 CHKERR moab_post_proc_mesh.set_coords(
601 &*map_gauss_pts.begin(), map_gauss_pts.size(), &*x_l2.data().begin());
602 CHKERR moab_post_proc_mesh.tag_set_data(
603 contactTreePtr->thSmallXH1, &*map_gauss_pts.begin(), map_gauss_pts.size(),
604 &*x_h1.data().begin());
605
607}
608
610 boost::shared_ptr<ContactTree> contact_tree_ptr,
611 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
612
613 moab::Interface *post_proc_mesh_ptr,
614 std::vector<EntityHandle> *map_gauss_pts_ptr
615
616 )
617 : UOP(NOSPACE, UOP::OPSPACE), contactTreePtr(contact_tree_ptr),
618 commonDataPtr(common_data_ptr), postProcMeshPtr(post_proc_mesh_ptr),
619 mapGaussPtsPtr(map_gauss_pts_ptr) {}
620
621MoFEMErrorCode OpTreeSearch::doWork(int side, EntityType type,
622 EntitiesFieldData::EntData &data) {
624
625 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
626 auto fe_id = id_from_handle(fe_ent);
627
629
630 const auto nb_gauss_pts = getGaussPts().size2();
631
632 auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
633 auto t_coords = getFTensor1CoordsAtGaussPts();
634
635 auto next = [&]() {
636 ++t_disp;
637 ++t_coords;
638 };
639
640 auto get_ele_centre = [i, this](auto t_ele_coords) {
641 FTensor::Tensor1<double, 3> t_ele_center;
642 t_ele_center(i) = 0;
643 for (int nn = 0; nn != 3; nn++) {
644 t_ele_center(i) += t_ele_coords(i);
645 ++t_ele_coords;
646 }
647 t_ele_center(i) /= 3;
648 return t_ele_center;
649 };
650
651 auto get_ele_radius = [this, i](auto t_ele_cenetr, auto t_ele_coords) {
653 t_n0(i) = t_ele_cenetr(i) - t_ele_coords(i);
654 return t_n0.l2();
655 };
656
657 auto get_face_conn = [this](auto face) {
658 const EntityHandle *conn;
659 int num_nodes;
660 CHK_MOAB_THROW(contactTreePtr->getPostProcMesh().get_connectivity(
661 face, conn, num_nodes, true),
662 "get conn");
663 return conn;
664 };
665
666 auto get_face_coords = [this](auto conn) {
667 std::array<double, 9> coords;
668 CHKERR contactTreePtr->getPostProcMesh().get_coords(conn, 3, coords.data());
669 return coords;
670 };
671
672 auto get_closet_face = [this](auto *point_ptr) {
673 FTensor::Tensor1<double, 3> t_point_out;
674 EntityHandle face_out;
675 CHK_MOAB_THROW(contactTreePtr->getTreeSurfPtr()->closest_to_location(
676 point_ptr, contactTreePtr->getRootSetSurf(),
677 &t_point_out(0), face_out),
678 "get cloasest faces");
679 return std::make_pair(face_out, t_point_out);
680 };
681
682 auto get_faces_out = [this](auto *point_ptr, auto *unit_ray_ptr, auto radius,
683 auto eps) {
684 std::vector<double> distances_out;
685 std::vector<EntityHandle> faces_out;
687
688 contactTreePtr->getTreeSurfPtr()->ray_intersect_triangles(
689 distances_out, faces_out, contactTreePtr->getRootSetSurf(), eps,
690 point_ptr, unit_ray_ptr, &radius),
691
692 "get cloasest faces");
693 return std::make_pair(faces_out, distances_out);
694 };
695
696 auto get_normal = [](auto &ele_coords) {
698 Tools::getTriNormal(ele_coords.data(), &t_normal(0));
699 return t_normal;
700 };
701
702 auto make_map = [](auto &face_out, auto &face_dist) {
703 std::map<double, EntityHandle> m;
704 for (auto ii = 0; ii != face_out.size(); ++ii) {
705 m[face_dist[ii]] = face_out[ii];
706 }
707 return m;
708 };
709
710 auto get_tag_data = [this](auto tag, auto face, auto &vec) {
712 int tag_length;
713 CHKERR contactTreePtr->getPostProcMesh().tag_get_length(tag, tag_length);
714 vec.resize(tag_length);
715 CHKERR contactTreePtr->getPostProcMesh().tag_get_data(tag, &face, 1,
716 &*vec.begin());
718 };
719
720 auto create_tag = [this](const std::string tag_name, const int size) {
721 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
722 Tag th;
723 if (postProcMeshPtr) {
724 CHKERR postProcMeshPtr->tag_get_handle(
725 tag_name.c_str(), size, MB_TYPE_DOUBLE, th,
726 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
727 }
728 return th;
729 };
730
731 auto set_float_precision = [](const double x) {
732 if (std::abs(x) < std::numeric_limits<float>::epsilon())
733 return 0.;
734 else
735 return x;
736 };
737
738 // scalars
739 auto save_scal_tag = [&](auto &th, auto v, const int gg) {
741 if (postProcMeshPtr) {
742 v = set_float_precision(v);
743 CHKERR postProcMeshPtr->tag_set_data(th, &(*mapGaussPtsPtr)[gg], 1, &v);
744 }
746 };
747
748 // vectors
749 VectorDouble3 v(3);
751 auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
753 if (postProcMeshPtr) {
754 t_v(i) = t_d(i);
755 for (auto &a : v.data())
756 a = set_float_precision(a);
757 CHKERR postProcMeshPtr->tag_set_data(th, &(*mapGaussPtsPtr)[gg], 1,
758 &*v.data().begin());
759 }
761 };
762
763 Tag th_gap = create_tag("contact_gap", 1);
764 Tag th_dist = create_tag("contact_dip", 3);
765 Tag th_test = create_tag("contact_test", 3);
766
767 auto t_centre = get_ele_centre(getFTensor1Coords());
768 auto radius = get_ele_radius(t_centre, getFTensor1Coords());
769
770 // MOFEM_LOG("EP", Sev::inform) << "radius " << radius;
771 // MOFEM_LOG("EP", Sev::inform) << "t_centre " << t_centre;
772 // MOFEM_LOG("EP", Sev::inform) << "elem coords " << getCoords();
773
774 {
776 auto t_coords = getFTensor1Coords();
777 t_center_a(i) = 0;
778 for (int nn = 0; nn != 3; nn++) {
779 t_center_a(i) += t_coords(i);
780 ++t_coords;
781 }
782 t_center_a(i) /= 3;
783 // MOFEM_LOG("EP", Sev::inform) << "t_centre_a " << t_center_a;
784 }
785
786 auto [face_close, t_close_global_coords] = get_closet_face(&t_centre(0));
787 auto face_close_conn = get_face_conn(face_close);
788 auto face_close_coords = get_face_coords(face_close_conn);
789 auto t_normal_face_close = get_normal(face_close_coords);
790 t_normal_face_close.normalize();
791
792 FTensor::Tensor1<double, 3> t_ray_point;
793 t_ray_point(i) = t_close_global_coords(i) - t_normal_face_close(i) * radius;
794
795 auto &shadow_vec = contactTreePtr->shadowDataMap[fe_ent];
796 shadow_vec.resize(nb_gauss_pts);
797 for (auto &shadow : shadow_vec) {
798 shadow.gap = radius;
799 shadow.tN(i) = 0;
800 }
801
802 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
803
804 FTensor::Tensor1<double, 3> t_spatial_coords;
805 t_spatial_coords(i) = t_coords(i) + t_disp(i);
806
808 t_unit_ray(i) = t_spatial_coords(i) - t_ray_point(i);
809 t_unit_ray.normalize();
810
811 constexpr double eps = 1e-6;
812 auto [faces_out, faces_dist] =
813 get_faces_out(&t_ray_point(0), &t_unit_ray(0), radius, eps * radius);
814
815 auto m = make_map(faces_out, faces_dist);
816 for (auto m_it = m.begin(); m_it != m.end(); ++m_it) {
817 auto face = m_it->second;
818 if(face != face_close) {
819 FTensor::Tensor1<double, 3> t_global_coords;
820 t_global_coords(i) = t_ray_point(i) + m_it->first * t_unit_ray(i);
821 t_global_coords(i) -= t_spatial_coords(i);
822 shadow_vec[gg].gap = t_global_coords(i) * t_unit_ray(i);
823 auto face_conn = get_face_conn(face);
824 auto face_coords = get_face_coords(face_conn);
825 auto t_normal_face = get_normal(face_coords);
826 t_normal_face.normalize();
827 shadow_vec[gg].tN(i) = t_normal_face(i);
828
829 CHKERR get_tag_data(contactTreePtr->thIds, face,
830 shadow_vec[gg].dofsIds);
831 CHKERR get_tag_data(contactTreePtr->thCoeffs, face,
832 shadow_vec[gg].dofsCoeffs);
833 CHKERR get_tag_data(contactTreePtr->thBases, face,
834 shadow_vec[gg].baseFuncs);
835
836 if (postProcMeshPtr) {
837 CHKERR save_scal_tag(th_gap, shadow_vec[gg].gap, gg);
838 CHKERR save_vec_tag(th_dist, t_global_coords, gg);
839 }
840
841 break;
842 }
843 }
844
845 next();
846 }
847
849}
850
851} // namespace EshelbianPlasticity
void tricircumcenter3d_tp(double a[3], double b[3], double c[3], double circumcenter[3], double *xi, double *eta)
constexpr double a
static const double eps
Kronecker Delta class.
Tensor1< T, Tensor_Dim > normalize()
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
@ NOSPACE
Definition: definitions.h:83
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#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_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'm', SPACE_DIM > m
static const double face_coords[4][9]
constexpr auto t_kd
double eta
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double constrain(double sdf, double tn)
constrain function
Definition: ContactOps.hpp:572
Definition: sdf.py:1
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
PostProcBrokenMeshInMoabBase< FaceElementForcesAndSourcesCore > Base
MoFEMErrorCode buildTree(Range &ents)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode postProcess()
function is run at the end of loop
ContactTree(MoFEM::Interface &m_field, boost::shared_ptr< moab::Core > core_mesh_ptr, int max_order, std::map< int, Range > &&body_map)
std::map< EntityHandle, std::vector< FaceData > > shadowDataMap
ContactOps::SurfaceDistanceFunction surfaceDistanceFunction
ContactOps::GradSurfaceDistanceFunction gradSurfaceDistanceFunction
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpConstrainBoundaryLhs_dP(const std::string row_field_name, boost::shared_ptr< ContactOps::CommonData > common_data_ptr)
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
ContactOps::SurfaceDistanceFunction surfaceDistanceFunction
ContactOps::GradSurfaceDistanceFunction gradSurfaceDistanceFunction
boost::shared_ptr< EntitiesFieldData::EntData > colDataPtr
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
OpConstrainBoundaryLhs_dU(const std::string row_field_name, boost::shared_ptr< ContactOps::CommonData > common_data_ptr, boost::shared_ptr< EntitiesFieldData::EntData > col_data_ptr)
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data)
ContactOps::HessSurfaceDistanceFunction hessSurfaceDistanceFunction
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
ContactOps::SurfaceDistanceFunction surfaceDistanceFunction
ContactOps::GradSurfaceDistanceFunction gradSurfaceDistanceFunction
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
OpConstrainBoundaryRhs(const std::string row_field_name, boost::shared_ptr< ContactOps::CommonData > common_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
ForcesAndSourcesCore::UserDataOperator UDO
OpLoopSideGetDataForSideEle(const std::string field_name, boost::shared_ptr< EntitiesFieldData::EntData > col_data_ptr)
boost::shared_ptr< EntitiesFieldData::EntData > colDataPtr
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
boost::shared_ptr< ContactTree > contactTreePtr
OpMoveNode(boost::shared_ptr< ContactTree > contact_tree_ptr, boost::shared_ptr< ContactOps::CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > u_h1_ptr, boost::shared_ptr< EntitiesFieldData::EntData > col_data_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< EntitiesFieldData::EntData > colDataPtr
FaceElementForcesAndSourcesCore::UserDataOperator UOP
boost::shared_ptr< MatrixDouble > uH1Ptr
OpTreeSearch(boost::shared_ptr< ContactTree > contact_tree_ptr, boost::shared_ptr< ContactOps::CommonData > common_data_ptr, moab::Interface *post_proc_mesh_ptr, std::vector< EntityHandle > *map_gauss_pts_ptr)
boost::shared_ptr< ContactTree > contactTreePtr
std::vector< EntityHandle > * mapGaussPtsPtr
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
FaceElementForcesAndSourcesCore::UserDataOperator UOP
virtual int get_comm_size() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
MatrixDouble locMat
local entity block matrix
boost::function< MoFEMErrorCode(ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
std::string optionsPrefix
Prefix for options.
auto & getPostProcMesh()
Get postprocessing mesh.
std::map< EntityType, PostProcGenerateRefMeshPtr > refElementsMap