Merge edges.
Sort all edges, where sorting is by quality calculated as edge length times quality of tets adjacent to the edge. Edge is merged if quality if the mesh is improved.
2307 CoreInterface &m_field =
cOre;
2315 CoreInterface &mField;
2316 const bool onlyIfImproveQuality;
2320 MergeNodes(CoreInterface &m_field,
const bool only_if_improve_quality,
2321 Tag
th,
bool update_mehsets)
2322 : mField(m_field), onlyIfImproveQuality(only_if_improve_quality),
2323 tH(
th), updateMehsets(update_mehsets) {
2324 mField.getInterface(nodeMergerPtr);
2326 NodeMergerInterface *nodeMergerPtr;
2329 bool add_child =
true) {
2334 CHKERR moab.get_adjacencies(conn, 2, 3,
false, vert_tets,
2335 moab::Interface::UNION);
2336 vert_tets = intersect(vert_tets, proc_tets);
2338 CHKERR nodeMergerPtr->mergeNodes(father, mother, out_tets, &vert_tets,
2339 onlyIfImproveQuality, 0, line_search,
2342 if (add_child && nodeMergerPtr->getSuccessMerge()) {
2344 Range::iterator lo, hi = proc_tets.begin();
2345 for (
auto pt = vert_tets.pair_begin(); pt != vert_tets.pair_end();
2347 lo = proc_tets.lower_bound(hi, proc_tets.end(), pt->first);
2348 if (lo != proc_tets.end()) {
2349 hi = proc_tets.upper_bound(lo, proc_tets.end(), pt->second);
2350 proc_tets.erase(lo, hi);
2354 proc_tets.merge(out_tets);
2356 auto &parent_child_map = nodeMergerPtr->getParentChildMap();
2358 struct ChangeChild {
2360 ChangeChild(
const EntityHandle child) : child(child) {}
2361 void operator()(NodeMergerInterface::ParentChild &p) {
2366 std::vector<decltype(parentsChildMap.get<0>().begin())> it_vec;
2367 it_vec.reserve(parentsChildMap.size());
2369 for (
auto &p : parent_child_map) {
2372 for (
auto it = parentsChildMap.get<0>().equal_range(p.pArent);
2373 it.first != it.second; ++it.first)
2374 it_vec.emplace_back(it.first);
2376 for (
auto it = parentsChildMap.get<1>().equal_range(p.pArent);
2377 it.first != it.second; ++it.first)
2378 it_vec.emplace_back(parentsChildMap.project<0>(it.first));
2380 for (
auto &it : it_vec)
2381 parentsChildMap.modify(it, ChangeChild(p.cHild));
2384 parentsChildMap.insert(parent_child_map.begin(),
2385 parent_child_map.end());
2391 Range ¬_merged_edges,
2397 std::vector<EntityHandle> parents_ents_vec(parentsChildMap.size());
2398 for (
auto &it : parentsChildMap)
2399 parents_ents_vec.emplace_back(it.pArent);
2401 parent_ents.insert_list(parents_ents_vec.begin(),
2402 parents_ents_vec.end());
2404 Range surf_parent_ents = intersect(new_surf, parent_ents);
2405 new_surf = subtract(new_surf, surf_parent_ents);
2406 Range child_surf_ents;
2407 CHKERR updateRangeByChilds(parentsChildMap, surf_parent_ents,
2409 new_surf.merge(child_surf_ents);
2411 Range edges_to_merge_parent_ents =
2412 intersect(edges_to_merge, parent_ents);
2413 edges_to_merge = subtract(edges_to_merge, edges_to_merge_parent_ents);
2414 Range merged_child_edge_ents;
2415 CHKERR updateRangeByChilds(parentsChildMap, edges_to_merge_parent_ents,
2416 merged_child_edge_ents);
2418 Range not_merged_edges_child_ents =
2419 intersect(not_merged_edges, parent_ents);
2421 subtract(not_merged_edges, not_merged_edges_child_ents);
2422 Range not_merged_child_edge_ents;
2423 CHKERR updateRangeByChilds(parentsChildMap, not_merged_edges_child_ents,
2424 not_merged_child_edge_ents);
2426 merged_child_edge_ents =
2427 subtract(merged_child_edge_ents, not_merged_child_edge_ents);
2428 edges_to_merge.merge(merged_child_edge_ents);
2429 not_merged_edges.merge(not_merged_child_edge_ents);
2431 if (updateMehsets) {
2433 (*mField.getInterface<MeshsetsManager>()), cubit_it)) {
2436 CHKERR moab.get_entities_by_handle(cubit_meshset, meshset_ents,
2439 CHKERR updateRangeByChilds(parentsChildMap, meshset_ents,
2441 CHKERR moab.add_entities(cubit_meshset, child_ents);
2451 std::vector<EntityHandle> childsVec;
2458 childsVec.reserve(parents.size());
2459 for (
auto pit = parents.pair_begin(); pit != parents.pair_end(); pit++) {
2460 auto it = parent_child_map.lower_bound(pit->first);
2461 if (it != parent_child_map.end()) {
2462 for (
auto hi_it = parent_child_map.upper_bound(pit->second);
2464 childsVec.emplace_back(it->cHild);
2467 childs.insert_list(childsVec.begin(), childsVec.end());
2477 CoreInterface &mField;
2480 LengthMap(CoreInterface &m_field, Tag
th,
double max_length)
2481 : tH(
th), mField(m_field), moab(m_field.get_moab()),
2485 acceptedThrasholdMergeQuality = 0.5;
2492 double acceptedThrasholdMergeQuality;
2497 double &ave)
const {
2500 std::array<double, 12> coords;
2502 VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
2503 VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
2506 struct NodeQuality {
2509 NodeQuality(
const EntityHandle node) : node(node), quality(1) {}
2512 typedef multi_index_container<
2513 NodeQuality, indexed_by<
2517 hashed_non_unique<tag<Ent_mi_tag>,
2519 &NodeQuality::node>>
2522 NodeQuality_sequence;
2524 NodeQuality_sequence node_quality_sequence;
2527 CHKERR moab.get_connectivity(tets, edges_nodes,
false);
2529 CHKERR moab.get_adjacencies(edges, 3,
false, edges_tets,
2530 moab::Interface::UNION);
2531 edges_tets = intersect(edges_tets, tets);
2533 for (
auto node : edges_nodes)
2534 node_quality_sequence.emplace_back(node);
2536 for (
auto tet : edges_tets) {
2538 CHKERR moab.get_connectivity(tet, conn, num_nodes,
true);
2540 CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2542 CHKERR moab.get_coords(conn, num_nodes, coords.data());
2546 for (
auto n : {0, 1, 2, 3}) {
2547 auto it = node_quality_sequence.get<1>().find(conn[
n]);
2548 if (it != node_quality_sequence.get<1>().end())
2549 const_cast<double &
>(it->quality) = std::min(q, it->quality);
2553 for (
auto edge : edges) {
2555 CHKERR moab.get_connectivity(edge, conn, num_nodes,
true);
2558 CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2560 CHKERR moab.get_coords(conn, num_nodes, coords.data());
2563 for (
auto n : {0, 1}) {
2564 auto it = node_quality_sequence.get<1>().find(conn[
n]);
2565 if (it != node_quality_sequence.get<1>().end())
2566 q = std::min(q, it->quality);
2569 if (q < acceptedThrasholdMergeQuality) {
2572 double val = pow(q, gammaQ) * pow(dot, 0.5 * gammaL);
2573 length_map.insert(LengthMapData(val, q, edge));
2578 for (LengthMapData_multi_index::nth_index<0>::type::iterator mit =
2579 length_map.get<0>().begin();
2580 mit != length_map.get<0>().end(); mit++) {
2581 ave += mit->qUality;
2583 ave /= length_map.size();
2593 CoreInterface &mField;
2595 Toplogy(CoreInterface &m_field, Tag
th) : mField(m_field), tH(
th) {}
2601 SURFACE_SKIN = 1 << 2,
2602 FRONT_ENDS = 1 << 3,
2604 FIX_CORNERS = 1 << 5
2607 typedef map<int, Range> SetsMap;
2610 const Range &fixed_edges,
2611 const Range &corner_nodes,
2612 const Range &constrain_surface,
2613 SetsMap &sets_map)
const {
2615 Skinner skin(&moab);
2618 sets_map[FIX_CORNERS].merge(corner_nodes);
2620 CHKERR moab.get_connectivity(fixed_edges, fixed_verts,
true);
2621 sets_map[FIX_EDGES].swap(fixed_verts);
2624 CHKERR skin.find_skin(0, tets,
false, tets_skin);
2625 Range tets_skin_edges;
2626 CHKERR moab.get_adjacencies(tets_skin, 1,
false, tets_skin_edges,
2627 moab::Interface::UNION);
2630 Range constrain_surface_verts;
2631 CHKERR moab.get_connectivity(constrain_surface, constrain_surface_verts,
2633 Range constrain_surface_edges;
2634 CHKERR moab.get_adjacencies(constrain_surface, 1,
false,
2635 constrain_surface_edges,
2636 moab::Interface::UNION);
2641 Range front_in_the_body;
2642 front_in_the_body = subtract(surface_skin, tets_skin_edges);
2643 Range front_in_the_body_verts;
2644 CHKERR moab.get_connectivity(front_in_the_body, front_in_the_body_verts,
2647 CHKERR skin.find_skin(0, front_in_the_body,
false, front_ends);
2649 intersect(front_in_the_body_verts, constrain_surface_verts));
2650 sets_map[FRONT_ENDS].swap(front_ends);
2652 Range surface_skin_verts;
2653 CHKERR moab.get_connectivity(surface_skin, surface_skin_verts,
true);
2654 sets_map[SURFACE_SKIN].swap(surface_skin_verts);
2657 Range surface_verts;
2659 sets_map[SURFACE_SKIN].merge(
2660 intersect(constrain_surface_verts, surface_verts));
2661 sets_map[SURFACE].swap(surface_verts);
2664 Range tets_skin_verts;
2665 CHKERR moab.get_connectivity(tets_skin, tets_skin_verts,
true);
2666 sets_map[SKIN].swap(tets_skin_verts);
2667 sets_map[SKIN].merge(constrain_surface_verts);
2670 CHKERR moab.get_connectivity(tets, tets_verts,
true);
2671 sets_map[FREE].swap(tets_verts);
2677 Range &proc_tets)
const {
2680 Range edges_to_merge_verts;
2681 CHKERR moab.get_connectivity(edges_to_merge, edges_to_merge_verts,
true);
2682 Range edges_to_merge_verts_tets;
2683 CHKERR moab.get_adjacencies(edges_to_merge_verts, 3,
false,
2684 edges_to_merge_verts_tets,
2685 moab::Interface::UNION);
2686 edges_to_merge_verts_tets = intersect(edges_to_merge_verts_tets, tets);
2687 proc_tets.swap(edges_to_merge_verts_tets);
2692 const Range &fixed_edges,
2693 const Range &corner_nodes,
2694 const Range &constrain_surface,
2695 Range &edges_to_merge,
2696 Range ¬_merged_edges) {
2701 Skinner skin(&moab);
2703 CHKERR skin.find_skin(0, tets,
false, tets_skin);
2708 Range constrain_surface_verts;
2709 CHKERR moab.get_connectivity(constrain_surface, constrain_surface_verts,
2711 Range constrain_surface_edges;
2712 CHKERR moab.get_adjacencies(constrain_surface, 1,
false,
2713 constrain_surface_edges,
2714 moab::Interface::UNION);
2717 Range tets_skin_edges;
2718 CHKERR moab.get_adjacencies(tets_skin, 1,
false, tets_skin_edges,
2719 moab::Interface::UNION);
2721 Range surface_front;
2722 surface_front = subtract(surface_skin, tets_skin_edges);
2723 Range surface_front_nodes;
2724 CHKERR moab.get_connectivity(surface_front, surface_front_nodes,
true);
2727 CHKERR skin.find_skin(0, surface_front,
false, ends_nodes);
2728 ends_nodes.merge(intersect(surface_front_nodes, constrain_surface_verts));
2731 surface_skin.merge(constrain_surface);
2732 tets_skin_edges.merge(constrain_surface_edges);
2735 Range surface_edges;
2736 CHKERR moab.get_adjacencies(
surface, 1,
false, surface_edges,
2737 moab::Interface::UNION);
2739 Range surface_edges_verts;
2740 CHKERR moab.get_connectivity(surface_edges, surface_edges_verts,
true);
2742 Range tets_skin_edges_verts;
2743 CHKERR moab.get_connectivity(tets_skin_edges, tets_skin_edges_verts,
2746 Range edges_to_remove;
2750 Range ents_nodes_and_edges;
2751 ents_nodes_and_edges.merge(tets_skin_edges_verts);
2752 ents_nodes_and_edges.merge(tets_skin_edges);
2753 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2756 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2757 not_merged_edges.merge(edges_to_remove);
2761 Range ents_nodes_and_edges;
2762 ents_nodes_and_edges.merge(surface_edges_verts);
2763 ents_nodes_and_edges.merge(surface_edges);
2764 ents_nodes_and_edges.merge(tets_skin_edges_verts);
2765 ents_nodes_and_edges.merge(tets_skin_edges);
2766 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2769 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2770 not_merged_edges.merge(edges_to_remove);
2773 Range fixed_edges_nodes;
2774 CHKERR moab.get_connectivity(fixed_edges, fixed_edges_nodes,
true);
2776 Range ents_nodes_and_edges;
2777 ents_nodes_and_edges.merge(fixed_edges_nodes);
2778 ents_nodes_and_edges.merge(ends_nodes);
2779 ents_nodes_and_edges.merge(corner_nodes);
2780 ents_nodes_and_edges.merge(fixed_edges);
2781 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2784 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2785 not_merged_edges.merge(edges_to_remove);
2788 CHKERR removeSelfConectingEdges(surface_edges, edges_to_remove,
false);
2789 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2790 not_merged_edges.merge(edges_to_remove);
2794 Range ents_nodes_and_edges;
2795 ents_nodes_and_edges.merge(surface_skin);
2796 ents_nodes_and_edges.merge(fixed_edges_nodes);
2797 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2800 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2801 not_merged_edges.merge(edges_to_remove);
2805 Range ents_nodes_and_edges;
2806 ents_nodes_and_edges.merge(surface_skin.subset_by_type(MBEDGE));
2807 ents_nodes_and_edges.merge(fixed_edges.subset_by_type(MBEDGE));
2808 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2811 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2812 not_merged_edges.merge(edges_to_remove);
2816 Range ents_nodes_and_edges;
2817 ents_nodes_and_edges.merge(surface_front_nodes);
2818 ents_nodes_and_edges.merge(surface_front);
2819 ents_nodes_and_edges.merge(tets_skin_edges_verts);
2820 ents_nodes_and_edges.merge(tets_skin_edges);
2821 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2824 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2825 not_merged_edges.merge(edges_to_remove);
2832 Range &edges_to_remove,
2837 Range ents_nodes = ents.subset_by_type(MBVERTEX);
2838 if (ents_nodes.empty()) {
2839 CHKERR moab.get_connectivity(ents, ents_nodes,
true);
2842 Range ents_nodes_edges;
2843 CHKERR moab.get_adjacencies(ents_nodes, 1,
false, ents_nodes_edges,
2844 moab::Interface::UNION);
2846 Range ents_nodes_edges_nodes;
2847 CHKERR moab.get_connectivity(ents_nodes_edges, ents_nodes_edges_nodes,
2850 ents_nodes_edges_nodes = subtract(ents_nodes_edges_nodes, ents_nodes);
2851 Range ents_nodes_edges_nodes_edges;
2852 CHKERR moab.get_adjacencies(ents_nodes_edges_nodes, 1,
false,
2853 ents_nodes_edges_nodes_edges,
2854 moab::Interface::UNION);
2857 subtract(ents_nodes_edges, ents_nodes_edges_nodes_edges);
2859 subtract(ents_nodes_edges, ents.subset_by_type(MBEDGE));
2861 edges_to_remove.swap(ents_nodes_edges);
2863 CHKERR SaveData(moab)(
"edges_to_remove.vtk", edges_to_remove);
2869 Range not_merged_edges;
2871 .removeBadEdges(
surface, tets, fixed_edges, corner_nodes,
2873 Toplogy::SetsMap sets_map;
2878 for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2879 sit != sets_map.rend(); sit++) {
2880 std::string name =
"classification_verts_" +
2881 boost::lexical_cast<std::string>(sit->first) +
".vtk";
2882 if (!sit->second.empty())
2883 CHKERR SaveData(moab)(name, sit->second);
2887 CHKERR Toplogy(m_field,
th).getProcTets(tets, edges_to_merge, proc_tets);
2888 out_tets = subtract(tets, proc_tets);
2891 Range all_out_ents = out_tets;
2892 for (
int dd = 2;
dd >= 0;
dd--) {
2893 CHKERR moab.get_adjacencies(out_tets,
dd,
false, all_out_ents,
2894 moab::Interface::UNION);
2896 CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(all_out_ents,
2900 int nb_nodes_merged = 0;
2904 auto save_merge_step = [&](
const int pp,
const Range collapsed_edges) {
2907 CHKERR moab.get_adjacencies(proc_tets, 2,
false, adj_faces,
2908 moab::Interface::UNION);
2910 name =
"node_merge_step_" + boost::lexical_cast<std::string>(pp) +
".vtk";
2912 name, unite(intersect(new_surf, adj_faces), collapsed_edges));
2914 "edges_to_merge_step_" + boost::lexical_cast<std::string>(pp) +
".vtk";
2916 name, unite(intersect(new_surf, adj_faces), edges_to_merge));
2923 double ave0 = 0, ave = 0, min = 0, min_p = 0, min_pp;
2926 int nb_nodes_merged_p = nb_nodes_merged;
2933 if(!length_map.empty())
2934 min = length_map.get<2>().begin()->qUality;
2940 Range collapsed_edges;
2941 MergeNodes merge_nodes(m_field,
true,
th, update_meshsets);
2943 for (
auto mit = length_map.get<0>().begin();
2944 mit != length_map.get<0>().end(); mit++, nn++) {
2948 auto get_father_and_mother =
2953 CHKERR moab.get_connectivity(mit->eDge, conn, num_nodes,
true);
2954 std::array<int, 2> conn_type = {0, 0};
2955 for (
int nn = 0; nn != 2; nn++) {
2956 for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2957 sit != sets_map.rend(); sit++) {
2958 if (sit->second.find(conn[nn]) != sit->second.end()) {
2959 conn_type[nn] |= sit->first;
2963 int type_father, type_mother;
2964 if (conn_type[0] > conn_type[1]) {
2967 type_father = conn_type[0];
2968 type_mother = conn_type[1];
2972 type_father = conn_type[1];
2973 type_mother = conn_type[0];
2975 if (type_father == type_mother) {
2981 int line_search = 0;
2983 CHKERR get_father_and_mother(father, mother, line_search);
2984 CHKERR merge_nodes.mergeNodes(line_search, father, mother, proc_tets);
2985 if (m_field.getInterface<NodeMergerInterface>()->getSuccessMerge()) {
2986 const EntityHandle father_and_mother[] = {father, mother};
2988 CHKERR moab.get_adjacencies(father_and_mother, 1, 3,
false, adj_tets);
2989 Range adj_tets_nodes;
2990 CHKERR moab.get_connectivity(adj_tets, adj_tets_nodes,
true);
2992 CHKERR moab.get_adjacencies(adj_tets_nodes, 1,
false, adj_edges,
2993 moab::Interface::UNION);
2994 for (
auto ait : adj_edges) {
2995 auto miit = length_map.get<1>().find(ait);
2996 if (miit != length_map.get<1>().end())
2997 (
const_cast<LengthMapData &
>(*miit)).skip =
true;
3000 collapsed_edges.insert(mit->eDge);
3003 if (nn >
static_cast<int>(length_map.size() / fraction_level))
3005 if (mit->qUality > ave)
3010 CHKERR merge_nodes.updateRangeByChilds(new_surf, edges_to_merge,
3011 not_merged_edges,
true);
3014 "(%d) Number of nodes merged %d ave q %3.4e min q %3.4e", pp,
3015 nb_nodes_merged, ave, min);
3018 CHKERR save_merge_step(pp + 1, collapsed_edges);
3020 if (nb_nodes_merged == nb_nodes_merged_p)
3022 if (min > 1e-2 && min == min_pp)
3028 CHKERR moab.get_adjacencies(proc_tets, 1,
false, adj_edges,
3029 moab::Interface::UNION);
3030 edges_to_merge = intersect(edges_to_merge, adj_edges);
3032 .removeBadEdges(new_surf, proc_tets, fixed_edges, corner_nodes,
3036 auto reconstruct_refined_ents = [&]() {
3047 CHKERR reconstruct_refined_ents();
3050 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(proc_tets,
3053 out_tets.merge(proc_tets);
3055 CHKERR moab.get_adjacencies(out_tets, 2,
false, adj_faces,
3056 moab::Interface::UNION);
3057 new_surf = intersect(new_surf, adj_faces);
Public Member Functions inherited from MoFEM::UnknownInterface