359 Tag th_griffith_force;
360 Tag th_griffith_force_projected;
362 Tag th_front_positions;
364 Tag th_material_force;
366 boost::shared_ptr<OrientedBoxTreeTool> tree_body_skin_ptr;
367 tree_body_skin_ptr = boost::shared_ptr<OrientedBoxTreeTool>(
372 std::map<EntityHandle, EntityHandle> verts_map;
377 struct SetFrontPositionsParameters {
379 double tolRelSnapToFixedEdges;
381 double tolAbsSnapToFixedEdges;
386 double endNodeEdgeDeltaFactor;
388 SetFrontPositionsParameters() {
389 tolRelSnapToFixedEdges = 0.1;
390 tolAbsSnapToFixedEdges = 0.;
393 endNodeEdgeDeltaFactor = 0.2;
395 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
401 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"cutting_",
402 "Set front positions",
"none");
404 CHKERR PetscOptionsScalar(
405 "-surf_corner_factor",
406 "extension of the edge beyond corner when node on corer edge",
"",
407 cornerFactor, &cornerFactor, PETSC_NULL);
408 CHKERR PetscOptionsScalar(
410 "extension of cutting surface when node on the skin",
"", skinFactor,
411 &skinFactor, PETSC_NULL);
412 CHKERR PetscOptionsScalar(
413 "-end_node_edge_delta_factor",
414 "extension of the front edge at the node on the surface",
"",
415 endNodeEdgeDeltaFactor, &endNodeEdgeDeltaFactor, PETSC_NULL);
416 CHKERR PetscOptionsScalar(
"-snap_to_fixed_edge_rtol",
417 "snap to fixed edges (relative tolerances)",
"",
418 tolRelSnapToFixedEdges, &tolRelSnapToFixedEdges,
420 CHKERR PetscOptionsScalar(
"-snap_to_fixed_edge_atol",
421 "snap to fixed edges (absolute tolerances)",
"",
422 tolAbsSnapToFixedEdges, &tolAbsSnapToFixedEdges,
426 "### Input parameter: -cutting_surf_corner_factor %6.4e",
429 "### Input parameter: -cutting_surf_skin_factor %6.4e",
432 "CPMeshCutWorld", Sev::inform,
433 "### Input parameter: -cutting_end_node_edge_delta_factor %6.4e",
434 endNodeEdgeDeltaFactor);
436 "### Input parameter: -cutting_snap_to_fixed_edge_rtol %6.4e",
437 tolRelSnapToFixedEdges);
439 "### Input parameter: -cutting_snap_to_fixed_edge_atol %6.4e",
440 tolAbsSnapToFixedEdges);
442 ierr = PetscOptionsEnd();
449 auto save_entities = [&](
const std::string name,
Range ents) {
462 auto get_tags = [&]() {
467 th_griffith_force_projected);
474 if (
rval == MB_SUCCESS)
477 double def_val[] = {0, 0, 0};
479 "FRONT_POSITIONS", 3, MB_TYPE_DOUBLE, th_front_positions,
480 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
483 if (
rval == MB_SUCCESS)
487 "AREA_CHANGE", 3, MB_TYPE_DOUBLE, th_area_change,
488 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
492 auto get_surface = [
this, &
bit]() {
499 return intersect(
surface, bit_faces);
502 auto get_crack_front = [&]() {
509 return intersect(crack_front, bit_edges);
513 auto get_cracked_body_blockset = [&]() {
514 Range cracked_body_tets;
517 return cracked_body_tets;
521 auto build_body_skin_tree =
522 [&](boost::shared_ptr<OrientedBoxTreeTool> &tree_body_skin_ptr,
530 last_bit_tets = intersect(last_bit_tets, get_cracked_body_blockset());
531 Range last_bit_tets_skin;
532 CHKERR skin.find_skin(0, last_bit_tets,
false, last_bit_tets_skin);
533 CHKERR tree_body_skin_ptr->build(last_bit_tets_skin,
534 rootset_tree_body_skin);
539 auto set_positions = [&](
const Range nodes, Tag *
th =
nullptr) {
551 auto get_nodes = [
this](
Range edges) {
557 auto get_surface_skin = [
this, &skin](
auto &
surface) {
563 auto get_crack_front_ends = [
this, &skin](
Range crack_front) {
565 CHKERR skin.find_skin(0, crack_front,
false, end_nodes);
576 if (conn[1] == node) {
587 auto get_node_vec = [](std::vector<double> &
v,
const int nn) {
591 auto check_if_node_is_in_range = [](
const EntityHandle node,
593 return (ents.find(node) != ents.end());
596 auto check_for_nan = [&](
auto &
v, std::string msg =
"") {
599 if (std::isnan(
d) || std::isinf(
d))
602 "(%s) Not number %3.4e", msg.c_str(),
d);
605 "Not number %3.4e",
d);
611 auto set_tag = [&](
const EntityHandle node,
auto coords) {
613 CHKERR check_for_nan(coords,
"set_tag");
619 auto get_str = [](
auto v) {
return boost::lexical_cast<std::string>(
v); };
621 auto get_tag_data = [
this](Tag
th,
const Range &ents) {
622 std::vector<double>
v(3 * ents.size());
627 auto get_coords_vec = [
this](
const Range &ents) {
628 std::vector<double>
v(3 * ents.size());
633 auto get_prisms = [&]() {
636 moab::Interface::UNION);
637 prisms = prisms.subset_by_type(MBPRISM);
641 auto get_prims_surface = [&](
Range prisms) {
642 Range prisms_surface;
643 for (
auto prism : prisms) {
646 prisms_surface.insert(face);
648 return prisms_surface;
651 auto get_random_double = []() {
652 return static_cast<double>(rand()) /
static_cast<double>(RAND_MAX);
656 auto is_point_is_outside = [&](
const auto &coords) {
659 get_random_double() + std::numeric_limits<double>::epsilon();
660 random_ray[1] = get_random_double();
661 random_ray[2] = get_random_double();
662 random_ray /= norm_2(random_ray);
663 std::vector<double> distances_out;
664 std::vector<EntityHandle> facets_out;
665 CHKERR tree_body_skin_ptr->ray_intersect_triangles(
666 distances_out, facets_out, rootset_tree_body_skin, 1e-12, &coords[0],
668 return ((facets_out.size() % 2) == 0);
672 auto is_point_is_outside_with_check = [&](
const auto &coords) {
673 bool is_out = is_point_is_outside(coords);
677 is_out = is_point_is_outside(coords);
678 }
while (test_out != is_out);
682 auto get_proj_front_disp = [&](
auto front_disp,
683 auto griffith_force_projected) {
684 const double mgn_force_projected = norm_2(griffith_force_projected);
687 if (mgn_force_projected > std::numeric_limits<double>::epsilon()) {
688 griffith_force_projected /= mgn_force_projected;
689 noalias(disp) = factor * griffith_force_projected *
690 fmax(0, inner_prod(front_disp, griffith_force_projected));
699 auto get_ave_a_change = [&](
Range &crack_front_nodes,
700 std::vector<double> &area_change) {
701 double ave_a_change = 0;
703 for (
auto node : crack_front_nodes) {
704 auto a_change = get_node_vec(area_change, nn);
705 ave_a_change += norm_2(a_change);
707 ave_a_change /= crack_front_nodes.size();
711 auto get_ray_pout = [&](
const auto coords,
const auto ray) {
712 std::vector<double> distances_out;
713 std::vector<EntityHandle> facets_out;
714 double ray_length = norm_2(ray);
717 CHKERR tree_body_skin_ptr->ray_intersect_triangles(
718 distances_out, facets_out, rootset_tree_body_skin, 1e-12, &coords[0],
719 &unit_ray[0], &ray_length);
720 if (distances_out.empty())
723 return *std::min_element(distances_out.begin(), distances_out.end());
729 SetFrontPositionsParameters front_params;
732 sUrface = get_prims_surface(get_prisms());
735 Range crack_front = intersect(get_crack_front(), surface_skin);
736 Range crack_front_nodes = get_nodes(crack_front);
737 Range end_nodes = get_crack_front_ends(crack_front);
738 Range surface_skin_nodes =
739 subtract(subtract(get_nodes(surface_skin), crack_front_nodes),
740 get_nodes(intersect(
fixedEdges, surface_skin)));
741 surface_skin_nodes.merge(intersect(
cornerNodes, get_nodes(surface_skin)));
743 CHKERR GriffithForceElement::SaveGriffithForceOnTags(
745 CHKERR set_positions(get_nodes(get_prisms()),
nullptr);
746 CHKERR build_body_skin_tree(tree_body_skin_ptr, rootset_tree_body_skin);
749 CHKERR save_entities(
"prisms_surface.vtk", get_prims_surface(get_prisms()));
751 CHKERR save_entities(
"surface_skin.vtk", surface_skin);
753 CHKERR save_entities(
"crack_front.vtk", crack_front);
755 CHKERR save_entities(
"crack_front_nodes.vtk", crack_front_nodes);
757 CHKERR save_entities(
"end_nodes.vtk", end_nodes);
759 auto surface_snap = [&]() {
762 surface_skin,
fixedEdges, front_params.tolRelSnapToFixedEdges,
763 front_params.tolAbsSnapToFixedEdges,
nullptr,
debug);
767 auto move_front_nodes = [&]() {
770 auto w = get_tag_data(th_w, crack_front_nodes);
771 auto griffith_forces_projected =
772 get_tag_data(th_griffith_force_projected, crack_front_nodes);
773 auto area_change = get_tag_data(th_area_change, crack_front_nodes);
774 auto mesh_node_positions = get_coords_vec(crack_front_nodes);
776 const double ave_a_change =
777 get_ave_a_change(crack_front_nodes, area_change);
780 for (
auto node : crack_front_nodes) {
782 auto front_disp = get_node_vec(
w, nn);
783 auto force_projected = get_node_vec(griffith_forces_projected, nn);
784 auto a_change = get_node_vec(area_change, nn);
785 auto position = get_node_vec(mesh_node_positions, nn);
787 CHKERR check_for_nan(front_disp,
"front_disp");
788 CHKERR check_for_nan(force_projected,
"force_projected");
789 CHKERR check_for_nan(a_change,
"a_change");
791 VectorDouble3 disp = get_proj_front_disp(front_disp, force_projected);
794 CHKERR check_for_nan(position);
795 CHKERR check_for_nan(disp);
801 Range adj_crack_front_nodes_edges;
803 adj_crack_front_nodes_edges,
804 moab::Interface::UNION);
805 Range adj_crack_front_nodes_on_fixed_edges;
806 adj_crack_front_nodes_on_fixed_edges =
807 intersect(adj_crack_front_nodes_edges,
fixedEdges);
808 adj_crack_front_nodes_on_fixed_edges =
809 intersect(adj_crack_front_nodes_on_fixed_edges, bit_edges);
811 CHKERR check_for_nan(front_disp,
"front_disp");
812 CHKERR check_for_nan(coords,
"coords");
813 CHKERR check_for_nan(disp,
"disp");
819 FIXED_EDGE_KIND0 = 1 << 2,
820 FIXED_EDGE_KIND1 = 1 << 3
823 unsigned int kind = FREE;
824 if (check_if_node_is_in_range(node, get_nodes(
fixedEdges)))
826 if (check_if_node_is_in_range(node, end_nodes))
833 ave_a_change * a_change /
834 (norm_2(a_change) + std::numeric_limits<double>::epsilon());
835 if (is_point_is_outside_with_check(coords)) {
837 coords += front_params.skinFactor * scaled_a_change;
838 CHKERR check_for_nan(coords,
839 "Pushed through surface (out by factor)");
841 MOFEM_LOG(
"CPMeshCutWorld", Sev::verbose)
842 <<
"Pushed through surface (out by factor): ";
847 coords + front_params.skinFactor * scaled_a_change;
848 if (is_point_is_outside_with_check(tmp_coords)) {
850 scaled_a_change / (norm_2(scaled_a_change) +
851 std::numeric_limits<double>::epsilon());
852 const double dist = get_ray_pout(coords, scaled_a_change);
854 coords +=
dist * ray;
856 coords += front_params.skinFactor * scaled_a_change;
858 coords,
"Pushing through surface (out by skin factor)");
860 MOFEM_LOG(
"CPMeshCutWorld", Sev::verbose)
861 <<
"Pushing through surface (out by skin factor): ";
866 if (kind & END_NODE) {
868 double min_dist = front_params.skinFactor * ave_a_change;
870 if (kind & FIXED_EDGE) {
874 &coords[0], 1,
fixedEdges, &min_dist, &pt_out[0]);
875 if (min_dist < front_params.skinFactor * ave_a_change)
878 CHKERR check_for_nan(coords,
"End node on the fixed edge");
880 MOFEM_LOG(
"CPMeshCutWorld", Sev::verbose)
881 <<
"End node on the fixed edge: ";
888 &coords[0], 1,
fixedEdges, &min_dist, &pt_out[0]);
889 if (min_dist < front_params.skinFactor * ave_a_change) {
891 if (inner_prod(a_change, ray) > 0) {
893 front_params.skinFactor * ave_a_change * ray / norm_2(ray);
894 coords += ray + scaled_ray;
896 coords,
"Pushing through surface dist (close to edge)");
897 MOFEM_LOG(
"CPMeshCutWorld", Sev::verbose)
898 <<
"Pushing through surface dist (close to edge): ";
904 if (!adj_crack_front_nodes_on_fixed_edges.empty()) {
906 const int nb_edges_on_fix_edges =
907 adj_crack_front_nodes_on_fixed_edges.size() - 1;
908 adj_crack_front_nodes_on_fixed_edges =
909 intersect(adj_crack_front_nodes_on_fixed_edges, surface_skin);
912 if (!adj_crack_front_nodes_on_fixed_edges.empty()) {
914 const int nb_edges_on_fix_edges_and_surface_skin =
915 adj_crack_front_nodes_on_fixed_edges.size();
916 const int diff_nb_edges =
917 nb_edges_on_fix_edges - nb_edges_on_fix_edges_and_surface_skin;
922 if (diff_nb_edges > 1) {
925 CHKERR get_edge_delta(adj_crack_front_nodes_on_fixed_edges[0], node,
927 coords += front_params.cornerFactor * edge_delta;
928 MOFEM_LOG(
"CPMeshCutWorld", Sev::verbose)
929 <<
"Node on fixed edge (crack meeting fixed edge >|): ";
931 coords,
"Node on fixed edge (crack meeting fixed edge >|)");
933 kind |= FIXED_EDGE_KIND0;
936 }
else if ((kind & END_NODE) && (kind & FIXED_EDGE)) {
941 Range adj_edges_on_skin =
942 intersect(adj_crack_front_nodes_edges, surface_skin);
943 if (!adj_edges_on_skin.empty()) {
946 CHKERR get_edge_delta(adj_edges_on_skin[0], node, edge_delta);
947 coords += front_params.cornerFactor * edge_delta;
949 coords,
"Node on fixed edge (crack crossing fixed edge =|>)");
951 MOFEM_LOG(
"CPMeshCutWorld", Sev::verbose)
952 <<
"Node on fixed edge (crack crossing fixed edge =|>): ";
954 kind |= FIXED_EDGE_KIND1;
963 !(kind & FIXED_EDGE_KIND0) &&
965 !(kind & FIXED_EDGE_KIND1)) {
969 Range adj_edges = intersect(adj_crack_front_nodes_edges, crack_front);
971 CHKERR get_edge_delta(adj_edges[0], node, edge_delta);
972 coords += edge_delta * front_params.endNodeEdgeDeltaFactor;
973 CHKERR check_for_nan(coords,
"End node on the skin");
976 MOFEM_LOG(
"CPMeshCutWorld", Sev::verbose) <<
"End node on the skin: ";
979 MOFEM_LOG_C(
"CPMeshCutWorld", Sev::verbose,
"Disp front %s at node %lu",
980 get_str(front_disp).c_str(), node);
982 CHKERR set_tag(node, coords);
990 auto smooth_front = [&]() {
992 Range smoothed_nodes = subtract(crack_front_nodes, end_nodes);
993 std::vector<double> new_positions(3 * smoothed_nodes.size());
995 for (
auto node : smoothed_nodes) {
998 edges = intersect(edges, crack_front);
1001 edges_nodes.erase(node);
1002 if (edges_nodes.size() != 2)
1004 "Expected two nodes but have %d", edges_nodes.size());
1016 const double l0 = norm_2(delta0);
1017 const double l2 = norm_2(delta2);
1019 if ((l0 / (l0 + l2)) < 0.5)
1020 node1 += 2 * (0.5 - (l0 / (l0 + l2))) * (delta2 / l2) * l0;
1022 node1 += 2 * (0.5 - (l2 / (l0 + l2))) * (delta0 / l0) * l2;
1029 &*new_positions.begin());
1034 auto delete_mofem_meshsets = [&]() {
1045 auto build_surface = [&]() {
1048 Tag th_interface_side;
1063 if (!(conn[0] != conn[1] && conn[0] != conn[2] && conn[1] != conn[2]))
1065 "Imposible trinagle");
1071 CHKERR check_for_nan(coords.data());
1074 for (
auto nn : {0, 1, 2}) {
1075 auto node_map = verts_map.find(conn[nn]);
1076 if (node_map == verts_map.end()) {
1079 verts_map[conn[nn]] = new_node;
1081 new_verts[nn] = verts_map[conn[nn]];
1083 for (
auto nn : {0, 1, 2})
1084 if (new_verts[nn] == 0)
1086 "Imposible trinagle");
1089 new_verts[0] = new_verts[1];
1094 new_surface.insert(ele);
1100 auto build_front = [&]() {
1102 Range new_surface_skin;
1106 for (
auto f : crack_front) {
1111 for (
int nn = 0; nn != 2; nn++) {
1112 auto node_map = verts_map.find(conn[nn]);
1113 if (node_map != verts_map.end())
1114 new_verts[nn] = node_map->second;
1117 "Expected that vertex should be already created");
1121 moab::Interface::INTERSECT);
1122 edges = intersect(edges, new_surface_skin);
1123 new_front.merge(edges);
1125 fRont.swap(new_front);
1129 auto edges_snap = [&]() {
1131 Range new_surface_skin;
1134 new_surface_skin,
fixedEdges, front_params.tolRelSnapToFixedEdges,
1135 front_params.tolAbsSnapToFixedEdges,
nullptr,
debug);
1139 auto delete_meshsets_with_crack_surfaces = [&]() {
1142 for (
int ll = 0; ll != 19; ++ll)
1149 CHKERR move_front_nodes();
1151 CHKERR delete_mofem_meshsets();
1159 CHKERR delete_meshsets_with_crack_surfaces();
1167 if (!save_mesh.empty())
Public Member Functions inherited from MoFEM::UnknownInterface