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