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
971 {
973
974 constexpr bool debug =
false;
976 constexpr auto sev = Sev::verbose;
977
980 auto body_skin = get_skin(
mField, body_ents);
981 Range body_skin_edges;
983 moab::Interface::UNION);
984 Range boundary_skin_verts;
986 boundary_skin_verts, true);
987
989 Range geometry_edges_verts;
991 geometry_edges_verts, true);
992 Range crack_faces_verts;
994 true);
995 Range crack_faces_edges;
997 *
crackFaces, 1,
true, crack_faces_edges, moab::Interface::UNION);
998 Range crack_faces_tets;
1000 *
crackFaces, 3,
true, crack_faces_tets, moab::Interface::UNION);
1001
1006 moab::Interface::UNION);
1007 Range front_verts_edges;
1009 front_verts, 1, true, front_verts_edges, moab::Interface::UNION);
1010
1011 auto get_tags_vec = [&](auto tag_name, int dim) {
1012 std::vector<Tag> tags(1);
1013
1014 if (dim > 3)
1016
1017 auto create_and_clean = [&]() {
1020 auto rval = moab.tag_get_handle(tag_name, tags[0]);
1021 if (rval == MB_SUCCESS) {
1022 moab.tag_delete(tags[0]);
1023 }
1024 double def_val[] = {0., 0., 0.};
1025 CHKERR moab.tag_get_handle(tag_name, dim, MB_TYPE_DOUBLE, tags[0],
1026 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
1028 };
1029
1031
1032 return tags;
1033 };
1034
1035 auto get_adj_front = [&](bool subtract_crack) {
1038 adj_front, moab::Interface::UNION);
1039 if (subtract_crack)
1040 adj_front = subtract(adj_front, *
crackFaces);
1041 return adj_front;
1042 };
1043
1045
1046 auto th_front_position = get_tags_vec("FrontPosition", 3);
1047 auto th_max_face_energy = get_tags_vec("MaxFaceEnergy", 1);
1048
1050
1051 auto get_crack_adj_tets = [&](
auto r) {
1052 Range crack_faces_conn;
1054 Range crack_faces_conn_tets;
1056 true, crack_faces_conn_tets,
1057 moab::Interface::UNION);
1058 return crack_faces_conn_tets;
1059 };
1060
1061 auto get_layers_for_sides = [&](auto &side) {
1062 std::vector<Range> layers;
1063 auto get = [&]() {
1065
1066 auto get_adj = [&](
auto &
r,
int dim) {
1069 moab::Interface::UNION);
1070 return adj;
1071 };
1072
1073 auto get_tets = [&](
auto r) {
return get_adj(r,
SPACE_DIM); };
1074
1077 true);
1078 Range front_faces = get_adj(front_nodes, 2);
1079 front_faces = subtract(front_faces, *
crackFaces);
1080 auto front_tets = get_tets(front_nodes);
1081 auto front_side = intersect(side, front_tets);
1082 layers.push_back(front_side);
1083 for (;;) {
1084 auto adj_faces = get_skin(
mField, layers.back());
1085 adj_faces = intersect(adj_faces, front_faces);
1086 auto adj_faces_tets = get_tets(adj_faces);
1087 adj_faces_tets = intersect(adj_faces_tets, front_tets);
1088 layers.push_back(unite(layers.back(), adj_faces_tets));
1089 if (layers.back().size() == layers[layers.size() - 2].size()) {
1090 break;
1091 }
1092 }
1094 };
1096 return layers;
1097 };
1098
1100 auto layers_top = get_layers_for_sides(sides_pair.first);
1101 auto layers_bottom = get_layers_for_sides(sides_pair.second);
1102
1103#ifndef NDEBUG
1107 "crack_tets_" +
1112 sides_pair.second);
1113 MOFEM_LOG(
"EP", sev) <<
"Nb. layers " << layers_top.size();
1115 for (auto &r : layers_top) {
1116 MOFEM_LOG(
"EP", sev) <<
"Layer " <<
l <<
" size " <<
r.size();
1119 "layers_top_" + boost::lexical_cast<std::string>(
l) +
".vtk", r);
1121 }
1122
1124 for (auto &r : layers_bottom) {
1125 MOFEM_LOG(
"EP", sev) <<
"Layer " <<
l <<
" size " <<
r.size();
1128 "layers_bottom_" + boost::lexical_cast<std::string>(
l) +
".vtk", r);
1130 }
1131 }
1132#endif
1133
1134 auto get_cross = [&](
auto t_dir,
auto f) {
1143 return t_cross;
1144 };
1145
1146 auto get_sense = [&](
auto f,
auto e) {
1147 int side, sense, offset;
1149 "get sense");
1150 return std::make_tuple(side, sense, offset);
1151 };
1152
1153 auto calculate_edge_direction = [&](auto e, auto normalize = true) {
1154 const EntityHandle *conn;
1155 int num_nodes;
1157 std::array<double, 6> coords;
1160 &coords[0], &coords[1], &coords[2]};
1162 &coords[3], &coords[4], &coords[5]};
1165 t_dir(
i) = t_p1(
i) - t_p0(
i);
1166 if (normalize)
1168 return t_dir;
1169 };
1170
1171 auto evaluate_face_energy_and_set_orientation = [&](auto front_edges,
1172 auto front_faces,
1173 auto &sides_pair,
1174 auto th_position) {
1176
1178 Tag th_material_force;
1183 th_face_energy);
1184
1185
1187 th_material_force);
1188
1189 break;
1190 default:
1191 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1192 "Unknown energy release selector");
1193 };
1194
1195
1196
1197
1198
1199
1200 auto find_maximal_face_energy = [&](auto front_edges, auto front_faces,
1201 auto &edge_face_max_energy_map) {
1203
1206 auto body_skin = get_skin(
mField, body_ents);
1207
1209
1210 for (auto e : front_edges) {
1211
1212 double griffith_force;
1214 &griffith_force);
1215
1218 faces = subtract(intersect(faces, front_faces), body_skin);
1219 std::vector<double> face_energy(faces.size());
1221 face_energy.data());
1222 auto max_energy_it =
1223 std::max_element(face_energy.begin(), face_energy.end());
1224 double max_energy =
1225 max_energy_it != face_energy.end() ? *max_energy_it : 0;
1226
1227 edge_face_max_energy_map[e] =
1228 std::make_tuple(faces[max_energy_it - face_energy.begin()],
1229 griffith_force, static_cast<double>(0));
1231 << "Edge " << e << " griffith force " << griffith_force
1232 << " max face energy " << max_energy << " factor "
1233 << max_energy / griffith_force;
1234
1235 max_faces.insert(faces[max_energy_it - face_energy.begin()]);
1236 }
1237
1238#ifndef NDEBUG
1242 "max_faces_" +
1244 ".vtk",
1245 max_faces);
1246 }
1247#endif
1248
1250 };
1251
1252
1253
1254
1255
1256
1257 auto calculate_face_orientation = [&](auto &edge_face_max_energy_map) {
1259
1260 auto up_down_face = [&](
1261
1262 auto &face_angle_map_up,
1263 auto &face_angle_map_down
1264
1265 ) {
1267
1268 for (
auto &
m : edge_face_max_energy_map) {
1270 auto [max_face, energy, opt_angle] =
m.second;
1271
1274 faces = intersect(faces, front_faces);
1277 false, adj_tets,
1278 moab::Interface::UNION);
1279 if (adj_tets.size()) {
1280
1283 false, adj_tets,
1284 moab::Interface::UNION);
1285 if (adj_tets.size()) {
1286
1287 Range adj_tets_faces;
1288
1290 adj_tets,
SPACE_DIM - 1,
false, adj_tets_faces,
1291 moab::Interface::UNION);
1292 adj_tets_faces = intersect(adj_tets_faces, faces);
1294
1295
1296 auto t_cross_max =
1297 get_cross(calculate_edge_direction(e, true), max_face);
1298 auto [side_max, sense_max, offset_max] = get_sense(max_face, e);
1299 t_cross_max(
i) *= sense_max;
1300
1301 for (
auto t : adj_tets) {
1302 Range adj_tets_faces;
1304 &
t, 1,
SPACE_DIM - 1,
false, adj_tets_faces);
1305 adj_tets_faces = intersect(adj_tets_faces, faces);
1306 adj_tets_faces =
1307 subtract(adj_tets_faces,
Range(max_face, max_face));
1308
1309 if (adj_tets_faces.size() == 1) {
1310
1311
1312
1313 auto t_cross = get_cross(calculate_edge_direction(e, true),
1314 adj_tets_faces[0]);
1315 auto [side, sense, offset] =
1316 get_sense(adj_tets_faces[0], e);
1317 t_cross(
i) *= sense;
1318 double dot = t_cross(
i) * t_cross_max(
i);
1319 auto angle = std::acos(dot);
1320
1321 double face_energy;
1323 th_face_energy, adj_tets_faces, &face_energy);
1324
1325 auto [side_face, sense_face, offset_face] =
1326 get_sense(
t, max_face);
1327
1328 if (sense_face > 0) {
1329 face_angle_map_up[e] = std::make_tuple(face_energy, angle,
1330 adj_tets_faces[0]);
1331
1332 } else {
1333 face_angle_map_down[e] = std::make_tuple(
1334 face_energy, -angle, adj_tets_faces[0]);
1335 }
1336 }
1337 }
1338 }
1339 }
1340 }
1341
1343 };
1344
1345 auto calc_optimal_angle = [&](
1346
1347 auto &face_angle_map_up,
1348 auto &face_angle_map_down
1349
1350 ) {
1352
1353 for (
auto &
m : edge_face_max_energy_map) {
1355 auto &[max_face, e0,
a0] =
m.second;
1356
1357 if (std::abs(e0) > std::numeric_limits<double>::epsilon()) {
1358
1359 if (face_angle_map_up.find(e) == face_angle_map_up.end() ||
1360 face_angle_map_down.find(e) == face_angle_map_down.end()) {
1361
1362 } else {
1363
1367
1368 Tag th_material_force;
1370 th_material_force);
1373 th_material_force, &e, 1, &t_material_force(0));
1374 auto material_force_magnitude = t_material_force.
l2();
1375 if (material_force_magnitude <
1376 std::numeric_limits<double>::epsilon()) {
1378
1379 } else {
1380
1381 auto t_edge_dir = calculate_edge_direction(e, true);
1382 auto t_cross_max = get_cross(t_edge_dir, max_face);
1383 auto [side, sense, offset] = get_sense(max_face, e);
1384 t_cross_max(sense) *= sense;
1385
1389
1391 t_cross_max.normalize();
1394 t_material_force(
J) * t_cross_max(K);
1395 a0 = -std::asin(t_cross(
I) * t_edge_dir(
I));
1396
1398 <<
"Optimal angle " <<
a0 <<
" energy " << e0;
1399 }
1400 break;
1401 }
1402 default: {
1403
1404 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
1405 "Unknown energy release selector");
1406 }
1407 }
1408 }
1409 }
1410 }
1411
1413 };
1414
1415 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
1416 face_angle_map_up;
1417 std::map<EntityHandle, std::tuple<double, double, EntityHandle>>
1418 face_angle_map_down;
1419 CHKERR up_down_face(face_angle_map_up, face_angle_map_down);
1420 CHKERR calc_optimal_angle(face_angle_map_up, face_angle_map_down);
1421
1422#ifndef NDEBUG
1424 auto th_angle = get_tags_vec("Angle", 1);
1426 for (
auto &
m : face_angle_map_up) {
1427 auto [e,
a, face] =
m.second;
1428 up.insert(face);
1430 }
1432 for (
auto &
m : face_angle_map_down) {
1433 auto [e,
a, face] =
m.second;
1434 down.insert(face);
1436 }
1437
1438 Range max_energy_faces;
1439 for (
auto &
m : edge_face_max_energy_map) {
1440 auto [face, e, angle] =
m.second;
1441 max_energy_faces.insert(face);
1443 &angle);
1444 }
1449 max_energy_faces);
1450 }
1451 }
1452#endif
1453
1455 };
1456
1457 auto get_conn = [&](auto e) {
1460 "get conn");
1461 return conn;
1462 };
1463
1464 auto get_adj = [&](auto e, auto dim) {
1467 e, dim, false, adj, moab::Interface::UNION),
1468 "get adj");
1469 return adj;
1470 };
1471
1472 auto get_coords = [&](
auto v) {
1475 "get coords");
1476 return t_coords;
1477 };
1478
1479
1480 auto get_rotated_normal = [&](
auto e,
auto f,
auto angle) {
1483 auto t_edge_dir = calculate_edge_direction(e, true);
1484 auto [side, sense, offset] = get_sense(
f, e);
1485 t_edge_dir(
i) *= sense;
1486 t_edge_dir.normalize();
1487 t_edge_dir(
i) *= angle;
1492 t_rotated_normal(
i) = t_R(
i,
j) * t_normal(
j);
1493 return std::make_tuple(t_normal, t_rotated_normal);
1494 };
1495
1496 auto set_coord = [&](
auto v,
auto &adj_vertex_tets_verts,
auto &coords,
1497 auto &t_move, auto gamma) {
1498 auto index = adj_vertex_tets_verts.index(
v);
1499 if (index >= 0) {
1500 for (auto ii : {0, 1, 2}) {
1501 coords[3 * index + ii] += gamma * t_move(ii);
1502 }
1503 return true;
1504 }
1505 return false;
1506 };
1507
1508 auto tets_quality = [&](auto quality, auto &adj_vertex_tets_verts,
1509 auto &adj_vertex_tets, auto &coords) {
1510 for (
auto t : adj_vertex_tets) {
1511 const EntityHandle *conn;
1512 int num_nodes;
1514 std::array<double, 12> tet_coords;
1515 for (
auto n = 0;
n != 4; ++
n) {
1516 auto index = adj_vertex_tets_verts.index(conn[
n]);
1517 if (index < 0) {
1519 }
1520 for (auto ii = 0; ii != 3; ++ii) {
1521 tet_coords[3 *
n + ii] = coords[3 * index + ii];
1522 }
1523 }
1524 double q = Tools::volumeLengthQuality(tet_coords.data());
1525 if (!std::isnormal(q))
1526 q = -2;
1527 quality = std::min(quality, q);
1528 };
1529
1530 return quality;
1531 };
1532
1533 auto calculate_free_face_node_displacement =
1534 [&](auto &edge_face_max_energy_map) {
1535
1536 auto get_vertex_edges = [&](auto vertex) {
1538
1539 auto impl = [&]() {
1542 vertex_edges);
1543 vertex_edges = subtract(vertex_edges, front_verts_edges);
1544
1545 if (boundary_skin_verts.size() &&
1546 boundary_skin_verts.find(vertex[0]) !=
1547 boundary_skin_verts.end()) {
1548 MOFEM_LOG(
"EP", sev) <<
"Boundary vertex";
1549 vertex_edges = intersect(vertex_edges, body_skin_edges);
1550 }
1551 if (geometry_edges_verts.size() &&
1552 geometry_edges_verts.find(vertex[0]) !=
1553 geometry_edges_verts.end()) {
1554 MOFEM_LOG(
"EP", sev) <<
"Geometry edge vertex";
1555 vertex_edges = intersect(vertex_edges, geometry_edges);
1556 }
1557 if (crack_faces_verts.size() &&
1558 crack_faces_verts.find(vertex[0]) !=
1559 crack_faces_verts.end()) {
1560 MOFEM_LOG(
"EP", sev) <<
"Crack face vertex";
1561 vertex_edges = intersect(vertex_edges, crack_faces_edges);
1562 }
1564 };
1565
1567
1568 return vertex_edges;
1569 };
1570
1571
1572
1573 using Bundle = std::vector<
1574
1575 std::tuple<EntityHandle, EntityHandle, EntityHandle,
1577
1578 >;
1579 std::map<EntityHandle, Bundle> edge_bundle_map;
1580
1581 for (
auto &
m : edge_face_max_energy_map) {
1582
1583 auto edge =
m.first;
1584 auto &[max_face, energy, opt_angle] =
m.second;
1585
1586
1587 auto [t_normal, t_rotated_normal] =
1588 get_rotated_normal(edge, max_face, opt_angle);
1589
1590 auto front_vertex = get_conn(
Range(
m.first,
m.first));
1591 auto adj_tets = get_adj(
Range(max_face, max_face), 3);
1592 auto adj_tets_faces = get_adj(adj_tets, 2);
1593 auto adj_front_faces = subtract(
1594 intersect(get_adj(
Range(edge, edge), 2), adj_tets_faces),
1596 if (adj_front_faces.size() > 3)
1598 "adj_front_faces.size()>3");
1599
1602 &t_material_force(0));
1603 std::vector<double> griffith_energy(adj_front_faces.size());
1605 th_face_energy, adj_front_faces, griffith_energy.data());
1606
1607 auto set_edge_bundle = [&](auto min_gamma) {
1608 for (auto rotated_f : adj_front_faces) {
1609
1610 double rotated_face_energy =
1611 griffith_energy[adj_front_faces.index(rotated_f)];
1612
1613 auto vertex = subtract(get_conn(
Range(rotated_f, rotated_f)),
1614 front_vertex);
1615 if (vertex.size() != 1) {
1617 "Wrong number of vertex to move");
1618 }
1619 auto front_vertex_edges_vertex = get_conn(
1620 intersect(get_adj(front_vertex, 1), crack_faces_edges));
1621 vertex = subtract(
1622 vertex, front_vertex_edges_vertex);
1623 if (vertex.empty()) {
1624 continue;
1625 }
1626
1627 auto face_cardinality = [&](
auto f,
auto &seen_front_edges) {
1628 auto whole_front =
1630 subtract(body_skin_edges, crack_faces_edges));
1633 for (;
c < 10; ++
c) {
1634 auto front_edges =
1635 subtract(get_adj(faces, 1), seen_front_edges);
1636 if (front_edges.size() == 0) {
1637 return 0;
1638 }
1639 auto front_connected_edges =
1640 intersect(front_edges, whole_front);
1641 if (front_connected_edges.size()) {
1642 seen_front_edges.merge(front_connected_edges);
1644 }
1645 faces.merge(get_adj(front_edges, 2));
1647 }
1649 };
1650
1652 double rotated_face_cardinality = face_cardinality(
1653 rotated_f,
1654 seen_edges);
1655
1656
1657
1658 rotated_face_cardinality = std::max(rotated_face_cardinality,
1659 1.);
1660
1661 auto t_vertex_coords = get_coords(vertex);
1662 auto vertex_edges = get_vertex_edges(vertex);
1663
1664 EntityHandle f0 = front_vertex[0];
1665 EntityHandle f1 = front_vertex[1];
1669
1671 for (auto e_used_to_move_detection : vertex_edges) {
1672 auto edge_conn = get_conn(
Range(e_used_to_move_detection,
1673 e_used_to_move_detection));
1674 edge_conn = subtract(edge_conn, vertex);
1675
1676
1677
1678
1679
1680
1681
1682
1684 t_v0(
i) = (t_v_e0(
i) + t_v_e1(
i)) / 2;
1688 (t_v0(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
1689 auto b =
1690 (t_v3(
i) - t_vertex_coords(
i)) * t_rotated_normal(
i);
1692
1693 constexpr double eps =
1694 std::numeric_limits<double>::epsilon();
1695 if (std::isnormal(gamma) && gamma < 1.0 -
eps &&
1696 gamma > -0.1) {
1698 t_move(
i) = gamma * (t_v3(
i) - t_vertex_coords(
i));
1699
1700 auto check_rotated_face_directoon = [&]() {
1702 t_delta(
i) = t_vertex_coords(
i) + t_move(
i) - t_v0(
i);
1704 auto dot =
1705 (t_material_force(
i) / t_material_force.
l2()) *
1707 return -dot > 0 ? true : false;
1708 };
1709
1710 if (check_rotated_face_directoon()) {
1711
1713 << "Crack edge " << edge << " moved face "
1714 << rotated_f
1715 << " edge: " << e_used_to_move_detection
1716 << " face direction/energy " << rotated_face_energy
1717 << " face cardinality " << rotated_face_cardinality
1718 << " gamma: " << gamma;
1719
1720 auto &bundle = edge_bundle_map[edge];
1721 bundle.emplace_back(rotated_f, e_used_to_move_detection,
1722 vertex[0], t_move, 1,
1723 rotated_face_cardinality, gamma);
1724 }
1725 }
1726 }
1727 }
1728 };
1729
1730 set_edge_bundle(std::numeric_limits<double>::epsilon());
1731 if (edge_bundle_map[edge].empty()) {
1732 set_edge_bundle(-1.);
1733 }
1734 }
1735
1736 return edge_bundle_map;
1737 };
1738
1739 auto get_sort_by_energy = [&](auto &edge_face_max_energy_map) {
1740 std::map<double, std::tuple<EntityHandle, EntityHandle, double>>
1741 sort_by_energy;
1742
1743 for (
auto &
m : edge_face_max_energy_map) {
1745 auto &[max_face, energy, opt_angle] =
m.second;
1746 auto abs_energy = std::abs(energy);
1747 sort_by_energy[abs_energy] = std::make_tuple(e, max_face, opt_angle);
1748 }
1749
1750 return sort_by_energy;
1751 };
1752
1753 auto set_tag = [&](auto &&adj_edges_map, auto &&sort_by_energy) {
1755
1756 Tag th_face_pressure;
1758 mField.
get_moab().tag_get_handle(
"FacePressure", th_face_pressure),
1759 "get tag");
1760 auto get_face_pressure = [&](auto face) {
1761 double pressure;
1763 1, &pressure),
1764 "get rag data");
1765 return pressure;
1766 };
1767
1769 << "Number of edges to check " << sort_by_energy.size();
1770
1771 enum face_energy { POSITIVE, NEGATIVE };
1772 constexpr bool skip_negative = true;
1773
1774 for (auto fe : {face_energy::POSITIVE, face_energy::NEGATIVE}) {
1775
1776 std::vector<double> energies;
1777 double max_pressure = -1;
1778
1779
1780 for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
1781 ++it) {
1782 auto energy = it->first;
1783 auto [max_edge, max_face, opt_angle] = it->second;
1784
1785 auto face_pressure = get_face_pressure(max_face);
1786 energies.push_back(energy);
1787
1788 max_pressure = std::max(max_pressure, face_pressure);
1789 }
1790
1791 double average_energy = 0;
1792 if (!energies.empty()) {
1793 average_energy =
1794 std::accumulate(energies.begin(), energies.end(), 0.) /
1795 energies.size();
1796 }
1797
1799 << "Average energy Griffiths energy of crack front "
1800 << average_energy;
1801
1802
1803
1804 for (auto it = sort_by_energy.rbegin(); it != sort_by_energy.rend();
1805 ++it) {
1806
1807 auto energy = it->first;
1808 auto [max_edge, max_face, opt_angle] = it->second;
1809
1810 auto face_pressure = get_face_pressure(max_face);
1811 if (skip_negative) {
1812 if (fe == face_energy::POSITIVE) {
1813 if (face_pressure <
1816 << "Skip negative face " << max_face << " with energy "
1817 << energy << " and pressure " << face_pressure;
1818 continue;
1819 }
1820 }
1821 }
1822
1824 << "Check face " << max_face << " edge " << max_edge
1825 << " energy " << energy << " optimal angle " << opt_angle
1826 << " face pressure " << face_pressure;
1827
1828
1829 if (!average_energy) {
1831 << "Average energy is zero, setting max Griffiths energy to "
1832 "current energy "
1833 << energy;
1834 average_energy = energy;
1835 }
1837 auto jt = adj_edges_map.find(max_edge);
1838 if (jt == adj_edges_map.end()) {
1840 << "Edge " << max_edge << " not found in adj_edges_map";
1841 continue;
1842 }
1843 auto &bundle = jt->second;
1844
1845 auto find_max_in_bundle_impl = [&](auto edge, auto &bundle,
1846 auto gamma) {
1848
1849 EntityHandle vertex_max = 0;
1850 EntityHandle face_max = 0;
1851 EntityHandle move_edge_max = 0;
1852 double max_quality = -2;
1853 double max_quality_evaluated = -2;
1854 double min_cardinality = std::numeric_limits<double>::max();
1855
1857
1858 for (auto &b : bundle) {
1859 auto &[face, move_edge, vertex, t_move, quality, cardinality,
1860 edge_gamma] = b;
1861
1862 auto adj_vertex_tets = get_adj(
Range(vertex, vertex), 3);
1863 auto adj_vertex_tets_verts = get_conn(adj_vertex_tets);
1864 std::vector<double> coords(3 * adj_vertex_tets_verts.size());
1866 adj_vertex_tets_verts, coords.data()),
1867 "get coords");
1868
1869 set_coord(vertex, adj_vertex_tets_verts, coords, t_move, gamma);
1870 quality = tets_quality(quality, adj_vertex_tets_verts,
1871 adj_vertex_tets, coords);
1872
1873 auto eval_quality = [](
auto q,
auto c,
auto edge_gamma) {
1874 if (q < 0) {
1875 return q;
1876 } else {
1877 return ((edge_gamma < 0) ? (q / 2) : q) / pow(
c, 2);
1878 }
1879 };
1880
1881 if (eval_quality(quality, cardinality, edge_gamma) >=
1882 max_quality_evaluated) {
1883 max_quality = quality;
1884 min_cardinality = cardinality;
1885 vertex_max = vertex;
1886 face_max = face;
1887 move_edge_max = move_edge;
1888 t_move_last(
i) = t_move(
i);
1889 max_quality_evaluated =
1890 eval_quality(max_quality, min_cardinality, edge_gamma);
1891 }
1892 }
1893
1894 return std::make_tuple(vertex_max, face_max, t_move_last,
1895 max_quality, min_cardinality);
1896 };
1897
1898 auto find_max_in_bundle = [&](auto edge, auto &bundle) {
1899 auto b_org_bundle = bundle;
1900 auto r = find_max_in_bundle_impl(edge, bundle, 1.);
1901 auto &[vertex_max, face_max, t_move_last, max_quality,
1903 if (max_quality < 0) {
1904 for (double gamma = 0.95; gamma >= 0.45; gamma -= 0.05) {
1905 bundle = b_org_bundle;
1906 r = find_max_in_bundle_impl(edge, bundle, gamma);
1907 auto &[vertex_max, face_max, t_move_last, max_quality,
1910 << "Back tracking: gamma " << gamma << " edge " << edge
1911 << " quality " << max_quality << " cardinality "
1912 << cardinality;
1913 if (max_quality > 0.01) {
1915 t_move_last(
I) *= gamma;
1917 }
1918 }
1921 }
1923 };
1924
1925
1926 auto set_tag_to_vertex_and_face = [&](
auto &&
r,
auto &quality) {
1928 auto &[
v,
f, t_move, q, cardinality] =
r;
1929
1930 if ((q > 0 && std::isnormal(q)) && energy > 0) {
1931
1933 <<
"Set tag: vertex " <<
v <<
" face " <<
f <<
" "
1934 << max_edge << " move " << t_move << " energy " << energy
1935 << " quality " << q << " cardinality " << cardinality;
1937 &t_move(0));
1939 1, &energy);
1940 }
1941
1942 quality = q;
1944 };
1945
1946 double quality = -2;
1947 CHKERR set_tag_to_vertex_and_face(
1948
1949 find_max_in_bundle(max_edge, bundle),
1950
1951 quality
1952
1953 );
1954
1955 if (quality > 0 && std::isnormal(quality) && energy > 0) {
1957 << "Crack face set with quality: " << quality;
1959 }
1960 }
1961
1962 if (!skip_negative)
1963 break;
1964 }
1965
1967 };
1968
1969
1970 MOFEM_LOG(
"EP", sev) <<
"Calculate orientation";
1971 std::map<EntityHandle, std::tuple<EntityHandle, double, double>>
1972 edge_face_max_energy_map;
1973 CHKERR find_maximal_face_energy(front_edges, front_faces,
1974 edge_face_max_energy_map);
1975 CHKERR calculate_face_orientation(edge_face_max_energy_map);
1976
1977 MOFEM_LOG(
"EP", sev) <<
"Calculate node positions";
1979
1980 calculate_free_face_node_displacement(edge_face_max_energy_map),
1981 get_sort_by_energy(edge_face_max_energy_map)
1982
1983 );
1984
1986 };
1987
1988 auto get_max_griffith_force = [&](
auto r) {
1990 std::vector<double> gc(
r.size());
1992 CHKERR moab.tag_get_handle(
"GriffithForce", th_gc);
1993 CHKERR moab.tag_get_data(th_gc, r, gc.data());
1994 double max_griffith_force = 0;
1995 for (
size_t i = 0;
i <
r.size(); ++
i) {
1996 max_griffith_force = std::max(max_griffith_force, std::abs(gc[
i]));
1997 }
1998 return max_griffith_force;
1999 };
2000
2002 if (std::abs(get_max_griffith_force(get_adj_front(true))) >
2003 std::numeric_limits<double>::epsilon()) {
2004 CHKERR evaluate_face_energy_and_set_orientation(
2005 *
frontEdges, get_adj_front(
true), sides_pair, th_front_position);
2006 } else {
2007 auto adj_front = get_adj_front(true);
2008 double zero[] = {0., 0., 0.};
2010 zero);
2011 }
2012 }
2013
2014
2017 SCATTER_FORWARD);
2019 SCATTER_FORWARD);
2024 SCATTER_FORWARD);
2028
2029 auto get_max_moved_faces = [&]() {
2030 Range max_moved_faces;
2031 auto adj_front = get_adj_front(false);
2032 std::vector<double> face_energy(adj_front.size());
2034 face_energy.data());
2035 for (
int i = 0;
i != adj_front.size(); ++
i) {
2036 if (face_energy[
i] > std::numeric_limits<double>::epsilon()) {
2037 max_moved_faces.insert(adj_front[
i]);
2038 }
2039 }
2040
2041 return boost::make_shared<Range>(max_moved_faces);
2042 };
2043
2044
2047
2048#ifndef NDEBUG
2052 "max_moved_faces_" +
2055 }
2056#endif
2057
2059}
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 enum EnergyReleaseSelector energyReleaseSelector
static auto exp(A &&t_w_vee, B &&theta)
Public Member Functions inherited from MoFEM::UnknownInterface