419 std::vector<size_t>& bnd_dom_ie_list,
423 using namespace details;
425 level_set_epsilon = opt.
epsilon;
432 std::vector<size_type> dis_idof;
437 std::array<list<pair<element_type,size_type> >,
439 list<point_basic<T> > gamma_node_list;
440 std::vector<size_type> j(
d+1);
441 std::vector<point_basic<T> > x(
d+1);
442 const size_type not_marked = numeric_limits<size_t>::max();
443 const size_type undef = numeric_limits<size_t>::max();
449 set<size_type> ext_marked_node_idx;
450 list<to_solve> node_to_solve;
456 extern_edge (edge_ownership, std::make_pair(not_marked,
point_basic<T>()));
457 set<size_type> ext_marked_edge_idx;
458 list<to_solve> edge_to_solve;
467 bnd_dom_ie_list.resize(0);
472 for (
size_type ie = 0, ne =
lambda.size(), bnd_ie = 0; ie < ne; ie++) {
477 Xh.dis_idof (K, dis_idof);
478 f.resize (dis_idof.size());
479 for (
size_type loc_idof = 0, loc_ndof = dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
480 f [loc_idof] = fh.
dis_dof (dis_idof[loc_idof]);
482 bool do_intersect =
false;
485 do_intersect = belongs_to_band_t (
f);
488 do_intersect = belongs_to_band_T (
f);
492 error_macro(
"level set: element type `" << K.
name() <<
"' not yet supported");
494 if (! do_intersect)
continue;
495 bnd_dom_ie_list.push_back (ie);
499 x.resize (dis_idof.size());
500 for (
size_type loc_idof = 0, loc_ndof = dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
503 x [loc_idof] =
lambda.dis_node(dis_inod);
505 element_type S (alloc);
513 subcompute_matrix_t (x,
f, j, a, b, length);
514 if (is_zero(
f[j[1]]) && is_zero(
f[j[2]])) {
520 if (node_ownership.
is_owned (dis_inod)) {
521 size_type inod = dis_inod - first_dis_inod;
522 if (marked_node [inod] == not_marked) {
523 marked_node [inod] = gamma_node_list.size();
524 gamma_node_list.push_back (
lambda.node(inod));
532 extern_node.dis_entry (dis_inod) = my_proc;
535 size_type loc_iedg = edge_t_iloc (j[1], j[2]);
537 if (edge_ownership.
is_owned (dis_iedg)) {
538 size_type iedg = dis_iedg - first_dis_iedg;
539 if (marked_edge [iedg] == not_marked) {
541 marked_edge [iedg] = gamma_side_list [S.variant()].size();
543 for (
size_type k = 0; k < S.n_node(); k++) {
545 if (node_ownership.
is_owned (dis_inod)) {
546 size_type inod = dis_inod - first_dis_inod;
547 S[k] = marked_node [inod];
550 node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
551 ext_marked_node_idx.insert (dis_inod);
554 gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
567 size_type S_ige = gamma_side_list [S.variant()].size ();
570 if (! is_zero(
f[j[k+1]]) && ! is_zero(
f[j[0]])) {
572 size_type loc_iedg = edge_t_iloc (j[0], j[k+1]);
574 if (edge_ownership.
is_owned (dis_iedg)) {
575 size_type iedg = dis_iedg - first_dis_iedg;
576 if (marked_edge [iedg] == not_marked) {
577 marked_edge [iedg] = gamma_node_list.size();
578 gamma_node_list.push_back (xx[k]);
580 S[k] = marked_edge [iedg];
583 edge_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_iedg));
584 ext_marked_edge_idx.insert (dis_iedg);
585 extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
588 size_type dis_inod = (!is_zero(
f[j[0]])) ? K[j[k+1]] : K[j[0]];
589 if (node_ownership.
is_owned (dis_inod)) {
590 size_type inod = dis_inod - first_dis_inod;
591 if (marked_node [inod] == not_marked) {
592 marked_node [inod] = gamma_node_list.size();
593 gamma_node_list.push_back (
lambda.node(inod));
595 S[k] = marked_node [inod];
598 node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
599 ext_marked_node_idx.insert (dis_inod);
600 extern_node.dis_entry (dis_inod) = my_proc;
605 check_macro (S[0] != S[1] || S[0] == undef,
"degenerate 2d intersection");
606 gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
618 subcompute_matrix_triangle_T (x,
f, j, a, b, c, aire);
619 if (is_zero(
f[j[1]]) && is_zero(
f[j[2]]) && is_zero(
f[j[3]])) {
623 if (node_ownership.
is_owned (dis_inod)) {
624 size_type inod = dis_inod - first_dis_inod;
625 if (marked_node [inod] == not_marked) {
626 marked_node [inod] = gamma_node_list.size();
627 gamma_node_list.push_back (
lambda.node(inod));
635 extern_node.dis_entry (dis_inod) = my_proc;
638 size_type loc_ifac = face_T_iloc (j[1], j[2], j[3]);
640 if (face_ownership.
is_owned (dis_ifac)) {
641 size_type ifac = dis_ifac - first_dis_ifac;
642 if (marked_face [ifac] == not_marked) {
644 marked_face [ifac] = gamma_side_list [S.variant()].size();
646 for (
size_type k = 0; k < S.n_node(); k++) {
648 if (node_ownership.
is_owned (dis_inod)) {
649 size_type inod = dis_inod - first_dis_inod;
650 S[k] = marked_node [inod];
653 node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
654 ext_marked_node_idx.insert (dis_inod);
657 gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
671 size_type S_ige = gamma_side_list [S.variant()].size();
674 if (! is_zero(
f[j[k+1]]) && ! is_zero(
f[j[0]])) {
676 size_type loc_iedg = edge_T_iloc (j[0], j[k+1]);
678 if (edge_ownership.
is_owned (dis_iedg)) {
679 size_type iedg = dis_iedg - first_dis_iedg;
680 if (marked_edge [iedg] == not_marked) {
681 marked_edge [iedg] = gamma_node_list.size();
682 gamma_node_list.push_back (xx[k]);
684 S[k] = marked_edge [iedg];
687 edge_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_iedg));
688 ext_marked_edge_idx.insert (dis_iedg);
689 extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
692 size_type dis_inod = (!is_zero(
f[j[0]])) ? K[j[k+1]] : K[j[0]];
693 if (node_ownership.
is_owned (dis_inod)) {
694 size_type inod = dis_inod - first_dis_inod;
695 if (marked_node [inod] == not_marked) {
696 marked_node [inod] = gamma_node_list.size();
697 gamma_node_list.push_back (
lambda.node(inod));
699 S[k] = marked_node [inod];
702 node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
703 ext_marked_node_idx.insert (dis_inod);
704 extern_node.dis_entry (dis_inod) = my_proc;
710 (S[1] != S[2] || S[1] == undef) &&
711 (S[2] != S[0] || S[2] == undef),
"degenerate 3d intersection");
712 gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
716 subcompute_matrix_quadrilateral_T (x,
f, q, a, b, c,
d, aire);
727 if (! is_zero(
f[s[k]]) && ! is_zero(
f[s[k1]])) {
729 size_type loc_iedg = edge_T_iloc (s[k], s[k1]);
731 if (edge_ownership.
is_owned (dis_iedg)) {
732 size_type iedg = dis_iedg - first_dis_iedg;
733 if (marked_edge [iedg] == not_marked) {
734 marked_edge [iedg] = gamma_node_list.size();
735 gamma_node_list.push_back (xx[k]);
737 S[k] = marked_edge [iedg];
743 if (iloc_S1[k] != undef) edge_to_solve.push_back (to_solve (
reference_element::t, S1_ige, iloc_S1[k], dis_iedg));
744 if (iloc_S2[k] != undef) edge_to_solve.push_back (to_solve (
reference_element::t, S2_ige, iloc_S2[k], dis_iedg));
746 ext_marked_edge_idx.insert (dis_iedg);
747 extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
750 size_type dis_inod = is_zero(
f[s[k]]) ? K[s[k]] : K[s[k1]];
751 if (node_ownership.
is_owned (dis_inod)) {
752 size_type inod = dis_inod - first_dis_inod;
753 if (marked_node [inod] == not_marked) {
754 marked_node [inod] = gamma_node_list.size();
755 gamma_node_list.push_back (
lambda.node(inod));
757 S[k] = marked_node [inod];
763 if (iloc_S1[k] != undef) node_to_solve.push_back (to_solve (
reference_element::t, S1_ige, iloc_S1[k], dis_inod));
764 if (iloc_S2[k] != undef) node_to_solve.push_back (to_solve (
reference_element::t, S2_ige, iloc_S2[k], dis_inod));
766 ext_marked_node_idx.insert (dis_inod);
772 (S[1] != S[2] || S[1] == undef) &&
773 (S[2] != S[3] || S[2] == undef) &&
774 (S[3] != S[0] || S[3] == undef),
"degenerate 3d intersection");
776 gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
781 element_type S1 (alloc);
784 if (iloc_S1[k] != undef) S1 [iloc_S1[k]] = S[k];
787 (S1[1] != S1[2] || S1[1] == undef) &&
788 (S1[2] != S1[0] || S1[2] == undef),
"degenerate 3d intersection");
789 gamma_side_list [S1.variant()].push_back (make_pair(S1,bnd_ie));
791 element_type S2 (alloc);
794 if (iloc_S2[k] != undef) S2 [iloc_S2[k]] = S[k];
797 (S2[1] != S2[2] || S2[1] == undef) &&
798 (S2[2] != S2[0] || S2[2] == undef),
"degenerate 3d intersection");
799 gamma_side_list [S2.variant()].push_back (make_pair(S2,bnd_ie));
805 error_macro(
"level set intersection: element not yet implemented: " << K.
name());
810 extern_node.dis_entry_assembly();
811 extern_edge.dis_entry_assembly();
818 for (
size_type inod = 0, nnod = node_ownership.
size(); inod < nnod; inod++) {
819 if (!(extern_node [inod] != not_marked && marked_node [inod] == not_marked))
continue;
823 size_type dis_inod = first_dis_inod + inod;
825 dis_inod_set.
insert (dis_inod);
826 orphan_node.dis_entry (iproc) += dis_inod_set;
828 orphan_node.dis_entry_assembly();
834 for (
size_type iedg = 0, nedg = edge_ownership.
size(); iedg < nedg; iedg++) {
835 if (!(extern_edge [iedg].first != not_marked && marked_edge [iedg] == not_marked))
continue;
838 size_type iproc = extern_edge[iedg].first;
839 size_type dis_iedg = first_dis_iedg + iedg;
841 dis_iedg_set.
insert (dis_iedg);
842 orphan_edge.dis_entry (iproc) += dis_iedg_set;
844 orphan_edge.dis_entry_assembly ();
845 check_macro (orphan_edge[0].size() == 0,
"unexpected orphan edges");
848 size_type orphan_gamma_nnod = orphan_node[0].size();
849 size_type regular_gamma_nnod = gamma_node_list.size();
850 size_type gamma_nnod = regular_gamma_nnod + orphan_gamma_nnod;
854 for (
size_type inod = 0, nnod = node_ownership.
size(); inod < nnod; inod++) {
855 if (marked_node [inod] == not_marked)
continue;
856 marked_node [inod] += gamma_first_dis_inod;
858 for (
size_type iedg = 0, nedg = edge_ownership.
size(); iedg < nedg; iedg++) {
859 if (marked_edge [iedg] == not_marked)
continue;
860 marked_edge [iedg] += gamma_first_dis_inod;
863 for (index_set::const_iterator
864 iter = orphan_node[0].begin(),
865 last = orphan_node[0].end(); iter != last; ++iter) {
867 marked_node.dis_entry(dis_inod) = gamma_first_dis_inod + gamma_node_list.size();
868 gamma_node_list.push_back (
lambda.dis_node(dis_inod));
870 marked_node.dis_entry_assembly();
871 marked_edge.dis_entry_assembly();
877 gamma_list2disarray (gamma_node_list, gamma_side_list, comm,
d, gamma_node, gamma_side, sid_ie2bnd_ie);
883 for (
size_type ige = 0, nge = gamma_side[variant].size(); ige < nge; ige++) {
884 element_type& S = gamma_side[variant][ige];
885 for (
size_type loc_inod = 0, loc_nnod = S.n_node(); loc_inod < loc_nnod; loc_inod++) {
886 if (S[loc_inod] == undef)
continue;
887 S[loc_inod] += gamma_first_dis_inod;
894 marked_node.set_dis_indexes (ext_marked_node_idx);
895 for (list<to_solve>::const_iterator iter = node_to_solve.begin(),
896 last = node_to_solve.end(); iter != last; iter++) {
897 const to_solve& x = *iter;
898 element_type& S = gamma_side[x.variant][x.S_ige];
899 check_macro (S[x.k] == undef,
"external index already solved");
902 size_type gamma_dis_inod = marked_node.dis_at (dis_inod);
903 S[x.k] = gamma_dis_inod;
906 marked_edge.set_dis_indexes (ext_marked_edge_idx);
907 for (list<to_solve>::const_iterator iter = edge_to_solve.begin(),
908 last = edge_to_solve.end(); iter != last; iter++) {
909 const to_solve& x = *iter;
910 element_type& S = gamma_side[x.variant][x.S_ige];
911 check_macro (S[x.k] == undef,
"external index already solved");
914 size_type gamma_dis_inod = marked_edge.dis_at (dis_iedg);
915 S[x.k] = gamma_dis_inod;