v0.15.4
Loading...
Searching...
No Matches
EshelbianFracture.cpp
Go to the documentation of this file.
1/**
2 * @file EshelbianFracture.cpp
3 * @author your name (you@domain.com)
4 * @brief
5 * @version 0.1
6 * @date 2025-12-21
7 *
8 * @copyright Copyright (c) 2025
9 *
10 */
11
12namespace EshelbianPlasticity {
13
14MoFEMErrorCode EshelbianCore::calculateFaceMaterialForce(const int tag, TS ts) {
15 MoFEMFunctionBegin;
16
17 constexpr bool debug = false;
18
19 auto get_tags_vec = [&](std::vector<std::pair<std::string, int>> names) {
20 std::vector<Tag> tags;
21 tags.reserve(names.size());
22 auto create_and_clean = [&]() {
23 MoFEMFunctionBegin;
24 for (auto n : names) {
25 tags.push_back(Tag());
26 auto &tag = tags.back();
27 auto &moab = mField.get_moab();
28 auto rval = moab.tag_get_handle(n.first.c_str(), tag);
29 if (rval == MB_SUCCESS) {
30 moab.tag_delete(tag);
31 }
32 double def_val[] = {0., 0., 0.};
33 CHKERR moab.tag_get_handle(n.first.c_str(), n.second, MB_TYPE_DOUBLE,
34 tag, MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
35 }
36 MoFEMFunctionReturn(0);
37 };
38 CHK_THROW_MESSAGE(create_and_clean(), "create_and_clean");
39 return tags;
40 };
41
42 enum ExhangeTags { MATERIALFORCE, AREAGROWTH, GRIFFITHFORCE, FACEPRESSURE };
43
44 auto tags = get_tags_vec({{"MaterialForce", 3},
45 {"AreaGrowth", 3},
46 {"GriffithForce", 1},
47 {"FacePressure", 1}});
48
49 auto calculate_material_forces = [&]() {
50 MoFEMFunctionBegin;
51
52 /**
53 * @brief Create element to integration faces energies
54 */
55 auto get_face_material_force_fe = [&]() {
56 using FaceEle = FaceElementForcesAndSourcesCore;
57 auto fe_ptr = boost::make_shared<FaceEle>(mField);
58 fe_ptr->getRuleHook = [](int, int, int) { return -1; };
59 fe_ptr->setRuleHook =
60 SetIntegrationAtFrontFace(frontVertices, frontAdjEdges);
61
62 // hybrid disp, evalated on face first
63 EshelbianPlasticity::AddHOOps<2, 2, 3>::add(
64 fe_ptr->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
65 fe_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
66 hybridSpatialDisp, dataAtPts->getHybridDispAtPts()));
67 fe_ptr->getOpPtrVector().push_back(
68 new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
69 hybridSpatialDisp, dataAtPts->getGradHybridDispAtPts()));
70 auto op_loop_domain_side =
71 new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
72 mField, elementVolumeName, SPACE_DIM, Sev::noisy);
73 fe_ptr->getOpPtrVector().push_back(op_loop_domain_side);
74 fe_ptr->getOpPtrVector().push_back(new OpFaceMaterialForce(dataAtPts));
75
76 // evaluated in side domain, that is op_loop_domain_side
77 op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
78 boost::make_shared<CGGUserPolynomialBase>();
79
80 EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
81 op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
82 materialH1Positions, frontAdjEdges, nullptr, nullptr, nullptr);
83 op_loop_domain_side->getOpPtrVector().push_back(
84 new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
85 piolaStress, dataAtPts->getApproxPAtPts()));
86 op_loop_domain_side->getOpPtrVector().push_back(
87 new OpCalculateHTensorTensorField<SPACE_DIM, SPACE_DIM>(
88 bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
89
90 op_loop_domain_side->getOpPtrVector().push_back(
91 new OpCalculateVectorFieldValues<SPACE_DIM>(
92 rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
93 if (noStretch || materialModel == MaterialModel::Hencky) {
94 // We have to use actual strains to evaluate J integral and energy,
95 // in this case. Note actual stresses, and actual energy can only drive
96 // crack growth
97
98 op_loop_domain_side->getOpPtrVector().push_back(
99 physicalEquations->returnOpCalculateStretchFromStress(
100 dataAtPts, physicalEquations));
101 } else {
102 // That will not work for problem with internal stress or strain, since
103 // we approximate mechanical stretch, not actual stretch. At some point
104 // in time we can change formulation so that actual stretch is
105 // approximated. However, the way how to do it is not clear.
106
107 op_loop_domain_side->getOpPtrVector().push_back(
108 new OpCalculateTensor2SymmetricFieldValues<SPACE_DIM>(
109 stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
110 }
111
112 op_loop_domain_side->getOpPtrVector().push_back(
113 new OpFaceSideMaterialForce(dataAtPts));
114
115 return fe_ptr;
116 };
117
118 auto integrate_face_material_force_fe = [&](auto &&face_energy_fe) {
119 MoFEMFunctionBegin;
120 CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(
121 dM, skeletonElement, face_energy_fe, 0, mField.get_comm_size());
122
123 auto face_exchange = CommInterface::createEntitiesPetscVector(
124 mField.get_comm(), mField.get_moab(), 2, 3, Sev::inform);
125
126 auto print_loc_size = [this](auto v, auto str, auto sev) {
127 MoFEMFunctionBegin;
128 int size;
129 CHKERR VecGetLocalSize(v.second, &size);
130 int low, high;
131 CHKERR VecGetOwnershipRange(v.second, &low, &high);
132 MOFEM_LOG("EPSYNC", sev) << str << " local size " << size << " ( "
133 << low << " " << high << " ) ";
134 MOFEM_LOG_SEVERITY_SYNC(mField.get_comm(), sev);
135 MoFEMFunctionReturn(0);
136 };
137 CHKERR print_loc_size(face_exchange, "material face_exchange",
138 Sev::verbose);
139
140 CHKERR CommInterface::updateEntitiesPetscVector(
141 mField.get_moab(), face_exchange, tags[ExhangeTags::MATERIALFORCE]);
142 CHKERR CommInterface::updateEntitiesPetscVector(
143 mField.get_moab(), faceExchange, tags[ExhangeTags::FACEPRESSURE]);
144
145 // #ifndef NDEBUG
146 if (debug) {
147 CHKERR save_range(mField.get_moab(),
148 "front_skin_faces_material_force_" +
149 std::to_string(mField.get_comm_rank()) + ".vtk",
150 *skeletonFaces);
151 }
152 // #endif
153
154 MoFEMFunctionReturn(0);
155 };
156
157 CHKERR integrate_face_material_force_fe(get_face_material_force_fe());
158
159 MoFEMFunctionReturn(0);
160 };
161
162 auto calculate_front_material_force = [&](auto nb_J_integral_contours) {
163 MoFEMFunctionBegin;
164 FTENSOR_INDEX(SPACE_DIM, I);
165
166 auto get_conn = [&](auto e) {
167 Range conn;
168 CHK_MOAB_THROW(mField.get_moab().get_connectivity(&e, 1, conn, true),
169 "get connectivity");
170 return conn;
171 };
172
173 auto get_conn_range = [&](auto e) {
174 Range conn;
175 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
176 "get connectivity");
177 return conn;
178 };
179
180 auto get_adj = [&](auto e, auto dim) {
181 Range adj;
182 CHK_MOAB_THROW(mField.get_moab().get_adjacencies(&e, 1, dim, true, adj),
183 "get adj");
184 return adj;
185 };
186
187 auto get_adj_range = [&](auto e, auto dim) {
188 Range adj;
189 CHK_MOAB_THROW(mField.get_moab().get_adjacencies(e, dim, true, adj,
190 moab::Interface::UNION),
191 "get adj");
192 return adj;
193 };
194
195 auto get_material_force = [&](auto r, auto th) {
196 MatrixDouble material_forces(r.size(), 3, false);
197 CHK_MOAB_THROW(
198 mField.get_moab().tag_get_data(th, r, material_forces.data().data()),
199 "get data");
200 return material_forces;
201 };
202
203 if (mField.get_comm_rank() == 0) {
204
205 auto crack_edges = get_adj_range(*crackFaces, 1);
206 auto front_nodes = get_conn_range(*frontEdges);
207 auto body_edges = get_range_from_block(mField, "EDGES", 1);
208 // auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
209 // front_block_edges = subtract(front_block_edges, crack_edges);
210 // auto front_block_edges_conn = get_conn_range(front_block_edges);
211
212 // #ifndef NDEBUG
213 Range all_skin_faces;
214 Range all_front_faces;
215 // #endif
216
217 auto calculate_edge_direction = [&](auto e) {
218 const EntityHandle *conn;
219 int num_nodes;
220 CHK_MOAB_THROW(
221 mField.get_moab().get_connectivity(e, conn, num_nodes, true),
222 "get connectivity");
223 std::array<double, 6> coords;
224 CHK_MOAB_THROW(
225 mField.get_moab().get_coords(conn, num_nodes, coords.data()),
226 "get coords");
227 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p0{
228 &coords[0], &coords[1], &coords[2]};
229 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p1{
230 &coords[3], &coords[4], &coords[5]};
231 FTensor::Tensor1<double, SPACE_DIM> t_dir;
232 FTENSOR_INDEX(SPACE_DIM, i);
233 t_dir(i) = t_p1(i) - t_p0(i);
234 return t_dir;
235 };
236
237 // take bubble tets at node, and then avarage over the edges
238 auto calculate_force_through_node = [&]() {
239 MoFEMFunctionBegin;
240
241 FTENSOR_INDEX(SPACE_DIM, I);
242 FTENSOR_INDEX(SPACE_DIM, J);
243 FTENSOR_INDEX(SPACE_DIM, K);
244 FTENSOR_INDEX(SPACE_DIM, L);
245
246 Range body_ents;
247 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
248 body_ents);
249 auto body_skin = get_skin(mField, body_ents);
250 auto body_skin_conn = get_conn_range(body_skin);
251
252 // calculate nodal material force
253 for (auto n : front_nodes) {
254 auto adj_tets = get_adj(n, 3);
255 for (int ll = 0; ll < nb_J_integral_contours; ++ll) {
256 auto conn = get_conn_range(adj_tets);
257 adj_tets = get_adj_range(conn, 3);
258 }
259 auto skin_faces = get_skin(mField, adj_tets);
260 auto material_forces = get_material_force(skin_faces, tags[0]);
261
262#ifndef NDEBUG
263 if (debug) {
264 all_skin_faces.merge(skin_faces);
265 }
266#endif
267
268 auto calculate_node_material_force = [&]() {
269 auto t_face_T =
270 getFTensor1FromPtr<SPACE_DIM>(material_forces.data().data());
271 FTensor::Tensor1<double, SPACE_DIM> t_node_force{0., 0., 0.};
272 for (auto face : skin_faces) {
273
274 FTensor::Tensor1<double, SPACE_DIM> t_face_force_tmp{0., 0., 0.};
275 t_face_force_tmp(I) = t_face_T(I);
276 ++t_face_T;
277
278 auto face_tets = intersect(get_adj(face, 3), adj_tets);
279
280 if (face_tets.empty()) {
281 continue;
282 }
283
284 if (face_tets.size() != 1) {
285 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
286 "face_tets.size() != 1");
287 }
288
289 int side_number, sense, offset;
290 CHK_MOAB_THROW(mField.get_moab().side_number(face_tets[0], face,
291 side_number, sense,
292 offset),
293 "moab side number");
294 t_face_force_tmp(I) *= sense;
295 t_node_force(I) += t_face_force_tmp(I);
296 }
297
298 return t_node_force;
299 };
300
301 auto calculate_crack_area_growth_direction = [&](auto n,
302 auto &t_node_force) {
303 // if skin is on body surface, project the direction on it
304 FTensor::Tensor1<double, SPACE_DIM> t_project{0., 0., 0.};
305 auto boundary_node = intersect(Range(n, n), body_skin_conn);
306 if (boundary_node.size()) {
307 auto faces = intersect(get_adj(n, 2), body_skin);
308 for (auto f : faces) {
309 FTensor::Tensor1<double, 3> t_normal_face;
310 CHKERR mField.getInterface<Tools>()->getTriNormal(
311 f, &t_normal_face(0));
312 t_project(I) += t_normal_face(I);
313 }
314 t_project.normalize();
315 }
316
317 // calculate surface projection matrix
318 FTensor::Tensor2<double, 3, 3> t_Q;
319 auto t_kd = FTensor::Kronecker_Delta<double>();
320 t_Q(I, J) = t_kd(I, J);
321 if (boundary_node.size()) {
322 t_Q(I, J) -= t_project(I) * t_project(J);
323 }
324
325 auto adj_faces = intersect(get_adj(n, 2), *crackFaces);
326 if (adj_faces.empty()) {
327 auto adj_edges =
328 intersect(get_adj(n, 1), unite(*frontEdges, body_edges));
329 double l = 0;
330 for (auto e : adj_edges) {
331 auto t_dir = calculate_edge_direction(e);
332 l += t_dir.l2();
333 }
334 l /= 2;
335 FTensor::Tensor1<double, SPACE_DIM> t_area_dir{0., 0., 0.};
336 FTensor::Tensor1<double, SPACE_DIM> t_node_force_tmp;
337 t_node_force_tmp(I) = t_node_force(I);
338 t_node_force_tmp.normalize();
339 t_area_dir(I) = -t_node_force_tmp(I);
340 t_area_dir(I) *= l / 2;
341 return t_area_dir;
342 }
343
344 // calculate direction
345 auto front_edges = get_adj(n, 1);
346 FTensor::Tensor1<double, 3> t_area_dir{0., 0., 0.};
347 for (auto f : adj_faces) {
348 int num_nodes;
349 const EntityHandle *conn;
350 CHKERR mField.get_moab().get_connectivity(f, conn, num_nodes,
351 true);
352 std::array<double, 9> coords;
353 CHKERR mField.get_moab().get_coords(conn, num_nodes,
354 coords.data());
355 FTensor::Tensor1<double, 3> t_face_normal;
356 FTensor::Tensor3<double, 3, 3, 3> t_d_normal;
357 CHKERR mField.getInterface<Tools>()->getTriNormal(
358 coords.data(), &t_face_normal(0), &t_d_normal(0, 0, 0));
359 auto n_it = std::find(conn, conn + num_nodes, n);
360 auto n_index = std::distance(conn, n_it);
361
362 FTensor::Tensor2<double, 3, 3> t_face_hessian{
363 t_d_normal(0, n_index, 0), t_d_normal(0, n_index, 1),
364 t_d_normal(0, n_index, 2),
365
366 t_d_normal(1, n_index, 0), t_d_normal(1, n_index, 1),
367 t_d_normal(1, n_index, 2),
368
369 t_d_normal(2, n_index, 0), t_d_normal(2, n_index, 1),
370 t_d_normal(2, n_index, 2)};
371
372 FTensor::Tensor2<double, 3, 3> t_projected_hessian;
373 t_projected_hessian(I, J) =
374 t_Q(I, K) * (t_face_hessian(K, L) * t_Q(L, J));
375 t_face_normal.normalize();
376 t_area_dir(K) +=
377 t_face_normal(I) * t_projected_hessian(I, K) / 2.;
378 }
379
380 return t_area_dir;
381 };
382
383 auto t_node_force = calculate_node_material_force();
384 t_node_force(I) /= griffithEnergy; // scale all by griffith energy
385 CHK_MOAB_THROW(
386 mField.get_moab().tag_set_data(tags[ExhangeTags::MATERIALFORCE],
387 &n, 1, &t_node_force(0)),
388 "set data");
389
390 auto get_area_dir = [&]() {
391 FTensor::Tensor1<double, SPACE_DIM> t_area_dir{0., 0., 0.};
392 auto adj_edges = intersect(get_adj_range(adj_tets, 1),
393 unite(*frontEdges, body_edges));
394 auto seed_n = get_conn_range(adj_edges);
395 auto skin_adj_edges = get_skin(mField, adj_edges);
396 skin_adj_edges = subtract(skin_adj_edges, body_skin_conn);
397 seed_n = subtract(seed_n, skin_adj_edges);
398 t_area_dir(I) = 0;
399 for (auto sn : seed_n) {
400 auto t_area_dir_sn =
401 calculate_crack_area_growth_direction(sn, t_node_force);
402 t_area_dir(I) += t_area_dir_sn(I);
403 }
404 for (auto sn : skin_adj_edges) {
405 auto t_area_dir_sn =
406 calculate_crack_area_growth_direction(sn, t_node_force);
407 t_area_dir(I) += t_area_dir_sn(I) / 2;
408 }
409 return t_area_dir;
410 };
411
412 auto t_area_dir = get_area_dir();
413
414 CHK_MOAB_THROW(
415 mField.get_moab().tag_set_data(tags[ExhangeTags::AREAGROWTH], &n,
416 1, &t_area_dir(0)),
417 "set data");
418 auto griffith = -t_node_force(I) * t_area_dir(I) /
419 (t_area_dir(K) * t_area_dir(K));
420 CHK_MOAB_THROW(
421 mField.get_moab().tag_set_data(tags[ExhangeTags::GRIFFITHFORCE],
422 &n, 1, &griffith),
423 "set data");
424 }
425
426 // iterate over edges, and calculate average edge material force
427 auto ave_node_force = [&](auto th) {
428 MoFEMFunctionBegin;
429
430 for (auto e : *frontEdges) {
431
432 auto conn = get_conn(e);
433 auto data = get_material_force(conn, th);
434 auto t_node = getFTensor1FromPtr<SPACE_DIM>(data.data().data());
435 FTensor::Tensor1<double, SPACE_DIM> t_edge{0., 0., 0.};
436 for (auto n : conn) {
437 NOT_USED(n);
438 t_edge(I) += t_node(I);
439 ++t_node;
440 }
441 t_edge(I) /= conn.size();
442
443 FTensor::Tensor1<double, SPACE_DIM> t_edge_direction =
444 calculate_edge_direction(e);
445 t_edge_direction.normalize();
446
447 FTENSOR_INDEX(SPACE_DIM, I);
448 FTENSOR_INDEX(SPACE_DIM, J);
449 FTENSOR_INDEX(SPACE_DIM, K);
450 FTensor::Tensor1<double, SPACE_DIM> t_cross;
451 t_cross(K) =
452 FTensor::levi_civita(I, J, K) * t_edge_direction(I) * t_edge(J);
453 t_edge(K) = FTensor::levi_civita(I, J, K) * t_edge_direction(J) *
454 t_cross(I);
455
456 CHKERR mField.get_moab().tag_set_data(th, &e, 1, &t_edge(0));
457 }
458 MoFEMFunctionReturn(0);
459 };
460
461 // iterate over edges, and calculate average edge griffith energy
462 auto ave_node_griffith_energy = [&](auto th) {
463 MoFEMFunctionBegin;
464 for (auto e : *frontEdges) {
465 FTensor::Tensor1<double, SPACE_DIM> t_edge_force;
466 CHKERR mField.get_moab().tag_get_data(
467 tags[ExhangeTags::MATERIALFORCE], &e, 1, &t_edge_force(0));
468 FTensor::Tensor1<double, SPACE_DIM> t_edge_area_dir;
469 CHKERR mField.get_moab().tag_get_data(tags[ExhangeTags::AREAGROWTH],
470 &e, 1, &t_edge_area_dir(0));
471 double griffith_energy = -t_edge_force(I) * t_edge_area_dir(I) /
472 (t_edge_area_dir(K) * t_edge_area_dir(K));
473 CHKERR mField.get_moab().tag_set_data(th, &e, 1, &griffith_energy);
474 }
475 MoFEMFunctionReturn(0);
476 };
477
478 CHKERR ave_node_force(tags[ExhangeTags::MATERIALFORCE]);
479 CHKERR ave_node_force(tags[ExhangeTags::AREAGROWTH]);
480 CHKERR ave_node_griffith_energy(tags[ExhangeTags::GRIFFITHFORCE]);
481
482 MoFEMFunctionReturn(0);
483 };
484
485 CHKERR calculate_force_through_node();
486
487 // calculate face cross
488 for (auto e : *frontEdges) {
489 auto adj_faces = get_adj(e, 2);
490 auto crack_face = intersect(get_adj(e, 2), *crackFaces);
491
492 // #ifndef NDEBUG
493 if (debug) {
494 all_front_faces.merge(adj_faces);
495 }
496 // #endif
497
498 FTensor::Tensor1<double, SPACE_DIM> t_edge_force;
499 CHKERR mField.get_moab().tag_get_data(tags[ExhangeTags::MATERIALFORCE],
500 &e, 1, &t_edge_force(0));
501 FTensor::Tensor1<double, SPACE_DIM> t_edge_direction =
502 calculate_edge_direction(e);
503 t_edge_direction.normalize();
504 FTENSOR_INDEX(SPACE_DIM, I);
505 FTENSOR_INDEX(SPACE_DIM, J);
506 FTENSOR_INDEX(SPACE_DIM, K);
507 FTensor::Tensor1<double, SPACE_DIM> t_cross;
508 t_cross(K) = FTensor::levi_civita(I, J, K) * t_edge_direction(I) *
509 t_edge_force(J);
510
511 for (auto f : adj_faces) {
512 FTensor::Tensor1<double, SPACE_DIM> t_normal;
513 CHKERR mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
514 t_normal.normalize();
515 int side_number, sense, offset;
516 CHKERR mField.get_moab().side_number(f, e, side_number, sense,
517 offset);
518 auto dot = -sense * t_cross(I) * t_normal(I);
519 CHK_MOAB_THROW(mField.get_moab().tag_set_data(
520 tags[ExhangeTags::GRIFFITHFORCE], &f, 1, &dot),
521 "set data");
522 }
523 }
524
525#ifndef NDEBUG
526 if (debug) {
527 int ts_step;
528 CHKERR TSGetStepNumber(ts, &ts_step);
529 CHKERR save_range(mField.get_moab(),
530 "front_edges_material_force_" +
531 std::to_string(ts_step) + ".vtk",
532 *frontEdges);
533 CHKERR save_range(mField.get_moab(),
534 "front_skin_faces_material_force_" +
535 std::to_string(ts_step) + ".vtk",
536 all_skin_faces);
537 CHKERR save_range(mField.get_moab(),
538 "front_faces_material_force_" +
539 std::to_string(ts_step) + ".vtk",
540 all_front_faces);
541 }
542#endif
543 }
544
545 auto edge_exchange = CommInterface::createEntitiesPetscVector(
546 mField.get_comm(), mField.get_moab(), 1, 3, Sev::inform);
547 CHKERR CommInterface::updateEntitiesPetscVector(
548 mField.get_moab(), edge_exchange, tags[ExhangeTags::MATERIALFORCE]);
549 CHKERR CommInterface::updateEntitiesPetscVector(
550 mField.get_moab(), edge_exchange, tags[ExhangeTags::AREAGROWTH]);
551 CHKERR CommInterface::updateEntitiesPetscVector(
552 mField.get_moab(), edgeExchange, tags[ExhangeTags::GRIFFITHFORCE]);
553
554 MoFEMFunctionReturn(0);
555 };
556
557 auto print_results = [&](auto nb_J_integral_conturs) {
558 MoFEMFunctionBegin;
559
560 auto get_conn_range = [&](auto e) {
561 Range conn;
562 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
563 "get connectivity");
564 return conn;
565 };
566
567 auto get_tag_data = [&](auto &ents, auto tag, auto dim) {
568 std::vector<double> data(ents.size() * dim);
569 CHK_MOAB_THROW(mField.get_moab().tag_get_data(tag, ents, data.data()),
570 "get data");
571 return data;
572 };
573
574 if (mField.get_comm_rank() == 0) {
575 auto at_nodes = [&]() {
576 MoFEMFunctionBegin;
577 auto conn = get_conn_range(*frontEdges);
578 auto material_force =
579 get_tag_data(conn, tags[ExhangeTags::MATERIALFORCE], 3);
580 auto area_growth = get_tag_data(conn, tags[ExhangeTags::AREAGROWTH], 3);
581 auto griffith_force =
582 get_tag_data(conn, tags[ExhangeTags::GRIFFITHFORCE], 1);
583 std::vector<double> coords(conn.size() * 3);
584 CHK_MOAB_THROW(mField.get_moab().get_coords(conn, coords.data()),
585 "get coords");
586 if (conn.size())
587 MOFEM_LOG("EPSELF", Sev::inform) << "Material force at nodes";
588 for (size_t i = 0; i < conn.size(); ++i) {
589 MOFEM_LOG("EPSELF", Sev::inform)
590 << "Node " << conn[i] << " coords " << coords[i * 3 + 0] << " "
591 << coords[i * 3 + 1] << " " << coords[i * 3 + 2]
592 << " material force " << material_force[i * 3 + 0] << " "
593 << material_force[i * 3 + 1] << " " << material_force[i * 3 + 2]
594 << " area growth " << area_growth[i * 3 + 0] << " "
595 << area_growth[i * 3 + 1] << " " << area_growth[i * 3 + 2]
596 << " griffith force " << std::setprecision(12)
597 << griffith_force[i] << " contour " << nb_J_integral_conturs;
598 }
599 MoFEMFunctionReturn(0);
600 };
601
602 at_nodes();
603 }
604 MoFEMFunctionReturn(0);
605 };
606
607 CHKERR calculate_material_forces();
608
609 PetscBool all_contours = PETSC_FALSE;
610 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "",
611 "-calculate_J_integral_all_levels", &all_contours,
612 PETSC_NULLPTR); // for backward compatibility
613 CHKERR PetscOptionsGetBool(
614 PETSC_NULLPTR, "", "-calculate_J_integral_all_contours", &all_contours,
615 PETSC_NULLPTR); // new name
616
617 if (all_contours == PETSC_TRUE) {
618 for (int l = 0; l < nbJIntegralContours; ++l) {
619 CHKERR calculate_front_material_force(l);
620 CHKERR print_results(l);
621 }
622 }
623
624 CHKERR calculate_front_material_force(nbJIntegralContours);
625 CHKERR print_results(nbJIntegralContours);
626
627 MoFEMFunctionReturn(0);
628}
629
630MoFEMErrorCode EshelbianCore::calculateOrientation(const int tag,
631 bool set_orientation) {
632 MoFEMFunctionBegin;
633
634 constexpr bool debug = false;
635 constexpr auto sev = Sev::verbose;
636
637 Range body_ents;
638 CHKERR mField.get_moab().get_entities_by_dimension(0, 3, body_ents);
639 auto body_skin = get_skin(mField, body_ents);
640 Range body_skin_edges;
641 CHKERR mField.get_moab().get_adjacencies(body_skin, 1, false, body_skin_edges,
642 moab::Interface::UNION);
643 Range boundary_skin_verts;
644 CHKERR mField.get_moab().get_connectivity(body_skin_edges,
645 boundary_skin_verts, true);
646
647 auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
648 Range geometry_edges_verts;
649 CHKERR mField.get_moab().get_connectivity(geometry_edges,
650 geometry_edges_verts, true);
651 Range crack_faces_verts;
652 CHKERR mField.get_moab().get_connectivity(*crackFaces, crack_faces_verts,
653 true);
654 Range crack_faces_edges;
655 CHKERR mField.get_moab().get_adjacencies(
656 *crackFaces, 1, true, crack_faces_edges, moab::Interface::UNION);
657 Range crack_faces_tets;
658 CHKERR mField.get_moab().get_adjacencies(
659 *crackFaces, 3, true, crack_faces_tets, moab::Interface::UNION);
660
661 Range front_verts;
662 CHKERR mField.get_moab().get_connectivity(*frontEdges, front_verts, true);
663 Range front_faces;
664 CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2, true, front_faces,
665 moab::Interface::UNION);
666 Range front_verts_edges;
667 CHKERR mField.get_moab().get_adjacencies(
668 front_verts, 1, true, front_verts_edges, moab::Interface::UNION);
669
670 auto get_tags_vec = [&](auto tag_name, int dim) {
671 std::vector<Tag> tags(1);
672
673 if (dim > 3)
674 CHK_THROW_MESSAGE(MOFEM_ATOM_TEST_INVALID, "dim>3");
675
676 auto create_and_clean = [&]() {
677 MoFEMFunctionBegin;
678 auto &moab = mField.get_moab();
679 auto rval = moab.tag_get_handle(tag_name, tags[0]);
680 if (rval == MB_SUCCESS) {
681 moab.tag_delete(tags[0]);
682 }
683 double def_val[] = {0., 0., 0.};
684 CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
685 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
686 MoFEMFunctionReturn(0);
687 };
688
689 CHK_THROW_MESSAGE(create_and_clean(), "create_and_clean");
690
691 return tags;
692 };
693
694 auto get_adj_front = [&](bool subtract_crack) {
695 Range adj_front;
696 CHKERR mField.get_moab().get_adjacencies(*frontEdges, SPACE_DIM - 1, true,
697 adj_front, moab::Interface::UNION);
698 if (subtract_crack)
699 adj_front = subtract(adj_front, *crackFaces);
700 return adj_front;
701 };
702
703 MOFEM_LOG_CHANNEL("SELF");
704
705 auto th_front_position = get_tags_vec("FrontPosition", 3);
706 auto th_max_face_energy = get_tags_vec("MaxFaceEnergy", 1);
707
708 if (mField.get_comm_rank() == 0) {
709
710 auto get_crack_adj_tets = [&](auto r) {
711 Range crack_faces_conn;
712 CHKERR mField.get_moab().get_connectivity(r, crack_faces_conn);
713 Range crack_faces_conn_tets;
714 CHKERR mField.get_moab().get_adjacencies(crack_faces_conn, SPACE_DIM,
715 true, crack_faces_conn_tets,
716 moab::Interface::UNION);
717 return crack_faces_conn_tets;
718 };
719
720 auto get_layers_for_sides = [&](auto &side) {
721 std::vector<Range> layers;
722 auto get = [&]() {
723 MoFEMFunctionBegin;
724
725 auto get_adj = [&](auto &r, int dim) {
726 Range adj;
727 CHKERR mField.get_moab().get_adjacencies(r, dim, true, adj,
728 moab::Interface::UNION);
729 return adj;
730 };
731
732 auto get_tets = [&](auto r) { return get_adj(r, SPACE_DIM); };
733
734 Range front_nodes;
735 CHKERR mField.get_moab().get_connectivity(*frontEdges, front_nodes,
736 true);
737 Range front_faces = get_adj(front_nodes, 2);
738 front_faces = subtract(front_faces, *crackFaces);
739 auto front_tets = get_tets(front_nodes);
740 auto front_side = intersect(side, front_tets);
741 layers.push_back(front_side);
742 for (;;) {
743 auto adj_faces = get_skin(mField, layers.back());
744 adj_faces = intersect(adj_faces, front_faces);
745 auto adj_faces_tets = get_tets(adj_faces);
746 adj_faces_tets = intersect(adj_faces_tets, front_tets);
747 layers.push_back(unite(layers.back(), adj_faces_tets));
748 if (layers.back().size() == layers[layers.size() - 2].size()) {
749 break;
750 }
751 }
752 MoFEMFunctionReturn(0);
753 };
754 CHK_THROW_MESSAGE(get(), "get_layers_for_sides");
755 return layers;
756 };
757
758 auto sides_pair = get_two_sides_of_crack_surface(mField, *crackFaces);
759 auto layers_top = get_layers_for_sides(sides_pair.first);
760 auto layers_bottom = get_layers_for_sides(sides_pair.second);
761
762#ifndef NDEBUG
763 if (debug) {
764 CHKERR save_range(
765 mField.get_moab(),
766 "crack_tets_" +
767 boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
768 get_crack_adj_tets(*crackFaces));
769 CHKERR save_range(mField.get_moab(), "sides_first.vtk", sides_pair.first);
770 CHKERR save_range(mField.get_moab(), "sides_second.vtk",
771 sides_pair.second);
772 MOFEM_LOG("EP", sev) << "Nb. layers " << layers_top.size();
773 int l = 0;
774 for (auto &r : layers_top) {
775 MOFEM_LOG("EP", sev) << "Layer " << l << " size " << r.size();
776 CHKERR save_range(
777 mField.get_moab(),
778 "layers_top_" + boost::lexical_cast<std::string>(l) + ".vtk", r);
779 ++l;
780 }
781
782 l = 0;
783 for (auto &r : layers_bottom) {
784 MOFEM_LOG("EP", sev) << "Layer " << l << " size " << r.size();
785 CHKERR save_range(
786 mField.get_moab(),
787 "layers_bottom_" + boost::lexical_cast<std::string>(l) + ".vtk", r);
788 ++l;
789 }
790 }
791#endif
792
793 auto get_cross = [&](auto t_dir, auto f) {
794 FTensor::Tensor1<double, SPACE_DIM> t_normal;
795 CHKERR mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
796 t_normal.normalize();
797 FTensor::Tensor1<double, SPACE_DIM> t_cross;
798 FTENSOR_INDEX(SPACE_DIM, i);
799 FTENSOR_INDEX(SPACE_DIM, j);
800 FTENSOR_INDEX(SPACE_DIM, k);
801 t_cross(i) = FTensor::levi_civita(i, j, k) * t_normal(j) * t_dir(k);
802 return t_cross;
803 };
804
805 auto get_sense = [&](auto f, auto e) {
806 int side, sense, offset;
807 CHK_MOAB_THROW(mField.get_moab().side_number(f, e, side, sense, offset),
808 "get sense");
809 return std::make_tuple(side, sense, offset);
810 };
811
812 auto calculate_edge_direction = [&](auto e, auto normalize = true) {
813 const EntityHandle *conn;
814 int num_nodes;
815 CHKERR mField.get_moab().get_connectivity(e, conn, num_nodes, true);
816 std::array<double, 6> coords;
817 CHKERR mField.get_moab().get_coords(conn, num_nodes, coords.data());
818 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p0{
819 &coords[0], &coords[1], &coords[2]};
820 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p1{
821 &coords[3], &coords[4], &coords[5]};
822 FTensor::Tensor1<double, SPACE_DIM> t_dir;
823 FTENSOR_INDEX(SPACE_DIM, i);
824 t_dir(i) = t_p1(i) - t_p0(i);
825 if (normalize)
826 t_dir.normalize();
827 return t_dir;
828 };
829
830 auto evaluate_face_energy_and_set_orientation = [&](auto front_edges,
831 auto front_faces,
832 auto &sides_pair,
833 auto th_position) {
834 MoFEMFunctionBegin;
835
836 Tag th_face_energy;
837 Tag th_material_force;
838 switch (energyReleaseSelector) {
839 case GRIFFITH_FORCE:
840 case GRIFFITH_SKELETON:
841 CHKERR mField.get_moab().tag_get_handle("GriffithForce",
842 th_face_energy);
843 CHKERR mField.get_moab().tag_get_handle("MaterialForce",
844 th_material_force);
845 break;
846 default:
847 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
848 "Unknown energy release selector");
849 };
850
851 /**
852 * Iterate over front edges, get adjacent faces, find maximal face energy.
853 * Maximal face energy is stored in the edge. Maximal face energy is
854 * magnitude of edge Griffith force.
855 */
856 auto find_maximal_face_energy = [&](auto front_edges, auto front_faces,
857 auto &edge_face_max_energy_map) {
858 MoFEMFunctionBegin;
859
860 Range body_ents;
861 CHKERR mField.get_moab().get_entities_by_dimension(0, 3, body_ents);
862 auto body_skin = get_skin(mField, body_ents);
863
864 Range max_faces;
865
866 for (auto e : front_edges) {
867
868 double griffith_force;
869 CHKERR mField.get_moab().tag_get_data(th_face_energy, &e, 1,
870 &griffith_force);
871
872 Range faces;
873 CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
874 faces = subtract(intersect(faces, front_faces), body_skin);
875 std::vector<double> face_energy(faces.size());
876 CHKERR mField.get_moab().tag_get_data(th_face_energy, faces,
877 face_energy.data());
878 auto max_energy_it =
879 std::max_element(face_energy.begin(), face_energy.end());
880 double max_energy =
881 max_energy_it != face_energy.end() ? *max_energy_it : 0;
882
883 edge_face_max_energy_map[e] =
884 std::make_tuple(faces[max_energy_it - face_energy.begin()],
885 griffith_force, static_cast<double>(0));
886 MOFEM_LOG("EP", Sev::inform)
887 << "Edge " << e << " griffith force " << griffith_force
888 << " max face energy " << max_energy << " factor "
889 << max_energy / griffith_force;
890
891 max_faces.insert(faces[max_energy_it - face_energy.begin()]);
892 }
893
894#ifndef NDEBUG
895 if (debug) {
896 CHKERR save_range(
897 mField.get_moab(),
898 "max_faces_" +
899 boost::lexical_cast<std::string>(mField.get_comm_rank()) +
900 ".vtk",
901 max_faces);
902 }
903#endif
904
905 MoFEMFunctionReturn(0);
906 };
907
908 /**
909 * For each front edge, find maximal face energy and orientation. This is
910 * by finding angle between edge material force and maximal face normal
911 *
912 */
913 auto calculate_face_orientation = [&](auto &edge_face_max_energy_map) {
914 MoFEMFunctionBegin;
915
916 auto up_down_face = [&](
917
918 auto &face_angle_map_up,
919 auto &face_angle_map_down
920
921 ) {
922 MoFEMFunctionBegin;
923
924 for (auto &m : edge_face_max_energy_map) {
925 auto e = m.first;
926 auto [max_face, energy, opt_angle] = m.second;
927
928 Range faces;
929 CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
930 faces = intersect(faces, front_faces);
931 Range adj_tets; // tetrahedrons adjacent to the face
932 CHKERR mField.get_moab().get_adjacencies(&max_face, 1, SPACE_DIM,
933 false, adj_tets,
934 moab::Interface::UNION);
935 if (adj_tets.size()) {
936
937 Range adj_tets; // tetrahedrons adjacent to the face
938 CHKERR mField.get_moab().get_adjacencies(&max_face, 1, SPACE_DIM,
939 false, adj_tets,
940 moab::Interface::UNION);
941 if (adj_tets.size()) {
942
943 Range adj_tets_faces;
944 // get faces
945 CHKERR mField.get_moab().get_adjacencies(
946 adj_tets, SPACE_DIM - 1, false, adj_tets_faces,
947 moab::Interface::UNION);
948 adj_tets_faces = intersect(adj_tets_faces, faces);
949 FTENSOR_INDEX(SPACE_DIM, i);
950
951 // cross product of face normal and edge direction
952 auto t_cross_max =
953 get_cross(calculate_edge_direction(e, true), max_face);
954 auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
955 t_cross_max(i) *= sense_max;
956
957 for (auto t : adj_tets) {
958 Range adj_tets_faces;
959 CHKERR mField.get_moab().get_adjacencies(
960 &t, 1, SPACE_DIM - 1, false, adj_tets_faces);
961 adj_tets_faces = intersect(adj_tets_faces, faces);
962 adj_tets_faces =
963 subtract(adj_tets_faces, Range(max_face, max_face));
964
965 if (adj_tets_faces.size() == 1) {
966
967 // cross product of adjacent face normal and edge
968 // direction
969 auto t_cross = get_cross(calculate_edge_direction(e, true),
970 adj_tets_faces[0]);
971 auto [side, sense, offset] =
972 get_sense(adj_tets_faces[0], e);
973 t_cross(i) *= sense;
974 double dot = t_cross(i) * t_cross_max(i);
975 auto angle = std::acos(dot);
976
977 double face_energy;
978 CHKERR mField.get_moab().tag_get_data(
979 th_face_energy, adj_tets_faces, &face_energy);
980
981 auto [side_face, sense_face, offset_face] =
982 get_sense(t, max_face);
983
984 if (sense_face > 0) {
985 face_angle_map_up[e] = std::make_tuple(face_energy, angle,
986 adj_tets_faces[0]);
987
988 } else {
989 face_angle_map_down[e] = std::make_tuple(
990 face_energy, -angle, adj_tets_faces[0]);
991 }
992 }
993 }
994 }
995 }
996 }
997
998 MoFEMFunctionReturn(0);
999 };
1000
1001 auto calc_optimal_angle = [&](
1002
1003 auto &face_angle_map_up,
1004 auto &face_angle_map_down
1005
1006 ) {
1007 MoFEMFunctionBegin;
1008
1009 for (auto &m : edge_face_max_energy_map) {
1010 auto e = m.first;
1011 auto &[max_face, e0, a0] = m.second;
1012
1013 if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
1014
1015 if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
1016 face_angle_map_down.find(e) == face_angle_map_down.end()) {
1017 // Do nothing
1018 } else {
1019
1020 switch (energyReleaseSelector) {
1021 case GRIFFITH_FORCE:
1022 case GRIFFITH_SKELETON: {
1023
1024 Tag th_material_force;
1025 CHKERR mField.get_moab().tag_get_handle("MaterialForce",
1026 th_material_force);
1027 FTensor::Tensor1<double, SPACE_DIM> t_material_force;
1028 CHKERR mField.get_moab().tag_get_data(
1029 th_material_force, &e, 1, &t_material_force(0));
1030 auto material_force_magnitude = t_material_force.l2();
1031 if (material_force_magnitude <
1032 std::numeric_limits<double>::epsilon()) {
1033 a0 = 0;
1034
1035 } else {
1036
1037 auto t_edge_dir = calculate_edge_direction(e, true);
1038 auto t_cross_max = get_cross(t_edge_dir, max_face);
1039 auto [side, sense, offset] = get_sense(max_face, e);
1040 t_cross_max(sense) *= sense;
1041
1042 FTENSOR_INDEX(SPACE_DIM, I);
1043 FTENSOR_INDEX(SPACE_DIM, J);
1044 FTENSOR_INDEX(SPACE_DIM, K);
1045
1046 t_material_force.normalize();
1047 t_cross_max.normalize();
1048 FTensor::Tensor1<double, SPACE_DIM> t_cross;
1049 t_cross(I) = FTensor::levi_civita(I, J, K) *
1050 t_material_force(J) * t_cross_max(K);
1051 a0 = -std::asin(t_cross(I) * t_edge_dir(I));
1052
1053 MOFEM_LOG("EP", sev)
1054 << "Optimal angle " << a0 << " energy " << e0;
1055 }
1056 break;
1057 }
1058 default: {
1059
1060 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1061 "Unknown energy release selector");
1062 }
1063 }
1064 }
1065 }
1066 }
1067
1068 MoFEMFunctionReturn(0);
1069 };
1070
1071 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
1072 face_angle_map_up;
1073 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
1074 face_angle_map_down;
1075 CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
1076 CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
1077
1078#ifndef NDEBUG
1079 if (debug) {
1080 auto th_angle = get_tags_vec("Angle", 1);
1081 Range up;
1082 for (auto &m : face_angle_map_up) {
1083 auto [e, a, face] = m.second;
1084 up.insert(face);
1085 CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1, &a);
1086 }
1087 Range down;
1088 for (auto &m : face_angle_map_down) {
1089 auto [e, a, face] = m.second;
1090 down.insert(face);
1091 CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1, &a);
1092 }
1093
1094 Range max_energy_faces;
1095 for (auto &m : edge_face_max_energy_map) {
1096 auto [face, e, angle] = m.second;
1097 max_energy_faces.insert(face);
1098 CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1,
1099 &angle);
1100 }
1101 if (mField.get_comm_rank() == 0) {
1102 CHKERR save_range(mField.get_moab(), "up_faces.vtk", up);
1103 CHKERR save_range(mField.get_moab(), "down_faces.vtk", down);
1104 CHKERR save_range(mField.get_moab(), "max_energy_faces.vtk",
1105 max_energy_faces);
1106 }
1107 }
1108#endif // NDEBUG
1109
1110 MoFEMFunctionReturn(0);
1111 };
1112
1113 auto get_conn = [&](auto e) {
1114 Range conn;
1115 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
1116 "get conn");
1117 return conn;
1118 };
1119
1120 auto get_adj = [&](auto e, auto dim) {
1121 Range adj;
1122 CHK_MOAB_THROW(mField.get_moab().get_adjacencies(
1123 e, dim, false, adj, moab::Interface::UNION),
1124 "get adj");
1125 return adj;
1126 };
1127
1128 auto get_coords = [&](auto v) {
1129 FTensor::Tensor1<double, 3> t_coords;
1130 CHK_MOAB_THROW(mField.get_moab().get_coords(v, &t_coords(0)),
1131 "get coords");
1132 return t_coords;
1133 };
1134
1135 // calulate normal of the max energy face
1136 auto get_rotated_normal = [&](auto e, auto f, auto angle) {
1137 FTENSOR_INDEX(SPACE_DIM, i);
1138 FTENSOR_INDEX(SPACE_DIM, j);
1139 auto t_edge_dir = calculate_edge_direction(e, true);
1140 auto [side, sense, offset] = get_sense(f, e);
1141 t_edge_dir(i) *= sense;
1142 t_edge_dir.normalize();
1143 t_edge_dir(i) *= angle;
1144 auto t_R = LieGroups::SO3::exp(t_edge_dir, angle);
1145 FTensor::Tensor1<double, SPACE_DIM> t_normal;
1146 mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
1147 FTensor::Tensor1<double, SPACE_DIM> t_rotated_normal;
1148 t_rotated_normal(i) = t_R(i, j) * t_normal(j);
1149 return std::make_tuple(t_normal, t_rotated_normal);
1150 };
1151
1152 auto set_coord = [&](auto v, auto &adj_vertex_tets_verts, auto &coords,
1153 auto &t_move, auto gamma) {
1154 auto index = adj_vertex_tets_verts.index(v);
1155 if (index >= 0) {
1156 for (auto ii : {0, 1, 2}) {
1157 coords[3 * index + ii] += gamma * t_move(ii);
1158 }
1159 return true;
1160 }
1161 return false;
1162 };
1163
1164 auto tets_quality = [&](auto quality, auto &adj_vertex_tets_verts,
1165 auto &adj_vertex_tets, auto &coords) {
1166 for (auto t : adj_vertex_tets) {
1167 const EntityHandle *conn;
1168 int num_nodes;
1169 CHKERR mField.get_moab().get_connectivity(t, conn, num_nodes, true);
1170 std::array<double, 12> tet_coords;
1171 for (auto n = 0; n != 4; ++n) {
1172 auto index = adj_vertex_tets_verts.index(conn[n]);
1173 if (index < 0) {
1174 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong index");
1175 }
1176 for (auto ii = 0; ii != 3; ++ii) {
1177 tet_coords[3 * n + ii] = coords[3 * index + ii];
1178 }
1179 }
1180 double q = Tools::volumeLengthQuality(tet_coords.data());
1181 if (!std::isnormal(q))
1182 q = -2;
1183 quality = std::min(quality, q);
1184 };
1185
1186 return quality;
1187 };
1188
1189 auto calculate_free_face_node_displacement =
1190 [&](auto &edge_face_max_energy_map) {
1191 // get edges adjacent to vertex along which nodes are moving
1192 auto get_vertex_edges = [&](auto vertex) {
1193 Range vertex_edges; // edges adjacent to vertex
1194
1195 auto impl = [&]() {
1196 MoFEMFunctionBegin;
1197 CHKERR mField.get_moab().get_adjacencies(vertex, 1, false,
1198 vertex_edges);
1199 vertex_edges = subtract(vertex_edges, front_verts_edges);
1200
1201 if (boundary_skin_verts.size() &&
1202 boundary_skin_verts.find(vertex[0]) !=
1203 boundary_skin_verts.end()) {
1204 MOFEM_LOG("EP", sev) << "Boundary vertex";
1205 vertex_edges = intersect(vertex_edges, body_skin_edges);
1206 }
1207 if (geometry_edges_verts.size() &&
1208 geometry_edges_verts.find(vertex[0]) !=
1209 geometry_edges_verts.end()) {
1210 MOFEM_LOG("EP", sev) << "Geometry edge vertex";
1211 vertex_edges = intersect(vertex_edges, geometry_edges);
1212 }
1213 if (crack_faces_verts.size() &&
1214 crack_faces_verts.find(vertex[0]) !=
1215 crack_faces_verts.end()) {
1216 MOFEM_LOG("EP", sev) << "Crack face vertex";
1217 vertex_edges = intersect(vertex_edges, crack_faces_edges);
1218 }
1219 MoFEMFunctionReturn(0);
1220 };
1221
1222 CHK_THROW_MESSAGE(impl(), "get_vertex_edges");
1223
1224 return vertex_edges;
1225 };
1226
1227 // vector of rotated faces, edge along node is moved, moved edge,
1228 // moved displacement, quality, cardinality, gamma
1229 using Bundle = std::vector<
1230
1231 std::tuple<EntityHandle, EntityHandle, EntityHandle,
1232 FTensor::Tensor1<double, 3>, double, double, double>
1233
1234 >;
1235 std::map<EntityHandle, Bundle> edge_bundle_map;
1236
1237 for (auto &m : edge_face_max_energy_map) {
1238
1239 auto edge = m.first;
1240 auto &[max_face, energy, opt_angle] = m.second;
1241
1242 // calculate rotation of max energy face
1243 auto [t_normal, t_rotated_normal] =
1244 get_rotated_normal(edge, max_face, opt_angle);
1245
1246 auto front_vertex = get_conn(Range(m.first, m.first));
1247 auto adj_tets = get_adj(Range(max_face, max_face), 3);
1248 auto adj_tets_faces = get_adj(adj_tets, 2);
1249 auto adj_front_faces = subtract(
1250 intersect(get_adj(Range(edge, edge), 2), adj_tets_faces),
1251 *crackFaces);
1252 if (adj_front_faces.size() > 3)
1253 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
1254 "adj_front_faces.size()>3");
1255
1256 FTensor::Tensor1<double, SPACE_DIM> t_material_force;
1257 CHKERR mField.get_moab().tag_get_data(th_material_force, &edge, 1,
1258 &t_material_force(0));
1259 std::vector<double> griffith_energy(adj_front_faces.size());
1260 CHKERR mField.get_moab().tag_get_data(
1261 th_face_energy, adj_front_faces, griffith_energy.data());
1262
1263 auto set_edge_bundle = [&](auto min_gamma) {
1264 for (auto rotated_f : adj_front_faces) {
1265
1266 double rotated_face_energy =
1267 griffith_energy[adj_front_faces.index(rotated_f)];
1268
1269 auto vertex = subtract(get_conn(Range(rotated_f, rotated_f)),
1270 front_vertex);
1271 if (vertex.size() != 1) {
1272 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
1273 "Wrong number of vertex to move");
1274 }
1275 auto front_vertex_edges_vertex = get_conn(
1276 intersect(get_adj(front_vertex, 1), crack_faces_edges));
1277 vertex = subtract(
1278 vertex, front_vertex_edges_vertex); // vertex free to move
1279 if (vertex.empty()) {
1280 continue;
1281 }
1282
1283 auto face_cardinality = [&](auto f, auto &seen_front_edges) {
1284 auto whole_front =
1285 unite(*frontEdges,
1286 subtract(body_skin_edges, crack_faces_edges));
1287 auto faces = Range(f, f);
1288 int c = 0;
1289 for (; c < 10; ++c) {
1290 auto front_edges =
1291 subtract(get_adj(faces, 1), seen_front_edges);
1292 if (front_edges.size() == 0) {
1293 return 0;
1294 }
1295 auto front_connected_edges =
1296 intersect(front_edges, whole_front);
1297 if (front_connected_edges.size()) {
1298 seen_front_edges.merge(front_connected_edges);
1299 return c;
1300 }
1301 faces.merge(get_adj(front_edges, 2));
1302 ++c;
1303 }
1304 return c;
1305 };
1306
1307 Range seen_edges = Range(edge, edge);
1308 double rotated_face_cardinality = face_cardinality(
1309 rotated_f,
1310 seen_edges); // add cardinality of max energy
1311 // face to rotated face cardinality
1312 // rotated_face_cardinality +=
1313 // face_cardinality(max_face, seen_edges);
1314 rotated_face_cardinality = std::max(rotated_face_cardinality,
1315 1.); // at least one edge
1316
1317 auto t_vertex_coords = get_coords(vertex);
1318 auto vertex_edges = get_vertex_edges(vertex);
1319
1320 EntityHandle f0 = front_vertex[0];
1321 EntityHandle f1 = front_vertex[1];
1322 FTensor::Tensor1<double, 3> t_v_e0, t_v_e1;
1323 CHKERR mField.get_moab().get_coords(&f0, 1, &t_v_e0(0));
1324 CHKERR mField.get_moab().get_coords(&f1, 1, &t_v_e1(0));
1325
1326 FTENSOR_INDEX(SPACE_DIM, i);
1327 for (auto e_used_to_move_detection : vertex_edges) {
1328 auto edge_conn = get_conn(Range(e_used_to_move_detection,
1329 e_used_to_move_detection));
1330 edge_conn = subtract(edge_conn, vertex);
1331 // Find displacement of the edge such that dot porduct with
1332 // normal is zero.
1333 //
1334 // { (t_v0 - t_vertex_coords) + gamma * (t_v3 -
1335 // t_vertex_coords) } * n = 0
1336 // where t_v0 is the edge vertex, t_v3 is the edge end
1337 // point, n is the rotated normal of the face gamma is the
1338 // factor by which the edge is moved
1339 FTensor::Tensor1<double, 3> t_v0;
1340 t_v0(i) = (t_v_e0(i) + t_v_e1(i)) / 2;
1341 FTensor::Tensor1<double, 3> t_v3;
1342 CHKERR mField.get_moab().get_coords(edge_conn, &t_v3(0));
1343 auto a =
1344 (t_v0(i) - t_vertex_coords(i)) * t_rotated_normal(i);
1345 auto b =
1346 (t_v3(i) - t_vertex_coords(i)) * t_rotated_normal(i);
1347 auto gamma = a / b;
1348
1349 constexpr double eps =
1350 std::numeric_limits<double>::epsilon();
1351 if (std::isnormal(gamma) && gamma < 1.0 - eps &&
1352 gamma > -0.1) {
1353 FTensor::Tensor1<double, 3> t_move;
1354 t_move(i) = gamma * (t_v3(i) - t_vertex_coords(i));
1355
1356 auto check_rotated_face_directoon = [&]() {
1357 FTensor::Tensor1<double, SPACE_DIM> t_delta;
1358 t_delta(i) = t_vertex_coords(i) + t_move(i) - t_v0(i);
1359 t_delta.normalize();
1360 auto dot =
1361 (t_material_force(i) / t_material_force.l2()) *
1362 t_delta(i);
1363 return -dot > 0 ? true : false;
1364 };
1365
1366 if (check_rotated_face_directoon()) {
1367
1368 MOFEM_LOG("EP", Sev::inform)
1369 << "Crack edge " << edge << " moved face "
1370 << rotated_f
1371 << " edge: " << e_used_to_move_detection
1372 << " face direction/energy " << rotated_face_energy
1373 << " face cardinality " << rotated_face_cardinality
1374 << " gamma: " << gamma;
1375
1376 auto &bundle = edge_bundle_map[edge];
1377 bundle.emplace_back(rotated_f, e_used_to_move_detection,
1378 vertex[0], t_move, 1,
1379 rotated_face_cardinality, gamma);
1380 }
1381 }
1382 }
1383 }
1384 };
1385
1386 set_edge_bundle(std::numeric_limits<double>::epsilon());
1387 if (edge_bundle_map[edge].empty()) {
1388 set_edge_bundle(-1.);
1389 }
1390 }
1391
1392 return edge_bundle_map;
1393 };
1394
1395 auto get_sort_by_energy = [&](auto &edge_face_max_energy_map) {
1396 std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
1397 sort_by_energy;
1398
1399 for (auto &m : edge_face_max_energy_map) {
1400 auto e = m.first;
1401 auto &[max_face, energy, opt_angle] = m.second;
1402 sort_by_energy[energy] = std::make_tuple(e, max_face, opt_angle);
1403 }
1404
1405 return sort_by_energy;
1406 };
1407
1408 auto set_tag = [&](auto &&adj_edges_map, auto &&sort_by_energy) {
1409 MoFEMFunctionBegin;
1410
1411 Tag th_face_pressure;
1412 CHK_MOAB_THROW(
1413 mField.get_moab().tag_get_handle("FacePressure", th_face_pressure),
1414 "get tag");
1415 auto get_face_pressure = [&](auto face) {
1416 double pressure;
1417 CHK_MOAB_THROW(mField.get_moab().tag_get_data(th_face_pressure, &face,
1418 1, &pressure),
1419 "get rag data");
1420 return pressure;
1421 };
1422
1423 MOFEM_LOG("EPSELF", Sev::inform)
1424 << "Number of edges to check " << sort_by_energy.size();
1425
1426 enum face_energy { POSITIVE, NEGATIVE };
1427 constexpr bool skip_negative = true;
1428
1429 for (auto fe : {face_energy::POSITIVE, face_energy::NEGATIVE}) {
1430
1431 // iterate edges wih maximal energy, and make them seed. Such edges,
1432 // will most likely will have also smallest node displacement
1433 for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
1434 ++it) {
1435
1436 auto energy = it->first;
1437 auto [max_edge, max_face, opt_angle] = it->second;
1438
1439 auto face_pressure = get_face_pressure(max_face);
1440 if (skip_negative) {
1441 if (fe == face_energy::POSITIVE) {
1442 if (face_pressure < 0) {
1443 MOFEM_LOG("EPSELF", Sev::inform)
1444 << "Skip negative face " << max_face << " with energy "
1445 << energy << " and pressure " << face_pressure;
1446 continue;
1447 }
1448 }
1449 }
1450
1451 MOFEM_LOG("EPSELF", Sev::inform)
1452 << "Check face " << max_face << " edge " << max_edge
1453 << " energy " << energy << " optimal angle " << opt_angle
1454 << " face pressure " << face_pressure;
1455
1456 auto jt = adj_edges_map.find(max_edge);
1457 if (jt == adj_edges_map.end()) {
1458 MOFEM_LOG("EPSELF", Sev::warning)
1459 << "Edge " << max_edge << " not found in adj_edges_map";
1460 continue;
1461 }
1462 auto &bundle = jt->second;
1463
1464 auto find_max_in_bundle_impl = [&](auto edge, auto &bundle,
1465 auto gamma) {
1466 FTENSOR_INDEX(SPACE_DIM, i);
1467
1468 EntityHandle vertex_max = 0;
1469 EntityHandle face_max = 0;
1470 EntityHandle move_edge_max = 0;
1471 double max_quality = -2;
1472 double max_quality_evaluated = -2;
1473 double min_cardinality = std::numeric_limits<double>::max();
1474
1475 FTensor::Tensor1<double, SPACE_DIM> t_move_last{0., 0., 0.};
1476
1477 for (auto &b : bundle) {
1478 auto &[face, move_edge, vertex, t_move, quality, cardinality,
1479 edge_gamma] = b;
1480
1481 auto adj_vertex_tets = get_adj(Range(vertex, vertex), 3);
1482 auto adj_vertex_tets_verts = get_conn(adj_vertex_tets);
1483 std::vector<double> coords(3 * adj_vertex_tets_verts.size());
1484 CHK_MOAB_THROW(mField.get_moab().get_coords(
1485 adj_vertex_tets_verts, coords.data()),
1486 "get coords");
1487
1488 set_coord(vertex, adj_vertex_tets_verts, coords, t_move, gamma);
1489 quality = tets_quality(quality, adj_vertex_tets_verts,
1490 adj_vertex_tets, coords);
1491
1492 auto eval_quality = [](auto q, auto c, auto edge_gamma) {
1493 if (q < 0) {
1494 return q;
1495 } else {
1496 return ((edge_gamma < 0) ? (q / 2) : q) / pow(c, 2);
1497 }
1498 };
1499
1500 if (eval_quality(quality, cardinality, edge_gamma) >=
1501 max_quality_evaluated) {
1502 max_quality = quality;
1503 min_cardinality = cardinality;
1504 vertex_max = vertex;
1505 face_max = face;
1506 move_edge_max = move_edge;
1507 t_move_last(i) = t_move(i);
1508 max_quality_evaluated =
1509 eval_quality(max_quality, min_cardinality, edge_gamma);
1510 }
1511 }
1512
1513 return std::make_tuple(vertex_max, face_max, t_move_last,
1514 max_quality, min_cardinality);
1515 };
1516
1517 auto find_max_in_bundle = [&](auto edge, auto &bundle) {
1518 auto b_org_bundle = bundle;
1519 auto r = find_max_in_bundle_impl(edge, bundle, 1.);
1520 auto &[vertex_max, face_max, t_move_last, max_quality,
1521 cardinality] = r;
1522 if (max_quality < 0) {
1523 for (double gamma = 0.95; gamma >= 0.45; gamma -= 0.05) {
1524 bundle = b_org_bundle;
1525 r = find_max_in_bundle_impl(edge, bundle, gamma);
1526 auto &[vertex_max, face_max, t_move_last, max_quality,
1527 cardinality] = r;
1528 MOFEM_LOG("EPSELF", Sev::warning)
1529 << "Back tracking: gamma " << gamma << " edge " << edge
1530 << " quality " << max_quality << " cardinality "
1531 << cardinality;
1532 if (max_quality > 0.01) {
1533 FTENSOR_INDEX(SPACE_DIM, I);
1534 t_move_last(I) *= gamma;
1535 return r;
1536 }
1537 }
1538 FTENSOR_INDEX(SPACE_DIM, I);
1539 t_move_last(I) = 0;
1540 }
1541 return r;
1542 };
1543
1544 // set tags with displacement of node and face energy
1545 auto set_tag_to_vertex_and_face = [&](auto &&r, auto &quality) {
1546 MoFEMFunctionBegin;
1547 auto &[v, f, t_move, q, cardinality] = r;
1548
1549 if ((q > 0 && std::isnormal(q)) && energy > 0) {
1550
1551 MOFEM_LOG("EPSELF", Sev::inform)
1552 << "Set tag: vertex " << v << " face " << f << " "
1553 << max_edge << " move " << t_move << " energy " << energy
1554 << " quality " << q << " cardinality " << cardinality;
1555 CHKERR mField.get_moab().tag_set_data(th_position[0], &v, 1,
1556 &t_move(0));
1557 CHKERR mField.get_moab().tag_set_data(th_max_face_energy[0], &f,
1558 1, &energy);
1559 }
1560
1561 quality = q;
1562 MoFEMFunctionReturn(0);
1563 };
1564
1565 double quality = -2;
1566 CHKERR set_tag_to_vertex_and_face(
1567
1568 find_max_in_bundle(max_edge, bundle),
1569
1570 quality
1571
1572 );
1573
1574 if (quality > 0 && std::isnormal(quality) && energy > 0) {
1575 MOFEM_LOG("EPSELF", Sev::inform)
1576 << "Crack face set with quality: " << quality;
1577 MoFEMFunctionReturnHot(0);
1578 }
1579 }
1580
1581 if (!skip_negative)
1582 break;
1583 }
1584
1585 MoFEMFunctionReturn(0);
1586 };
1587
1588 // map: {edge, {face, energy, optimal_angle}}
1589 MOFEM_LOG("EP", sev) << "Calculate orientation";
1590 std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
1591 edge_face_max_energy_map;
1592 CHKERR find_maximal_face_energy(front_edges, front_faces,
1593 edge_face_max_energy_map);
1594 CHKERR calculate_face_orientation(edge_face_max_energy_map);
1595
1596 MOFEM_LOG("EP", sev) << "Calculate node positions";
1597 CHKERR set_tag(
1598
1599 calculate_free_face_node_displacement(edge_face_max_energy_map),
1600 get_sort_by_energy(edge_face_max_energy_map)
1601
1602 );
1603
1604 MoFEMFunctionReturn(0);
1605 };
1606
1607 MOFEM_LOG("EP", sev) << "Front edges " << frontEdges->size();
1608 CHKERR evaluate_face_energy_and_set_orientation(
1609 *frontEdges, get_adj_front(false), sides_pair, th_front_position);
1610 }
1611
1612 // exchange positions and energies from processor zero to all other
1613 CHKERR VecZeroEntries(vertexExchange.second);
1614 CHKERR VecGhostUpdateBegin(vertexExchange.second, INSERT_VALUES,
1615 SCATTER_FORWARD);
1616 CHKERR VecGhostUpdateEnd(vertexExchange.second, INSERT_VALUES,
1617 SCATTER_FORWARD);
1618 CHKERR mField.getInterface<CommInterface>()->updateEntitiesPetscVector(
1619 mField.get_moab(), vertexExchange, th_front_position[0]);
1620 CHKERR VecZeroEntries(faceExchange.second);
1621 CHKERR VecGhostUpdateBegin(faceExchange.second, INSERT_VALUES,
1622 SCATTER_FORWARD);
1623 CHKERR VecGhostUpdateEnd(faceExchange.second, INSERT_VALUES, SCATTER_FORWARD);
1624 CHKERR mField.getInterface<CommInterface>()->updateEntitiesPetscVector(
1625 mField.get_moab(), faceExchange, th_max_face_energy[0]);
1626
1627 auto get_max_moved_faces = [&]() {
1628 Range max_moved_faces;
1629 auto adj_front = get_adj_front(false);
1630 std::vector<double> face_energy(adj_front.size());
1631 CHKERR mField.get_moab().tag_get_data(th_max_face_energy[0], adj_front,
1632 face_energy.data());
1633 for (int i = 0; i != adj_front.size(); ++i) {
1634 if (face_energy[i] > std::numeric_limits<double>::epsilon()) {
1635 max_moved_faces.insert(adj_front[i]);
1636 }
1637 }
1638
1639 return boost::make_shared<Range>(max_moved_faces);
1640 };
1641
1642 // move all faces with energy larger than 0
1643 maxMovedFaces = get_max_moved_faces();
1644 MOFEM_LOG("EP", sev) << "Number of of moved faces: " << maxMovedFaces->size();
1645
1646#ifndef NDEBUG
1647 if (debug) {
1648 CHKERR save_range(
1649 mField.get_moab(),
1650 "max_moved_faces_" +
1651 boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
1652 *maxMovedFaces);
1653 }
1654#endif
1655
1656 MoFEMFunctionReturn(0);
1657}
1658
1659MoFEMErrorCode EshelbianCore::setNewFrontCoordinates() {
1660 MoFEMFunctionBegin;
1661
1662 if (!maxMovedFaces)
1663 MoFEMFunctionReturnHot(0);
1664
1665 Tag th_front_position;
1666 auto rval =
1667 mField.get_moab().tag_get_handle("FrontPosition", th_front_position);
1668 if (rval == MB_SUCCESS && maxMovedFaces) {
1669 Range verts;
1670 CHKERR mField.get_moab().get_connectivity(*maxMovedFaces, verts, true);
1671 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(verts);
1672 std::vector<double> coords(3 * verts.size());
1673 CHKERR mField.get_moab().get_coords(verts, coords.data());
1674 std::vector<double> pos(3 * verts.size());
1675 CHKERR mField.get_moab().tag_get_data(th_front_position, verts, pos.data());
1676 for (int i = 0; i != 3 * verts.size(); ++i) {
1677 coords[i] += pos[i];
1678 }
1679 CHKERR mField.get_moab().set_coords(verts, coords.data());
1680 double zero[] = {0., 0., 0.};
1681 CHKERR mField.get_moab().tag_clear_data(th_front_position, verts, zero);
1682 }
1683
1684#ifndef NDEBUG
1685 constexpr bool debug = false;
1686 if (debug) {
1687
1688 CHKERR save_range(
1689 mField.get_moab(),
1690 "set_coords_faces_" +
1691 boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
1692 *maxMovedFaces);
1693 }
1694#endif
1695 MoFEMFunctionReturn(0);
1696}
1697
1698MoFEMErrorCode EshelbianCore::addCrackSurfaces(const bool debug) {
1699 MoFEMFunctionBegin;
1700
1701 constexpr bool potential_crack_debug = false;
1702 if constexpr (potential_crack_debug) {
1703
1704 auto add_ents = get_range_from_block(mField, "POTENTIAL", SPACE_DIM - 1);
1705 Range crack_front_verts;
1706 CHKERR mField.get_moab().get_connectivity(*frontEdges, crack_front_verts,
1707 true);
1708 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1709 crack_front_verts);
1710 Range crack_front_faces;
1711 CHKERR mField.get_moab().get_adjacencies(crack_front_verts, SPACE_DIM - 1,
1712 true, crack_front_faces,
1713 moab::Interface::UNION);
1714 crack_front_faces = intersect(crack_front_faces, add_ents);
1715 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1716 crack_front_faces);
1717 CHKERR mField.getInterface<MeshsetsManager>()->addEntitiesToMeshset(
1718 BLOCKSET, addCrackMeshsetId, crack_front_faces);
1719 }
1720
1721 auto get_crack_faces = [&]() {
1722 if (maxMovedFaces) {
1723 return unite(*crackFaces, *maxMovedFaces);
1724 } else {
1725 return *crackFaces;
1726 }
1727 };
1728
1729 auto get_extended_crack_faces = [&]() {
1730 auto get_faces_of_crack_front_verts = [&](auto crack_faces_org) {
1731 ParallelComm *pcomm =
1732 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
1733
1734 Range crack_faces;
1735
1736 if (!pcomm->rank()) {
1737
1738 auto get_nodes = [&](auto &&e) {
1739 Range nodes;
1740 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, nodes, true),
1741 "get connectivity");
1742 return nodes;
1743 };
1744
1745 auto get_adj = [&](auto &&e, auto dim,
1746 auto t = moab::Interface::UNION) {
1747 Range adj;
1748 CHK_MOAB_THROW(
1749 mField.get_moab().get_adjacencies(e, dim, true, adj, t),
1750 "get adj");
1751 return adj;
1752 };
1753
1754 Range body_ents;
1755 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
1756 body_ents);
1757 auto body_skin = get_skin(mField, body_ents);
1758 auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
1759 auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
1760 auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
1761 auto front_block_nodes = get_nodes(front_block_edges);
1762
1763 size_t s;
1764 do {
1765 s = crack_faces.size();
1766
1767 auto crack_face_nodes = get_nodes(crack_faces_org);
1768 auto crack_faces_edges =
1769 get_adj(crack_faces_org, 1, moab::Interface::UNION);
1770
1771 auto crack_skin = get_skin(mField, crack_faces_org);
1772 front_block_edges = subtract(front_block_edges, crack_skin);
1773 auto crack_skin_nodes = get_nodes(crack_skin);
1774 crack_skin_nodes.merge(front_block_nodes);
1775
1776 auto crack_skin_faces =
1777 get_adj(crack_skin, 2, moab::Interface::UNION);
1778 crack_skin_faces =
1779 subtract(subtract(crack_skin_faces, crack_faces_org), body_skin);
1780
1781 crack_faces = crack_faces_org;
1782 for (auto f : crack_skin_faces) {
1783 auto edges = intersect(
1784 get_adj(Range(f, f), 1, moab::Interface::UNION), crack_skin);
1785
1786 // if other edge is part of body skin, e.g. crack punching through
1787 // body surface
1788 if (edges.size() == 2) {
1789 edges.merge(
1790 intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
1791 body_skin_edges));
1792 }
1793
1794 if (edges.size() == 2) {
1795 auto edge_conn = get_nodes(Range(edges));
1796 auto faces = intersect(get_adj(edges, 2, moab::Interface::UNION),
1797 crack_faces_org);
1798 if (faces.size() == 2) {
1799 auto edge0_conn = get_nodes(Range(edges[0], edges[0]));
1800 auto edge1_conn = get_nodes(Range(edges[1], edges[1]));
1801 auto edges_conn = intersect(intersect(edge0_conn, edge1_conn),
1802 crack_skin_nodes); // node at apex
1803 if (edges_conn.size() == 1) {
1804
1805 auto node_edges =
1806 subtract(intersect(get_adj(edges_conn, 1,
1807 moab::Interface::INTERSECT),
1808 crack_faces_edges),
1809 crack_skin); // nodes on crack surface, but not
1810 // at the skin
1811
1812 if (node_edges.size()) {
1813 FTENSOR_INDEX(SPACE_DIM, i);
1814 FTensor::Tensor1<double, SPACE_DIM> t_v0;
1815 CHKERR mField.get_moab().get_coords(edges_conn, &t_v0(0));
1816
1817 auto get_t_dir = [&](auto e_conn) {
1818 auto other_node = subtract(e_conn, edges_conn);
1819 FTensor::Tensor1<double, SPACE_DIM> t_dir;
1820 CHKERR mField.get_moab().get_coords(other_node,
1821 &t_dir(0));
1822 t_dir(i) -= t_v0(i);
1823 return t_dir;
1824 };
1825
1826 FTensor::Tensor1<double, SPACE_DIM> t_ave_dir;
1827 t_ave_dir(i) =
1828 get_t_dir(edge0_conn)(i) + get_t_dir(edge1_conn)(i);
1829
1830 FTensor::Tensor1<double, SPACE_DIM> t_crack_surface_ave_dir;
1831 t_crack_surface_ave_dir(i) = 0;
1832 for (auto e : node_edges) {
1833 auto e_conn = get_nodes(Range(e, e));
1834 auto t_dir = get_t_dir(e_conn);
1835 t_crack_surface_ave_dir(i) += t_dir(i);
1836 }
1837
1838 auto dot = t_ave_dir(i) * t_crack_surface_ave_dir(i);
1839 // ave edges is in opposite direction to crack surface, so
1840 // thus crack is not turning back
1841 if (dot < -std::numeric_limits<double>::epsilon()) {
1842 crack_faces.insert(f);
1843 }
1844 } else {
1845 crack_faces.insert(f);
1846 }
1847 }
1848 }
1849 } else if (edges.size() == 3) {
1850 crack_faces.insert(f);
1851 }
1852
1853 // if other edge is part of geometry edge, e.g. keyway
1854 if (edges.size() == 1) {
1855 edges.merge(
1856 intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
1857 geometry_edges));
1858 edges.merge(
1859 intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
1860 front_block_edges));
1861 if (edges.size() == 2) {
1862 crack_faces.insert(f);
1863 continue;
1864 }
1865 }
1866 }
1867
1868 crack_faces_org = crack_faces;
1869
1870 } while (s != crack_faces.size());
1871 };
1872
1873 return crack_faces; // send_type(mField, crack_faces, MBTRI);
1874 };
1875
1876 return get_faces_of_crack_front_verts(get_crack_faces());
1877 };
1878
1879 if (debug) {
1880 CHKERR save_range(mField.get_moab(), "new_crack_surface_debug.vtk",
1881 get_extended_crack_faces());
1882 }
1883
1884 auto reconstruct_crack_faces = [&](auto crack_faces) {
1885 ParallelComm *pcomm =
1886 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
1887
1888 auto impl = [&]() {
1889 MoFEMFunctionBegin;
1890
1891 Range new_crack_faces;
1892 if (!pcomm->rank()) {
1893
1894 auto get_nodes = [&](auto &&e) {
1895 Range nodes;
1896 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, nodes, true),
1897 "get connectivity");
1898 return nodes;
1899 };
1900
1901 auto get_adj = [&](auto &&e, auto dim,
1902 auto t = moab::Interface::UNION) {
1903 Range adj;
1904 CHK_MOAB_THROW(
1905 mField.get_moab().get_adjacencies(e, dim, true, adj, t),
1906 "get adj");
1907 return adj;
1908 };
1909
1910 auto get_test_on_crack_surface = [&]() {
1911 auto crack_faces_nodes =
1912 get_nodes(crack_faces); // nodes on crac faces
1913 auto crack_faces_tets =
1914 get_adj(crack_faces_nodes, 3,
1915 moab::Interface::UNION); // adjacent
1916 // tets to
1917 // crack
1918 // faces throug nodes
1919 auto crack_faces_tets_nodes =
1920 get_nodes(crack_faces_tets); // nodes on crack faces tets
1921 crack_faces_tets_nodes =
1922 subtract(crack_faces_tets_nodes, crack_faces_nodes);
1923 crack_faces_tets =
1924 subtract(crack_faces_tets, get_adj(crack_faces_tets_nodes, 3,
1925 moab::Interface::UNION));
1926 new_crack_faces =
1927 get_adj(crack_faces_tets, 2,
1928 moab::Interface::UNION); // adjacency faces to crack
1929 // faces through tets
1930 new_crack_faces.merge(crack_faces); // add original crack faces
1931
1932 return std::make_tuple(new_crack_faces, crack_faces_tets);
1933 };
1934
1935 auto carck_faces_test_edges = [&](auto faces, auto tets) {
1936 auto adj_tets_faces = get_adj(tets, 2, moab::Interface::UNION);
1937 auto adj_faces_edges = get_adj(subtract(faces, adj_tets_faces), 1,
1938 moab::Interface::UNION);
1939 auto adj_tets_edges = get_adj(tets, 1, moab::Interface::UNION);
1940 auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
1941 auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
1942 adj_faces_edges.merge(geometry_edges); // geometry edges
1943 adj_faces_edges.merge(front_block_edges); // front block edges
1944
1945 auto boundary_tets_edges = intersect(adj_tets_edges, adj_faces_edges);
1946 auto boundary_test_nodes = get_nodes(boundary_tets_edges);
1947 auto boundary_test_nodes_edges =
1948 get_adj(boundary_test_nodes, 1, moab::Interface::UNION);
1949 auto boundary_test_nodes_edges_nodes = subtract(
1950 get_nodes(boundary_test_nodes_edges), boundary_test_nodes);
1951
1952 boundary_tets_edges =
1953 subtract(boundary_test_nodes_edges,
1954 get_adj(boundary_test_nodes_edges_nodes, 1,
1955 moab::Interface::UNION));
1956
1957 Range body_ents;
1958 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
1959 body_ents);
1960 auto body_skin = get_skin(mField, body_ents);
1961
1962 auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
1963 body_skin_edges = intersect(get_adj(tets, 1, moab::Interface::UNION),
1964 body_skin_edges);
1965 body_skin = intersect(body_skin, adj_tets_faces);
1966 body_skin_edges = subtract(
1967 body_skin_edges, get_adj(body_skin, 1, moab::Interface::UNION));
1968
1969 save_range(mField.get_moab(), "body_skin_edges.vtk", body_skin_edges);
1970 for (auto e : body_skin_edges) {
1971 auto adj_tet = intersect(
1972 get_adj(Range(e, e), 3, moab::Interface::INTERSECT), tets);
1973 if (adj_tet.size() == 1) {
1974 boundary_tets_edges.insert(e);
1975 }
1976 }
1977
1978 return boundary_tets_edges;
1979 };
1980
1981 auto p = get_test_on_crack_surface();
1982 auto &[new_crack_faces, crack_faces_tets] = p;
1983
1984 if (debug) {
1985 CHKERR save_range(mField.get_moab(), "hole_crack_faces_debug.vtk",
1986 crack_faces);
1987 CHKERR save_range(mField.get_moab(), "new_crack_faces_debug.vtk",
1988 new_crack_faces);
1989 CHKERR save_range(mField.get_moab(), "new_crack_tets_debug.vtk",
1990 crack_faces_tets);
1991 }
1992
1993 auto boundary_tets_edges =
1994 carck_faces_test_edges(new_crack_faces, crack_faces_tets);
1995 CHKERR save_range(mField.get_moab(), "boundary_tets_edges.vtk",
1996 boundary_tets_edges);
1997
1998 auto resolve_surface = [&](auto boundary_tets_edges,
1999 auto crack_faces_tets) {
2000 auto boundary_tets_edges_nodes = get_nodes(boundary_tets_edges);
2001 auto crack_faces_tets_faces =
2002 get_adj(crack_faces_tets, 2, moab::Interface::UNION);
2003
2004 Range all_removed_faces;
2005 Range all_removed_tets;
2006 int counter = 0;
2007
2008 int size = 0;
2009 while (size != crack_faces_tets.size()) {
2010 auto tets_faces =
2011 get_adj(crack_faces_tets, 2, moab::Interface::UNION);
2012 auto skin_tets = get_skin(mField, crack_faces_tets);
2013 auto skin_skin =
2014 get_skin(mField, subtract(crack_faces_tets_faces, tets_faces));
2015 auto skin_skin_nodes = get_nodes(skin_skin);
2016
2017 size = crack_faces_tets.size();
2018 MOFEM_LOG("SELF", Sev::inform)
2019 << "Crack faces tets size " << crack_faces_tets.size()
2020 << " crack faces size " << crack_faces_tets_faces.size();
2021 auto skin_tets_nodes = subtract(
2022 get_nodes(skin_tets),
2023 boundary_tets_edges_nodes); // not remove tets which are
2024 // adjagasent to crack faces nodes
2025 skin_tets_nodes = subtract(skin_tets_nodes, skin_skin_nodes);
2026
2027 Range removed_nodes;
2028 Range tets_to_remove;
2029 Range faces_to_remove;
2030 for (auto n : skin_tets_nodes) {
2031 auto tets =
2032 intersect(get_adj(Range(n, n), 3, moab::Interface::INTERSECT),
2033 crack_faces_tets);
2034 if (tets.size() == 0) {
2035 continue;
2036 }
2037
2038 auto hole_detetction = [&]() {
2039 auto adj_tets =
2040 get_adj(Range(n, n), 3, moab::Interface::INTERSECT);
2041 adj_tets =
2042 subtract(adj_tets,
2043 crack_faces_tets); // tetst adjacent to the node
2044 // but not part of crack surface
2045 if (adj_tets.size() == 0) {
2046 return std::make_pair(
2047 intersect(
2048 get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
2049 tets_faces),
2050 tets);
2051 }
2052
2053 std::vector<Range> tets_groups;
2054 auto test_adj_tets = adj_tets;
2055 while (test_adj_tets.size()) {
2056 auto seed_size = 0;
2057 Range seed = Range(test_adj_tets[0], test_adj_tets[0]);
2058 while (seed.size() != seed_size) {
2059 auto adj_faces =
2060 subtract(get_adj(seed, 2, moab::Interface::UNION),
2061 tets_faces); // edges which are not
2062 // part of the node
2063 seed_size = seed.size();
2064 seed.merge(
2065 intersect(get_adj(adj_faces, 3, moab::Interface::UNION),
2066 test_adj_tets));
2067 }
2068 tets_groups.push_back(seed);
2069 test_adj_tets = subtract(test_adj_tets, seed);
2070 }
2071 if (tets_groups.size() == 1) {
2072
2073 return std::make_pair(
2074 intersect(
2075 get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
2076 tets_faces),
2077 tets);
2078 }
2079
2080 Range tets_to_remove;
2081 Range faces_to_remove;
2082 for (auto &r : tets_groups) {
2083 auto f = get_adj(r, 2, moab::Interface::UNION);
2084 auto t = intersect(get_adj(f, 3, moab::Interface::UNION),
2085 crack_faces_tets); // tets
2086
2087 if (f.size() > faces_to_remove.size() ||
2088 faces_to_remove.size() == 0) {
2089 faces_to_remove = f;
2090 tets_to_remove = t; // largest group of tets
2091 }
2092 }
2093 MOFEM_LOG("EPSELF", Sev::inform)
2094 << "Hole detection: faces to remove "
2095 << faces_to_remove.size() << " tets to remove "
2096 << tets_to_remove.size();
2097 return std::make_pair(faces_to_remove, tets_to_remove);
2098 };
2099
2100 if (tets.size() < tets_to_remove.size() ||
2101 tets_to_remove.size() == 0) {
2102 removed_nodes = Range(n, n);
2103 auto [h_faces_to_remove, h_tets_to_remove] =
2104 hole_detetction(); // find faces and tets to remove
2105 faces_to_remove = h_faces_to_remove;
2106 tets_to_remove = h_tets_to_remove;
2107
2108 // intersect(
2109 // get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
2110 // tets_faces);
2111
2112 } // find tets which is largest adjacencty size, so that it is
2113 // removed first, and then faces are removed
2114 all_removed_faces.merge(faces_to_remove);
2115 all_removed_tets.merge(tets_to_remove);
2116 }
2117
2118 crack_faces_tets = subtract(crack_faces_tets, tets_to_remove);
2119 crack_faces_tets_faces =
2120 subtract(crack_faces_tets_faces, faces_to_remove);
2121
2122 if (debug) {
2123 save_range(mField.get_moab(),
2124 "removed_nodes_" +
2125 boost::lexical_cast<std::string>(counter) + ".vtk",
2126 removed_nodes);
2127 save_range(mField.get_moab(),
2128 "faces_to_remove_" +
2129 boost::lexical_cast<std::string>(counter) + ".vtk",
2130 faces_to_remove);
2131 save_range(mField.get_moab(),
2132 "tets_to_remove_" +
2133 boost::lexical_cast<std::string>(counter) + ".vtk",
2134 tets_to_remove);
2135 save_range(mField.get_moab(),
2136 "crack_faces_tets_faces_" +
2137 boost::lexical_cast<std::string>(counter) + ".vtk",
2138 crack_faces_tets_faces);
2139 save_range(mField.get_moab(),
2140 "crack_faces_tets_" +
2141 boost::lexical_cast<std::string>(counter) + ".vtk",
2142 crack_faces_tets);
2143 }
2144 counter++;
2145 }
2146
2147 auto cese_internal_faces = [&]() {
2148 MoFEMFunctionBegin;
2149 auto skin_tets = get_skin(mField, crack_faces_tets);
2150 auto adj_faces = get_adj(skin_tets, 2, moab::Interface::UNION);
2151 adj_faces =
2152 subtract(adj_faces, skin_tets); // remove skin tets faces
2153 auto adj_tets = get_adj(adj_faces, 3,
2154 moab::Interface::UNION); // tets which are
2155 // adjacent to skin
2156 crack_faces_tets =
2157 subtract(crack_faces_tets,
2158 adj_tets); // remove tets which are adjacent to
2159 // skin, so that they are not removed
2160 crack_faces_tets_faces =
2161 subtract(crack_faces_tets_faces, adj_faces);
2162
2163 all_removed_faces.merge(adj_faces);
2164 all_removed_tets.merge(adj_tets);
2165
2166 MOFEM_LOG("EPSELF", Sev::inform)
2167 << "Remove internal faces size " << adj_faces.size()
2168 << " tets size " << adj_tets.size();
2169 MoFEMFunctionReturn(0);
2170 };
2171
2172 auto case_only_one_free_edge = [&]() {
2173 MoFEMFunctionBegin;
2174
2175 for (auto t : Range(crack_faces_tets)) {
2176
2177 auto adj_faces = get_adj(
2178 Range(t, t), 2,
2179 moab::Interface::UNION); // faces of tet which can be removed
2180 auto crack_surface_edges =
2181 get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
2182 adj_faces),
2183 1,
2184 moab::Interface::UNION); // edges not on the tet but
2185 // on crack surface
2186 auto adj_edges =
2187 subtract(get_adj(Range(t, t), 1, moab::Interface::INTERSECT),
2188 crack_surface_edges); // free edges
2189 adj_edges = subtract(
2190 adj_edges,
2191 boundary_tets_edges); // edges which are not part of gemetry
2192
2193 if (adj_edges.size() == 1) {
2194 crack_faces_tets =
2195 subtract(crack_faces_tets,
2196 Range(t, t)); // remove tets which are adjacent to
2197 // skin, so that they are not removed
2198
2199 auto faces_to_remove =
2200 get_adj(adj_edges, 2, moab::Interface::UNION); // faces
2201 // which can
2202 // be removed
2203 crack_faces_tets_faces =
2204 subtract(crack_faces_tets_faces, faces_to_remove);
2205
2206 all_removed_faces.merge(faces_to_remove);
2207 all_removed_tets.merge(Range(t, t));
2208
2209 MOFEM_LOG("EPSELF", Sev::inform) << "Remove free one edges ";
2210 }
2211 }
2212
2213 crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
2214 crack_faces_tets_faces =
2215 subtract(crack_faces_tets_faces, all_removed_faces);
2216
2217 MoFEMFunctionReturn(0);
2218 };
2219
2220 auto cese_flat_tet = [&](auto max_adj_edges) {
2221 MoFEMFunctionBegin;
2222
2223 Range body_ents;
2224 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
2225 body_ents);
2226 auto body_skin = get_skin(mField, body_ents);
2227 auto body_skin_edges =
2228 get_adj(body_skin, 1, moab::Interface::UNION);
2229
2230 for (auto t : Range(crack_faces_tets)) {
2231
2232 auto adj_faces = get_adj(
2233 Range(t, t), 2,
2234 moab::Interface::UNION); // faces of tet which can be removed
2235 auto crack_surface_edges =
2236 get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
2237 adj_faces),
2238 1,
2239 moab::Interface::UNION); // edges not on the tet but
2240 // on crack surface
2241 auto adj_edges =
2242 subtract(get_adj(Range(t, t), 1, moab::Interface::INTERSECT),
2243 crack_surface_edges); // free edges
2244 adj_edges = subtract(adj_edges, body_skin_edges);
2245
2246 auto tet_edges = get_adj(Range(t, t), 1,
2247 moab::Interface::UNION); // edges of
2248 // tet
2249 tet_edges = subtract(tet_edges, adj_edges);
2250
2251 for (auto e : tet_edges) {
2252 constexpr int opposite_edge[] = {5, 3, 4, 1, 2, 0};
2253 auto get_side = [&](auto e) {
2254 int side, sense, offset;
2255 CHK_MOAB_THROW(
2256 mField.get_moab().side_number(t, e, side, sense, offset),
2257 "get side number failed");
2258 return side;
2259 };
2260 auto get_side_ent = [&](auto side) {
2261 EntityHandle side_edge;
2262 CHK_MOAB_THROW(
2263 mField.get_moab().side_element(t, 1, side, side_edge),
2264 "get side failed");
2265 return side_edge;
2266 };
2267 adj_edges.erase(get_side_ent(opposite_edge[get_side(e)]));
2268 }
2269
2270 if (adj_edges.size() <= max_adj_edges) {
2271
2272 double dot = 1;
2273 Range faces_to_remove;
2274 for (auto e : adj_edges) {
2275 auto edge_adj_faces =
2276 get_adj(Range(e, e), 2, moab::Interface::UNION);
2277 edge_adj_faces = intersect(edge_adj_faces, adj_faces);
2278 if (edge_adj_faces.size() != 2) {
2279 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
2280 "Adj faces size is not 2 for edge " +
2281 boost::lexical_cast<std::string>(e));
2282 }
2283
2284 auto get_normal = [&](auto f) {
2285 FTensor::Tensor1<double, SPACE_DIM> t_n;
2286 CHK_THROW_MESSAGE(
2287 mField.getInterface<Tools>()->getTriNormal(f, &t_n(0)),
2288 "get tri normal failed");
2289 return t_n;
2290 };
2291 auto t_n0 = get_normal(edge_adj_faces[0]);
2292 auto t_n1 = get_normal(edge_adj_faces[1]);
2293 auto get_sense = [&](auto f) {
2294 int side, sense, offset;
2295 CHK_MOAB_THROW(mField.get_moab().side_number(t, f, side,
2296 sense, offset),
2297 "get side number failed");
2298 return sense;
2299 };
2300 auto sense0 = get_sense(edge_adj_faces[0]);
2301 auto sense1 = get_sense(edge_adj_faces[1]);
2302 t_n0.normalize();
2303 t_n1.normalize();
2304
2305 FTENSOR_INDEX(SPACE_DIM, i);
2306 auto dot_e = (sense0 * sense1) * t_n0(i) * t_n1(i);
2307 if (dot_e < dot || e == adj_edges[0]) {
2308 dot = dot_e;
2309 faces_to_remove = edge_adj_faces;
2310 }
2311 }
2312
2313 all_removed_faces.merge(faces_to_remove);
2314 all_removed_tets.merge(Range(t, t));
2315
2316 MOFEM_LOG("EPSELF", Sev::inform)
2317 << "Remove free edges on flat tet, with considered nb. of "
2318 "edges "
2319 << adj_edges.size();
2320 }
2321 }
2322
2323 crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
2324 crack_faces_tets_faces =
2325 subtract(crack_faces_tets_faces, all_removed_faces);
2326
2327 MoFEMFunctionReturn(0);
2328 };
2329
2330 CHK_THROW_MESSAGE(case_only_one_free_edge(),
2331 "Case only one free edge failed");
2332 for (auto max_adj_edges : {0, 1, 2, 3}) {
2333 CHK_THROW_MESSAGE(cese_flat_tet(max_adj_edges),
2334 "Case only one free edge failed");
2335 }
2336 CHK_THROW_MESSAGE(cese_internal_faces(),
2337 "Case internal faces failed");
2338
2339 if (debug) {
2340 save_range(mField.get_moab(),
2341 "crack_faces_tets_faces_" +
2342 boost::lexical_cast<std::string>(counter) + ".vtk",
2343 crack_faces_tets_faces);
2344 save_range(mField.get_moab(),
2345 "crack_faces_tets_" +
2346 boost::lexical_cast<std::string>(counter) + ".vtk",
2347 crack_faces_tets);
2348 }
2349
2350 return std::make_tuple(crack_faces_tets_faces, crack_faces_tets,
2351 all_removed_faces, all_removed_tets);
2352 };
2353
2354 auto [resolved_faces, resolved_tets, all_removed_faces,
2355 all_removed_tets] =
2356 resolve_surface(boundary_tets_edges, crack_faces_tets);
2357 resolved_faces.merge(subtract(crack_faces, all_removed_faces));
2358 if (debug) {
2359 CHKERR save_range(mField.get_moab(), "resolved_faces.vtk",
2360 resolved_faces);
2361 CHKERR save_range(mField.get_moab(), "resolved_tets.vtk",
2362 resolved_tets);
2363 }
2364
2365 crack_faces = resolved_faces;
2366 }
2367
2368 MoFEMFunctionReturn(0);
2369 };
2370
2371 CHK_THROW_MESSAGE(impl(), "resolve new crack surfaces");
2372
2373 return crack_faces; // send_type(mField, crack_faces, MBTRI);
2374 };
2375
2376 auto resolve_consisten_crack_extension = [&]() {
2377 MoFEMFunctionBegin;
2378 auto crack_meshset =
2379 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
2380 addCrackMeshsetId, BLOCKSET);
2381 auto meshset = crack_meshset->getMeshset();
2382
2383 if (!mField.get_comm_rank()) {
2384 Range old_crack_faces;
2385 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTRI,
2386 old_crack_faces);
2387 auto extendeded_crack_faces = get_extended_crack_faces();
2388 auto reconstructed_crack_faces =
2389 subtract(reconstruct_crack_faces(extendeded_crack_faces),
2390 subtract(*crackFaces, old_crack_faces));
2391 if (nbCrackFaces >= reconstructed_crack_faces.size()) {
2392 MOFEM_LOG("EPSELF", Sev::warning)
2393 << "No new crack faces to add, skipping adding to meshset";
2394 extendeded_crack_faces = subtract(
2395 extendeded_crack_faces, subtract(*crackFaces, old_crack_faces));
2396 MOFEM_LOG("EPSELF", Sev::inform)
2397 << "Number crack faces size (extended) "
2398 << extendeded_crack_faces.size();
2399 CHKERR mField.get_moab().clear_meshset(&meshset, 1);
2400 CHKERR mField.get_moab().add_entities(meshset, extendeded_crack_faces);
2401 } else {
2402 CHKERR mField.get_moab().clear_meshset(&meshset, 1);
2403 CHKERR mField.get_moab().add_entities(meshset,
2404 reconstructed_crack_faces);
2405 MOFEM_LOG("EPSELF", Sev::inform)
2406 << "Number crack faces size (reconstructed) "
2407 << reconstructed_crack_faces.size();
2408 nbCrackFaces = reconstructed_crack_faces.size();
2409 }
2410 }
2411
2412 Range crack_faces;
2413 if (!mField.get_comm_rank()) {
2414 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTRI,
2415 crack_faces);
2416 }
2417 crack_faces = send_type(mField, crack_faces, MBTRI);
2418 if (mField.get_comm_rank()) {
2419 CHKERR mField.get_moab().clear_meshset(&meshset, 1);
2420 CHKERR mField.get_moab().add_entities(meshset, crack_faces);
2421 }
2422
2423 MoFEMFunctionReturn(0);
2424 };
2425
2426 CHKERR resolve_consisten_crack_extension();
2427
2428 MoFEMFunctionReturn(0);
2429};
2430
2431MoFEMErrorCode EshelbianCore::saveOrgCoords() {
2432 MoFEMFunctionBegin;
2433 auto crack_faces =
2434 get_range_from_block(mField, "CRACK_COMPUTED", SPACE_DIM - 1);
2435 Range conn;
2436 CHKERR mField.get_moab().get_connectivity(crack_faces, conn, true);
2437 Range verts;
2438 CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, verts);
2439 verts = subtract(verts, conn);
2440 std::vector<double> coords(3 * verts.size());
2441 CHKERR mField.get_moab().get_coords(verts, coords.data());
2442 double def_coords[] = {0., 0., 0.};
2443 Tag th_org_coords;
2444 CHKERR mField.get_moab().tag_get_handle(
2445 "ORG_COORDS", 3, MB_TYPE_DOUBLE, th_org_coords,
2446 MB_TAG_CREAT | MB_TAG_DENSE, def_coords);
2447 CHKERR mField.get_moab().tag_set_data(th_org_coords, verts, coords.data());
2448 MoFEMFunctionReturn(0);
2449}
2450
2451MoFEMErrorCode EshelbianCore::createCrackSurfaceMeshset() {
2452 MoFEMFunctionBegin;
2453 auto meshset_mng = mField.getInterface<MeshsetsManager>();
2454 while (meshset_mng->checkMeshset(addCrackMeshsetId, BLOCKSET))
2455 ++addCrackMeshsetId;
2456 MOFEM_LOG("EP", Sev::inform)
2457 << "Crack added surface meshset " << addCrackMeshsetId;
2458 CHKERR meshset_mng->addMeshset(BLOCKSET, addCrackMeshsetId, "CRACK_COMPUTED");
2459 MoFEMFunctionReturn(0);
2460};
2461
2462MoFEMErrorCode
2463EshelbianCore::calculateCrackArea(boost::shared_ptr<double> area_ptr) {
2464 MoFEMFunctionBegin;
2465
2466 if (!area_ptr) {
2467 CHK_THROW_MESSAGE(MOFEM_INVALID_DATA, "area_ptr is null");
2468 }
2469
2470 int success;
2471 *area_ptr = 0;
2472 if (mField.get_comm_rank() == 0) {
2473 MOFEM_LOG("EP", Sev::inform) << "Calculate crack area";
2474 auto crack_faces = get_range_from_block(mField, "CRACK", SPACE_DIM - 1);
2475 for (auto f : crack_faces) {
2476 *area_ptr += mField.getInterface<Tools>()->getTriArea(f);
2477 }
2478 success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0, mField.get_comm());
2479 } else {
2480 success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0, mField.get_comm());
2481 }
2482 if (success != MPI_SUCCESS) {
2483 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "bcast failed");
2484 }
2485 MoFEMFunctionReturn(0);
2486}
2487
2488} // namespace EshelbianPlasticity
send_type
static auto send_type(MoFEM::Interface &m_field, Range r, const EntityType type)
Definition EshelbianPlasticity.cpp:64
get_two_sides_of_crack_surface
static auto get_two_sides_of_crack_surface(MoFEM::Interface &m_field, Range crack_faces)
Definition EshelbianPlasticity.cpp:244
MOFEM_LOG_SEVERITY_SYNC
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
Definition LogManager.hpp:360
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition Templates.hpp:2110
a
constexpr double a
Definition approx_sphere.cpp:30
eps
static const double eps
Definition check_base_functions_derivatives_on_tet.cpp:11
SPACE_DIM
constexpr int SPACE_DIM
Definition child_and_parent.cpp:17
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition Kronecker_Delta.hpp:15
FTensor::Tensor1::normalize
Tensor1< T, Tensor_Dim > normalize()
Definition Tensor1_value.hpp:77
NOT_USED
#define NOT_USED(x)
Definition definitions.h:255
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition definitions.h:612
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition definitions.h:592
BLOCKSET
@ BLOCKSET
Definition definitions.h:161
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition definitions.h:36
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
debug
static const bool debug
Definition dm_create_subdm.cpp:12
t_kd
constexpr auto t_kd
Definition free_surface.cpp:163
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition LogManager.hpp:292
a0
constexpr double a0
Definition hcurl_check_approx_in_2d.cpp:37
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
c
const double c
speed of light (cm/ns)
Definition initial_diffusion.cpp:39
v
const double v
phase velocity of light in medium (cm/ns)
Definition initial_diffusion.cpp:40
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38
J
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
EshelbianPlasticity::get_range_from_block
static auto get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
Definition EshelbianContact.cpp:1889
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition Levi_Civita.hpp:617
I
constexpr IntegrationType I
Definition operators_tests.cpp:31
t
constexpr double t
plate stiffness
Definition plate.cpp:58
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
EshelbianCore::createCrackSurfaceMeshset
MoFEMErrorCode createCrackSurfaceMeshset()
Definition EshelbianFracture.cpp:2451
EshelbianCore::calculateFaceMaterialForce
MoFEMErrorCode calculateFaceMaterialForce(const int tag, TS ts)
Definition EshelbianFracture.cpp:14
EshelbianCore::calculateCrackArea
MoFEMErrorCode calculateCrackArea(boost::shared_ptr< double > area_ptr)
Definition EshelbianFracture.cpp:2463
EshelbianCore::calculateOrientation
MoFEMErrorCode calculateOrientation(const int tag, bool set_orientation)
Definition EshelbianFracture.cpp:630
EshelbianCore::setNewFrontCoordinates
MoFEMErrorCode setNewFrontCoordinates()
Definition EshelbianFracture.cpp:1659
EshelbianCore::addCrackSurfaces
MoFEMErrorCode addCrackSurfaces(const bool debug=false)
Definition EshelbianFracture.cpp:1698
EshelbianCore::saveOrgCoords
MoFEMErrorCode saveOrgCoords()
Definition EshelbianFracture.cpp:2431
LieGroups::SO3::exp
static auto exp(A &&t_w_vee, B &&theta)
Definition Lie.hpp:48
save_range
auto save_range
Definition thermo_elastic.cpp:121
Generated by Doxygen 1.9.8 and hosted at University of Glasgow