1 #include "rheolef/level_set.h"
34 edge_t_iloc (
size_t l,
size_t m)
36 static const int edge_t_iloc_table [3][3] = {
40 return size_t(edge_t_iloc_table [l][
m]);
44 edge_T_iloc (
size_t l,
size_t m)
46 static const int edge_T_iloc_table [4][4] = {
51 return size_t(edge_T_iloc_table [l][
m]);
55 face_T_iloc (
size_t l,
size_t m,
size_t n)
57 static size_t face_T_iloc_table [4][4][4];
58 bool static initialized =
false;
60 for (
size_t i = 0; i < 4; i++)
61 for (
size_t j = 0; j < 4; j++)
62 for (
size_t k = 0; k < 4; k++)
63 face_T_iloc_table [i][j][k] =
size_t(-1);
65 for (
size_t i_face = 0; i_face < hat_K.n_face(); i_face++) {
67 for (
size_t k = 0; k < 3; k++)
p[k] = hat_K.subgeo_local_vertex(2,i_face,k);
68 face_T_iloc_table [
p[0]][
p[1]][
p[2]] = i_face;
69 face_T_iloc_table [
p[0]][
p[2]][
p[1]] = i_face;
70 face_T_iloc_table [
p[1]][
p[0]][
p[2]] = i_face;
71 face_T_iloc_table [
p[1]][
p[2]][
p[0]] = i_face;
72 face_T_iloc_table [
p[2]][
p[0]][
p[1]] = i_face;
73 face_T_iloc_table [
p[2]][
p[1]][
p[0]] = i_face;
76 return face_T_iloc_table [l][
m][
n];
86 is_zero (
const T& x) {
93 have_same_sign (
const T& x,
const T& y) {
94 using namespace details;
95 return !is_zero(x) && !is_zero(y) && x*y > 0;
100 have_opposite_sign (
const T& x,
const T& y) {
101 using namespace details;
102 return !is_zero(x) && !is_zero(y) && x*y < 0;
111 belongs_to_band_t (
const std::vector<T>&
f) {
112 using namespace details;
113 if (have_same_sign(
f[0],
f[1]) && have_same_sign(
f[0],
f[2]))
return false;
115 if (is_zero(
f[0]) && is_zero(
f[1]) && is_zero(
f[2]))
return false;
117 if (is_zero(
f[0]) && have_same_sign(
f[1],
f[2]))
return false;
118 if (is_zero(
f[1]) && have_same_sign(
f[0],
f[2]))
return false;
119 if (is_zero(
f[2]) && have_same_sign(
f[0],
f[1]))
return false;
126 isolated_vertex_t (
const std::vector<T>&
f)
128 using namespace details;
130 if (have_same_sign (
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]))
return 2;
131 if (have_opposite_sign(
f[0],
f[1]) && have_same_sign (
f[0],
f[2]))
return 1;
132 if (have_opposite_sign(
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]))
return 0;
133 if (is_zero(
f[0]) && have_opposite_sign(
f[1],
f[2]))
return 1;
135 if (is_zero(
f[1]) && have_opposite_sign(
f[0],
f[2]))
return 0;
136 if (is_zero(
f[2]) && have_opposite_sign(
f[0],
f[1]))
return 0;
137 if (is_zero(
f[0]) && is_zero(
f[1]) && ! is_zero(
f[2]))
return 2;
138 if (is_zero(
f[0]) && is_zero(
f[2]) && ! is_zero(
f[1]))
return 1;
144 subcompute_matrix_t (
146 const std::vector<T>&
f,
147 std::vector<size_t>& j,
152 using namespace details;
154 j[0] = isolated_vertex_t (
f);
158 if (! is_zero(
f[j[0]]) &&
f[j[0]] < 0) std::swap (j[1], j[2]);
159 T theta_1=
f[j[1]]/(
f[j[1]]-
f[j[0]]);
160 T theta_2=
f[j[2]]/(
f[j[2]]-
f[j[0]]);
162 a = theta_1*x[j[0]]+(1-theta_1)*x[j[1]];
163 b = theta_2*x[j[0]]+(1-theta_2)*x[j[2]];
165 if (is_zero(
f[j[1]]) && is_zero(
f[j[2]])) {
175 quadruplet (
size_t a=0,
size_t b=0,
size_t c=0,
size_t d=0) {
181 size_t operator[] (
size_t i)
const {
184 size_t& operator[] (
size_t i) {
187 friend std::ostream&
operator<< (std::ostream&
out,
const quadruplet& q) {
188 out <<
"((" << q[0] <<
"," << q[1] <<
"), (" << q[2] <<
"," << q[3] <<
"))";
198 belongs_to_band_T (
const std::vector<T>&
f)
200 using namespace details;
201 if (have_same_sign(
f[0],
f[1]) && have_same_sign(
f[0],
f[2]) && have_same_sign(
f[2],
f[3]))
return false;
203 if (is_zero(
f[0]) && is_zero(
f[1]) && is_zero(
f[2]) && is_zero(
f[3]))
return false;
205 if (is_zero(
f[0]) && have_same_sign(
f[1],
f[2]) && have_same_sign(
f[1],
f[3]))
return false;
206 if (is_zero(
f[1]) && have_same_sign(
f[0],
f[2]) && have_same_sign(
f[0],
f[3]))
return false;
207 if (is_zero(
f[2]) && have_same_sign(
f[0],
f[1]) && have_same_sign(
f[0],
f[3]))
return false;
208 if (is_zero(
f[3]) && have_same_sign(
f[0],
f[1]) && have_same_sign(
f[0],
f[2]))
return false;
210 if (is_zero(
f[0]) && is_zero(
f[1]) && have_same_sign(
f[2],
f[3]))
return false;
211 if (is_zero(
f[0]) && is_zero(
f[2]) && have_same_sign(
f[1],
f[3]))
return false;
212 if (is_zero(
f[0]) && is_zero(
f[3]) && have_same_sign(
f[1],
f[2]))
return false;
213 if (is_zero(
f[1]) && is_zero(
f[2]) && have_same_sign(
f[0],
f[3]))
return false;
214 if (is_zero(
f[1]) && is_zero(
f[3]) && have_same_sign(
f[0],
f[2]))
return false;
215 if (is_zero(
f[2]) && is_zero(
f[3]) && have_same_sign(
f[1],
f[0]))
return false;
223 if (have_same_sign(
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]) && have_same_sign(
f[2],
f[3])) {
224 if (
f[0] > 0) q = quadruplet(0,1, 2,3);
225 else q = quadruplet(2,3, 0,1);
228 if (have_opposite_sign(
f[0],
f[1]) && have_same_sign(
f[0],
f[2]) && have_opposite_sign(
f[2],
f[3])) {
229 if (
f[0] < 0) q = quadruplet(0,2, 1,3);
230 else q = quadruplet(1,3, 0,2);
233 if (have_opposite_sign(
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]) && have_opposite_sign(
f[2],
f[3])) {
234 if (
f[0] > 0) q = quadruplet(0,3, 1,2);
235 else q = quadruplet(1,2, 0,3);
244 isolated_vertex_T (
const std::vector<T>&
f)
246 using namespace details;
248 if (have_opposite_sign(
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]) && have_same_sign (
f[2],
f[3]))
return 0;
249 if (have_same_sign (
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]) && have_opposite_sign(
f[2],
f[3]))
return 2;
250 if (have_same_sign (
f[0],
f[1]) && have_same_sign (
f[0],
f[2]) && have_opposite_sign(
f[2],
f[3]))
return 3;
252 if (have_opposite_sign(
f[0],
f[1]) && have_same_sign(
f[0],
f[2]) && have_same_sign(
f[2],
f[3]))
return 1;
253 if (is_zero(
f[0]) && have_same_sign (
f[1],
f[2]) && have_opposite_sign(
f[1],
f[3]))
return 3;
254 if (is_zero(
f[0]) && have_opposite_sign(
f[1],
f[2]) && have_same_sign (
f[1],
f[3]))
return 2;
255 if (is_zero(
f[0]) && have_opposite_sign(
f[1],
f[2]) && have_opposite_sign(
f[1],
f[3]))
return 1;
256 if (is_zero(
f[1]) && have_opposite_sign(
f[0],
f[2]) && have_same_sign (
f[0],
f[3]))
return 2;
257 if (is_zero(
f[1]) && have_same_sign (
f[0],
f[2]) && have_opposite_sign(
f[0],
f[3]))
return 3;
258 if (is_zero(
f[1]) && have_opposite_sign(
f[0],
f[2]) && have_opposite_sign(
f[0],
f[3]))
return 0;
259 if (is_zero(
f[2]) && have_opposite_sign(
f[0],
f[1]) && have_same_sign (
f[0],
f[3]))
return 1;
260 if (is_zero(
f[2]) && have_same_sign (
f[0],
f[1]) && have_opposite_sign(
f[0],
f[3]))
return 3;
261 if (is_zero(
f[2]) && have_opposite_sign(
f[0],
f[1]) && have_opposite_sign(
f[0],
f[3]))
return 0;
262 if (is_zero(
f[3]) && have_opposite_sign(
f[0],
f[1]) && have_same_sign (
f[0],
f[2]))
return 1;
263 if (is_zero(
f[3]) && have_same_sign (
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]))
return 2;
264 if (is_zero(
f[3]) && have_opposite_sign(
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]))
return 0;
266 if (is_zero(
f[0]) && is_zero(
f[1]) && have_opposite_sign(
f[2],
f[3]))
return 2;
267 if (is_zero(
f[0]) && is_zero(
f[2]) && have_opposite_sign(
f[1],
f[3]))
return 1;
268 if (is_zero(
f[0]) && is_zero(
f[3]) && have_opposite_sign(
f[1],
f[2]))
return 1;
269 if (is_zero(
f[1]) && is_zero(
f[2]) && have_opposite_sign(
f[0],
f[3]))
return 0;
270 if (is_zero(
f[1]) && is_zero(
f[3]) && have_opposite_sign(
f[0],
f[2]))
return 0;
271 if (is_zero(
f[2]) && is_zero(
f[3]) && have_opposite_sign(
f[0],
f[1]))
return 0;
274 if (is_zero(
f[0]) && is_zero(
f[1]) && is_zero(
f[2]) && ! is_zero(
f[3]))
return 3;
275 if (is_zero(
f[0]) && is_zero(
f[1]) && is_zero(
f[3]) && ! is_zero(
f[2]))
return 2;
276 if (is_zero(
f[1]) && is_zero(
f[2]) && is_zero(
f[3]) && ! is_zero(
f[0]))
return 0;
282 subcompute_matrix_triangle_T (
284 const std::vector<T>&
f,
285 std::vector<size_t>& j,
291 using namespace details;
293 j[0] = isolated_vertex_T (
f);
298 if (! is_zero(
f[j[0]]) && ((
f[j[0]] > 0 && j[0] % 2 == 0) || (
f[j[0]] < 0 && j[0] % 2 == 1)))
299 std::swap (j[1], j[2]);
300 T theta_1 =
f[j[1]]/(
f[j[1]]-
f[j[0]]);
301 T theta_2 =
f[j[2]]/(
f[j[2]]-
f[j[0]]);
302 T theta_3 =
f[j[3]]/(
f[j[3]]-
f[j[0]]);
305 a = theta_1*x[j[0]]+(1-theta_1)*x[j[1]];
306 b = theta_2*x[j[0]]+(1-theta_2)*x[j[2]];
307 c = theta_3*x[j[0]]+(1-theta_3)*x[j[3]];
309 if (is_zero(
f[j[1]]) && is_zero(
f[j[2]]) && is_zero(
f[j[3]])) {
316 subcompute_matrix_quadrilateral_T (
318 const std::vector<T>&
f,
332 T theta_1 =
f[q[2]]/(
f[q[2]]-
f[q[0]]);
333 T theta_2 =
f[q[2]]/(
f[q[2]]-
f[q[1]]);
334 T theta_3 =
f[q[3]]/(
f[q[3]]-
f[q[0]]);
335 T theta_4 =
f[q[3]]/(
f[q[3]]-
f[q[1]]);
337 a = theta_1*x[q[0]]+(1-theta_1)*x[q[2]];
338 b = theta_2*x[q[1]]+(1-theta_2)*x[q[2]];
339 c = theta_3*x[q[0]]+(1-theta_3)*x[q[3]];
340 d = theta_4*x[q[1]]+(1-theta_4)*x[q[3]];
348 template <
class T,
class M>
352 std::array<std::list<std::pair<element_type,size_t> >,
354 const communicator& comm,
365 gamma_node.resize (gamma_node_ownership);
368 iter = gamma_node_list.begin(),
369 last = gamma_node_list.end();
370 iter != last; iter++, node_iter++) {
384 gamma_side[
variant].resize (gamma_sidv_ownership, element_init);
386 for (
typename std::list<std::pair<element_type,size_type> >::
const_iterator
387 iter = gamma_side_list[
variant].begin(),
388 last = gamma_side_list[
variant].end();
389 iter != last; iter++, side_iter++) {
390 *side_iter = (*iter).first;
395 const size_type undef = std::numeric_limits<size_t>::max();
396 sid_ie2bnd_ie.resize (gamma_sid_ownership, undef);
400 for (
typename std::list<std::pair<element_type,size_type> >::
const_iterator
401 iter = gamma_side_list[
variant].begin(),
402 last = gamma_side_list[
variant].end();
403 iter != last; iter++, idx_iter++) {
404 *idx_iter = (*iter).second;
409 size_t variant, S_ige, k, dis_i;
410 to_solve (
size_t variant1,
size_t S_ige1=0,
size_t k1=0,
size_t dis_i1=0)
411 :
variant(variant1), S_ige(S_ige1), k(k1), dis_i(dis_i1) {}
413 template <
class T,
class M>
418 std::vector<size_t>& bnd_dom_ie_list,
422 using namespace details;
436 std::array<list<pair<element_type,size_type> >,
438 list<point_basic<T> > gamma_node_list;
439 std::vector<size_type> j(
d+1);
440 std::vector<point_basic<T> > x(
d+1);
441 const size_type not_marked = numeric_limits<size_t>::max();
442 const size_type undef = numeric_limits<size_t>::max();
448 set<size_type> ext_marked_node_idx;
449 list<to_solve> node_to_solve;
455 extern_edge (edge_ownership, std::make_pair(not_marked,
point_basic<T>()));
456 set<size_type> ext_marked_edge_idx;
457 list<to_solve> edge_to_solve;
463 communicator comm = node_ownership.
comm();
466 bnd_dom_ie_list.resize(0);
471 for (
size_type ie = 0, ne =
lambda.size(), bnd_ie = 0; ie < ne; ie++) {
478 for (
size_type loc_idof = 0, loc_ndof =
dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
481 bool do_intersect =
false;
484 do_intersect = belongs_to_band_t (
f);
487 do_intersect = belongs_to_band_T (
f);
491 error_macro(
"level set: element type `" << K.
name() <<
"' not yet supported");
493 if (! do_intersect)
continue;
494 bnd_dom_ie_list.push_back (ie);
499 for (
size_type loc_idof = 0, loc_ndof =
dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
512 subcompute_matrix_t (x,
f, j,
a,
b, length);
513 if (is_zero(
f[j[1]]) && is_zero(
f[j[2]])) {
521 if (marked_node [inod] == not_marked) {
522 marked_node [inod] = gamma_node_list.size();
523 gamma_node_list.push_back (
lambda.node(inod));
531 extern_node.dis_entry (
dis_inod) = my_proc;
534 size_type loc_iedg = edge_t_iloc (j[1], j[2]);
536 if (edge_ownership.
is_owned (dis_iedg)) {
537 size_type iedg = dis_iedg - first_dis_iedg;
538 if (marked_edge [iedg] == not_marked) {
540 marked_edge [iedg] = gamma_side_list [S.variant()].size();
542 for (
size_type k = 0; k < S.n_node(); k++) {
546 S[k] = marked_node [inod];
549 node_to_solve.push_back (to_solve (S.variant(), S_ige, k,
dis_inod));
550 ext_marked_node_idx.insert (
dis_inod);
553 gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
566 size_type S_ige = gamma_side_list [S.variant()].size ();
569 if (! is_zero(
f[j[k+1]]) && ! is_zero(
f[j[0]])) {
571 size_type loc_iedg = edge_t_iloc (j[0], j[k+1]);
573 if (edge_ownership.
is_owned (dis_iedg)) {
574 size_type iedg = dis_iedg - first_dis_iedg;
575 if (marked_edge [iedg] == not_marked) {
576 marked_edge [iedg] = gamma_node_list.size();
577 gamma_node_list.push_back (xx[k]);
579 S[k] = marked_edge [iedg];
582 edge_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_iedg));
583 ext_marked_edge_idx.insert (dis_iedg);
584 extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
590 if (marked_node [inod] == not_marked) {
591 marked_node [inod] = gamma_node_list.size();
592 gamma_node_list.push_back (
lambda.node(inod));
594 S[k] = marked_node [inod];
597 node_to_solve.push_back (to_solve (S.variant(), S_ige, k,
dis_inod));
598 ext_marked_node_idx.insert (
dis_inod);
599 extern_node.dis_entry (
dis_inod) = my_proc;
604 check_macro (S[0] != S[1] || S[0] == undef,
"degenerate 2d intersection");
605 gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
617 subcompute_matrix_triangle_T (x,
f, j,
a,
b,
c, aire);
618 if (is_zero(
f[j[1]]) && is_zero(
f[j[2]]) && is_zero(
f[j[3]])) {
624 if (marked_node [inod] == not_marked) {
625 marked_node [inod] = gamma_node_list.size();
626 gamma_node_list.push_back (
lambda.node(inod));
634 extern_node.dis_entry (
dis_inod) = my_proc;
637 size_type loc_ifac = face_T_iloc (j[1], j[2], j[3]);
639 if (face_ownership.
is_owned (dis_ifac)) {
640 size_type ifac = dis_ifac - first_dis_ifac;
641 if (marked_face [ifac] == not_marked) {
643 marked_face [ifac] = gamma_side_list [S.variant()].size();
645 for (
size_type k = 0; k < S.n_node(); k++) {
649 S[k] = marked_node [inod];
652 node_to_solve.push_back (to_solve (S.variant(), S_ige, k,
dis_inod));
653 ext_marked_node_idx.insert (
dis_inod);
656 gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
670 size_type S_ige = gamma_side_list [S.variant()].size();
673 if (! is_zero(
f[j[k+1]]) && ! is_zero(
f[j[0]])) {
675 size_type loc_iedg = edge_T_iloc (j[0], j[k+1]);
677 if (edge_ownership.
is_owned (dis_iedg)) {
678 size_type iedg = dis_iedg - first_dis_iedg;
679 if (marked_edge [iedg] == not_marked) {
680 marked_edge [iedg] = gamma_node_list.size();
681 gamma_node_list.push_back (xx[k]);
683 S[k] = marked_edge [iedg];
686 edge_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_iedg));
687 ext_marked_edge_idx.insert (dis_iedg);
688 extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
694 if (marked_node [inod] == not_marked) {
695 marked_node [inod] = gamma_node_list.size();
696 gamma_node_list.push_back (
lambda.node(inod));
698 S[k] = marked_node [inod];
701 node_to_solve.push_back (to_solve (S.variant(), S_ige, k,
dis_inod));
702 ext_marked_node_idx.insert (
dis_inod);
703 extern_node.dis_entry (
dis_inod) = my_proc;
709 (S[1] != S[2] || S[1] == undef) &&
710 (S[2] != S[0] || S[2] == undef),
"degenerate 3d intersection");
711 gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
715 subcompute_matrix_quadrilateral_T (x,
f, q,
a,
b,
c,
d, aire);
726 if (! is_zero(
f[s[k]]) && ! is_zero(
f[s[k1]])) {
728 size_type loc_iedg = edge_T_iloc (s[k], s[k1]);
730 if (edge_ownership.
is_owned (dis_iedg)) {
731 size_type iedg = dis_iedg - first_dis_iedg;
732 if (marked_edge [iedg] == not_marked) {
733 marked_edge [iedg] = gamma_node_list.size();
734 gamma_node_list.push_back (xx[k]);
736 S[k] = marked_edge [iedg];
742 if (iloc_S1[k] != undef) edge_to_solve.push_back (to_solve (
reference_element::t, S1_ige, iloc_S1[k], dis_iedg));
743 if (iloc_S2[k] != undef) edge_to_solve.push_back (to_solve (
reference_element::t, S2_ige, iloc_S2[k], dis_iedg));
745 ext_marked_edge_idx.insert (dis_iedg);
746 extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
752 if (marked_node [inod] == not_marked) {
753 marked_node [inod] = gamma_node_list.size();
754 gamma_node_list.push_back (
lambda.node(inod));
756 S[k] = marked_node [inod];
765 ext_marked_node_idx.insert (
dis_inod);
771 (S[1] != S[2] || S[1] == undef) &&
772 (S[2] != S[3] || S[2] == undef) &&
773 (S[3] != S[0] || S[3] == undef),
"degenerate 3d intersection");
775 gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
783 if (iloc_S1[k] != undef) S1 [iloc_S1[k]] = S[k];
786 (S1[1] != S1[2] || S1[1] == undef) &&
787 (S1[2] != S1[0] || S1[2] == undef),
"degenerate 3d intersection");
788 gamma_side_list [S1.variant()].push_back (make_pair(S1,bnd_ie));
793 if (iloc_S2[k] != undef) S2 [iloc_S2[k]] = S[k];
796 (S2[1] != S2[2] || S2[1] == undef) &&
797 (S2[2] != S2[0] || S2[2] == undef),
"degenerate 3d intersection");
798 gamma_side_list [S2.variant()].push_back (make_pair(S2,bnd_ie));
804 error_macro(
"level set intersection: element not yet implemented: " << K.
name());
809 extern_node.dis_entry_assembly();
810 extern_edge.dis_entry_assembly();
818 if (!(extern_node [inod] != not_marked && marked_node [inod] == not_marked))
continue;
825 orphan_node.dis_entry (iproc) += dis_inod_set;
827 orphan_node.dis_entry_assembly();
833 for (
size_type iedg = 0, nedg = edge_ownership.
size(); iedg < nedg; iedg++) {
834 if (!(extern_edge [iedg].first != not_marked && marked_edge [iedg] == not_marked))
continue;
837 size_type iproc = extern_edge[iedg].first;
838 size_type dis_iedg = first_dis_iedg + iedg;
840 dis_iedg_set.
insert (dis_iedg);
841 orphan_edge.dis_entry (iproc) += dis_iedg_set;
843 orphan_edge.dis_entry_assembly ();
844 check_macro (orphan_edge[0].size() == 0,
"unexpected orphan edges");
847 size_type orphan_gamma_nnod = orphan_node[0].size();
848 size_type regular_gamma_nnod = gamma_node_list.size();
849 size_type gamma_nnod = regular_gamma_nnod + orphan_gamma_nnod;
854 if (marked_node [inod] == not_marked)
continue;
855 marked_node [inod] += gamma_first_dis_inod;
857 for (
size_type iedg = 0, nedg = edge_ownership.
size(); iedg < nedg; iedg++) {
858 if (marked_edge [iedg] == not_marked)
continue;
859 marked_edge [iedg] += gamma_first_dis_inod;
862 for (index_set::const_iterator
863 iter = orphan_node[0].begin(),
864 last = orphan_node[0].end(); iter != last; ++iter) {
866 marked_node.dis_entry(
dis_inod) = gamma_first_dis_inod + gamma_node_list.size();
869 marked_node.dis_entry_assembly();
870 marked_edge.dis_entry_assembly();
876 gamma_list2disarray (gamma_node_list, gamma_side_list, comm,
d, gamma_node, gamma_side, sid_ie2bnd_ie);
882 for (
size_type ige = 0, nge = gamma_side[
variant].size(); ige < nge; ige++) {
884 for (
size_type loc_inod = 0, loc_nnod = S.n_node(); loc_inod < loc_nnod; loc_inod++) {
885 if (S[loc_inod] == undef)
continue;
886 S[loc_inod] += gamma_first_dis_inod;
893 marked_node.set_dis_indexes (ext_marked_node_idx);
894 for (list<to_solve>::const_iterator iter = node_to_solve.begin(),
895 last = node_to_solve.end(); iter != last; iter++) {
896 const to_solve& x = *iter;
898 check_macro (S[x.k] == undef,
"external index already solved");
902 S[x.k] = gamma_dis_inod;
905 marked_edge.set_dis_indexes (ext_marked_edge_idx);
906 for (list<to_solve>::const_iterator iter = edge_to_solve.begin(),
907 last = edge_to_solve.end(); iter != last; iter++) {
908 const to_solve& x = *iter;
910 check_macro (S[x.k] == undef,
"external index already solved");
913 size_type gamma_dis_inod = marked_edge.dis_at (dis_iedg);
914 S[x.k] = gamma_dis_inod;
922 template <
class T,
class M>
929 std::vector<size_type> bnd_dom_ie_list;
936 #define _RHEOLEF_instanciation(T,M) \
937 template geo_basic<T,M> level_set_internal ( \
938 const field_basic<T,M>&, \
939 const level_set_option&, \
940 std::vector<size_t>&, \
941 disarray<size_t,M>&); \
942 template geo_basic<T,M> level_set ( \
943 const field_basic<T,M>&, \
944 const level_set_option&);
947 #ifdef _RHEOLEF_HAVE_MPI
field::size_type size_type
see the Float page for the full documentation
see the disarray page for the full documentation
rep::base::iterator iterator
see the distributor page for the full documentation
size_type find_owner(size_type dis_i) const
find iproc associated to a global index dis_i: CPU=log(nproc)
size_type size(size_type iproc) const
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
bool is_owned(size_type dis_i, size_type iproc) const
true when dis_i in [first_index(iproc):last_index(iproc)[
static const size_type decide
const communicator_type & comm() const
const space_type & get_space() const
const T & dis_dof(size_type dis_idof) const
std::string get_approx() const
const geo_type & get_geo() const
void dis_dof_update() const
see the geo_element page for the full documentation
size_type edge(size_type i) const
size_type face(size_type i) const
reference_element::size_type size_type
variant_type variant() const
void insert(size_type dis_i)
static const variant_type q
static const variant_type e
static const variant_type max_variant
static variant_type last_variant_by_dimension(size_type dim)
static variant_type first_variant_by_dimension(size_type dim)
static const variant_type T
static const variant_type t
geo_element_auto< heap_allocator< geo_element::size_type > > element_type
static Float level_set_epsilon
#define error_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
This file is part of Rheolef.
bool intersection_is_quadrilateral_T(const std::vector< T > &f, quadruplet &q)
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
void gamma_list2disarray(const std::list< point_basic< T > > &gamma_node_list, std::array< std::list< std::pair< element_type, size_t > >, reference_element::max_variant > gamma_side_list, const communicator &comm, size_t d, disarray< point_basic< T >, M > &gamma_node, std::array< disarray< element_type, M >, reference_element::max_variant > &gamma_side, disarray< size_t, M > &sid_ie2bnd_ie)
std::ostream & operator<<(std::ostream &os, const catchmark &m)
geo_basic< T, M > level_set_internal(const field_basic< T, M > &, const level_set_option &, std::vector< size_t > &, disarray< size_t, M > &)
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
geo_basic< T, M > level_set(const field_basic< T, M > &fh, const level_set_option &opt)