209 auto refined_ents_ptr = m_field.get_ref_ents();
210 auto refined_finite_elements_ptr = m_field.get_ref_finite_elements();
211 ReadUtilIface *read_util;
214 CHKERR m_field.get_moab().query_interface(read_util);
218 map<EntityHandle, EntityHandle> entParentMap;
225 for (Range::iterator eit = ref_edges->begin(); eit != ref_edges->end();
227 RefEntity_multiIndex::index<
228 Composite_ParentEnt_And_EntType_mi_tag>::type::iterator vit =
229 ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
230 boost::make_tuple(MBVERTEX, *eit));
232 ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
233 RefEntity_multiIndex::iterator e_eit = ref_ents_ptr->find(*eit);
234 if (e_eit == ref_ents_ptr->end()) {
236 "Edge not found %ld", *eit);
238 cerr <<
"Parent edge" << endl << **e_eit << endl;
239 if (entParentMap.find(*eit) != entParentMap.end()) {
240 RefEntity_multiIndex::iterator v_eit =
241 ref_ents_ptr->find(entParentMap[*eit]);
242 if (v_eit == ref_ents_ptr->end()) {
246 cerr <<
"Vertex " << **v_eit << endl;
248 RefEntity_multiIndex::index<Ent_Ent_mi_tag>::type::iterator ee_it,
250 ee_it = ref_ents_ptr->get<Ent_Ent_mi_tag>().lower_bound(*eit);
251 ee_hi_it = ref_ents_ptr->get<Ent_Ent_mi_tag>().upper_bound(*eit);
252 for (; ee_it != ee_hi_it; ++ee_it) {
253 cerr <<
"Ent having edge parent by parent " << **ee_it << endl;
256 tmp_index.insert(ref_ents_ptr->begin(), ref_ents_ptr->end());
257 RefEntity_multiIndex::index<
258 Composite_ParentEnt_And_EntType_mi_tag>::type::iterator vvit =
259 tmp_index.get<Composite_ParentEnt_And_EntType_mi_tag>().find(
260 boost::make_tuple(MBVERTEX, *eit));
262 tmp_index.get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
263 cerr <<
"Tmp idx Vertex " << **vvit << endl;
266 "No vertex on trim edges, that make no sense");
268 entParentMap[vit->get()->getParentEnt()] = vit->get()->getEnt();
276 map<EntityHandle, EntityHandle> parentsToChange;
282 RefEntity_multiIndex::iterator it = ref_ents_ptr->find(ent);
283 if (it != ref_ents_ptr->end()) {
284 if (it->get()->getParentEnt() != parent && ent != parent) {
285 parentsToChange[ent] = parent;
297 for (map<EntityHandle, EntityHandle>::iterator mit =
298 parentsToChange.begin();
299 mit != parentsToChange.end(); ++mit) {
300 RefEntity_multiIndex::iterator it = ref_ents_ptr->find(mit->first);
302 ->modify(it, RefEntity_change_parent(mit->second));
305 "impossible to set parent");
311 SetParent set_parent;
313 Range ents_to_set_bit;
317 CHKERR check(ref_edges_ptr, refined_ents_ptr,
cOre);
323 typedef const RefEntity_multiIndex::index<
324 Composite_EntType_and_ParentEntType_mi_tag>
::type RefEntsByComposite;
325 RefEntsByComposite &ref_ents =
326 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
328 ref_parent_ents_view.insert(
329 ref_ents.lower_bound(boost::make_tuple(MBVERTEX, MBEDGE)),
330 ref_ents.upper_bound(boost::make_tuple(MBVERTEX, MBEDGE)));
332 RefElementByEnt &ref_finite_element =
333 refined_finite_elements_ptr->get<Ent_mi_tag>();
334 typedef const RefElement_multiIndex_parents_view::index<
335 Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>
::type
336 RefEntByParentAndRefEdges;
338 ref_ele_parent_view.insert(
339 refined_finite_elements_ptr->get<Ent_mi_tag>().lower_bound(
340 get_id_for_min_type<MBTET>()),
341 refined_finite_elements_ptr->get<Ent_mi_tag>().upper_bound(
342 get_id_for_max_type<MBTET>()));
343 RefEntByParentAndRefEdges &ref_ele_by_parent_and_ref_edges =
345 .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
347 if (respect_interface) {
349 "not implemented, set last parameter in refine_TET to false");
352 Range tets = _tets.subset_by_type(MBTET);
354 std::vector<EntityHandle> parent_ents_refined_and_created;
355 std::vector<EntityHandle> vertices_of_created_tets;
356 std::vector<BitRefEdges> parent_edges_bit_vec;
357 std::vector<int> nb_new_tets_vec;
358 std::vector<int> sub_type_vec;
360 parent_ents_refined_and_created.reserve(tets.size());
361 vertices_of_created_tets.reserve(4 * tets.size());
362 parent_edges_bit_vec.reserve(tets.size());
363 nb_new_tets_vec.reserve(tets.size());
364 sub_type_vec.reserve(tets.size());
367 for (
auto tit : tets) {
369 if (ref_finite_element.find(tit) == ref_finite_element.end()) {
371 "this tet is not in refinedFiniteElements");
377 CHKERR moab.get_connectivity(tit, conn, num_nodes,
true);
382 int split_edges[] = {-1, -1, -1, -1, -1, -1};
384 for (
int ee = 0; ee != 6; ++ee) {
386 CHKERR moab.side_element(tit, 1, ee, edge);
387 RefEntity_multiIndex_view_by_hashed_parent_entity::iterator eit =
388 ref_parent_ents_view.find(edge);
389 if (eit != ref_parent_ents_view.end()) {
390 if (((*eit)->getBitRefLevel() & bit).any()) {
391 edge_new_nodes[ee] = (*eit)->getEnt();
395 moab.get_connectivity(edge, conn_edge, num_nodes,
true);
396 if (conn_edge[0] == edge_new_nodes[ee]) {
398 "node 0 on the edges is mid node, that make no sense");
400 if (conn_edge[1] == edge_new_nodes[ee]) {
402 "node 1 on the edges is mid node, that make no sense");
405 split_edges[parent_edges_bit.count()] = ee;
406 parent_edges_bit.set(ee, 1);
413 for (
int ee = 0; ee != 6; ee++) {
416 for (
int nn = 0; nn != 4; nn++) {
417 if (conn[nn] == edge_new_nodes[ee]) {
418 std::cerr <<
"problem on edge: " << ee << endl;
419 std::cerr <<
"tet conn : " << conn[0] <<
" " << conn[1] <<
" "
420 << conn[2] <<
" " << conn[3] <<
" : "
421 << edge_new_nodes[ee] << std::endl;
422 for (
int eee = 0; eee != 6; ++eee) {
424 CHKERR moab.side_element(tit, 1, eee, edge);
427 CHKERR moab.get_connectivity(edge, conn_edge, num_nodes,
true);
428 std::cerr << eee <<
" : " << conn_edge[0] <<
" " << conn_edge[1]
432 "nodes used to refine are not part of tet");
440 std::copy(&conn[0], &conn[4], &_conn_[0]);
444 std::fill(&new_tets_conns[0], &new_tets_conns[8 * 4],
no_handle);
446 int sub_type = -1, nb_new_tets = 0;
447 switch (parent_edges_bit.count()) {
449 RefEntity_multiIndex::iterator tit_miit;
450 tit_miit = refined_ents_ptr->find(tit);
451 if (tit_miit == refined_ents_ptr->end())
453 "can not find this tet");
454 ents_to_set_bit.insert(tit);
459 tet_type_1(_conn_, split_edges[0], edge_new_nodes[split_edges[0]],
465 tet_type_2(_conn_, split_edges, edge_new_nodes, new_tets_conns);
466 if (sub_type & (4 | 8 | 16)) {
469 }
else if (sub_type == 1) {
477 tet_type_3(_conn_, split_edges, edge_new_nodes, new_tets_conns);
481 }
else if (sub_type <= 7) {
488 tet_type_4(_conn_, split_edges, edge_new_nodes, new_tets_conns);
492 }
else if (sub_type <= 7) {
498 sub_type =
tet_type_5(moab, _conn_, edge_new_nodes, new_tets_conns);
503 tet_type_6(moab, _conn_, edge_new_nodes, new_tets_conns);
511 auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.equal_range(
512 boost::make_tuple(tit, parent_edges_bit.to_ulong()));
514 std::vector<EntityHandle> ents_to_set_bit_vec;
515 if (std::distance(it_by_ref_edges.first, it_by_ref_edges.second) ==
516 static_cast<size_t>(nb_new_tets)) {
517 ents_to_set_bit_vec.reserve(nb_new_tets);
518 for (
int tt = 0; it_by_ref_edges.first != it_by_ref_edges.second;
519 it_by_ref_edges.first++, tt++) {
520 auto tet = (*it_by_ref_edges.first)->getEnt();
521 ents_to_set_bit_vec.emplace_back(tet);
525 RefEntity_multiIndex::iterator ref_tet_it;
526 ref_tet_it = refined_ents_ptr->find(tet);
527 if (ref_tet_it == refined_ents_ptr->end()) {
529 "Tetrahedron should be in database");
533 ents_to_set_bit.insert_list(ents_to_set_bit_vec.begin(),
534 ents_to_set_bit_vec.end());
540 parent_ents_refined_and_created.emplace_back(tit);
541 for (
int tt = 0; tt != nb_new_tets; ++tt) {
542 for (
auto nn : {0, 1, 2, 3})
543 vertices_of_created_tets.emplace_back(new_tets_conns[4 * tt + nn]);
545 parent_edges_bit_vec.emplace_back(parent_edges_bit);
546 nb_new_tets_vec.emplace_back(nb_new_tets);
547 sub_type_vec.emplace_back(sub_type);
554 const int nb_new_tets = vertices_of_created_tets.size() / 4;
555 read_util->get_element_connect(nb_new_tets, 4, MBTET, 0, start_e, conn_moab);
556 std::copy(vertices_of_created_tets.begin(), vertices_of_created_tets.end(),
558 CHKERR read_util->update_adjacencies(start_e, nb_new_tets, 4, conn_moab);
559 ents_to_set_bit.insert(start_e, start_e + nb_new_tets - 1);
562 for (
auto d : {1, 2}) {
564 CHKERR moab.get_adjacencies(ents_to_set_bit,
d,
true, ents,
565 moab::Interface::UNION);
569 for (
int idx = 0; idx != parent_ents_refined_and_created.size(); ++idx) {
571 const EntityHandle tit = parent_ents_refined_and_created[idx];
572 const BitRefEdges &parent_edges_bit = parent_edges_bit_vec[idx];
573 const int nb_new_tets = nb_new_tets_vec[idx];
574 const int sub_type = sub_type_vec[idx];
576 std::array<EntityHandle, 8> ref_tets;
577 for (
int tt = 0; tt != nb_new_tets; ++tt, ++start_e)
578 ref_tets[tt] = start_e;
583 ref_type[0] = parent_edges_bit.count();
584 ref_type[1] = sub_type;
585 for (
int tt = 0; tt != nb_new_tets; ++tt) {
596 map_ref_nodes_by_edges;
599 CHKERR moab.get_adjacencies(&tit, 1, 1,
false, tet_edges);
600 for (
auto edge : tet_edges) {
601 RefEntity_multiIndex_view_by_hashed_parent_entity::iterator eit =
602 ref_parent_ents_view.find(edge);
603 if (eit != ref_parent_ents_view.end()) {
604 if (((*eit)->getBitRefLevel() & bit).any()) {
605 map_ref_nodes_by_edges[(*eit)->getParentEnt()] =
606 eit->get()->getEnt();
613 Range tit_edges, tit_faces;
614 CHKERR moab.get_adjacencies(&tit, 1, 1,
false, tit_edges);
615 CHKERR moab.get_adjacencies(&tit, 1, 2,
false, tit_faces);
616 Range edges_nodes[6], faces_nodes[4];
620 Range::iterator eit = tit_edges.begin();
621 for (
int ee = 0; eit != tit_edges.end(); eit++, ee++) {
622 CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee],
true);
623 std::map<EntityHandle, EntityHandle>::iterator map_miit =
624 map_ref_nodes_by_edges.find(*eit);
625 if (map_miit != map_ref_nodes_by_edges.end()) {
626 edges_nodes[ee].insert(map_miit->second);
632 Range::iterator fit = tit_faces.begin();
633 for (
int ff = 0; fit != tit_faces.end(); fit++, ff++) {
634 CHKERR moab.get_connectivity(&*fit, 1, faces_nodes[ff],
true);
637 CHKERR moab.get_adjacencies(&*fit, 1, 1,
false, fit_edges);
638 for (Range::iterator eit2 = fit_edges.begin(); eit2 != fit_edges.end();
640 std::map<EntityHandle, EntityHandle>::iterator map_miit =
641 map_ref_nodes_by_edges.find(*eit2);
642 if (map_miit != map_ref_nodes_by_edges.end()) {
643 faces_nodes[ff].insert(map_miit->second);
650 CHKERR moab.get_connectivity(&tit, 1, tet_nodes,
true);
651 for (std::map<EntityHandle, EntityHandle>::iterator map_miit =
652 map_ref_nodes_by_edges.begin();
653 map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
654 tet_nodes.insert(map_miit->second);
658 CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 1,
false,
659 ref_edges, moab::Interface::UNION);
661 for (Range::iterator reit = ref_edges.begin(); reit != ref_edges.end();
663 Range ref_edges_nodes;
664 CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes,
true);
665 if (ref_edges_nodes.size() != 2) {
667 "data inconsistency, edge should have 2 nodes");
671 for (; ee < 6; ee++) {
675 if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
677 CHKERR set_parent(*reit, edge, refined_ents_ptr,
cOre);
685 for (; ff < 4; ff++) {
688 if (intersect(faces_nodes[ff], ref_edges_nodes).size() == 2) {
690 CHKERR set_parent(*reit, face, refined_ents_ptr,
cOre);
699 if (intersect(tet_nodes, ref_edges_nodes).size() == 2) {
700 CHKERR set_parent(*reit, tit, refined_ents_ptr,
cOre);
707 "impossible refined edge");
711 CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 2,
false,
712 ref_faces, moab::Interface::UNION);
713 Tag th_interface_side;
714 const int def_side[] = {0};
715 CHKERR moab.tag_get_handle(
"INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
717 MB_TAG_CREAT | MB_TAG_SPARSE, def_side);
719 for (
auto rfit : ref_faces) {
720 Range ref_faces_nodes;
721 CHKERR moab.get_connectivity(&rfit, 1, ref_faces_nodes,
true);
724 for (; ff < 4; ff++) {
726 if (intersect(faces_nodes[ff], ref_faces_nodes).size() == 3) {
728 CHKERR set_parent(rfit, face, refined_ents_ptr,
cOre);
731 CHKERR moab.tag_get_data(th_interface_side, &face, 1, &side);
732 CHKERR moab.tag_set_data(th_interface_side, &rfit, 1, &side);
740 if (intersect(tet_nodes, ref_faces_nodes).size() == 3) {
741 CHKERR set_parent(rfit, tit, refined_ents_ptr,
cOre);
745 "impossible refined face");
751 CHKERR check(ref_edges_ptr, refined_ents_ptr,
cOre);
752 CHKERR set_parent(refined_ents_ptr);
754 CHKERR check(ref_edges_ptr, refined_ents_ptr,
cOre);
758 CHKERR check(ref_edges_ptr, refined_ents_ptr,
cOre);