Iterate over front edges, get adjacent faces, find maximal face energy. Maximal face energy is stored in the edge. Maximal face energy is magnitude of edge Griffith force.
For each front edge, find maximal face energy and orientation. This is by finding angle between edge material force and maximal face normal
958 {
960
961 constexpr bool debug =
false;
963 constexpr auto sev = Sev::verbose;
964
967 auto body_skin = get_skin(
mField, body_ents);
968 Range body_skin_edges;
970 moab::Interface::UNION);
971 Range boundary_skin_verts;
973 boundary_skin_verts, true);
974
976 Range geometry_edges_verts;
978 geometry_edges_verts, true);
979 Range crack_faces_verts;
981 true);
982 Range crack_faces_edges;
984 *
crackFaces, 1,
true, crack_faces_edges, moab::Interface::UNION);
985 Range crack_faces_tets;
987 *
crackFaces, 3,
true, crack_faces_tets, moab::Interface::UNION);
988
993 moab::Interface::UNION);
994 Range front_verts_edges;
996 front_verts, 1, true, front_verts_edges, moab::Interface::UNION);
997
998 auto get_tags_vec = [&](auto tag_name, int dim) {
999 std::vector<Tag> tags(1);
1000
1001 if (dim > 3)
1003
1004 auto create_and_clean = [&]() {
1007 auto rval = moab.tag_get_handle(tag_name, tags[0]);
1008 if (rval == MB_SUCCESS) {
1009 moab.tag_delete(tags[0]);
1010 }
1011 double def_val[] = {0., 0., 0.};
1012 CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
1013 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
1015 };
1016
1018
1019 return tags;
1020 };
1021
1022 auto get_adj_front = [&](bool subtract_crack) {
1025 adj_front, moab::Interface::UNION);
1026 if (subtract_crack)
1027 adj_front = subtract(adj_front, *
crackFaces);
1028 return adj_front;
1029 };
1030
1032
1033 auto th_front_position = get_tags_vec("FrontPosition", 3);
1034 auto th_max_face_energy = get_tags_vec("MaxFaceEnergy", 1);
1035
1037
1038 auto get_layers_for_sides = [&](auto &side) {
1039 std::vector<Range> layers;
1040 auto get = [&]() {
1042
1043 auto get_adj = [&](
auto &
r,
int dim) {
1046 moab::Interface::UNION);
1047 return adj;
1048 };
1049
1050 auto get_tets = [&](
auto r) {
return get_adj(r,
SPACE_DIM); };
1051
1054 true);
1055 Range front_faces = get_adj(front_nodes, 2);
1056 front_faces = subtract(front_faces, *
crackFaces);
1057 auto front_tets = get_tets(front_nodes);
1058 auto front_side = intersect(side, front_tets);
1059 layers.push_back(front_side);
1060 for (;;) {
1061 auto adj_faces = get_skin(
mField, layers.back());
1062 adj_faces = intersect(adj_faces, front_faces);
1063 auto adj_faces_tets = get_tets(adj_faces);
1064 adj_faces_tets = intersect(adj_faces_tets, front_tets);
1065 layers.push_back(unite(layers.back(), adj_faces_tets));
1066 if (layers.back().size() == layers[layers.size() - 2].size()) {
1067 break;
1068 }
1069 }
1071 };
1073 return layers;
1074 };
1075
1077 auto layers_top = get_layers_for_sides(sides_pair.first);
1078 auto layers_bottom = get_layers_for_sides(sides_pair.second);
1079
1080#ifndef NDEBUG
1082 auto get_crack_adj_tets = [&](
auto r) {
1083 Range crack_faces_conn;
1085 Range crack_faces_conn_tets;
1087 crack_faces_conn,
SPACE_DIM,
true, crack_faces_conn_tets,
1088 moab::Interface::UNION);
1089 return crack_faces_conn_tets;
1090 };
1093 "crack_tets_" +
1098 sides_pair.second);
1099 MOFEM_LOG(
"EP", sev) <<
"Nb. layers " << layers_top.size();
1101 for (auto &r : layers_top) {
1102 MOFEM_LOG(
"EP", sev) <<
"Layer " <<
l <<
" size " <<
r.size();
1105 "layers_top_" + boost::lexical_cast<std::string>(
l) +
".vtk", r);
1107 }
1108
1110 for (auto &r : layers_bottom) {
1111 MOFEM_LOG(
"EP", sev) <<
"Layer " <<
l <<
" size " <<
r.size();
1114 "layers_bottom_" + boost::lexical_cast<std::string>(
l) +
".vtk", r);
1116 }
1117 }
1118#endif
1119
1120 auto get_cross = [&](
auto t_dir,
auto f) {
1129 return t_cross;
1130 };
1131
1132 auto get_sense = [&](
auto f,
auto e) {
1133 int side, sense, offset;
1135 "get sense");
1136 return std::make_tuple(side, sense, offset);
1137 };
1138
1139 auto calculate_edge_direction = [&](auto e, auto normalize = true) {
1140 const EntityHandle *conn;
1141 int num_nodes;
1143 std::array<double, 6> coords;
1146 &coords[0], &coords[1], &coords[2]};
1148 &coords[3], &coords[4], &coords[5]};
1151 t_dir(
i) = t_p1(
i) - t_p0(
i);
1152 if (normalize)
1154 return t_dir;
1155 };
1156
1157 auto evaluate_face_energy_and_set_orientation = [&](auto front_edges,
1158 auto front_faces,
1159 auto &sides_pair,
1160 auto th_position) {
1162
1164 Tag th_material_force;
1169 th_face_energy);
1170
1171
1173 th_material_force);
1174
1175 break;
1176 default:
1177 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1178 "Unknown energy release selector");
1179 };
1180
1181
1182
1183
1184
1185
1186 auto find_maximal_face_energy = [&](auto front_edges, auto front_faces,
1187 auto &edge_face_max_energy_map) {
1189
1192 auto body_skin = get_skin(
mField, body_ents);
1193
1195
1196 for (auto e : front_edges) {
1197
1198 double griffith_force;
1200 &griffith_force);
1201
1204 faces = subtract(intersect(faces, front_faces), body_skin);
1205 std::vector<double> face_energy(faces.size());
1207 face_energy.data());
1208 auto max_energy_it =
1209 std::max_element(face_energy.begin(), face_energy.end());
1210 double max_energy =
1211 max_energy_it != face_energy.end() ? *max_energy_it : 0;
1212
1213 edge_face_max_energy_map[e] =
1214 std::make_tuple(faces[max_energy_it - face_energy.begin()],
1215 griffith_force, static_cast<double>(0));
1217 << "Edge " << e << " griffith force " << griffith_force
1218 << " max face energy " << max_energy << " factor "
1219 << max_energy / griffith_force;
1220
1221 max_faces.insert(faces[max_energy_it - face_energy.begin()]);
1222 }
1223
1224#ifndef NDEBUG
1228 "max_faces_" +
1230 ".vtk",
1231 max_faces);
1232 }
1233#endif
1234
1236 };
1237
1238
1239
1240
1241
1242
1243 auto calculate_face_orientation = [&](auto &edge_face_max_energy_map) {
1245
1246 auto up_down_face = [&](
1247
1248 auto &face_angle_map_up,
1249 auto &face_angle_map_down
1250
1251 ) {
1253
1254 for (
auto &
m : edge_face_max_energy_map) {
1256 auto [max_face, energy, opt_angle] =
m.second;
1257
1260 faces = intersect(faces, front_faces);
1263 false, adj_tets,
1264 moab::Interface::UNION);
1265 if (adj_tets.size()) {
1266
1269 false, adj_tets,
1270 moab::Interface::UNION);
1271 if (adj_tets.size()) {
1272
1273 Range adj_tets_faces;
1274
1276 adj_tets,
SPACE_DIM - 1,
false, adj_tets_faces,
1277 moab::Interface::UNION);
1278 adj_tets_faces = intersect(adj_tets_faces, faces);
1280
1281
1282 auto t_cross_max =
1283 get_cross(calculate_edge_direction(e, true), max_face);
1284 auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
1285 t_cross_max(
i) *= sense_max;
1286
1287 for (
auto t : adj_tets) {
1288 Range adj_tets_faces;
1290 &
t, 1,
SPACE_DIM - 1,
false, adj_tets_faces);
1291 adj_tets_faces = intersect(adj_tets_faces, faces);
1292 adj_tets_faces =
1293 subtract(adj_tets_faces,
Range(max_face, max_face));
1294
1295 if (adj_tets_faces.size() == 1) {
1296
1297
1298
1299 auto t_cross = get_cross(calculate_edge_direction(e, true),
1300 adj_tets_faces[0]);
1301 auto [side, sense, offset] =
1302 get_sense(adj_tets_faces[0], e);
1303 t_cross(
i) *= sense;
1304 double dot = t_cross(
i) * t_cross_max(
i);
1305 auto angle = std::acos(dot);
1306
1307 double face_energy;
1309 th_face_energy, adj_tets_faces, &face_energy);
1310
1311 auto [side_face, sense_face, offset_face] =
1312 get_sense(
t, max_face);
1313
1314 if (sense_face > 0) {
1315 face_angle_map_up[e] = std::make_tuple(face_energy, angle,
1316 adj_tets_faces[0]);
1317
1318 } else {
1319 face_angle_map_down[e] = std::make_tuple(
1320 face_energy, -angle, adj_tets_faces[0]);
1321 }
1322 }
1323 }
1324 }
1325 }
1326 }
1327
1329 };
1330
1331 auto calc_optimal_angle = [&](
1332
1333 auto &face_angle_map_up,
1334 auto &face_angle_map_down
1335
1336 ) {
1338
1339 for (
auto &
m : edge_face_max_energy_map) {
1341 auto &[max_face, e0,
a0] =
m.second;
1342
1343 if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
1344
1345 if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
1346 face_angle_map_down.find(e) == face_angle_map_down.end()) {
1347
1348 } else {
1349
1353
1354 Tag th_material_force;
1356 th_material_force);
1359 th_material_force, &e, 1, &t_material_force(0));
1360 auto material_force_magnitude = t_material_force.
l2();
1361 if (material_force_magnitude <
1362 std::numeric_limits<double>::epsilon()) {
1364
1365 } else {
1366
1367 auto t_edge_dir = calculate_edge_direction(e, true);
1368 auto t_cross_max = get_cross(t_edge_dir, max_face);
1369 auto [side, sense, offset] = get_sense(max_face, e);
1370 t_cross_max(sense) *= sense;
1371
1375
1377 t_cross_max.normalize();
1380 t_material_force(
J) * t_cross_max(K);
1381 a0 = -std::asin(t_cross(
I) * t_edge_dir(
I));
1382
1384 <<
"Optimal angle " <<
a0 <<
" energy " << e0;
1385 }
1386 break;
1387 }
1388 default: {
1389
1390 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1391 "Unknown energy release selector");
1392 }
1393 }
1394 }
1395 }
1396 }
1397
1399 };
1400
1401 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
1402 face_angle_map_up;
1403 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
1404 face_angle_map_down;
1405 CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
1406 CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
1407
1408#ifndef NDEBUG
1410 auto th_angle = get_tags_vec("Angle", 1);
1412 for (
auto &
m : face_angle_map_up) {
1413 auto [e,
a, face] =
m.second;
1414 up.insert(face);
1416 }
1418 for (
auto &
m : face_angle_map_down) {
1419 auto [e,
a, face] =
m.second;
1420 down.insert(face);
1422 }
1423
1424 Range max_energy_faces;
1425 for (
auto &
m : edge_face_max_energy_map) {
1426 auto [face, e, angle] =
m.second;
1427 max_energy_faces.insert(face);
1429 &angle);
1430 }
1435 max_energy_faces);
1436 }
1437 }
1438#endif
1439
1441 };
1442
1443 auto get_conn = [&](auto e) {
1446 "get conn");
1447 return conn;
1448 };
1449
1450 auto get_adj = [&](auto e, auto dim) {
1453 e, dim, false, adj, moab::Interface::UNION),
1454 "get adj");
1455 return adj;
1456 };
1457
1458 auto get_coords = [&](
auto v) {
1461 "get coords");
1462 return t_coords;
1463 };
1464
1465
1466 auto get_rotated_normal = [&](
auto e,
auto f,
auto angle) {
1469 auto t_edge_dir = calculate_edge_direction(e, true);
1470 auto [side, sense, offset] = get_sense(
f, e);
1471 t_edge_dir(
i) *= sense;
1472 t_edge_dir.normalize();
1473 t_edge_dir(
i) *= angle;
1478 t_rotated_normal(
i) = t_R(
i,
j) * t_normal(
j);
1479 return std::make_tuple(t_normal, t_rotated_normal);
1480 };
1481
1482 auto set_coord = [&](
auto v,
auto &adj_vertex_tets_verts,
auto &coords,
1483 auto &t_move, auto gamma) {
1484 auto index = adj_vertex_tets_verts.index(
v);
1485 if (index >= 0) {
1486 for (auto ii : {0, 1, 2}) {
1487 coords[3 * index + ii] += gamma * t_move(ii);
1488 }
1489 return true;
1490 }
1491 return false;
1492 };
1493
1494 auto tets_quality = [&](auto quality, auto &adj_vertex_tets_verts,
1495 auto &adj_vertex_tets, auto &coords) {
1496 for (
auto t : adj_vertex_tets) {
1497 const EntityHandle *conn;
1498 int num_nodes;
1500 std::array<double, 12> tet_coords;
1501 for (
auto n = 0;
n != 4; ++
n) {
1502 auto index = adj_vertex_tets_verts.index(conn[
n]);
1503 if (index < 0) {
1505 }
1506 for (auto ii = 0; ii != 3; ++ii) {
1507 tet_coords[3 *
n + ii] = coords[3 * index + ii];
1508 }
1509 }
1510 double q = Tools::volumeLengthQuality(tet_coords.data());
1511 if (!std::isnormal(
q))
1513 quality = std::min(quality,
q);
1514 };
1515
1516 return quality;
1517 };
1518
1519 auto calculate_free_face_node_displacement =
1520 [&](auto &edge_face_max_energy_map) {
1521
1522 auto get_vertex_edges = [&](auto vertex) {
1524
1525 auto impl = [&]() {
1528 vertex_edges);
1529 vertex_edges = subtract(vertex_edges, front_verts_edges);
1530
1531 if (boundary_skin_verts.size() &&
1532 boundary_skin_verts.find(vertex[0]) !=
1533 boundary_skin_verts.end()) {
1534 MOFEM_LOG(
"EP", sev) <<
"Boundary vertex";
1535 vertex_edges = intersect(vertex_edges, body_skin_edges);
1536 }
1537 if (geometry_edges_verts.size() &&
1538 geometry_edges_verts.find(vertex[0]) !=
1539 geometry_edges_verts.end()) {
1540 MOFEM_LOG(
"EP", sev) <<
"Geometry edge vertex";
1541 vertex_edges = intersect(vertex_edges, geometry_edges);
1542 }
1543 if (crack_faces_verts.size() &&
1544 crack_faces_verts.find(vertex[0]) !=
1545 crack_faces_verts.end()) {
1546 MOFEM_LOG(
"EP", sev) <<
"Crack face vertex";
1547 vertex_edges = intersect(vertex_edges, crack_faces_edges);
1548 }
1550 };
1551
1553
1554 return vertex_edges;
1555 };
1556
1557
1558
1559 using Bundle = std::vector<
1560
1561 std::tuple<EntityHandle, EntityHandle, EntityHandle,
1563
1564 >;
1565 std::map<EntityHandle, Bundle> edge_bundle_map;
1566
1567 for (
auto &
m : edge_face_max_energy_map) {
1568
1569 auto edge =
m.first;
1570 auto &[max_face, energy, opt_angle] =
m.second;
1571
1572
1573 auto [t_normal, t_rotated_normal] =
1574 get_rotated_normal(edge, max_face, opt_angle);
1575
1576 auto front_vertex = get_conn(
Range(
m.first,
m.first));
1577 auto adj_tets = get_adj(
Range(max_face, max_face), 3);
1578 auto adj_tets_faces = get_adj(adj_tets, 2);
1579 auto adj_front_faces = subtract(
1580 intersect(get_adj(
Range(edge, edge), 2), adj_tets_faces),
1582 if (adj_front_faces.size() > 3)
1584 "adj_front_faces.size()>3");
1585
1588 &t_material_force(0));
1589 std::vector<double> griffith_energy(adj_front_faces.size());
1591 th_face_energy, adj_front_faces, griffith_energy.data());
1592
1593 auto set_edge_bundle = [&](auto min_gamma) {
1594 for (auto rotated_f : adj_front_faces) {
1595
1596 double rotated_face_energy =
1597 griffith_energy[adj_front_faces.index(rotated_f)];
1598
1599 auto vertex = subtract(get_conn(
Range(rotated_f, rotated_f)),
1600 front_vertex);
1601 if (vertex.size() != 1) {
1603 "Wrong number of vertex to move");
1604 }
1605 auto front_vertex_edges_vertex = get_conn(
1606 intersect(get_adj(front_vertex, 1), crack_faces_edges));
1607 vertex = subtract(
1608 vertex, front_vertex_edges_vertex);
1609 if (vertex.empty()) {
1610 continue;
1611 }
1612
1613 auto face_cardinality = [&](
auto f,
auto &seen_front_edges) {
1614 auto whole_front =
1616 subtract(body_skin_edges, crack_faces_edges));
1619 for (;
c < 10; ++
c) {
1620 auto front_edges =
1621 subtract(get_adj(faces, 1), seen_front_edges);
1622 if (front_edges.size() == 0) {
1623 return 0;
1624 }
1625 auto front_connected_edges =
1626 intersect(front_edges, whole_front);
1627 if (front_connected_edges.size()) {
1628 seen_front_edges.merge(front_connected_edges);
1630 }
1631 faces.merge(get_adj(front_edges, 2));
1633 }
1635 };
1636
1638 double rotated_face_cardinality = face_cardinality(
1639 rotated_f,
1640 seen_edges);
1641
1642
1643
1644 rotated_face_cardinality = std::max(rotated_face_cardinality,
1645 1.);
1646
1647 auto t_vertex_coords = get_coords(vertex);
1648 auto vertex_edges = get_vertex_edges(vertex);
1649
1650 EntityHandle f0 = front_vertex[0];
1651 EntityHandle f1 = front_vertex[1];
1655
1657 for (auto e_used_to_move_detection : vertex_edges) {
1658 auto edge_conn = get_conn(
Range(e_used_to_move_detection,
1659 e_used_to_move_detection));
1660 edge_conn = subtract(edge_conn, vertex);
1661
1662
1663
1664
1665
1666
1667
1668
1670 t_v0(
i) = (t_v_e0(
i) + t_v_e1(
i)) / 2;
1674 (t_v0(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
1675 auto b =
1676 (t_v3(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
1678
1679 constexpr double eps =
1680 std::numeric_limits<double>::epsilon();
1681 if (std::isnormal(gamma) && gamma < 1.0 -
eps &&
1682 gamma > -0.1) {
1684 t_move(
i) = gamma * (t_v3(
i) - t_vertex_coords(
i));
1685
1686 auto check_rotated_face_directoon = [&]() {
1688 t_delta(
i) = t_vertex_coords(
i) + t_move(
i) - t_v0(
i);
1690 auto dot =
1691 (t_material_force(
i) / t_material_force.
l2()) *
1693 return -dot > 0 ? true : false;
1694 };
1695
1696 if (check_rotated_face_directoon()) {
1697
1699 << "Crack edge " << edge << " moved face "
1700 << rotated_f
1701 << " edge: " << e_used_to_move_detection
1702 << " face direction/energy " << rotated_face_energy
1703 << " face cardinality " << rotated_face_cardinality
1704 << " gamma: " << gamma;
1705
1706 auto &bundle = edge_bundle_map[edge];
1707 bundle.emplace_back(rotated_f, e_used_to_move_detection,
1708 vertex[0], t_move, 1,
1709 rotated_face_cardinality, gamma);
1710 }
1711 }
1712 }
1713 }
1714 };
1715
1716 set_edge_bundle(std::numeric_limits<double>::epsilon());
1717 if (edge_bundle_map[edge].empty()) {
1718 set_edge_bundle(-1.);
1719 }
1720 }
1721
1722 return edge_bundle_map;
1723 };
1724
1725 auto get_sort_by_energy = [&](auto &edge_face_max_energy_map) {
1726 std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
1727 sort_by_energy;
1728
1729 for (
auto &
m : edge_face_max_energy_map) {
1731 auto &[max_face, energy, opt_angle] =
m.second;
1732 auto abs_energy = std::abs(energy);
1733 sort_by_energy[abs_energy] = std::make_tuple(e, max_face, opt_angle);
1734 }
1735
1736 return sort_by_energy;
1737 };
1738
1739 auto set_tag = [&](auto &&adj_edges_map, auto &&sort_by_energy) {
1741
1742 Tag th_face_pressure;
1744 mField.
get_moab().tag_get_handle(
"FacePressure", th_face_pressure),
1745 "get tag");
1746 auto get_face_pressure = [&](auto face) {
1747 double pressure;
1749 1, &pressure),
1750 "get rag data");
1751 return pressure;
1752 };
1753
1755 << "Number of edges to check " << sort_by_energy.size();
1756
1757 enum face_energy { POSITIVE, NEGATIVE };
1758 constexpr bool skip_negative = true;
1759
1760 for (auto fe : {face_energy::POSITIVE, face_energy::NEGATIVE}) {
1761
1762 std::vector<double> energies;
1763 double max_pressure = -1;
1764
1765
1766 for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
1767 ++it) {
1768 auto energy = it->first;
1769 auto [max_edge, max_face, opt_angle] = it->second;
1770
1771 auto face_pressure = get_face_pressure(max_face);
1773 << "Faces to check: " << max_face << " energy " << energy
1774 << " face pressure " << face_pressure;
1775
1776 const bool pressure_check =
1778 if (energy > 0 && pressure_check) {
1779 energies.push_back(energy);
1780 }
1781 max_pressure = std::max(max_pressure, face_pressure);
1782 }
1783
1784 double average_energy = 0;
1785 if (!energies.empty()) {
1786 average_energy =
1787 std::accumulate(energies.begin(), energies.end(), 0.) /
1788 energies.size();
1789 }
1790
1792 << "Average energy Griffiths energy of crack front "
1793 << average_energy;
1794
1795 bool positive_pressure_face_found = false;
1796
1797
1798
1799 for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
1800 ++it) {
1801
1802 auto energy = it->first;
1803 auto [max_edge, max_face, opt_angle] = it->second;
1804
1805 auto face_pressure = get_face_pressure(max_face);
1806 if (skip_negative) {
1807 if (fe == face_energy::POSITIVE) {
1808 if (face_pressure <
1811 << "Skip negative face " << max_face << " with energy "
1812 << energy << " and pressure " << face_pressure;
1813 continue;
1814 }
1815 }
1816 }
1817
1818 if (fe == face_energy::POSITIVE)
1819 positive_pressure_face_found = true;
1820
1822 << "Check face " << max_face << " edge " << max_edge
1823 << " energy " << energy << " optimal angle " << opt_angle
1824 << " face pressure " << face_pressure;
1825
1826
1827 if (!average_energy) {
1829 << "Average energy is zero, setting max Griffiths energy to "
1830 "current energy "
1831 << energy;
1832 average_energy = energy;
1833 }
1835 auto jt = adj_edges_map.find(max_edge);
1836 if (jt == adj_edges_map.end()) {
1838 << "Edge " << max_edge << " not found in adj_edges_map";
1839 continue;
1840 }
1841 auto &bundle = jt->second;
1842
1843 auto find_max_in_bundle_impl = [&](auto edge, auto &bundle,
1844 auto gamma) {
1846
1847 EntityHandle vertex_max = 0;
1848 EntityHandle face_max = 0;
1849 EntityHandle move_edge_max = 0;
1850 double max_quality = -2;
1851 double max_quality_evaluated = -2;
1852 double min_cardinality = std::numeric_limits<double>::max();
1853
1855
1856 for (auto &b : bundle) {
1857 auto &[face, move_edge, vertex, t_move, quality, cardinality,
1858 edge_gamma] = b;
1859
1860 auto adj_vertex_tets = get_adj(
Range(vertex, vertex), 3);
1861 auto adj_vertex_tets_verts = get_conn(adj_vertex_tets);
1862 std::vector<double> coords(3 * adj_vertex_tets_verts.size());
1864 adj_vertex_tets_verts, coords.data()),
1865 "get coords");
1866
1867 set_coord(vertex, adj_vertex_tets_verts, coords, t_move, gamma);
1868 quality = tets_quality(quality, adj_vertex_tets_verts,
1869 adj_vertex_tets, coords);
1870
1871 auto eval_quality = [](
auto q,
auto c,
auto edge_gamma) {
1874 } else {
1875 return ((edge_gamma < 0) ? (
q / 2) :
q) / pow(
c, 2);
1876 }
1877 };
1878
1879 if (eval_quality(quality, cardinality, edge_gamma) >=
1880 max_quality_evaluated) {
1881 max_quality = quality;
1882 min_cardinality = cardinality;
1883 vertex_max = vertex;
1884 face_max = face;
1885 move_edge_max = move_edge;
1886 t_move_last(
i) = t_move(
i);
1887 max_quality_evaluated =
1888 eval_quality(max_quality, min_cardinality, edge_gamma);
1889 }
1890 }
1891
1892 return std::make_tuple(vertex_max, face_max, t_move_last,
1893 max_quality, min_cardinality);
1894 };
1895
1896 auto find_max_in_bundle = [&](auto edge, auto &bundle) {
1897 auto b_org_bundle = bundle;
1898 auto r = find_max_in_bundle_impl(edge, bundle, 1.);
1899 auto &[vertex_max, face_max, t_move_last, max_quality,
1901 if (max_quality < 0) {
1902 for (double gamma = 0.95; gamma >= 0.45; gamma -= 0.05) {
1903 bundle = b_org_bundle;
1904 r = find_max_in_bundle_impl(edge, bundle, gamma);
1905 auto &[vertex_max, face_max, t_move_last, max_quality,
1908 << "Back tracking: gamma " << gamma << " edge " << edge
1909 << " quality " << max_quality << " cardinality "
1910 << cardinality;
1911 if (max_quality > 0.01) {
1913 t_move_last(
I) *= gamma;
1915 }
1916 }
1919 }
1921 };
1922
1923
1924 auto set_tag_to_vertex_and_face = [&](
auto &&
r,
auto &quality) {
1926 auto &[
v,
f, t_move,
q, cardinality] =
r;
1927
1928 if ((
q > 0 && std::isnormal(
q)) && energy > 0) {
1929
1931 <<
"Set tag: vertex " <<
v <<
" face " <<
f <<
" "
1932 << max_edge << " move " << t_move << " energy " << energy
1933 <<
" quality " <<
q <<
" cardinality " << cardinality;
1935 &t_move(0));
1937 1, &energy);
1938 }
1939
1942 };
1943
1944 double quality = -2;
1945 CHKERR set_tag_to_vertex_and_face(
1946
1947 find_max_in_bundle(max_edge, bundle),
1948
1949 quality
1950
1951 );
1952
1953 if (quality > 0 && std::isnormal(quality) && energy > 0) {
1955 << "Crack face set with quality: " << quality;
1957 }
1958 }
1959
1960 if (fe == face_energy::POSITIVE && !positive_pressure_face_found) {
1964 << "POTENTIAL ARREST: No suitable face found with positive "
1965 "face pressure to propagate crack";
1966 } else {
1968 << "POTENTIAL ARREST: No suitable face found with positive "
1969 "face pressure to propagate crack; continuing because "
1970 "propagation under compression is enabled";
1971 }
1972 }
1973
1974 if (!skip_negative)
1975 break;
1976 }
1977
1979 };
1980
1981
1982 MOFEM_LOG(
"EP", sev) <<
"Calculate orientation";
1983 std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
1984 edge_face_max_energy_map;
1985 CHKERR find_maximal_face_energy(front_edges, front_faces,
1986 edge_face_max_energy_map);
1987 CHKERR calculate_face_orientation(edge_face_max_energy_map);
1988
1989 MOFEM_LOG(
"EP", sev) <<
"Calculate node positions";
1991
1992 calculate_free_face_node_displacement(edge_face_max_energy_map),
1993 get_sort_by_energy(edge_face_max_energy_map)
1994
1995 );
1996
1998 };
1999
2000 auto get_max_griffith_force = [&](
auto r) {
2002 std::vector<double> gc(
r.size());
2004 CHKERR moab.tag_get_handle(
"GriffithForce", th_gc);
2005 CHKERR moab.tag_get_data(th_gc, r, gc.data());
2006 double max_griffith_force = 0;
2007 for (
size_t i = 0;
i <
r.size(); ++
i) {
2008 max_griffith_force = std::max(max_griffith_force, std::abs(gc[
i]));
2009 }
2010 return max_griffith_force;
2011 };
2012
2014 if (std::abs(get_max_griffith_force(get_adj_front(true))) >
2015 std::numeric_limits<double>::epsilon()) {
2016 CHKERR evaluate_face_energy_and_set_orientation(
2017 *
frontEdges, get_adj_front(
true), sides_pair, th_front_position);
2018 } else {
2019 auto adj_front = get_adj_front(true);
2020 double zero[] = {0., 0., 0.};
2022 zero);
2023 }
2024 }
2025
2026
2029 SCATTER_FORWARD);
2031 SCATTER_FORWARD);
2036 SCATTER_FORWARD);
2040
2041 auto get_max_moved_faces = [&]() {
2042 Range max_moved_faces;
2043 auto adj_front = get_adj_front(false);
2044 std::vector<double> face_energy(adj_front.size());
2046 face_energy.data());
2047 for (
int i = 0;
i != adj_front.size(); ++
i) {
2048 if (face_energy[
i] > std::numeric_limits<double>::epsilon()) {
2049 max_moved_faces.insert(adj_front[
i]);
2050 }
2051 }
2052
2053 return boost::make_shared<Range>(max_moved_faces);
2054 };
2055
2056
2059
2060#ifndef NDEBUG
2064 "max_moved_faces_" +
2067 }
2068#endif
2069
2071}
static auto get_two_sides_of_crack_surface(MoFEM::Interface &m_field, Range crack_faces)
@ MOFEM_ATOM_TEST_INVALID
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
const double c
speed of light (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FTensor::Index< 'm', 3 > m
static double crackingAtol
Cracking absolute tolerance.
static double crackingRtol
Cracking relative tolerance.
double avgGriffithsEnergy
static bool potentialCrackArrest
static PetscBool propagateUnderCompression
static enum EnergyReleaseSelector energyReleaseSelector
static auto exp(A &&t_w_vee, B &&theta)
Public Member Functions inherited from MoFEM::UnknownInterface