v0.15.5
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 (void)debug;
636 constexpr auto sev = Sev::verbose;
637
638 Range body_ents;
639 CHKERR mField.get_moab().get_entities_by_dimension(0, 3, body_ents);
640 auto body_skin = get_skin(mField, body_ents);
641 Range body_skin_edges;
642 CHKERR mField.get_moab().get_adjacencies(body_skin, 1, false, body_skin_edges,
643 moab::Interface::UNION);
644 Range boundary_skin_verts;
645 CHKERR mField.get_moab().get_connectivity(body_skin_edges,
646 boundary_skin_verts, true);
647
648 auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
649 Range geometry_edges_verts;
650 CHKERR mField.get_moab().get_connectivity(geometry_edges,
651 geometry_edges_verts, true);
652 Range crack_faces_verts;
653 CHKERR mField.get_moab().get_connectivity(*crackFaces, crack_faces_verts,
654 true);
655 Range crack_faces_edges;
656 CHKERR mField.get_moab().get_adjacencies(
657 *crackFaces, 1, true, crack_faces_edges, moab::Interface::UNION);
658 Range crack_faces_tets;
659 CHKERR mField.get_moab().get_adjacencies(
660 *crackFaces, 3, true, crack_faces_tets, moab::Interface::UNION);
661
662 Range front_verts;
663 CHKERR mField.get_moab().get_connectivity(*frontEdges, front_verts, true);
664 Range front_faces;
665 CHKERR mField.get_moab().get_adjacencies(*frontEdges, 2, true, front_faces,
666 moab::Interface::UNION);
667 Range front_verts_edges;
668 CHKERR mField.get_moab().get_adjacencies(
669 front_verts, 1, true, front_verts_edges, moab::Interface::UNION);
670
671 auto get_tags_vec = [&](auto tag_name, int dim) {
672 std::vector<Tag> tags(1);
673
674 if (dim > 3)
675 CHK_THROW_MESSAGE(MOFEM_ATOM_TEST_INVALID, "dim>3");
676
677 auto create_and_clean = [&]() {
678 MoFEMFunctionBegin;
679 auto &moab = mField.get_moab();
680 auto rval = moab.tag_get_handle(tag_name, tags[0]);
681 if (rval == MB_SUCCESS) {
682 moab.tag_delete(tags[0]);
683 }
684 double def_val[] = {0., 0., 0.};
685 CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
686 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
687 MoFEMFunctionReturn(0);
688 };
689
690 CHK_THROW_MESSAGE(create_and_clean(), "create_and_clean");
691
692 return tags;
693 };
694
695 auto get_adj_front = [&](bool subtract_crack) {
696 Range adj_front;
697 CHKERR mField.get_moab().get_adjacencies(*frontEdges, SPACE_DIM - 1, true,
698 adj_front, moab::Interface::UNION);
699 if (subtract_crack)
700 adj_front = subtract(adj_front, *crackFaces);
701 return adj_front;
702 };
703
704 MOFEM_LOG_CHANNEL("SELF");
705
706 auto th_front_position = get_tags_vec("FrontPosition", 3);
707 auto th_max_face_energy = get_tags_vec("MaxFaceEnergy", 1);
708
709 if (mField.get_comm_rank() == 0) {
710
711 auto get_crack_adj_tets = [&](auto r) {
712 Range crack_faces_conn;
713 CHKERR mField.get_moab().get_connectivity(r, crack_faces_conn);
714 Range crack_faces_conn_tets;
715 CHKERR mField.get_moab().get_adjacencies(crack_faces_conn, SPACE_DIM,
716 true, crack_faces_conn_tets,
717 moab::Interface::UNION);
718 return crack_faces_conn_tets;
719 };
720
721 auto get_layers_for_sides = [&](auto &side) {
722 std::vector<Range> layers;
723 auto get = [&]() {
724 MoFEMFunctionBegin;
725
726 auto get_adj = [&](auto &r, int dim) {
727 Range adj;
728 CHKERR mField.get_moab().get_adjacencies(r, dim, true, adj,
729 moab::Interface::UNION);
730 return adj;
731 };
732
733 auto get_tets = [&](auto r) { return get_adj(r, SPACE_DIM); };
734
735 Range front_nodes;
736 CHKERR mField.get_moab().get_connectivity(*frontEdges, front_nodes,
737 true);
738 Range front_faces = get_adj(front_nodes, 2);
739 front_faces = subtract(front_faces, *crackFaces);
740 auto front_tets = get_tets(front_nodes);
741 auto front_side = intersect(side, front_tets);
742 layers.push_back(front_side);
743 for (;;) {
744 auto adj_faces = get_skin(mField, layers.back());
745 adj_faces = intersect(adj_faces, front_faces);
746 auto adj_faces_tets = get_tets(adj_faces);
747 adj_faces_tets = intersect(adj_faces_tets, front_tets);
748 layers.push_back(unite(layers.back(), adj_faces_tets));
749 if (layers.back().size() == layers[layers.size() - 2].size()) {
750 break;
751 }
752 }
753 MoFEMFunctionReturn(0);
754 };
755 CHK_THROW_MESSAGE(get(), "get_layers_for_sides");
756 return layers;
757 };
758
759 auto sides_pair = get_two_sides_of_crack_surface(mField, *crackFaces);
760 auto layers_top = get_layers_for_sides(sides_pair.first);
761 auto layers_bottom = get_layers_for_sides(sides_pair.second);
762
763#ifndef NDEBUG
764 if (debug) {
765 CHKERR save_range(
766 mField.get_moab(),
767 "crack_tets_" +
768 boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
769 get_crack_adj_tets(*crackFaces));
770 CHKERR save_range(mField.get_moab(), "sides_first.vtk", sides_pair.first);
771 CHKERR save_range(mField.get_moab(), "sides_second.vtk",
772 sides_pair.second);
773 MOFEM_LOG("EP", sev) << "Nb. layers " << layers_top.size();
774 int l = 0;
775 for (auto &r : layers_top) {
776 MOFEM_LOG("EP", sev) << "Layer " << l << " size " << r.size();
777 CHKERR save_range(
778 mField.get_moab(),
779 "layers_top_" + boost::lexical_cast<std::string>(l) + ".vtk", r);
780 ++l;
781 }
782
783 l = 0;
784 for (auto &r : layers_bottom) {
785 MOFEM_LOG("EP", sev) << "Layer " << l << " size " << r.size();
786 CHKERR save_range(
787 mField.get_moab(),
788 "layers_bottom_" + boost::lexical_cast<std::string>(l) + ".vtk", r);
789 ++l;
790 }
791 }
792#endif
793
794 auto get_cross = [&](auto t_dir, auto f) {
795 FTensor::Tensor1<double, SPACE_DIM> t_normal;
796 CHKERR mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
797 t_normal.normalize();
798 FTensor::Tensor1<double, SPACE_DIM> t_cross;
799 FTENSOR_INDEX(SPACE_DIM, i);
800 FTENSOR_INDEX(SPACE_DIM, j);
801 FTENSOR_INDEX(SPACE_DIM, k);
802 t_cross(i) = FTensor::levi_civita(i, j, k) * t_normal(j) * t_dir(k);
803 return t_cross;
804 };
805
806 auto get_sense = [&](auto f, auto e) {
807 int side, sense, offset;
808 CHK_MOAB_THROW(mField.get_moab().side_number(f, e, side, sense, offset),
809 "get sense");
810 return std::make_tuple(side, sense, offset);
811 };
812
813 auto calculate_edge_direction = [&](auto e, auto normalize = true) {
814 const EntityHandle *conn;
815 int num_nodes;
816 CHKERR mField.get_moab().get_connectivity(e, conn, num_nodes, true);
817 std::array<double, 6> coords;
818 CHKERR mField.get_moab().get_coords(conn, num_nodes, coords.data());
819 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p0{
820 &coords[0], &coords[1], &coords[2]};
821 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, SPACE_DIM> t_p1{
822 &coords[3], &coords[4], &coords[5]};
823 FTensor::Tensor1<double, SPACE_DIM> t_dir;
824 FTENSOR_INDEX(SPACE_DIM, i);
825 t_dir(i) = t_p1(i) - t_p0(i);
826 if (normalize)
827 t_dir.normalize();
828 return t_dir;
829 };
830
831 auto evaluate_face_energy_and_set_orientation = [&](auto front_edges,
832 auto front_faces,
833 auto &sides_pair,
834 auto th_position) {
835 MoFEMFunctionBegin;
836
837 Tag th_face_energy;
838 Tag th_material_force;
839 switch (energyReleaseSelector) {
840 case GRIFFITH_FORCE:
841 case GRIFFITH_SKELETON:
842 CHKERR mField.get_moab().tag_get_handle("GriffithForce",
843 th_face_energy);
844 CHKERR mField.get_moab().tag_get_handle("MaterialForce",
845 th_material_force);
846 break;
847 default:
848 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
849 "Unknown energy release selector");
850 };
851
852 /**
853 * Iterate over front edges, get adjacent faces, find maximal face energy.
854 * Maximal face energy is stored in the edge. Maximal face energy is
855 * magnitude of edge Griffith force.
856 */
857 auto find_maximal_face_energy = [&](auto front_edges, auto front_faces,
858 auto &edge_face_max_energy_map) {
859 MoFEMFunctionBegin;
860
861 Range body_ents;
862 CHKERR mField.get_moab().get_entities_by_dimension(0, 3, body_ents);
863 auto body_skin = get_skin(mField, body_ents);
864
865 Range max_faces;
866
867 for (auto e : front_edges) {
868
869 double griffith_force;
870 CHKERR mField.get_moab().tag_get_data(th_face_energy, &e, 1,
871 &griffith_force);
872
873 Range faces;
874 CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
875 faces = subtract(intersect(faces, front_faces), body_skin);
876 std::vector<double> face_energy(faces.size());
877 CHKERR mField.get_moab().tag_get_data(th_face_energy, faces,
878 face_energy.data());
879 auto max_energy_it =
880 std::max_element(face_energy.begin(), face_energy.end());
881 double max_energy =
882 max_energy_it != face_energy.end() ? *max_energy_it : 0;
883
884 edge_face_max_energy_map[e] =
885 std::make_tuple(faces[max_energy_it - face_energy.begin()],
886 griffith_force, static_cast<double>(0));
887 MOFEM_LOG("EP", Sev::inform)
888 << "Edge " << e << " griffith force " << griffith_force
889 << " max face energy " << max_energy << " factor "
890 << max_energy / griffith_force;
891
892 max_faces.insert(faces[max_energy_it - face_energy.begin()]);
893 }
894
895#ifndef NDEBUG
896 if (debug) {
897 CHKERR save_range(
898 mField.get_moab(),
899 "max_faces_" +
900 boost::lexical_cast<std::string>(mField.get_comm_rank()) +
901 ".vtk",
902 max_faces);
903 }
904#endif
905
906 MoFEMFunctionReturn(0);
907 };
908
909 /**
910 * For each front edge, find maximal face energy and orientation. This is
911 * by finding angle between edge material force and maximal face normal
912 *
913 */
914 auto calculate_face_orientation = [&](auto &edge_face_max_energy_map) {
915 MoFEMFunctionBegin;
916
917 auto up_down_face = [&](
918
919 auto &face_angle_map_up,
920 auto &face_angle_map_down
921
922 ) {
923 MoFEMFunctionBegin;
924
925 for (auto &m : edge_face_max_energy_map) {
926 auto e = m.first;
927 auto [max_face, energy, opt_angle] = m.second;
928
929 Range faces;
930 CHKERR mField.get_moab().get_adjacencies(&e, 1, 2, false, faces);
931 faces = intersect(faces, front_faces);
932 Range adj_tets; // tetrahedrons adjacent to the face
933 CHKERR mField.get_moab().get_adjacencies(&max_face, 1, SPACE_DIM,
934 false, adj_tets,
935 moab::Interface::UNION);
936 if (adj_tets.size()) {
937
938 Range adj_tets; // tetrahedrons adjacent to the face
939 CHKERR mField.get_moab().get_adjacencies(&max_face, 1, SPACE_DIM,
940 false, adj_tets,
941 moab::Interface::UNION);
942 if (adj_tets.size()) {
943
944 Range adj_tets_faces;
945 // get faces
946 CHKERR mField.get_moab().get_adjacencies(
947 adj_tets, SPACE_DIM - 1, false, adj_tets_faces,
948 moab::Interface::UNION);
949 adj_tets_faces = intersect(adj_tets_faces, faces);
950 FTENSOR_INDEX(SPACE_DIM, i);
951
952 // cross product of face normal and edge direction
953 auto t_cross_max =
954 get_cross(calculate_edge_direction(e, true), max_face);
955 auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
956 t_cross_max(i) *= sense_max;
957
958 for (auto t : adj_tets) {
959 Range adj_tets_faces;
960 CHKERR mField.get_moab().get_adjacencies(
961 &t, 1, SPACE_DIM - 1, false, adj_tets_faces);
962 adj_tets_faces = intersect(adj_tets_faces, faces);
963 adj_tets_faces =
964 subtract(adj_tets_faces, Range(max_face, max_face));
965
966 if (adj_tets_faces.size() == 1) {
967
968 // cross product of adjacent face normal and edge
969 // direction
970 auto t_cross = get_cross(calculate_edge_direction(e, true),
971 adj_tets_faces[0]);
972 auto [side, sense, offset] =
973 get_sense(adj_tets_faces[0], e);
974 t_cross(i) *= sense;
975 double dot = t_cross(i) * t_cross_max(i);
976 auto angle = std::acos(dot);
977
978 double face_energy;
979 CHKERR mField.get_moab().tag_get_data(
980 th_face_energy, adj_tets_faces, &face_energy);
981
982 auto [side_face, sense_face, offset_face] =
983 get_sense(t, max_face);
984
985 if (sense_face > 0) {
986 face_angle_map_up[e] = std::make_tuple(face_energy, angle,
987 adj_tets_faces[0]);
988
989 } else {
990 face_angle_map_down[e] = std::make_tuple(
991 face_energy, -angle, adj_tets_faces[0]);
992 }
993 }
994 }
995 }
996 }
997 }
998
999 MoFEMFunctionReturn(0);
1000 };
1001
1002 auto calc_optimal_angle = [&](
1003
1004 auto &face_angle_map_up,
1005 auto &face_angle_map_down
1006
1007 ) {
1008 MoFEMFunctionBegin;
1009
1010 for (auto &m : edge_face_max_energy_map) {
1011 auto e = m.first;
1012 auto &[max_face, e0, a0] = m.second;
1013
1014 if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
1015
1016 if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
1017 face_angle_map_down.find(e) == face_angle_map_down.end()) {
1018 // Do nothing
1019 } else {
1020
1021 switch (energyReleaseSelector) {
1022 case GRIFFITH_FORCE:
1023 case GRIFFITH_SKELETON: {
1024
1025 Tag th_material_force;
1026 CHKERR mField.get_moab().tag_get_handle("MaterialForce",
1027 th_material_force);
1028 FTensor::Tensor1<double, SPACE_DIM> t_material_force;
1029 CHKERR mField.get_moab().tag_get_data(
1030 th_material_force, &e, 1, &t_material_force(0));
1031 auto material_force_magnitude = t_material_force.l2();
1032 if (material_force_magnitude <
1033 std::numeric_limits<double>::epsilon()) {
1034 a0 = 0;
1035
1036 } else {
1037
1038 auto t_edge_dir = calculate_edge_direction(e, true);
1039 auto t_cross_max = get_cross(t_edge_dir, max_face);
1040 auto [side, sense, offset] = get_sense(max_face, e);
1041 t_cross_max(sense) *= sense;
1042
1043 FTENSOR_INDEX(SPACE_DIM, I);
1044 FTENSOR_INDEX(SPACE_DIM, J);
1045 FTENSOR_INDEX(SPACE_DIM, K);
1046
1047 t_material_force.normalize();
1048 t_cross_max.normalize();
1049 FTensor::Tensor1<double, SPACE_DIM> t_cross;
1050 t_cross(I) = FTensor::levi_civita(I, J, K) *
1051 t_material_force(J) * t_cross_max(K);
1052 a0 = -std::asin(t_cross(I) * t_edge_dir(I));
1053
1054 MOFEM_LOG("EP", sev)
1055 << "Optimal angle " << a0 << " energy " << e0;
1056 }
1057 break;
1058 }
1059 default: {
1060
1061 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1062 "Unknown energy release selector");
1063 }
1064 }
1065 }
1066 }
1067 }
1068
1069 MoFEMFunctionReturn(0);
1070 };
1071
1072 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
1073 face_angle_map_up;
1074 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
1075 face_angle_map_down;
1076 CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
1077 CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
1078
1079#ifndef NDEBUG
1080 if (debug) {
1081 auto th_angle = get_tags_vec("Angle", 1);
1082 Range up;
1083 for (auto &m : face_angle_map_up) {
1084 auto [e, a, face] = m.second;
1085 up.insert(face);
1086 CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1, &a);
1087 }
1088 Range down;
1089 for (auto &m : face_angle_map_down) {
1090 auto [e, a, face] = m.second;
1091 down.insert(face);
1092 CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1, &a);
1093 }
1094
1095 Range max_energy_faces;
1096 for (auto &m : edge_face_max_energy_map) {
1097 auto [face, e, angle] = m.second;
1098 max_energy_faces.insert(face);
1099 CHKERR mField.get_moab().tag_set_data(th_angle[0], &face, 1,
1100 &angle);
1101 }
1102 if (mField.get_comm_rank() == 0) {
1103 CHKERR save_range(mField.get_moab(), "up_faces.vtk", up);
1104 CHKERR save_range(mField.get_moab(), "down_faces.vtk", down);
1105 CHKERR save_range(mField.get_moab(), "max_energy_faces.vtk",
1106 max_energy_faces);
1107 }
1108 }
1109#endif // NDEBUG
1110
1111 MoFEMFunctionReturn(0);
1112 };
1113
1114 auto get_conn = [&](auto e) {
1115 Range conn;
1116 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, conn, true),
1117 "get conn");
1118 return conn;
1119 };
1120
1121 auto get_adj = [&](auto e, auto dim) {
1122 Range adj;
1123 CHK_MOAB_THROW(mField.get_moab().get_adjacencies(
1124 e, dim, false, adj, moab::Interface::UNION),
1125 "get adj");
1126 return adj;
1127 };
1128
1129 auto get_coords = [&](auto v) {
1130 FTensor::Tensor1<double, 3> t_coords;
1131 CHK_MOAB_THROW(mField.get_moab().get_coords(v, &t_coords(0)),
1132 "get coords");
1133 return t_coords;
1134 };
1135
1136 // calulate normal of the max energy face
1137 auto get_rotated_normal = [&](auto e, auto f, auto angle) {
1138 FTENSOR_INDEX(SPACE_DIM, i);
1139 FTENSOR_INDEX(SPACE_DIM, j);
1140 auto t_edge_dir = calculate_edge_direction(e, true);
1141 auto [side, sense, offset] = get_sense(f, e);
1142 t_edge_dir(i) *= sense;
1143 t_edge_dir.normalize();
1144 t_edge_dir(i) *= angle;
1145 auto t_R = LieGroups::SO3::exp(t_edge_dir, angle);
1146 FTensor::Tensor1<double, SPACE_DIM> t_normal;
1147 mField.getInterface<Tools>()->getTriNormal(f, &t_normal(0));
1148 FTensor::Tensor1<double, SPACE_DIM> t_rotated_normal;
1149 t_rotated_normal(i) = t_R(i, j) * t_normal(j);
1150 return std::make_tuple(t_normal, t_rotated_normal);
1151 };
1152
1153 auto set_coord = [&](auto v, auto &adj_vertex_tets_verts, auto &coords,
1154 auto &t_move, auto gamma) {
1155 auto index = adj_vertex_tets_verts.index(v);
1156 if (index >= 0) {
1157 for (auto ii : {0, 1, 2}) {
1158 coords[3 * index + ii] += gamma * t_move(ii);
1159 }
1160 return true;
1161 }
1162 return false;
1163 };
1164
1165 auto tets_quality = [&](auto quality, auto &adj_vertex_tets_verts,
1166 auto &adj_vertex_tets, auto &coords) {
1167 for (auto t : adj_vertex_tets) {
1168 const EntityHandle *conn;
1169 int num_nodes;
1170 CHKERR mField.get_moab().get_connectivity(t, conn, num_nodes, true);
1171 std::array<double, 12> tet_coords;
1172 for (auto n = 0; n != 4; ++n) {
1173 auto index = adj_vertex_tets_verts.index(conn[n]);
1174 if (index < 0) {
1175 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong index");
1176 }
1177 for (auto ii = 0; ii != 3; ++ii) {
1178 tet_coords[3 * n + ii] = coords[3 * index + ii];
1179 }
1180 }
1181 double q = Tools::volumeLengthQuality(tet_coords.data());
1182 if (!std::isnormal(q))
1183 q = -2;
1184 quality = std::min(quality, q);
1185 };
1186
1187 return quality;
1188 };
1189
1190 auto calculate_free_face_node_displacement =
1191 [&](auto &edge_face_max_energy_map) {
1192 // get edges adjacent to vertex along which nodes are moving
1193 auto get_vertex_edges = [&](auto vertex) {
1194 Range vertex_edges; // edges adjacent to vertex
1195
1196 auto impl = [&]() {
1197 MoFEMFunctionBegin;
1198 CHKERR mField.get_moab().get_adjacencies(vertex, 1, false,
1199 vertex_edges);
1200 vertex_edges = subtract(vertex_edges, front_verts_edges);
1201
1202 if (boundary_skin_verts.size() &&
1203 boundary_skin_verts.find(vertex[0]) !=
1204 boundary_skin_verts.end()) {
1205 MOFEM_LOG("EP", sev) << "Boundary vertex";
1206 vertex_edges = intersect(vertex_edges, body_skin_edges);
1207 }
1208 if (geometry_edges_verts.size() &&
1209 geometry_edges_verts.find(vertex[0]) !=
1210 geometry_edges_verts.end()) {
1211 MOFEM_LOG("EP", sev) << "Geometry edge vertex";
1212 vertex_edges = intersect(vertex_edges, geometry_edges);
1213 }
1214 if (crack_faces_verts.size() &&
1215 crack_faces_verts.find(vertex[0]) !=
1216 crack_faces_verts.end()) {
1217 MOFEM_LOG("EP", sev) << "Crack face vertex";
1218 vertex_edges = intersect(vertex_edges, crack_faces_edges);
1219 }
1220 MoFEMFunctionReturn(0);
1221 };
1222
1223 CHK_THROW_MESSAGE(impl(), "get_vertex_edges");
1224
1225 return vertex_edges;
1226 };
1227
1228 // vector of rotated faces, edge along node is moved, moved edge,
1229 // moved displacement, quality, cardinality, gamma
1230 using Bundle = std::vector<
1231
1232 std::tuple<EntityHandle, EntityHandle, EntityHandle,
1233 FTensor::Tensor1<double, 3>, double, double, double>
1234
1235 >;
1236 std::map<EntityHandle, Bundle> edge_bundle_map;
1237
1238 for (auto &m : edge_face_max_energy_map) {
1239
1240 auto edge = m.first;
1241 auto &[max_face, energy, opt_angle] = m.second;
1242
1243 // calculate rotation of max energy face
1244 auto [t_normal, t_rotated_normal] =
1245 get_rotated_normal(edge, max_face, opt_angle);
1246
1247 auto front_vertex = get_conn(Range(m.first, m.first));
1248 auto adj_tets = get_adj(Range(max_face, max_face), 3);
1249 auto adj_tets_faces = get_adj(adj_tets, 2);
1250 auto adj_front_faces = subtract(
1251 intersect(get_adj(Range(edge, edge), 2), adj_tets_faces),
1252 *crackFaces);
1253 if (adj_front_faces.size() > 3)
1254 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
1255 "adj_front_faces.size()>3");
1256
1257 FTensor::Tensor1<double, SPACE_DIM> t_material_force;
1258 CHKERR mField.get_moab().tag_get_data(th_material_force, &edge, 1,
1259 &t_material_force(0));
1260 std::vector<double> griffith_energy(adj_front_faces.size());
1261 CHKERR mField.get_moab().tag_get_data(
1262 th_face_energy, adj_front_faces, griffith_energy.data());
1263
1264 auto set_edge_bundle = [&](auto min_gamma) {
1265 for (auto rotated_f : adj_front_faces) {
1266
1267 double rotated_face_energy =
1268 griffith_energy[adj_front_faces.index(rotated_f)];
1269
1270 auto vertex = subtract(get_conn(Range(rotated_f, rotated_f)),
1271 front_vertex);
1272 if (vertex.size() != 1) {
1273 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
1274 "Wrong number of vertex to move");
1275 }
1276 auto front_vertex_edges_vertex = get_conn(
1277 intersect(get_adj(front_vertex, 1), crack_faces_edges));
1278 vertex = subtract(
1279 vertex, front_vertex_edges_vertex); // vertex free to move
1280 if (vertex.empty()) {
1281 continue;
1282 }
1283
1284 auto face_cardinality = [&](auto f, auto &seen_front_edges) {
1285 auto whole_front =
1286 unite(*frontEdges,
1287 subtract(body_skin_edges, crack_faces_edges));
1288 auto faces = Range(f, f);
1289 int c = 0;
1290 for (; c < 10; ++c) {
1291 auto front_edges =
1292 subtract(get_adj(faces, 1), seen_front_edges);
1293 if (front_edges.size() == 0) {
1294 return 0;
1295 }
1296 auto front_connected_edges =
1297 intersect(front_edges, whole_front);
1298 if (front_connected_edges.size()) {
1299 seen_front_edges.merge(front_connected_edges);
1300 return c;
1301 }
1302 faces.merge(get_adj(front_edges, 2));
1303 ++c;
1304 }
1305 return c;
1306 };
1307
1308 Range seen_edges = Range(edge, edge);
1309 double rotated_face_cardinality = face_cardinality(
1310 rotated_f,
1311 seen_edges); // add cardinality of max energy
1312 // face to rotated face cardinality
1313 // rotated_face_cardinality +=
1314 // face_cardinality(max_face, seen_edges);
1315 rotated_face_cardinality = std::max(rotated_face_cardinality,
1316 1.); // at least one edge
1317
1318 auto t_vertex_coords = get_coords(vertex);
1319 auto vertex_edges = get_vertex_edges(vertex);
1320
1321 EntityHandle f0 = front_vertex[0];
1322 EntityHandle f1 = front_vertex[1];
1323 FTensor::Tensor1<double, 3> t_v_e0, t_v_e1;
1324 CHKERR mField.get_moab().get_coords(&f0, 1, &t_v_e0(0));
1325 CHKERR mField.get_moab().get_coords(&f1, 1, &t_v_e1(0));
1326
1327 FTENSOR_INDEX(SPACE_DIM, i);
1328 for (auto e_used_to_move_detection : vertex_edges) {
1329 auto edge_conn = get_conn(Range(e_used_to_move_detection,
1330 e_used_to_move_detection));
1331 edge_conn = subtract(edge_conn, vertex);
1332 // Find displacement of the edge such that dot porduct with
1333 // normal is zero.
1334 //
1335 // { (t_v0 - t_vertex_coords) + gamma * (t_v3 -
1336 // t_vertex_coords) } * n = 0
1337 // where t_v0 is the edge vertex, t_v3 is the edge end
1338 // point, n is the rotated normal of the face gamma is the
1339 // factor by which the edge is moved
1340 FTensor::Tensor1<double, 3> t_v0;
1341 t_v0(i) = (t_v_e0(i) + t_v_e1(i)) / 2;
1342 FTensor::Tensor1<double, 3> t_v3;
1343 CHKERR mField.get_moab().get_coords(edge_conn, &t_v3(0));
1344 auto a =
1345 (t_v0(i) - t_vertex_coords(i)) * t_rotated_normal(i);
1346 auto b =
1347 (t_v3(i) - t_vertex_coords(i)) * t_rotated_normal(i);
1348 auto gamma = a / b;
1349
1350 constexpr double eps =
1351 std::numeric_limits<double>::epsilon();
1352 if (std::isnormal(gamma) && gamma < 1.0 - eps &&
1353 gamma > -0.1) {
1354 FTensor::Tensor1<double, 3> t_move;
1355 t_move(i) = gamma * (t_v3(i) - t_vertex_coords(i));
1356
1357 auto check_rotated_face_directoon = [&]() {
1358 FTensor::Tensor1<double, SPACE_DIM> t_delta;
1359 t_delta(i) = t_vertex_coords(i) + t_move(i) - t_v0(i);
1360 t_delta.normalize();
1361 auto dot =
1362 (t_material_force(i) / t_material_force.l2()) *
1363 t_delta(i);
1364 return -dot > 0 ? true : false;
1365 };
1366
1367 if (check_rotated_face_directoon()) {
1368
1369 MOFEM_LOG("EP", Sev::inform)
1370 << "Crack edge " << edge << " moved face "
1371 << rotated_f
1372 << " edge: " << e_used_to_move_detection
1373 << " face direction/energy " << rotated_face_energy
1374 << " face cardinality " << rotated_face_cardinality
1375 << " gamma: " << gamma;
1376
1377 auto &bundle = edge_bundle_map[edge];
1378 bundle.emplace_back(rotated_f, e_used_to_move_detection,
1379 vertex[0], t_move, 1,
1380 rotated_face_cardinality, gamma);
1381 }
1382 }
1383 }
1384 }
1385 };
1386
1387 set_edge_bundle(std::numeric_limits<double>::epsilon());
1388 if (edge_bundle_map[edge].empty()) {
1389 set_edge_bundle(-1.);
1390 }
1391 }
1392
1393 return edge_bundle_map;
1394 };
1395
1396 auto get_sort_by_energy = [&](auto &edge_face_max_energy_map) {
1397 std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
1398 sort_by_energy;
1399
1400 for (auto &m : edge_face_max_energy_map) {
1401 auto e = m.first;
1402 auto &[max_face, energy, opt_angle] = m.second;
1403 sort_by_energy[energy] = std::make_tuple(e, max_face, opt_angle);
1404 }
1405
1406 return sort_by_energy;
1407 };
1408
1409 auto set_tag = [&](auto &&adj_edges_map, auto &&sort_by_energy) {
1410 MoFEMFunctionBegin;
1411
1412 Tag th_face_pressure;
1413 CHK_MOAB_THROW(
1414 mField.get_moab().tag_get_handle("FacePressure", th_face_pressure),
1415 "get tag");
1416 auto get_face_pressure = [&](auto face) {
1417 double pressure;
1418 CHK_MOAB_THROW(mField.get_moab().tag_get_data(th_face_pressure, &face,
1419 1, &pressure),
1420 "get rag data");
1421 return pressure;
1422 };
1423
1424 MOFEM_LOG("EPSELF", Sev::inform)
1425 << "Number of edges to check " << sort_by_energy.size();
1426
1427 enum face_energy { POSITIVE, NEGATIVE };
1428 constexpr bool skip_negative = true;
1429
1430 for (auto fe : {face_energy::POSITIVE, face_energy::NEGATIVE}) {
1431
1432 // iterate edges wih maximal energy, and make them seed. Such edges,
1433 // will most likely will have also smallest node displacement
1434 for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
1435 ++it) {
1436
1437 auto energy = it->first;
1438 auto [max_edge, max_face, opt_angle] = it->second;
1439
1440 auto face_pressure = get_face_pressure(max_face);
1441 if (skip_negative) {
1442 if (fe == face_energy::POSITIVE) {
1443 if (face_pressure < 0) {
1444 MOFEM_LOG("EPSELF", Sev::inform)
1445 << "Skip negative face " << max_face << " with energy "
1446 << energy << " and pressure " << face_pressure;
1447 continue;
1448 }
1449 }
1450 }
1451
1452 MOFEM_LOG("EPSELF", Sev::inform)
1453 << "Check face " << max_face << " edge " << max_edge
1454 << " energy " << energy << " optimal angle " << opt_angle
1455 << " face pressure " << face_pressure;
1456
1457 auto jt = adj_edges_map.find(max_edge);
1458 if (jt == adj_edges_map.end()) {
1459 MOFEM_LOG("EPSELF", Sev::warning)
1460 << "Edge " << max_edge << " not found in adj_edges_map";
1461 continue;
1462 }
1463 auto &bundle = jt->second;
1464
1465 auto find_max_in_bundle_impl = [&](auto edge, auto &bundle,
1466 auto gamma) {
1467 FTENSOR_INDEX(SPACE_DIM, i);
1468
1469 EntityHandle vertex_max = 0;
1470 EntityHandle face_max = 0;
1471 EntityHandle move_edge_max = 0;
1472 double max_quality = -2;
1473 double max_quality_evaluated = -2;
1474 double min_cardinality = std::numeric_limits<double>::max();
1475
1476 FTensor::Tensor1<double, SPACE_DIM> t_move_last{0., 0., 0.};
1477
1478 for (auto &b : bundle) {
1479 auto &[face, move_edge, vertex, t_move, quality, cardinality,
1480 edge_gamma] = b;
1481
1482 auto adj_vertex_tets = get_adj(Range(vertex, vertex), 3);
1483 auto adj_vertex_tets_verts = get_conn(adj_vertex_tets);
1484 std::vector<double> coords(3 * adj_vertex_tets_verts.size());
1485 CHK_MOAB_THROW(mField.get_moab().get_coords(
1486 adj_vertex_tets_verts, coords.data()),
1487 "get coords");
1488
1489 set_coord(vertex, adj_vertex_tets_verts, coords, t_move, gamma);
1490 quality = tets_quality(quality, adj_vertex_tets_verts,
1491 adj_vertex_tets, coords);
1492
1493 auto eval_quality = [](auto q, auto c, auto edge_gamma) {
1494 if (q < 0) {
1495 return q;
1496 } else {
1497 return ((edge_gamma < 0) ? (q / 2) : q) / pow(c, 2);
1498 }
1499 };
1500
1501 if (eval_quality(quality, cardinality, edge_gamma) >=
1502 max_quality_evaluated) {
1503 max_quality = quality;
1504 min_cardinality = cardinality;
1505 vertex_max = vertex;
1506 face_max = face;
1507 move_edge_max = move_edge;
1508 t_move_last(i) = t_move(i);
1509 max_quality_evaluated =
1510 eval_quality(max_quality, min_cardinality, edge_gamma);
1511 }
1512 }
1513
1514 return std::make_tuple(vertex_max, face_max, t_move_last,
1515 max_quality, min_cardinality);
1516 };
1517
1518 auto find_max_in_bundle = [&](auto edge, auto &bundle) {
1519 auto b_org_bundle = bundle;
1520 auto r = find_max_in_bundle_impl(edge, bundle, 1.);
1521 auto &[vertex_max, face_max, t_move_last, max_quality,
1522 cardinality] = r;
1523 if (max_quality < 0) {
1524 for (double gamma = 0.95; gamma >= 0.45; gamma -= 0.05) {
1525 bundle = b_org_bundle;
1526 r = find_max_in_bundle_impl(edge, bundle, gamma);
1527 auto &[vertex_max, face_max, t_move_last, max_quality,
1528 cardinality] = r;
1529 MOFEM_LOG("EPSELF", Sev::warning)
1530 << "Back tracking: gamma " << gamma << " edge " << edge
1531 << " quality " << max_quality << " cardinality "
1532 << cardinality;
1533 if (max_quality > 0.01) {
1534 FTENSOR_INDEX(SPACE_DIM, I);
1535 t_move_last(I) *= gamma;
1536 return r;
1537 }
1538 }
1539 FTENSOR_INDEX(SPACE_DIM, I);
1540 t_move_last(I) = 0;
1541 }
1542 return r;
1543 };
1544
1545 // set tags with displacement of node and face energy
1546 auto set_tag_to_vertex_and_face = [&](auto &&r, auto &quality) {
1547 MoFEMFunctionBegin;
1548 auto &[v, f, t_move, q, cardinality] = r;
1549
1550 if ((q > 0 && std::isnormal(q)) && energy > 0) {
1551
1552 MOFEM_LOG("EPSELF", Sev::inform)
1553 << "Set tag: vertex " << v << " face " << f << " "
1554 << max_edge << " move " << t_move << " energy " << energy
1555 << " quality " << q << " cardinality " << cardinality;
1556 CHKERR mField.get_moab().tag_set_data(th_position[0], &v, 1,
1557 &t_move(0));
1558 CHKERR mField.get_moab().tag_set_data(th_max_face_energy[0], &f,
1559 1, &energy);
1560 }
1561
1562 quality = q;
1563 MoFEMFunctionReturn(0);
1564 };
1565
1566 double quality = -2;
1567 CHKERR set_tag_to_vertex_and_face(
1568
1569 find_max_in_bundle(max_edge, bundle),
1570
1571 quality
1572
1573 );
1574
1575 if (quality > 0 && std::isnormal(quality) && energy > 0) {
1576 MOFEM_LOG("EPSELF", Sev::inform)
1577 << "Crack face set with quality: " << quality;
1578 MoFEMFunctionReturnHot(0);
1579 }
1580 }
1581
1582 if (!skip_negative)
1583 break;
1584 }
1585
1586 MoFEMFunctionReturn(0);
1587 };
1588
1589 // map: {edge, {face, energy, optimal_angle}}
1590 MOFEM_LOG("EP", sev) << "Calculate orientation";
1591 std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
1592 edge_face_max_energy_map;
1593 CHKERR find_maximal_face_energy(front_edges, front_faces,
1594 edge_face_max_energy_map);
1595 CHKERR calculate_face_orientation(edge_face_max_energy_map);
1596
1597 MOFEM_LOG("EP", sev) << "Calculate node positions";
1598 CHKERR set_tag(
1599
1600 calculate_free_face_node_displacement(edge_face_max_energy_map),
1601 get_sort_by_energy(edge_face_max_energy_map)
1602
1603 );
1604
1605 MoFEMFunctionReturn(0);
1606 };
1607
1608 MOFEM_LOG("EP", sev) << "Front edges " << frontEdges->size();
1609 CHKERR evaluate_face_energy_and_set_orientation(
1610 *frontEdges, get_adj_front(false), sides_pair, th_front_position);
1611 }
1612
1613 // exchange positions and energies from processor zero to all other
1614 CHKERR VecZeroEntries(vertexExchange.second);
1615 CHKERR VecGhostUpdateBegin(vertexExchange.second, INSERT_VALUES,
1616 SCATTER_FORWARD);
1617 CHKERR VecGhostUpdateEnd(vertexExchange.second, INSERT_VALUES,
1618 SCATTER_FORWARD);
1619 CHKERR mField.getInterface<CommInterface>()->updateEntitiesPetscVector(
1620 mField.get_moab(), vertexExchange, th_front_position[0]);
1621 CHKERR VecZeroEntries(faceExchange.second);
1622 CHKERR VecGhostUpdateBegin(faceExchange.second, INSERT_VALUES,
1623 SCATTER_FORWARD);
1624 CHKERR VecGhostUpdateEnd(faceExchange.second, INSERT_VALUES, SCATTER_FORWARD);
1625 CHKERR mField.getInterface<CommInterface>()->updateEntitiesPetscVector(
1626 mField.get_moab(), faceExchange, th_max_face_energy[0]);
1627
1628 auto get_max_moved_faces = [&]() {
1629 Range max_moved_faces;
1630 auto adj_front = get_adj_front(false);
1631 std::vector<double> face_energy(adj_front.size());
1632 CHKERR mField.get_moab().tag_get_data(th_max_face_energy[0], adj_front,
1633 face_energy.data());
1634 for (int i = 0; i != adj_front.size(); ++i) {
1635 if (face_energy[i] > std::numeric_limits<double>::epsilon()) {
1636 max_moved_faces.insert(adj_front[i]);
1637 }
1638 }
1639
1640 return boost::make_shared<Range>(max_moved_faces);
1641 };
1642
1643 // move all faces with energy larger than 0
1644 maxMovedFaces = get_max_moved_faces();
1645 MOFEM_LOG("EP", sev) << "Number of of moved faces: " << maxMovedFaces->size();
1646
1647#ifndef NDEBUG
1648 if (debug) {
1649 CHKERR save_range(
1650 mField.get_moab(),
1651 "max_moved_faces_" +
1652 boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
1653 *maxMovedFaces);
1654 }
1655#endif
1656
1657 MoFEMFunctionReturn(0);
1658}
1659
1660MoFEMErrorCode EshelbianCore::setNewFrontCoordinates() {
1661 MoFEMFunctionBegin;
1662
1663 if (!maxMovedFaces)
1664 MoFEMFunctionReturnHot(0);
1665
1666 Tag th_front_position;
1667 auto rval =
1668 mField.get_moab().tag_get_handle("FrontPosition", th_front_position);
1669 if (rval == MB_SUCCESS && maxMovedFaces) {
1670 Range verts;
1671 CHKERR mField.get_moab().get_connectivity(*maxMovedFaces, verts, true);
1672 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(verts);
1673 std::vector<double> coords(3 * verts.size());
1674 CHKERR mField.get_moab().get_coords(verts, coords.data());
1675 std::vector<double> pos(3 * verts.size());
1676 CHKERR mField.get_moab().tag_get_data(th_front_position, verts, pos.data());
1677 for (int i = 0; i != 3 * verts.size(); ++i) {
1678 coords[i] += pos[i];
1679 }
1680 CHKERR mField.get_moab().set_coords(verts, coords.data());
1681 double zero[] = {0., 0., 0.};
1682 CHKERR mField.get_moab().tag_clear_data(th_front_position, verts, zero);
1683 }
1684
1685#ifndef NDEBUG
1686 constexpr bool debug = false;
1687 if (debug) {
1688
1689 CHKERR save_range(
1690 mField.get_moab(),
1691 "set_coords_faces_" +
1692 boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
1693 *maxMovedFaces);
1694 }
1695#endif
1696 MoFEMFunctionReturn(0);
1697}
1698
1699MoFEMErrorCode EshelbianCore::addCrackSurfaces(const bool debug) {
1700 MoFEMFunctionBegin;
1701
1702 constexpr bool potential_crack_debug = false;
1703 if constexpr (potential_crack_debug) {
1704
1705 auto add_ents = get_range_from_block(mField, "POTENTIAL", SPACE_DIM - 1);
1706 Range crack_front_verts;
1707 CHKERR mField.get_moab().get_connectivity(*frontEdges, crack_front_verts,
1708 true);
1709 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1710 crack_front_verts);
1711 Range crack_front_faces;
1712 CHKERR mField.get_moab().get_adjacencies(crack_front_verts, SPACE_DIM - 1,
1713 true, crack_front_faces,
1714 moab::Interface::UNION);
1715 crack_front_faces = intersect(crack_front_faces, add_ents);
1716 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
1717 crack_front_faces);
1718 CHKERR mField.getInterface<MeshsetsManager>()->addEntitiesToMeshset(
1719 BLOCKSET, addCrackMeshsetId, crack_front_faces);
1720 }
1721
1722 auto get_crack_faces = [&]() {
1723 if (maxMovedFaces) {
1724 return unite(*crackFaces, *maxMovedFaces);
1725 } else {
1726 return *crackFaces;
1727 }
1728 };
1729
1730 auto get_extended_crack_faces = [&]() {
1731 auto get_faces_of_crack_front_verts = [&](auto crack_faces_org) {
1732 ParallelComm *pcomm =
1733 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
1734
1735 Range crack_faces;
1736
1737 if (!pcomm->rank()) {
1738
1739 auto get_nodes = [&](auto &&e) {
1740 Range nodes;
1741 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, nodes, true),
1742 "get connectivity");
1743 return nodes;
1744 };
1745
1746 auto get_adj = [&](auto &&e, auto dim,
1747 auto t = moab::Interface::UNION) {
1748 Range adj;
1749 CHK_MOAB_THROW(
1750 mField.get_moab().get_adjacencies(e, dim, true, adj, t),
1751 "get adj");
1752 return adj;
1753 };
1754
1755 Range body_ents;
1756 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
1757 body_ents);
1758 auto body_skin = get_skin(mField, body_ents);
1759 auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
1760 auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
1761 auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
1762 auto front_block_nodes = get_nodes(front_block_edges);
1763
1764 size_t s;
1765 do {
1766 s = crack_faces.size();
1767
1768 auto crack_face_nodes = get_nodes(crack_faces_org);
1769 auto crack_faces_edges =
1770 get_adj(crack_faces_org, 1, moab::Interface::UNION);
1771
1772 auto crack_skin = get_skin(mField, crack_faces_org);
1773 front_block_edges = subtract(front_block_edges, crack_skin);
1774 auto crack_skin_nodes = get_nodes(crack_skin);
1775 crack_skin_nodes.merge(front_block_nodes);
1776
1777 auto crack_skin_faces =
1778 get_adj(crack_skin, 2, moab::Interface::UNION);
1779 crack_skin_faces =
1780 subtract(subtract(crack_skin_faces, crack_faces_org), body_skin);
1781
1782 crack_faces = crack_faces_org;
1783 for (auto f : crack_skin_faces) {
1784 auto edges = intersect(
1785 get_adj(Range(f, f), 1, moab::Interface::UNION), crack_skin);
1786
1787 // if other edge is part of body skin, e.g. crack punching through
1788 // body surface
1789 if (edges.size() == 2) {
1790 edges.merge(
1791 intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
1792 body_skin_edges));
1793 }
1794
1795 if (edges.size() == 2) {
1796 auto edge_conn = get_nodes(Range(edges));
1797 auto faces = intersect(get_adj(edges, 2, moab::Interface::UNION),
1798 crack_faces_org);
1799 if (faces.size() == 2) {
1800 auto edge0_conn = get_nodes(Range(edges[0], edges[0]));
1801 auto edge1_conn = get_nodes(Range(edges[1], edges[1]));
1802 auto edges_conn = intersect(intersect(edge0_conn, edge1_conn),
1803 crack_skin_nodes); // node at apex
1804 if (edges_conn.size() == 1) {
1805
1806 auto node_edges =
1807 subtract(intersect(get_adj(edges_conn, 1,
1808 moab::Interface::INTERSECT),
1809 crack_faces_edges),
1810 crack_skin); // nodes on crack surface, but not
1811 // at the skin
1812
1813 if (node_edges.size()) {
1814 FTENSOR_INDEX(SPACE_DIM, i);
1815 FTensor::Tensor1<double, SPACE_DIM> t_v0;
1816 CHKERR mField.get_moab().get_coords(edges_conn, &t_v0(0));
1817
1818 auto get_t_dir = [&](auto e_conn) {
1819 auto other_node = subtract(e_conn, edges_conn);
1820 FTensor::Tensor1<double, SPACE_DIM> t_dir;
1821 CHKERR mField.get_moab().get_coords(other_node,
1822 &t_dir(0));
1823 t_dir(i) -= t_v0(i);
1824 return t_dir;
1825 };
1826
1827 FTensor::Tensor1<double, SPACE_DIM> t_ave_dir;
1828 t_ave_dir(i) =
1829 get_t_dir(edge0_conn)(i) + get_t_dir(edge1_conn)(i);
1830
1831 FTensor::Tensor1<double, SPACE_DIM> t_crack_surface_ave_dir;
1832 t_crack_surface_ave_dir(i) = 0;
1833 for (auto e : node_edges) {
1834 auto e_conn = get_nodes(Range(e, e));
1835 auto t_dir = get_t_dir(e_conn);
1836 t_crack_surface_ave_dir(i) += t_dir(i);
1837 }
1838
1839 auto dot = t_ave_dir(i) * t_crack_surface_ave_dir(i);
1840 // ave edges is in opposite direction to crack surface, so
1841 // thus crack is not turning back
1842 if (dot < -std::numeric_limits<double>::epsilon()) {
1843 crack_faces.insert(f);
1844 }
1845 } else {
1846 crack_faces.insert(f);
1847 }
1848 }
1849 }
1850 } else if (edges.size() == 3) {
1851 crack_faces.insert(f);
1852 }
1853
1854 // if other edge is part of geometry edge, e.g. keyway
1855 if (edges.size() == 1) {
1856 edges.merge(
1857 intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
1858 geometry_edges));
1859 edges.merge(
1860 intersect(get_adj(Range(f, f), 1, moab::Interface::UNION),
1861 front_block_edges));
1862 if (edges.size() == 2) {
1863 crack_faces.insert(f);
1864 continue;
1865 }
1866 }
1867 }
1868
1869 crack_faces_org = crack_faces;
1870
1871 } while (s != crack_faces.size());
1872 };
1873
1874 return crack_faces; // send_type(mField, crack_faces, MBTRI);
1875 };
1876
1877 return get_faces_of_crack_front_verts(get_crack_faces());
1878 };
1879
1880 if (debug) {
1881 CHKERR save_range(mField.get_moab(), "new_crack_surface_debug.vtk",
1882 get_extended_crack_faces());
1883 }
1884
1885 auto reconstruct_crack_faces = [&](auto crack_faces) {
1886 ParallelComm *pcomm =
1887 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
1888
1889 auto impl = [&]() {
1890 MoFEMFunctionBegin;
1891
1892 Range new_crack_faces;
1893 if (!pcomm->rank()) {
1894
1895 auto get_nodes = [&](auto &&e) {
1896 Range nodes;
1897 CHK_MOAB_THROW(mField.get_moab().get_connectivity(e, nodes, true),
1898 "get connectivity");
1899 return nodes;
1900 };
1901
1902 auto get_adj = [&](auto &&e, auto dim,
1903 auto t = moab::Interface::UNION) {
1904 Range adj;
1905 CHK_MOAB_THROW(
1906 mField.get_moab().get_adjacencies(e, dim, true, adj, t),
1907 "get adj");
1908 return adj;
1909 };
1910
1911 auto get_test_on_crack_surface = [&]() {
1912 auto crack_faces_nodes =
1913 get_nodes(crack_faces); // nodes on crac faces
1914 auto crack_faces_tets =
1915 get_adj(crack_faces_nodes, 3,
1916 moab::Interface::UNION); // adjacent
1917 // tets to
1918 // crack
1919 // faces throug nodes
1920 auto crack_faces_tets_nodes =
1921 get_nodes(crack_faces_tets); // nodes on crack faces tets
1922 crack_faces_tets_nodes =
1923 subtract(crack_faces_tets_nodes, crack_faces_nodes);
1924 crack_faces_tets =
1925 subtract(crack_faces_tets, get_adj(crack_faces_tets_nodes, 3,
1926 moab::Interface::UNION));
1927 new_crack_faces =
1928 get_adj(crack_faces_tets, 2,
1929 moab::Interface::UNION); // adjacency faces to crack
1930 // faces through tets
1931 new_crack_faces.merge(crack_faces); // add original crack faces
1932
1933 return std::make_tuple(new_crack_faces, crack_faces_tets);
1934 };
1935
1936 auto carck_faces_test_edges = [&](auto faces, auto tets) {
1937 auto adj_tets_faces = get_adj(tets, 2, moab::Interface::UNION);
1938 auto adj_faces_edges = get_adj(subtract(faces, adj_tets_faces), 1,
1939 moab::Interface::UNION);
1940 auto adj_tets_edges = get_adj(tets, 1, moab::Interface::UNION);
1941 auto geometry_edges = get_range_from_block(mField, "EDGES", 1);
1942 auto front_block_edges = get_range_from_block(mField, "FRONT", 1);
1943 adj_faces_edges.merge(geometry_edges); // geometry edges
1944 adj_faces_edges.merge(front_block_edges); // front block edges
1945
1946 auto boundary_tets_edges = intersect(adj_tets_edges, adj_faces_edges);
1947 auto boundary_test_nodes = get_nodes(boundary_tets_edges);
1948 auto boundary_test_nodes_edges =
1949 get_adj(boundary_test_nodes, 1, moab::Interface::UNION);
1950 auto boundary_test_nodes_edges_nodes = subtract(
1951 get_nodes(boundary_test_nodes_edges), boundary_test_nodes);
1952
1953 boundary_tets_edges =
1954 subtract(boundary_test_nodes_edges,
1955 get_adj(boundary_test_nodes_edges_nodes, 1,
1956 moab::Interface::UNION));
1957
1958 Range body_ents;
1959 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
1960 body_ents);
1961 auto body_skin = get_skin(mField, body_ents);
1962
1963 auto body_skin_edges = get_adj(body_skin, 1, moab::Interface::UNION);
1964 body_skin_edges = intersect(get_adj(tets, 1, moab::Interface::UNION),
1965 body_skin_edges);
1966 body_skin = intersect(body_skin, adj_tets_faces);
1967 body_skin_edges = subtract(
1968 body_skin_edges, get_adj(body_skin, 1, moab::Interface::UNION));
1969
1970 save_range(mField.get_moab(), "body_skin_edges.vtk", body_skin_edges);
1971 for (auto e : body_skin_edges) {
1972 auto adj_tet = intersect(
1973 get_adj(Range(e, e), 3, moab::Interface::INTERSECT), tets);
1974 if (adj_tet.size() == 1) {
1975 boundary_tets_edges.insert(e);
1976 }
1977 }
1978
1979 return boundary_tets_edges;
1980 };
1981
1982 auto p = get_test_on_crack_surface();
1983 auto &[new_crack_faces, crack_faces_tets] = p;
1984
1985 if (debug) {
1986 CHKERR save_range(mField.get_moab(), "hole_crack_faces_debug.vtk",
1987 crack_faces);
1988 CHKERR save_range(mField.get_moab(), "new_crack_faces_debug.vtk",
1989 new_crack_faces);
1990 CHKERR save_range(mField.get_moab(), "new_crack_tets_debug.vtk",
1991 crack_faces_tets);
1992 }
1993
1994 auto boundary_tets_edges =
1995 carck_faces_test_edges(new_crack_faces, crack_faces_tets);
1996 CHKERR save_range(mField.get_moab(), "boundary_tets_edges.vtk",
1997 boundary_tets_edges);
1998
1999 auto resolve_surface = [&](auto boundary_tets_edges,
2000 auto crack_faces_tets) {
2001 auto boundary_tets_edges_nodes = get_nodes(boundary_tets_edges);
2002 auto crack_faces_tets_faces =
2003 get_adj(crack_faces_tets, 2, moab::Interface::UNION);
2004
2005 Range all_removed_faces;
2006 Range all_removed_tets;
2007 int counter = 0;
2008
2009 int size = 0;
2010 while (size != crack_faces_tets.size()) {
2011 auto tets_faces =
2012 get_adj(crack_faces_tets, 2, moab::Interface::UNION);
2013 auto skin_tets = get_skin(mField, crack_faces_tets);
2014 auto skin_skin =
2015 get_skin(mField, subtract(crack_faces_tets_faces, tets_faces));
2016 auto skin_skin_nodes = get_nodes(skin_skin);
2017
2018 size = crack_faces_tets.size();
2019 MOFEM_LOG("SELF", Sev::inform)
2020 << "Crack faces tets size " << crack_faces_tets.size()
2021 << " crack faces size " << crack_faces_tets_faces.size();
2022 auto skin_tets_nodes = subtract(
2023 get_nodes(skin_tets),
2024 boundary_tets_edges_nodes); // not remove tets which are
2025 // adjagasent to crack faces nodes
2026 skin_tets_nodes = subtract(skin_tets_nodes, skin_skin_nodes);
2027
2028 Range removed_nodes;
2029 Range tets_to_remove;
2030 Range faces_to_remove;
2031 for (auto n : skin_tets_nodes) {
2032 auto tets =
2033 intersect(get_adj(Range(n, n), 3, moab::Interface::INTERSECT),
2034 crack_faces_tets);
2035 if (tets.size() == 0) {
2036 continue;
2037 }
2038
2039 auto hole_detetction = [&]() {
2040 auto adj_tets =
2041 get_adj(Range(n, n), 3, moab::Interface::INTERSECT);
2042 adj_tets =
2043 subtract(adj_tets,
2044 crack_faces_tets); // tetst adjacent to the node
2045 // but not part of crack surface
2046 if (adj_tets.size() == 0) {
2047 return std::make_pair(
2048 intersect(
2049 get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
2050 tets_faces),
2051 tets);
2052 }
2053
2054 std::vector<Range> tets_groups;
2055 auto test_adj_tets = adj_tets;
2056 while (test_adj_tets.size()) {
2057 auto seed_size = 0;
2058 Range seed = Range(test_adj_tets[0], test_adj_tets[0]);
2059 while (seed.size() != seed_size) {
2060 auto adj_faces =
2061 subtract(get_adj(seed, 2, moab::Interface::UNION),
2062 tets_faces); // edges which are not
2063 // part of the node
2064 seed_size = seed.size();
2065 seed.merge(
2066 intersect(get_adj(adj_faces, 3, moab::Interface::UNION),
2067 test_adj_tets));
2068 }
2069 tets_groups.push_back(seed);
2070 test_adj_tets = subtract(test_adj_tets, seed);
2071 }
2072 if (tets_groups.size() == 1) {
2073
2074 return std::make_pair(
2075 intersect(
2076 get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
2077 tets_faces),
2078 tets);
2079 }
2080
2081 Range tets_to_remove;
2082 Range faces_to_remove;
2083 for (auto &r : tets_groups) {
2084 auto f = get_adj(r, 2, moab::Interface::UNION);
2085 auto t = intersect(get_adj(f, 3, moab::Interface::UNION),
2086 crack_faces_tets); // tets
2087
2088 if (f.size() > faces_to_remove.size() ||
2089 faces_to_remove.size() == 0) {
2090 faces_to_remove = f;
2091 tets_to_remove = t; // largest group of tets
2092 }
2093 }
2094 MOFEM_LOG("EPSELF", Sev::inform)
2095 << "Hole detection: faces to remove "
2096 << faces_to_remove.size() << " tets to remove "
2097 << tets_to_remove.size();
2098 return std::make_pair(faces_to_remove, tets_to_remove);
2099 };
2100
2101 if (tets.size() < tets_to_remove.size() ||
2102 tets_to_remove.size() == 0) {
2103 removed_nodes = Range(n, n);
2104 auto [h_faces_to_remove, h_tets_to_remove] =
2105 hole_detetction(); // find faces and tets to remove
2106 faces_to_remove = h_faces_to_remove;
2107 tets_to_remove = h_tets_to_remove;
2108
2109 // intersect(
2110 // get_adj(Range(n, n), 2, moab::Interface::INTERSECT),
2111 // tets_faces);
2112
2113 } // find tets which is largest adjacencty size, so that it is
2114 // removed first, and then faces are removed
2115 all_removed_faces.merge(faces_to_remove);
2116 all_removed_tets.merge(tets_to_remove);
2117 }
2118
2119 crack_faces_tets = subtract(crack_faces_tets, tets_to_remove);
2120 crack_faces_tets_faces =
2121 subtract(crack_faces_tets_faces, faces_to_remove);
2122
2123 if (debug) {
2124 save_range(mField.get_moab(),
2125 "removed_nodes_" +
2126 boost::lexical_cast<std::string>(counter) + ".vtk",
2127 removed_nodes);
2128 save_range(mField.get_moab(),
2129 "faces_to_remove_" +
2130 boost::lexical_cast<std::string>(counter) + ".vtk",
2131 faces_to_remove);
2132 save_range(mField.get_moab(),
2133 "tets_to_remove_" +
2134 boost::lexical_cast<std::string>(counter) + ".vtk",
2135 tets_to_remove);
2136 save_range(mField.get_moab(),
2137 "crack_faces_tets_faces_" +
2138 boost::lexical_cast<std::string>(counter) + ".vtk",
2139 crack_faces_tets_faces);
2140 save_range(mField.get_moab(),
2141 "crack_faces_tets_" +
2142 boost::lexical_cast<std::string>(counter) + ".vtk",
2143 crack_faces_tets);
2144 }
2145 counter++;
2146 }
2147
2148 auto cese_internal_faces = [&]() {
2149 MoFEMFunctionBegin;
2150 auto skin_tets = get_skin(mField, crack_faces_tets);
2151 auto adj_faces = get_adj(skin_tets, 2, moab::Interface::UNION);
2152 adj_faces =
2153 subtract(adj_faces, skin_tets); // remove skin tets faces
2154 auto adj_tets = get_adj(adj_faces, 3,
2155 moab::Interface::UNION); // tets which are
2156 // adjacent to skin
2157 crack_faces_tets =
2158 subtract(crack_faces_tets,
2159 adj_tets); // remove tets which are adjacent to
2160 // skin, so that they are not removed
2161 crack_faces_tets_faces =
2162 subtract(crack_faces_tets_faces, adj_faces);
2163
2164 all_removed_faces.merge(adj_faces);
2165 all_removed_tets.merge(adj_tets);
2166
2167 MOFEM_LOG("EPSELF", Sev::inform)
2168 << "Remove internal faces size " << adj_faces.size()
2169 << " tets size " << adj_tets.size();
2170 MoFEMFunctionReturn(0);
2171 };
2172
2173 auto case_only_one_free_edge = [&]() {
2174 MoFEMFunctionBegin;
2175
2176 for (auto t : Range(crack_faces_tets)) {
2177
2178 auto adj_faces = get_adj(
2179 Range(t, t), 2,
2180 moab::Interface::UNION); // faces of tet which can be removed
2181 auto crack_surface_edges =
2182 get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
2183 adj_faces),
2184 1,
2185 moab::Interface::UNION); // edges not on the tet but
2186 // on crack surface
2187 auto adj_edges =
2188 subtract(get_adj(Range(t, t), 1, moab::Interface::INTERSECT),
2189 crack_surface_edges); // free edges
2190 adj_edges = subtract(
2191 adj_edges,
2192 boundary_tets_edges); // edges which are not part of gemetry
2193
2194 if (adj_edges.size() == 1) {
2195 crack_faces_tets =
2196 subtract(crack_faces_tets,
2197 Range(t, t)); // remove tets which are adjacent to
2198 // skin, so that they are not removed
2199
2200 auto faces_to_remove =
2201 get_adj(adj_edges, 2, moab::Interface::UNION); // faces
2202 // which can
2203 // be removed
2204 crack_faces_tets_faces =
2205 subtract(crack_faces_tets_faces, faces_to_remove);
2206
2207 all_removed_faces.merge(faces_to_remove);
2208 all_removed_tets.merge(Range(t, t));
2209
2210 MOFEM_LOG("EPSELF", Sev::inform) << "Remove free one edges ";
2211 }
2212 }
2213
2214 crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
2215 crack_faces_tets_faces =
2216 subtract(crack_faces_tets_faces, all_removed_faces);
2217
2218 MoFEMFunctionReturn(0);
2219 };
2220
2221 auto cese_flat_tet = [&](auto max_adj_edges) {
2222 MoFEMFunctionBegin;
2223
2224 Range body_ents;
2225 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
2226 body_ents);
2227 auto body_skin = get_skin(mField, body_ents);
2228 auto body_skin_edges =
2229 get_adj(body_skin, 1, moab::Interface::UNION);
2230
2231 for (auto t : Range(crack_faces_tets)) {
2232
2233 auto adj_faces = get_adj(
2234 Range(t, t), 2,
2235 moab::Interface::UNION); // faces of tet which can be removed
2236 auto crack_surface_edges =
2237 get_adj(subtract(unite(crack_faces_tets_faces, crack_faces),
2238 adj_faces),
2239 1,
2240 moab::Interface::UNION); // edges not on the tet but
2241 // on crack surface
2242 auto adj_edges =
2243 subtract(get_adj(Range(t, t), 1, moab::Interface::INTERSECT),
2244 crack_surface_edges); // free edges
2245 adj_edges = subtract(adj_edges, body_skin_edges);
2246
2247 auto tet_edges = get_adj(Range(t, t), 1,
2248 moab::Interface::UNION); // edges of
2249 // tet
2250 tet_edges = subtract(tet_edges, adj_edges);
2251
2252 for (auto e : tet_edges) {
2253 constexpr int opposite_edge[] = {5, 3, 4, 1, 2, 0};
2254 auto get_side = [&](auto e) {
2255 int side, sense, offset;
2256 CHK_MOAB_THROW(
2257 mField.get_moab().side_number(t, e, side, sense, offset),
2258 "get side number failed");
2259 return side;
2260 };
2261 auto get_side_ent = [&](auto side) {
2262 EntityHandle side_edge;
2263 CHK_MOAB_THROW(
2264 mField.get_moab().side_element(t, 1, side, side_edge),
2265 "get side failed");
2266 return side_edge;
2267 };
2268 adj_edges.erase(get_side_ent(opposite_edge[get_side(e)]));
2269 }
2270
2271 if (adj_edges.size() <= max_adj_edges) {
2272
2273 double dot = 1;
2274 Range faces_to_remove;
2275 for (auto e : adj_edges) {
2276 auto edge_adj_faces =
2277 get_adj(Range(e, e), 2, moab::Interface::UNION);
2278 edge_adj_faces = intersect(edge_adj_faces, adj_faces);
2279 if (edge_adj_faces.size() != 2) {
2280 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY,
2281 "Adj faces size is not 2 for edge " +
2282 boost::lexical_cast<std::string>(e));
2283 }
2284
2285 auto get_normal = [&](auto f) {
2286 FTensor::Tensor1<double, SPACE_DIM> t_n;
2287 CHK_THROW_MESSAGE(
2288 mField.getInterface<Tools>()->getTriNormal(f, &t_n(0)),
2289 "get tri normal failed");
2290 return t_n;
2291 };
2292 auto t_n0 = get_normal(edge_adj_faces[0]);
2293 auto t_n1 = get_normal(edge_adj_faces[1]);
2294 auto get_sense = [&](auto f) {
2295 int side, sense, offset;
2296 CHK_MOAB_THROW(mField.get_moab().side_number(t, f, side,
2297 sense, offset),
2298 "get side number failed");
2299 return sense;
2300 };
2301 auto sense0 = get_sense(edge_adj_faces[0]);
2302 auto sense1 = get_sense(edge_adj_faces[1]);
2303 t_n0.normalize();
2304 t_n1.normalize();
2305
2306 FTENSOR_INDEX(SPACE_DIM, i);
2307 auto dot_e = (sense0 * sense1) * t_n0(i) * t_n1(i);
2308 if (dot_e < dot || e == adj_edges[0]) {
2309 dot = dot_e;
2310 faces_to_remove = edge_adj_faces;
2311 }
2312 }
2313
2314 all_removed_faces.merge(faces_to_remove);
2315 all_removed_tets.merge(Range(t, t));
2316
2317 MOFEM_LOG("EPSELF", Sev::inform)
2318 << "Remove free edges on flat tet, with considered nb. of "
2319 "edges "
2320 << adj_edges.size();
2321 }
2322 }
2323
2324 crack_faces_tets = subtract(crack_faces_tets, all_removed_tets);
2325 crack_faces_tets_faces =
2326 subtract(crack_faces_tets_faces, all_removed_faces);
2327
2328 MoFEMFunctionReturn(0);
2329 };
2330
2331 CHK_THROW_MESSAGE(case_only_one_free_edge(),
2332 "Case only one free edge failed");
2333 for (auto max_adj_edges : {0, 1, 2, 3}) {
2334 CHK_THROW_MESSAGE(cese_flat_tet(max_adj_edges),
2335 "Case only one free edge failed");
2336 }
2337 CHK_THROW_MESSAGE(cese_internal_faces(),
2338 "Case internal faces failed");
2339
2340 if (debug) {
2341 save_range(mField.get_moab(),
2342 "crack_faces_tets_faces_" +
2343 boost::lexical_cast<std::string>(counter) + ".vtk",
2344 crack_faces_tets_faces);
2345 save_range(mField.get_moab(),
2346 "crack_faces_tets_" +
2347 boost::lexical_cast<std::string>(counter) + ".vtk",
2348 crack_faces_tets);
2349 }
2350
2351 return std::make_tuple(crack_faces_tets_faces, crack_faces_tets,
2352 all_removed_faces, all_removed_tets);
2353 };
2354
2355 auto [resolved_faces, resolved_tets, all_removed_faces,
2356 all_removed_tets] =
2357 resolve_surface(boundary_tets_edges, crack_faces_tets);
2358 resolved_faces.merge(subtract(crack_faces, all_removed_faces));
2359 if (debug) {
2360 CHKERR save_range(mField.get_moab(), "resolved_faces.vtk",
2361 resolved_faces);
2362 CHKERR save_range(mField.get_moab(), "resolved_tets.vtk",
2363 resolved_tets);
2364 }
2365
2366 crack_faces = resolved_faces;
2367 }
2368
2369 MoFEMFunctionReturn(0);
2370 };
2371
2372 CHK_THROW_MESSAGE(impl(), "resolve new crack surfaces");
2373
2374 return crack_faces; // send_type(mField, crack_faces, MBTRI);
2375 };
2376
2377 auto resolve_consisten_crack_extension = [&]() {
2378 MoFEMFunctionBegin;
2379 auto crack_meshset =
2380 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
2381 addCrackMeshsetId, BLOCKSET);
2382 auto meshset = crack_meshset->getMeshset();
2383
2384 if (!mField.get_comm_rank()) {
2385 Range old_crack_faces;
2386 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTRI,
2387 old_crack_faces);
2388 auto extendeded_crack_faces = get_extended_crack_faces();
2389 auto reconstructed_crack_faces =
2390 subtract(reconstruct_crack_faces(extendeded_crack_faces),
2391 subtract(*crackFaces, old_crack_faces));
2392 if (nbCrackFaces >= reconstructed_crack_faces.size()) {
2393 MOFEM_LOG("EPSELF", Sev::warning)
2394 << "No new crack faces to add, skipping adding to meshset";
2395 extendeded_crack_faces = subtract(
2396 extendeded_crack_faces, subtract(*crackFaces, old_crack_faces));
2397 MOFEM_LOG("EPSELF", Sev::inform)
2398 << "Number crack faces size (extended) "
2399 << extendeded_crack_faces.size();
2400 CHKERR mField.get_moab().clear_meshset(&meshset, 1);
2401 CHKERR mField.get_moab().add_entities(meshset, extendeded_crack_faces);
2402 } else {
2403 CHKERR mField.get_moab().clear_meshset(&meshset, 1);
2404 CHKERR mField.get_moab().add_entities(meshset,
2405 reconstructed_crack_faces);
2406 MOFEM_LOG("EPSELF", Sev::inform)
2407 << "Number crack faces size (reconstructed) "
2408 << reconstructed_crack_faces.size();
2409 nbCrackFaces = reconstructed_crack_faces.size();
2410 }
2411 }
2412
2413 Range crack_faces;
2414 if (!mField.get_comm_rank()) {
2415 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTRI,
2416 crack_faces);
2417 }
2418 crack_faces = send_type(mField, crack_faces, MBTRI);
2419 if (mField.get_comm_rank()) {
2420 CHKERR mField.get_moab().clear_meshset(&meshset, 1);
2421 CHKERR mField.get_moab().add_entities(meshset, crack_faces);
2422 }
2423
2424 MoFEMFunctionReturn(0);
2425 };
2426
2427 CHKERR resolve_consisten_crack_extension();
2428
2429 MoFEMFunctionReturn(0);
2430};
2431
2432MoFEMErrorCode EshelbianCore::saveOrgCoords() {
2433 MoFEMFunctionBegin;
2434 auto crack_faces =
2435 get_range_from_block(mField, "CRACK_COMPUTED", SPACE_DIM - 1);
2436 Range conn;
2437 CHKERR mField.get_moab().get_connectivity(crack_faces, conn, true);
2438 Range verts;
2439 CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, verts);
2440 verts = subtract(verts, conn);
2441 std::vector<double> coords(3 * verts.size());
2442 CHKERR mField.get_moab().get_coords(verts, coords.data());
2443 double def_coords[] = {0., 0., 0.};
2444 Tag th_org_coords;
2445 CHKERR mField.get_moab().tag_get_handle(
2446 "ORG_COORDS", 3, MB_TYPE_DOUBLE, th_org_coords,
2447 MB_TAG_CREAT | MB_TAG_DENSE, def_coords);
2448 CHKERR mField.get_moab().tag_set_data(th_org_coords, verts, coords.data());
2449 MoFEMFunctionReturn(0);
2450}
2451
2452MoFEMErrorCode EshelbianCore::createCrackSurfaceMeshset() {
2453 MoFEMFunctionBegin;
2454 auto meshset_mng = mField.getInterface<MeshsetsManager>();
2455 while (meshset_mng->checkMeshset(addCrackMeshsetId, BLOCKSET))
2456 ++addCrackMeshsetId;
2457 MOFEM_LOG("EP", Sev::inform)
2458 << "Crack added surface meshset " << addCrackMeshsetId;
2459 CHKERR meshset_mng->addMeshset(BLOCKSET, addCrackMeshsetId, "CRACK_COMPUTED");
2460 MoFEMFunctionReturn(0);
2461};
2462
2463MoFEMErrorCode
2464EshelbianCore::calculateCrackArea(boost::shared_ptr<double> area_ptr) {
2465 MoFEMFunctionBegin;
2466
2467 if (!area_ptr) {
2468 CHK_THROW_MESSAGE(MOFEM_INVALID_DATA, "area_ptr is null");
2469 }
2470
2471 int success;
2472 *area_ptr = 0;
2473 if (mField.get_comm_rank() == 0) {
2474 MOFEM_LOG("EP", Sev::inform) << "Calculate crack area";
2475 auto crack_faces = get_range_from_block(mField, "CRACK", SPACE_DIM - 1);
2476 for (auto f : crack_faces) {
2477 *area_ptr += mField.getInterface<Tools>()->getTriArea(f);
2478 }
2479 success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0, mField.get_comm());
2480 } else {
2481 success = MPI_Bcast(area_ptr.get(), 1, MPI_DOUBLE, 0, mField.get_comm());
2482 }
2483 if (success != MPI_SUCCESS) {
2484 CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "bcast failed");
2485 }
2486 MoFEMFunctionReturn(0);
2487}
2488
2489} // namespace EshelbianPlasticity
send_type
static auto send_type(MoFEM::Interface &m_field, Range r, const EntityType type)
Definition EshelbianPlasticity.cpp:57
get_two_sides_of_crack_surface
static auto get_two_sides_of_crack_surface(MoFEM::Interface &m_field, Range crack_faces)
Definition EshelbianPlasticity.cpp:237
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:1882
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:2452
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:2464
EshelbianCore::calculateOrientation
MoFEMErrorCode calculateOrientation(const int tag, bool set_orientation)
Definition EshelbianFracture.cpp:630
EshelbianCore::setNewFrontCoordinates
MoFEMErrorCode setNewFrontCoordinates()
Definition EshelbianFracture.cpp:1660
EshelbianCore::addCrackSurfaces
MoFEMErrorCode addCrackSurfaces(const bool debug=false)
Definition EshelbianFracture.cpp:1699
EshelbianCore::saveOrgCoords
MoFEMErrorCode saveOrgCoords()
Definition EshelbianFracture.cpp:2432
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