321 {
322
325 auto refined_ents_ptr = m_field.get_ref_ents();
326 auto refined_finite_elements_ptr = m_field.get_ref_finite_elements();
327 ReadUtilIface *read_util;
329
330 CHKERR m_field.get_moab().query_interface(read_util);
331
332 Range ents_to_set_bit;
333
334
335
336
337 auto get_parent_ents_view = [&](const auto type_child,
338 const auto type_parent) {
339 auto &ref_ents =
340 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
341 auto range =
342 ref_ents.equal_range(boost::make_tuple(type_child, type_parent));
343
345
346 using I =
decltype(range.first);
347
348 boost::function<
bool(
I it)> tester = [&](
I it) {
349 return ((*it)->getBitRefLevel() &
bit).any();
350 };
351
353 ref_parent_ents_view.insert(
f, s);
354 return 0;
355 };
356
358
359 return ref_parent_ents_view;
360 };
361
362 auto ref_parent_ents_view = get_parent_ents_view(MBVERTEX, MBEDGE);
363
364 auto tets = _tets.subset_by_type(MBTET);
365
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);
373 }
374 return ref_ele_parent_view;
375 };
376
377 auto ref_ele_parent_view = get_ele_parent_view();
378
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;
384
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());
390
391
392 for (auto tit : tets) {
393
394
395 const EntityHandle *conn;
396 int num_nodes;
397 CHKERR moab.get_connectivity(tit, conn, num_nodes,
true);
398
399
401 EntityHandle edge_new_nodes[] = {0, 0, 0, 0, 0, 0};
402 int split_edges[] = {-1, -1, -1, -1, -1, -1};
403
404 CHKERR set_edge_bits(moab, ref_parent_ents_view, tit, parent_edges_bit,
405 edge_new_nodes, split_edges);
406
407 if (parent_edges_bit.count()) {
408
409
410 EntityHandle _conn_[4];
411 std::copy(&conn[0], &conn[4], &_conn_[0]);
412
413
414 EntityHandle new_tets_conns[8 * 4];
415 std::fill(&new_tets_conns[0], &new_tets_conns[8 * 4],
no_handle);
416
417 int sub_type = -1, nb_new_tets = 0;
418 switch (parent_edges_bit.count()) {
419 case 1:
420 sub_type = 0;
421 tet_type_1(_conn_, split_edges[0], edge_new_nodes[split_edges[0]],
422 new_tets_conns);
423 nb_new_tets = 2;
424 break;
425 case 2:
426 sub_type =
427 tet_type_2(_conn_, split_edges, edge_new_nodes, new_tets_conns);
428 if (sub_type & (4 | 8 | 16)) {
429 nb_new_tets = 3;
430 break;
431 } else if (sub_type == 1) {
432 nb_new_tets = 4;
433 break;
434 };
436 break;
437 case 3:
438 sub_type =
439 tet_type_3(_conn_, split_edges, edge_new_nodes, new_tets_conns);
440 if (sub_type <= 4) {
441 nb_new_tets = 4;
442 break;
443 } else if (sub_type <= 7) {
444 nb_new_tets = 5;
445 break;
446 }
448 case 4:
449 sub_type =
450 tet_type_4(_conn_, split_edges, edge_new_nodes, new_tets_conns);
451 if (sub_type == 0) {
452 nb_new_tets = 5;
453 break;
454 } else if (sub_type <= 7) {
455 nb_new_tets = 6;
456 break;
457 }
459 case 5:
460 sub_type =
tet_type_5(moab, _conn_, edge_new_nodes, new_tets_conns);
461 nb_new_tets = 7;
462 break;
463 case 6:
464 sub_type = 0;
465 tet_type_6(moab, _conn_, edge_new_nodes, new_tets_conns);
466 nb_new_tets = 8;
467 break;
468 default:
470 }
471
472
473 auto &ref_ele_by_parent_and_ref_edges =
474 ref_ele_parent_view
475 .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
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()));
478
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);
487
489
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");
494 }
495 }
496 }
497 ents_to_set_bit.insert_list(ents_to_set_bit_vec.begin(),
498 ents_to_set_bit_vec.end());
499
500 } else {
501
502
503
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]);
508 }
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);
512 }
513 } else {
514
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");
521 }
522
523 ents_to_set_bit.insert(tit);
524 }
525 }
526
527
528 EntityHandle start_e = 0;
529 EntityHandle *conn_moab;
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(),
533 conn_moab);
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);
536
537
539 for (
auto d : {1, 2}) {
540 CHKERR moab.get_adjacencies(ents_to_set_bit,
d,
true, ents,
541 moab::Interface::UNION);
542 }
543
544 SetParent set_parent;
545
546
547 for (int idx = 0; idx != parent_tets_refinded.size(); ++idx) {
548
549 const EntityHandle tit = parent_tets_refinded[idx];
550 const BitRefEdges &parent_edges_bit = parent_edges_bit_vec[idx];
551 const int nb_new_tets = nb_new_tets_vec[idx];
552 const int sub_type = sub_type_vec[idx];
553
554 std::array<EntityHandle, 8> ref_tets;
555 for (int tt = 0; tt != nb_new_tets; ++tt, ++start_e)
556 ref_tets[tt] = start_e;
557
558 if (nb_new_tets) {
559
560 int ref_type[2];
561 ref_type[0] = parent_edges_bit.count();
562 ref_type[1] = sub_type;
563
565 nb_new_tets, ref_type);
567 nb_new_tets, &parent_edges_bit);
569 nb_new_tets, &tit);
570
571
572 std::map<EntityHandle , EntityHandle >
573 map_ref_nodes_by_edges;
574
576 CHKERR moab.get_adjacencies(&tit, 1, 1,
false, tit_edges);
577 for (auto edge : tit_edges) {
578 auto eit = ref_parent_ents_view.find(edge);
579 if (eit != ref_parent_ents_view.end()) {
580 map_ref_nodes_by_edges[(*eit)->getParentEnt()] = eit->get()->getEnt();
581 }
582 }
583
584
585
587 CHKERR moab.get_adjacencies(&tit, 1, 2,
false, tit_faces);
588 Range edges_nodes[6], faces_nodes[4];
589
590
591
592 Range::iterator eit = tit_edges.begin();
593 for (int ee = 0; eit != tit_edges.end(); eit++, ee++) {
594 CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee],
true);
595 auto map_miit = map_ref_nodes_by_edges.find(*eit);
596 if (map_miit != map_ref_nodes_by_edges.end()) {
597 edges_nodes[ee].insert(map_miit->second);
598 }
599 }
600
601
602
603 Range::iterator fit = tit_faces.begin();
604 for (int ff = 0; fit != tit_faces.end(); fit++, ff++) {
605 CHKERR moab.get_connectivity(&*fit, 1, faces_nodes[ff],
true);
606
608 CHKERR moab.get_adjacencies(&*fit, 1, 1,
false, fit_edges);
609 for (Range::iterator eit2 = fit_edges.begin(); eit2 != fit_edges.end();
610 eit2++) {
611 auto map_miit = map_ref_nodes_by_edges.find(*eit2);
612 if (map_miit != map_ref_nodes_by_edges.end()) {
613 faces_nodes[ff].insert(map_miit->second);
614 }
615 }
616 }
617
618
620 CHKERR moab.get_connectivity(&tit, 1, tet_nodes,
true);
621 for (std::map<EntityHandle, EntityHandle>::iterator map_miit =
622 map_ref_nodes_by_edges.begin();
623 map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
624 tet_nodes.insert(map_miit->second);
625 }
627
628 CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 1,
false,
629 ref_edges, moab::Interface::UNION);
630
631 for (Range::iterator reit = ref_edges.begin(); reit != ref_edges.end();
632 reit++) {
633 Range ref_edges_nodes;
634 CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes,
true);
635
636
637 int ee = 0;
638 for (; ee < 6; ee++) {
639
640
641
642 if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
643 EntityHandle edge = tit_edges[ee];
644 CHKERR set_parent(*reit, edge, refined_ents_ptr,
cOre);
645 break;
646 }
647 }
648 if (ee < 6)
649 continue;
650
651 int ff = 0;
652 for (; ff < 4; ff++) {
653
654
655 if (intersect(faces_nodes[ff], ref_edges_nodes).size() == 2) {
656 EntityHandle face = tit_faces[ff];
657 CHKERR set_parent(*reit, face, refined_ents_ptr,
cOre);
658 break;
659 }
660 }
661 if (ff < 4)
662 continue;
663
664
665
666 if (intersect(tet_nodes, ref_edges_nodes).size() == 2) {
667 CHKERR set_parent(*reit, tit, refined_ents_ptr,
cOre);
668 continue;
669 }
670
671
672
674 "impossible refined edge");
675 }
676
678 CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 2,
false,
679 ref_faces, moab::Interface::UNION);
680 Tag th_interface_side;
681 const int def_side[] = {0};
682 CHKERR moab.tag_get_handle(
"INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
683 th_interface_side,
684 MB_TAG_CREAT | MB_TAG_SPARSE, def_side);
685
686 for (auto rfit : ref_faces) {
687 Range ref_faces_nodes;
688 CHKERR moab.get_connectivity(&rfit, 1, ref_faces_nodes,
true);
689
690 int ff = 0;
691 for (; ff < 4; ff++) {
692
693 if (intersect(faces_nodes[ff], ref_faces_nodes).size() == 3) {
694 EntityHandle face = tit_faces[ff];
695 CHKERR set_parent(rfit, face, refined_ents_ptr,
cOre);
696 int side = 0;
697
698 CHKERR moab.tag_get_data(th_interface_side, &face, 1, &side);
699 CHKERR moab.tag_set_data(th_interface_side, &rfit, 1, &side);
700 break;
701 }
702 }
703 if (ff < 4)
704 continue;
705
706
707 if (intersect(tet_nodes, ref_faces_nodes).size() == 3) {
708 CHKERR set_parent(rfit, tit, refined_ents_ptr,
cOre);
709 continue;
710 }
712 "impossible refined face");
713 }
714 }
715 }
716
717 CHKERR set_parent(refined_ents_ptr);
718 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
719 ents_to_set_bit,
bit,
true, verb, &ents);
721}
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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)
int tet_type_2(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
int tet_type_5(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
int tet_type_4(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
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.
Tag get_th_RefType() const
Public Member Functions inherited from MoFEM::UnknownInterface