11typedef multi_index_container<
12 boost::shared_ptr<RefElement>,
15 ordered_unique<tag<Ent_mi_tag>,
23 tag<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>,
28 const_mem_fun<RefElement,
int,
40 : cOre(const_cast<
Core &>(core)) {}
46 moab::Interface &moab = m_field.
get_moab();
49 CHKERR moab.get_entities_by_type(meshset, MBEDGE, edges, recursive);
52 CHKERR moab.get_entities_by_type(meshset, MBTET, tets, recursive);
53 CHKERR moab.get_adjacencies(tets, 1,
true, edges, moab::Interface::UNION);
56 CHKERR moab.get_entities_by_type(meshset, MBPRISM, prisms, recursive);
57 for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
60 CHKERR moab.get_connectivity(*pit, conn, num_nodes,
true);
61 assert(num_nodes == 6);
64 CHKERR moab.get_adjacencies(&conn[0], 2, 1,
true, edge);
65 assert(edge.size() == 1);
66 edges.insert(edge[0]);
68 CHKERR moab.get_adjacencies(&conn[1], 2, 1,
true, edge);
69 assert(edge.size() == 1);
70 edges.insert(edge[0]);
73 CHKERR moab.get_adjacencies(conn_edge2, 2, 1,
true, edge);
74 assert(edge.size() == 1);
75 edges.insert(edge[0]);
78 CHKERR moab.get_adjacencies(&conn[3], 2, 1,
true, edge);
79 assert(edge.size() == 1);
80 edges.insert(edge[0]);
82 CHKERR moab.get_adjacencies(&conn[4], 2, 1,
true, edge);
83 assert(edge.size() == 1);
84 edges.insert(edge[0]);
87 CHKERR moab.get_adjacencies(conn_edge8, 2, 1,
true, edge);
88 assert(edge.size() == 1);
89 edges.insert(edge[0]);
99 moab::Interface &moab = m_field.
get_moab();
104 .lower_bound(boost::make_tuple(MBVERTEX, MBEDGE));
107 .upper_bound(boost::make_tuple(MBVERTEX, MBEDGE));
109 ref_parent_ents_view.insert(miit, hi_miit);
111 Range edges = ents.subset_by_type(MBEDGE);
114 <<
"ref level " <<
bit <<
" nb. edges to refine " << edges.size();
117 std::array<std::vector<double>, 3> vert_coords;
118 for (
auto &vc : vert_coords)
119 vc.reserve(edges.size());
121 std::vector<EntityHandle> parent_edge;
122 parent_edge.reserve(edges.size());
124 std::array<double, 6> coords;
126 &coords[0], &coords[1], &coords[2]};
128 &coords[3], &coords[4], &coords[5]};
132 for (
auto p_eit = edges.pair_begin(); p_eit != edges.pair_end(); ++p_eit) {
133 auto miit_view = ref_parent_ents_view.lower_bound(p_eit->first);
135 Range edge_having_parent_vertex;
136 if (miit_view != ref_parent_ents_view.end()) {
137 for (
auto hi_miit_view = ref_parent_ents_view.upper_bound(p_eit->second);
138 miit_view != hi_miit_view; ++miit_view) {
139 edge_having_parent_vertex.insert(edge_having_parent_vertex.end(),
140 miit_view->get()->getParentEnt());
141 add_bit.insert(add_bit.end(), miit_view->get()->getEnt());
145 Range add_vertex_edges =
146 subtract(
Range(p_eit->first, p_eit->second), edge_having_parent_vertex);
148 for (
auto e : add_vertex_edges)
149 parent_edge.emplace_back(e);
151 for (
auto e : add_vertex_edges) {
155 CHKERR moab.get_connectivity(e, conn, num_nodes,
true);
156 if (PetscUnlikely(num_nodes != 2)) {
158 "edge should have 2 edges");
160 CHKERR moab.get_coords(conn, num_nodes, coords.data());
164 for (
auto j : {0, 1, 2})
165 vert_coords[
j].emplace_back(t_0(
j));
171 if (!vert_coords[0].empty()) {
172 ReadUtilIface *read_util;
173 CHKERR moab.query_interface(read_util);
174 int num_nodes = vert_coords[0].size();
175 vector<double *> arrays_coord;
176 CHKERR read_util->get_node_coords(3, num_nodes, 0, start_v, arrays_coord);
177 Range verts(start_v, start_v + num_nodes - 1);
178 for (
auto dd : {0, 1, 2}) {
179 std::copy(vert_coords[dd].begin(), vert_coords[dd].end(),
183 &*parent_edge.begin());
194 moab::Interface &moab = m_field.
get_moab();
197 CHKERR moab.get_entities_by_type(meshset, MBTET, tets,
false);
207 auto it = ref_ents_ptr->find(ent);
208 if (it != ref_ents_ptr->end()) {
209 if (it->get()->getParentEnt() != parent && ent != parent) {
224 for (
auto mit = parentsToChange.begin(); mit != parentsToChange.end();
226 auto it = ref_ents_ptr->find(mit->first);
231 "impossible to set parent");
243 auto set_edge_bits = [](moab::Interface &moab,
245 &ref_parent_ents_view,
250 for (
int ee = 0; ee != 6; ++ee) {
252 CHKERR moab.side_element(tet, 1, ee, edge);
253 auto eit = ref_parent_ents_view.find(edge);
254 if (eit != ref_parent_ents_view.end()) {
255 edge_new_nodes[ee] = (*eit)->getEnt();
256 split_edges[parent_edges_bit.count()] = ee;
257 parent_edges_bit.set(ee, 1);
276 auto set_edge_bits = [](moab::Interface &moab,
278 &ref_parent_ents_view,
283 for (
int ee = 0; ee != 6; ++ee) {
285 CHKERR moab.side_element(tet, 1, ee, edge);
286 auto eit = ref_parent_ents_view.find(edge);
287 if (eit != ref_parent_ents_view.end()) {
288 edge_new_nodes[ee] = (*eit)->getEnt();
289 split_edges[parent_edges_bit.count()] = ee;
290 parent_edges_bit.set(ee, 1);
294 if(parent_edges_bit.count() < 6)
295 parent_edges_bit.reset();
310 moab::Interface &moab = m_field.
get_moab();
313 CHKERR moab.get_entities_by_type(meshset, MBTET, tets,
false);
321 int verb,
const bool debug) {
324 moab::Interface &moab = m_field.
get_moab();
327 ReadUtilIface *read_util;
332 Range ents_to_set_bit;
337 auto get_parent_ents_view = [&](
const auto type_child,
338 const auto type_parent) {
342 ref_ents.equal_range(boost::make_tuple(type_child, type_parent));
346 using I =
decltype(range.first);
348 boost::function<
bool(
I it)> tester = [&](
I it) {
349 return ((*it)->getBitRefLevel() &
bit).any();
353 ref_parent_ents_view.insert(f, s);
359 return ref_parent_ents_view;
362 auto ref_parent_ents_view = get_parent_ents_view(MBVERTEX, MBEDGE);
364 auto tets = _tets.subset_by_type(MBTET);
366 auto get_ele_parent_view = [&]() {
367 auto &ref_ents = refined_finite_elements_ptr->get<
Ent_mi_tag>();
369 for (
auto p = tets.pair_begin();
p != tets.pair_end(); ++
p) {
370 auto tmi = ref_ents.lower_bound(
p->first);
371 auto tme = ref_ents.upper_bound(
p->second);
372 ref_ele_parent_view.insert(tmi, tme);
374 return ref_ele_parent_view;
377 auto ref_ele_parent_view = get_ele_parent_view();
379 std::vector<EntityHandle> parent_tets_refinded;
380 std::vector<EntityHandle> vertices_of_created_tets;
381 std::vector<BitRefEdges> parent_edges_bit_vec;
382 std::vector<int> nb_new_tets_vec;
383 std::vector<int> sub_type_vec;
385 parent_tets_refinded.reserve(tets.size());
386 vertices_of_created_tets.reserve(4 * tets.size());
387 parent_edges_bit_vec.reserve(tets.size());
388 nb_new_tets_vec.reserve(tets.size());
389 sub_type_vec.reserve(tets.size());
392 for (
auto tit : tets) {
397 CHKERR moab.get_connectivity(tit, conn, num_nodes,
true);
402 int split_edges[] = {-1, -1, -1, -1, -1, -1};
404 CHKERR set_edge_bits(moab, ref_parent_ents_view, tit, parent_edges_bit,
405 edge_new_nodes, split_edges);
407 if (parent_edges_bit.count()) {
411 std::copy(&conn[0], &conn[4], &_conn_[0]);
415 std::fill(&new_tets_conns[0], &new_tets_conns[8 * 4],
no_handle);
417 int sub_type = -1, nb_new_tets = 0;
418 switch (parent_edges_bit.count()) {
421 tet_type_1(_conn_, split_edges[0], edge_new_nodes[split_edges[0]],
427 tet_type_2(_conn_, split_edges, edge_new_nodes, new_tets_conns);
428 if (sub_type & (4 | 8 | 16)) {
431 }
else if (sub_type == 1) {
439 tet_type_3(_conn_, split_edges, edge_new_nodes, new_tets_conns);
443 }
else if (sub_type <= 7) {
450 tet_type_4(_conn_, split_edges, edge_new_nodes, new_tets_conns);
454 }
else if (sub_type <= 7) {
460 sub_type =
tet_type_5(moab, _conn_, edge_new_nodes, new_tets_conns);
465 tet_type_6(moab, _conn_, edge_new_nodes, new_tets_conns);
473 auto &ref_ele_by_parent_and_ref_edges =
476 auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.equal_range(
477 boost::make_tuple(tit, parent_edges_bit.to_ulong()));
479 std::vector<EntityHandle> ents_to_set_bit_vec;
480 if (std::distance(it_by_ref_edges.first, it_by_ref_edges.second) ==
481 static_cast<size_t>(nb_new_tets)) {
482 ents_to_set_bit_vec.reserve(nb_new_tets);
483 for (
int tt = 0; it_by_ref_edges.first != it_by_ref_edges.second;
484 it_by_ref_edges.first++, tt++) {
485 auto tet = (*it_by_ref_edges.first)->getEnt();
486 ents_to_set_bit_vec.emplace_back(tet);
490 auto ref_tet_it = refined_ents_ptr->find(tet);
491 if (ref_tet_it == refined_ents_ptr->end()) {
493 "Tetrahedron should be in database");
497 ents_to_set_bit.insert_list(ents_to_set_bit_vec.begin(),
498 ents_to_set_bit_vec.end());
504 parent_tets_refinded.emplace_back(tit);
505 for (
int tt = 0; tt != nb_new_tets; ++tt) {
506 for (
auto nn : {0, 1, 2, 3})
507 vertices_of_created_tets.emplace_back(new_tets_conns[4 * tt + nn]);
509 parent_edges_bit_vec.emplace_back(parent_edges_bit);
510 nb_new_tets_vec.emplace_back(nb_new_tets);
511 sub_type_vec.emplace_back(sub_type);
516 RefEntity_multiIndex::iterator tit_miit;
517 tit_miit = refined_ents_ptr->find(tit);
518 if (tit_miit == refined_ents_ptr->end())
520 "can not find this tet");
523 ents_to_set_bit.insert(tit);
530 const int nb_new_tets = vertices_of_created_tets.size() / 4;
531 read_util->get_element_connect(nb_new_tets, 4, MBTET, 0, start_e, conn_moab);
532 std::copy(vertices_of_created_tets.begin(), vertices_of_created_tets.end(),
534 CHKERR read_util->update_adjacencies(start_e, nb_new_tets, 4, conn_moab);
535 ents_to_set_bit.insert(start_e, start_e + nb_new_tets - 1);
539 for (
auto d : {1, 2}) {
540 CHKERR moab.get_adjacencies(ents_to_set_bit, d,
true, ents,
541 moab::Interface::UNION);
547 for (
int idx = 0; idx != parent_tets_refinded.size(); ++idx) {
550 const BitRefEdges &parent_edges_bit = parent_edges_bit_vec[idx];
551 const int nb_new_tets = nb_new_tets_vec[idx];
553 std::array<EntityHandle, 8> ref_tets;
554 for (
int tt = 0; tt != nb_new_tets; ++tt, ++start_e)
555 ref_tets[tt] = start_e;
560 nb_new_tets, &parent_edges_bit);
566 map_ref_nodes_by_edges;
569 CHKERR moab.get_adjacencies(&tit, 1, 1,
false, tit_edges);
570 for (
auto edge : tit_edges) {
571 auto eit = ref_parent_ents_view.find(edge);
572 if (eit != ref_parent_ents_view.end()) {
573 map_ref_nodes_by_edges[(*eit)->getParentEnt()] = eit->get()->getEnt();
580 CHKERR moab.get_adjacencies(&tit, 1, 2,
false, tit_faces);
581 Range edges_nodes[6], faces_nodes[4];
585 Range::iterator eit = tit_edges.begin();
586 for (
int ee = 0; eit != tit_edges.end(); eit++, ee++) {
587 CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee],
true);
588 auto map_miit = map_ref_nodes_by_edges.find(*eit);
589 if (map_miit != map_ref_nodes_by_edges.end()) {
590 edges_nodes[ee].insert(map_miit->second);
596 Range::iterator fit = tit_faces.begin();
597 for (
int ff = 0; fit != tit_faces.end(); fit++, ff++) {
598 CHKERR moab.get_connectivity(&*fit, 1, faces_nodes[ff],
true);
601 CHKERR moab.get_adjacencies(&*fit, 1, 1,
false, fit_edges);
602 for (Range::iterator eit2 = fit_edges.begin(); eit2 != fit_edges.end();
604 auto map_miit = map_ref_nodes_by_edges.find(*eit2);
605 if (map_miit != map_ref_nodes_by_edges.end()) {
606 faces_nodes[ff].insert(map_miit->second);
613 CHKERR moab.get_connectivity(&tit, 1, tet_nodes,
true);
614 for (std::map<EntityHandle, EntityHandle>::iterator map_miit =
615 map_ref_nodes_by_edges.begin();
616 map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
617 tet_nodes.insert(map_miit->second);
621 CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 1,
false,
622 ref_edges, moab::Interface::UNION);
624 for (Range::iterator reit = ref_edges.begin(); reit != ref_edges.end();
626 Range ref_edges_nodes;
627 CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes,
true);
631 for (; ee < 6; ee++) {
635 if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
637 CHKERR set_parent(*reit, edge, refined_ents_ptr,
cOre);
645 for (; ff < 4; ff++) {
648 if (intersect(faces_nodes[ff], ref_edges_nodes).size() == 2) {
650 CHKERR set_parent(*reit, face, refined_ents_ptr,
cOre);
659 if (intersect(tet_nodes, ref_edges_nodes).size() == 2) {
660 CHKERR set_parent(*reit, tit, refined_ents_ptr,
cOre);
667 "impossible refined edge");
671 CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 2,
false,
672 ref_faces, moab::Interface::UNION);
673 Tag th_interface_side;
674 const int def_side[] = {0};
675 CHKERR moab.tag_get_handle(
"INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
677 MB_TAG_CREAT | MB_TAG_SPARSE, def_side);
679 for (
auto rfit : ref_faces) {
680 Range ref_faces_nodes;
681 CHKERR moab.get_connectivity(&rfit, 1, ref_faces_nodes,
true);
684 for (; ff < 4; ff++) {
686 if (intersect(faces_nodes[ff], ref_faces_nodes).size() == 3) {
688 CHKERR set_parent(rfit, face, refined_ents_ptr,
cOre);
691 CHKERR moab.tag_get_data(th_interface_side, &face, 1, &side);
692 CHKERR moab.tag_set_data(th_interface_side, &rfit, 1, &side);
700 if (intersect(tet_nodes, ref_faces_nodes).size() == 3) {
701 CHKERR set_parent(rfit, tit, refined_ents_ptr,
cOre);
705 "impossible refined face");
710 CHKERR set_parent(refined_ents_ptr);
712 ents_to_set_bit,
bit,
true, verb, &ents);
719 moab::Interface &moab = m_field.
get_moab();
729 ref_ele_parent_view.insert(
730 refined_finite_elements_ptr->get<
Ent_mi_tag>().lower_bound(
731 get_id_for_min_type<MBPRISM>()),
732 refined_finite_elements_ptr->get<
Ent_mi_tag>().upper_bound(
733 get_id_for_max_type<MBPRISM>()));
734 auto &ref_ele_by_parent_and_ref_edges =
738 auto &ref_ents_by_comp =
741 ref_parent_ents_view.insert(
742 ref_ents_by_comp.lower_bound(boost::make_tuple(MBVERTEX, MBEDGE)),
743 ref_ents_by_comp.upper_bound(boost::make_tuple(MBVERTEX, MBEDGE)));
745 CHKERR moab.get_entities_by_type(meshset, MBPRISM, prisms,
false);
746 Range::iterator pit = prisms.begin();
747 for (; pit != prisms.end(); pit++) {
748 auto miit_prism = refined_ents_ptr->get<
Ent_mi_tag>().find(*pit);
749 if (miit_prism == refined_ents_ptr->end()) {
751 "this prism is not in ref database");
755 <<
"ref prism " << **miit_prism;
760 CHKERR moab.get_connectivity(*pit, conn, num_nodes,
true);
761 assert(num_nodes == 6);
764 for (
int ee = 0; ee < 3; ee++) {
765 CHKERR moab.side_element(*pit, 1, ee, edges[ee]);
767 for (
int ee = 6; ee < 9; ee++) {
768 CHKERR moab.side_element(*pit, 1, ee, edges[ee - 3]);
773 std::fill(&edge_nodes[0], &edge_nodes[6],
no_handle);
774 for (
int ee = 0; ee < 6; ee++) {
775 RefEntity_multiIndex_view_by_hashed_parent_entity::iterator miit_view =
776 ref_parent_ents_view.find(edges[ee]);
777 if (miit_view != ref_parent_ents_view.end()) {
778 if (((*miit_view)->getBitRefLevel() &
bit).any()) {
779 edge_nodes[ee] = (*miit_view)->getEnt();
784 if (split_edges.count() == 0) {
785 *(
const_cast<RefEntity *
>(miit_prism->get())->getBitRefLevelPtr()) |=
bit;
793 std::ostringstream ss;
794 ss <<
"prism split edges " << split_edges <<
" count "
795 << split_edges.count() << std::endl;
796 PetscPrintf(m_field.
get_comm(), ss.str().c_str());
800 std::fill(&new_prism_conn[0], &new_prism_conn[4 * 6],
no_handle);
801 int nb_new_prisms = 0;
802 switch (split_edges.count()) {
818 std::ostringstream ss;
819 ss << split_edges <<
" : [ " << conn[0] <<
" " << conn[1] <<
" "
820 << conn[2] <<
" " << conn[3] <<
" " << conn[4] <<
" " << conn[5]
825 std::bitset<4> ref_prism_bit(0);
826 auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.lower_bound(
827 boost::make_tuple(*pit, split_edges.to_ulong()));
828 auto hi_it_by_ref_edges = ref_ele_by_parent_and_ref_edges.upper_bound(
829 boost::make_tuple(*pit, split_edges.to_ulong()));
830 auto it_by_ref_edges2 = it_by_ref_edges;
831 for (
int pp = 0; it_by_ref_edges2 != hi_it_by_ref_edges;
832 it_by_ref_edges2++, pp++) {
834 *(
const_cast<RefElement *
>(it_by_ref_edges2->get())
835 ->getBitRefLevelPtr()) |=
bit;
836 ref_prism_bit.set(pp, 1);
838 std::ostringstream ss;
839 ss <<
"is refined " << *it_by_ref_edges2 << std::endl;
840 PetscPrintf(m_field.
get_comm(), ss.str().c_str());
843 if (it_by_ref_edges != hi_it_by_ref_edges) {
844 if (ref_prism_bit.count() != (
unsigned int)nb_new_prisms)
846 "data inconsistency");
850 for (
int pp = 0; pp < nb_new_prisms; pp++) {
852 std::ostringstream ss;
853 ss <<
"ref prism " << ref_prism_bit << std::endl;
854 PetscPrintf(m_field.
get_comm(), ss.str().c_str());
856 if (!ref_prism_bit.test(pp)) {
857 CHKERR moab.create_element(MBPRISM, &new_prism_conn[6 * pp], 6,
860 &ref_prisms[pp], 1, &*pit);
867 ->insert(boost::shared_ptr<RefEntity>(
new RefEntity(
869 std::pair<RefElement_multiIndex::iterator, bool> p_fe;
873 ->insert(boost::shared_ptr<RefElement>(
878 ref_prism_bit.set(pp);
881 std::ostringstream ss;
882 ss <<
"add prism: " << **p_fe.first << std::endl;
884 for (
int nn = 0; nn < 6; nn++) {
885 ss << new_prism_conn[nn] <<
" ";
889 PetscPrintf(m_field.
get_comm(), ss.str().c_str());
899 const bool recursive,
int verb) {
903 auto miit = refined_ents_ptr->find(meshset);
904 if (miit == refined_ents_ptr->end()) {
906 "this meshset is not in ref database");
909 meshset,
bit, meshset, MBMAXTYPE, recursive, verb);
910 *(
const_cast<RefEntity *
>(miit->get())->getBitRefLevelPtr()) |=
bit;
918 moab::Interface &moab = m_field.
get_moab();
921 CHKERR moab.get_entities_by_type(meshset, MBTRI, tris,
false);
931 auto set_edge_bits = [](moab::Interface &moab,
933 &ref_parent_ents_view,
938 for (
int ee = 0; ee != 3; ++ee) {
940 CHKERR moab.side_element(tet, 1, ee, edge);
941 auto eit = ref_parent_ents_view.find(edge);
942 if (eit != ref_parent_ents_view.end()) {
943 edge_new_nodes[ee] = (*eit)->getEnt();
944 split_edges[parent_edges_bit.count()] = ee;
945 parent_edges_bit.set(ee, 1);
964 auto set_edge_bits = [](moab::Interface &moab,
966 &ref_parent_ents_view,
971 for (
int ee = 0; ee != 3; ++ee) {
973 CHKERR moab.side_element(tet, 1, ee, edge);
974 auto eit = ref_parent_ents_view.find(edge);
975 if (eit != ref_parent_ents_view.end()) {
976 edge_new_nodes[ee] = (*eit)->getEnt();
977 split_edges[parent_edges_bit.count()] = ee;
978 parent_edges_bit.set(ee, 1);
982 if (parent_edges_bit.count() < 3) {
983 parent_edges_bit.reset();
999 moab::Interface &moab = m_field.
get_moab();
1002 CHKERR moab.get_entities_by_type(meshset, MBTRI, tris,
false);
1010 int verb,
const bool debug) {
1012 moab::Interface &moab = m_field.
get_moab();
1015 ReadUtilIface *read_util;
1019 auto get_parent_ents_view = [&](
const auto type_child,
1020 const auto type_parent) {
1024 ref_ents.equal_range(boost::make_tuple(type_child, type_parent));
1028 using I =
decltype(range.first);
1030 boost::function<
bool(
I it)> tester = [&](
I it) {
1031 return ((*it)->getBitRefLevel() &
bit).any();
1035 ref_parent_ents_view.insert(f, s);
1041 return ref_parent_ents_view;
1044 auto ref_parent_ents_view = get_parent_ents_view(MBVERTEX, MBEDGE);
1046 auto tris = _tris.subset_by_type(MBTRI);
1048 auto get_ele_parent_view = [&]() {
1049 auto &ref_ents = refined_finite_elements_ptr->get<
Ent_mi_tag>();
1051 for (
auto p = tris.pair_begin();
p != tris.pair_end(); ++
p) {
1052 auto tmi = ref_ents.lower_bound(
p->first);
1053 auto tme = ref_ents.upper_bound(
p->second);
1054 ref_ele_parent_view.insert(tmi, tme);
1056 return ref_ele_parent_view;
1059 auto ref_ele_parent_view = get_ele_parent_view();
1061 Range ents_to_set_bit;
1063 std::vector<EntityHandle> parent_tris_refined;
1064 std::vector<EntityHandle> vertices_of_created_tris;
1065 std::vector<BitRefEdges> parent_edges_bit_vec;
1066 std::vector<int> nb_new_tris_vec;
1068 parent_tris_refined.reserve(tris.size());
1069 vertices_of_created_tris.reserve(3 * tris.size());
1070 parent_edges_bit_vec.reserve(tris.size());
1071 nb_new_tris_vec.reserve(tris.size());
1074 for (
auto tit : tris) {
1079 CHKERR moab.get_connectivity(tit, conn, num_nodes,
true);
1084 int split_edges[] = {-1, -1, -1};
1086 CHKERR set_edge_bits(moab, ref_parent_ents_view, tit, parent_edges_bit,
1087 edge_new_nodes, split_edges);
1089 if (parent_edges_bit.count()) {
1093 std::copy(&conn[0], &conn[3], &_conn_[0]);
1097 std::fill(&new_tets_conns[0], &new_tets_conns[4 * 3],
no_handle);
1099 int nb_new_tris = 0;
1100 switch (parent_edges_bit.count()) {
1119 auto &ref_ele_by_parent_and_ref_edges =
1122 auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.equal_range(
1123 boost::make_tuple(tit, parent_edges_bit.to_ulong()));
1125 std::vector<EntityHandle> ents_to_set_bit_vec;
1126 if (std::distance(it_by_ref_edges.first, it_by_ref_edges.second) ==
1127 static_cast<size_t>(nb_new_tris)) {
1128 ents_to_set_bit_vec.reserve(nb_new_tris);
1129 for (
int tt = 0; it_by_ref_edges.first != it_by_ref_edges.second;
1130 it_by_ref_edges.first++, tt++) {
1131 auto tet = (*it_by_ref_edges.first)->getEnt();
1132 ents_to_set_bit_vec.emplace_back(tet);
1136 auto ref_tet_it = refined_ents_ptr->find(tet);
1137 if (ref_tet_it == refined_ents_ptr->end()) {
1139 "Tetrahedron should be in database");
1143 ents_to_set_bit.insert_list(ents_to_set_bit_vec.begin(),
1144 ents_to_set_bit_vec.end());
1150 parent_tris_refined.emplace_back(tit);
1151 for (
int tt = 0; tt != nb_new_tris; ++tt) {
1152 for (
auto nn : {0, 1, 2})
1153 vertices_of_created_tris.emplace_back(new_tets_conns[3 * tt + nn]);
1155 parent_edges_bit_vec.emplace_back(parent_edges_bit);
1156 nb_new_tris_vec.emplace_back(nb_new_tris);
1161 RefEntity_multiIndex::iterator tit_miit;
1162 tit_miit = refined_ents_ptr->find(tit);
1163 if (tit_miit == refined_ents_ptr->end())
1165 "can not find this tet");
1168 ents_to_set_bit.insert(tit);
1172 const int nb_new_tris = vertices_of_created_tris.size() / 3;
1180 CHKERR read_util->get_element_connect(nb_new_tris, 3, MBTRI, 0, start_e,
1182 std::copy(vertices_of_created_tris.begin(), vertices_of_created_tris.end(),
1184 CHKERR read_util->update_adjacencies(start_e, nb_new_tris, 3, conn_moab);
1185 ents_to_set_bit.insert(start_e, start_e + nb_new_tris - 1);
1189 CHKERR moab.add_entities(*meshset_ptr, ents_to_set_bit);
1190 CHKERR moab.write_file(
"out_refinded_debug.vtk",
"VTK",
"",
1191 meshset_ptr->get_ptr(), 1);
1197 for (
auto d : {1}) {
1198 CHKERR moab.get_adjacencies(ents_to_set_bit, d,
true, ents,
1199 moab::Interface::UNION);
1205 for (
int idx = 0; idx != parent_tris_refined.size(); ++idx) {
1207 const auto tit = parent_tris_refined[idx];
1208 const auto &parent_edges_bit =
1209 parent_edges_bit_vec[idx];
1210 const auto nb_new_tris =
1211 nb_new_tris_vec[idx];
1214 std::array<EntityHandle, 4> ref_tris;
1215 for (
int tt = 0; tt != nb_new_tris; ++tt, ++start_e)
1216 ref_tris[tt] = start_e;
1221 nb_new_tris, &parent_edges_bit);
1227 map_ref_nodes_by_edges;
1232 CHKERR moab.get_adjacencies(&tit, 1, 1,
false, tit_edges);
1233 for (
auto edge : tit_edges) {
1234 auto eit = ref_parent_ents_view.find(edge);
1235 if (eit != ref_parent_ents_view.end()) {
1236 map_ref_nodes_by_edges[(*eit)->getParentEnt()] = eit->get()->getEnt();
1240 Range edges_nodes[3];
1244 Range::iterator eit = tit_edges.begin();
1245 for (
int ee = 0; eit != tit_edges.end(); eit++, ee++) {
1246 CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee],
true);
1247 auto map_miit = map_ref_nodes_by_edges.find(*eit);
1248 if (map_miit != map_ref_nodes_by_edges.end()) {
1249 edges_nodes[ee].insert(map_miit->second);
1256 CHKERR moab.get_connectivity(&tit, 1, tri_nodes,
true);
1257 for (
auto map_miit = map_ref_nodes_by_edges.begin();
1258 map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
1259 tri_nodes.insert(map_miit->second);
1264 CHKERR moab.get_adjacencies(ref_tris.data(), nb_new_tris, 1,
false,
1265 ref_edges, moab::Interface::UNION);
1268 for (
auto reit = ref_edges.begin(); reit != ref_edges.end(); ++reit) {
1269 Range ref_edges_nodes;
1270 CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes,
true);
1273 for (; ee < 3; ee++) {
1277 if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
1279 CHKERR set_parent(*reit, edge, refined_ents_ptr,
cOre);
1288 if (intersect(tri_nodes, ref_edges_nodes).size() == 2) {
1289 CHKERR set_parent(*reit, tit, refined_ents_ptr,
cOre);
1296 "impossible refined edge");
1301 CHKERR set_parent(refined_ents_ptr);
1303 ents_to_set_bit,
bit,
true, verb, &ents);
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< hashed_non_unique< const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, hashed_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< boost::shared_ptr< RefEntity >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex_view_by_hashed_parent_entity
multi-index view of RefEntity by parent entity
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
multi_index_container< boost::shared_ptr< RefElement >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getEnt > > > > RefElement_multiIndex
virtual const RefEntity_multiIndex * get_ref_ents() const =0
Get the ref ents object.
virtual const RefElement_multiIndex * get_ref_finite_elements() const =0
Get the ref finite elements object.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFEDGES_SIZE > BitRefEdges
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
implementation of Data Operators for Forces and Sources
int tet_type_3(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
void tet_type_6(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
MoFEMErrorCode tri_type_2(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tris_conn)
MoFEMErrorCode prism_type_2(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
MoFEMErrorCode tri_type_1(const EntityHandle *conn, const int split_edge, const EntityHandle edge_new_node, EntityHandle *new_tris_conn)
RefEntityTmp< 0 > RefEntity
MoFEMErrorCode prism_type_3(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
MoFEMErrorCode tri_type_3(const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tris_conn)
MoFEMErrorCode prism_type_1(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
int tet_type_2(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
multi_index_container< boost::shared_ptr< RefElement >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getEnt > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getParentEnt > >, ordered_non_unique< tag< Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag >, composite_key< RefElement, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getParentEnt >, const_mem_fun< RefElement, int, &RefElement::getBitRefEdgesUlong > > > > > RefElement_multiIndex_parents_view
int tet_type_5(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
const EntityHandle no_handle
No entity handle is indicated by zero handle, i.e. root meshset.
int tet_type_4(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
void tet_type_1(const EntityHandle *conn, const int split_edge, const EntityHandle edge_new_node, EntityHandle *new_tets_conn)
auto rangeInserter(const I f, const I s, boost::function< bool(I it)> tester, boost::function< MoFEMErrorCode(I f, I s)> inserter)
Insert ranges.
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_non_unique< const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, hashed_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< boost::shared_ptr< RefEntity >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex_view_by_ordered_parent_entity
constexpr IntegrationType I
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.
MoFEMErrorCode addPrismToDatabase(const EntityHandle prism, int verb=DEFAULT_VERBOSITY)
add prim element
Tag get_th_RefBitLevel() const
Tag get_th_RefParentHandle() const
Tag get_th_RefBitEdge() const
Deprecated interface functions.
MoFEMErrorCode operator()(const EntityHandle ent, const EntityHandle parent, const RefEntity_multiIndex *ref_ents_ptr, MoFEM::Core &cOre)
map< EntityHandle, EntityHandle > parentsToChange
Mesh refinement interface.
MoFEMErrorCode refineTetsHangingNodes(const Range &tets, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
MoFEMErrorCode refineTets(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
MoFEMErrorCode addVerticesInTheMiddleOfEdges(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
make vertices in the middle of edges in meshset and add them to refinement levels defined by bit
MoFEMErrorCode refinePrisms(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET)
refine PRISM in the meshset
MoFEMErrorCode refineTrisHangingNodes(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TRI in the meshset
MeshRefinement(const MoFEM::Core &core)
boost::function< MoFEMErrorCode(moab::Interface &moab, RefEntity_multiIndex_view_by_ordered_parent_entity &ref_parent_ents_view, EntityHandle tet, BitRefEdges &parent_edges_bit, EntityHandle *edge_new_nodes, int *split_edges)> SetEdgeBitsFun
Functions setting edges for refinemnt on enetity level.
MoFEMErrorCode refineTris(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine triangles in the meshset
MoFEMErrorCode refineMeshset(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET)
refine meshset, i.e. add child of refined entities to meshset
keeps data about abstract PRISM finite element
keeps data about abstract refined finite element
interface_RefEntity< RefEntity > interface_type_RefEntity
virtual int getBitRefEdgesUlong() const
Struct keeps handle to refined handle.
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
EntityHandle getParentEnt() const
EntityHandle getEnt() const