25#include "rheolef/level_set.h"
26#include "rheolef/field.h"
35edge_t_iloc (
size_t l,
size_t m)
37 static const int edge_t_iloc_table [3][3] = {
41 return size_t(edge_t_iloc_table [l][
m]);
45edge_T_iloc (
size_t l,
size_t m)
47 static const int edge_T_iloc_table [4][4] = {
52 return size_t(edge_T_iloc_table [l][
m]);
56face_T_iloc (
size_t l,
size_t m,
size_t n)
58 static size_t face_T_iloc_table [4][4][4];
59 bool static initialized =
false;
61 for (
size_t i = 0; i < 4; i++)
62 for (
size_t j = 0; j < 4; j++)
63 for (
size_t k = 0; k < 4; k++)
64 face_T_iloc_table [i][j][k] =
size_t(-1);
66 for (
size_t i_face = 0; i_face < hat_K.n_face(); i_face++) {
68 for (
size_t k = 0; k < 3; k++)
p[k] = hat_K.subgeo_local_vertex(2,i_face,k);
69 face_T_iloc_table [
p[0]][
p[1]][
p[2]] = i_face;
70 face_T_iloc_table [
p[0]][
p[2]][
p[1]] = i_face;
71 face_T_iloc_table [
p[1]][
p[0]][
p[2]] = i_face;
72 face_T_iloc_table [
p[1]][
p[2]][
p[0]] = i_face;
73 face_T_iloc_table [
p[2]][
p[0]][
p[1]] = i_face;
74 face_T_iloc_table [
p[2]][
p[1]][
p[0]] = i_face;
77 return face_T_iloc_table [l][
m][
n];
94have_same_sign (
const T& x,
const T& y) {
95 using namespace details;
96 return !is_zero(x) && !is_zero(y) && x*y > 0;
101have_opposite_sign (
const T& x,
const T& y) {
102 using namespace details;
103 return !is_zero(x) && !is_zero(y) && x*y < 0;
112belongs_to_band_t (
const std::vector<T>&
f) {
113 using namespace details;
114 if (have_same_sign(
f[0],
f[1]) && have_same_sign(
f[0],
f[2]))
return false;
116 if (is_zero(
f[0]) && is_zero(
f[1]) && is_zero(
f[2]))
return false;
118 if (is_zero(
f[0]) && have_same_sign(
f[1],
f[2]))
return false;
119 if (is_zero(
f[1]) && have_same_sign(
f[0],
f[2]))
return false;
120 if (is_zero(
f[2]) && have_same_sign(
f[0],
f[1]))
return false;
127isolated_vertex_t (
const std::vector<T>&
f)
129 using namespace details;
131 if (have_same_sign (
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]))
return 2;
132 if (have_opposite_sign(
f[0],
f[1]) && have_same_sign (
f[0],
f[2]))
return 1;
133 if (have_opposite_sign(
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]))
return 0;
134 if (is_zero(
f[0]) && have_opposite_sign(
f[1],
f[2]))
return 1;
136 if (is_zero(
f[1]) && have_opposite_sign(
f[0],
f[2]))
return 0;
137 if (is_zero(
f[2]) && have_opposite_sign(
f[0],
f[1]))
return 0;
138 if (is_zero(
f[0]) && is_zero(
f[1]) && ! is_zero(
f[2]))
return 2;
139 if (is_zero(
f[0]) && is_zero(
f[2]) && ! is_zero(
f[1]))
return 1;
147 const std::vector<T>&
f,
148 std::vector<size_t>& j,
153 using namespace details;
155 j[0] = isolated_vertex_t (
f);
159 if (! is_zero(
f[j[0]]) &&
f[j[0]] < 0) std::swap (j[1], j[2]);
160 T theta_1=
f[j[1]]/(
f[j[1]]-
f[j[0]]);
161 T theta_2=
f[j[2]]/(
f[j[2]]-
f[j[0]]);
163 a = theta_1*x[j[0]]+(1-theta_1)*x[j[1]];
164 b = theta_2*x[j[0]]+(1-theta_2)*x[j[2]];
166 if (is_zero(
f[j[1]]) && is_zero(
f[j[2]])) {
176 quadruplet (
size_t a=0,
size_t b=0,
size_t c=0,
size_t d=0) {
182 size_t operator[] (
size_t i)
const {
185 size_t& operator[] (
size_t i) {
188 friend std::ostream&
operator<< (std::ostream&
out,
const quadruplet& q) {
189 out <<
"((" << q[0] <<
"," << q[1] <<
"), (" << q[2] <<
"," << q[3] <<
"))";
199belongs_to_band_T (
const std::vector<T>&
f)
201 using namespace details;
202 if (have_same_sign(
f[0],
f[1]) && have_same_sign(
f[0],
f[2]) && have_same_sign(
f[2],
f[3]))
return false;
204 if (is_zero(
f[0]) && is_zero(
f[1]) && is_zero(
f[2]) && is_zero(
f[3]))
return false;
206 if (is_zero(
f[0]) && have_same_sign(
f[1],
f[2]) && have_same_sign(
f[1],
f[3]))
return false;
207 if (is_zero(
f[1]) && have_same_sign(
f[0],
f[2]) && have_same_sign(
f[0],
f[3]))
return false;
208 if (is_zero(
f[2]) && have_same_sign(
f[0],
f[1]) && have_same_sign(
f[0],
f[3]))
return false;
209 if (is_zero(
f[3]) && have_same_sign(
f[0],
f[1]) && have_same_sign(
f[0],
f[2]))
return false;
211 if (is_zero(
f[0]) && is_zero(
f[1]) && have_same_sign(
f[2],
f[3]))
return false;
212 if (is_zero(
f[0]) && is_zero(
f[2]) && have_same_sign(
f[1],
f[3]))
return false;
213 if (is_zero(
f[0]) && is_zero(
f[3]) && have_same_sign(
f[1],
f[2]))
return false;
214 if (is_zero(
f[1]) && is_zero(
f[2]) && have_same_sign(
f[0],
f[3]))
return false;
215 if (is_zero(
f[1]) && is_zero(
f[3]) && have_same_sign(
f[0],
f[2]))
return false;
216 if (is_zero(
f[2]) && is_zero(
f[3]) && have_same_sign(
f[1],
f[0]))
return false;
224 if (have_same_sign(
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]) && have_same_sign(
f[2],
f[3])) {
225 if (
f[0] > 0) q = quadruplet(0,1, 2,3);
226 else q = quadruplet(2,3, 0,1);
229 if (have_opposite_sign(
f[0],
f[1]) && have_same_sign(
f[0],
f[2]) && have_opposite_sign(
f[2],
f[3])) {
230 if (
f[0] < 0) q = quadruplet(0,2, 1,3);
231 else q = quadruplet(1,3, 0,2);
234 if (have_opposite_sign(
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]) && have_opposite_sign(
f[2],
f[3])) {
235 if (
f[0] > 0) q = quadruplet(0,3, 1,2);
236 else q = quadruplet(1,2, 0,3);
245isolated_vertex_T (
const std::vector<T>&
f)
247 using namespace details;
249 if (have_opposite_sign(
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]) && have_same_sign (
f[2],
f[3]))
return 0;
250 if (have_same_sign (
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]) && have_opposite_sign(
f[2],
f[3]))
return 2;
251 if (have_same_sign (
f[0],
f[1]) && have_same_sign (
f[0],
f[2]) && have_opposite_sign(
f[2],
f[3]))
return 3;
253 if (have_opposite_sign(
f[0],
f[1]) && have_same_sign(
f[0],
f[2]) && have_same_sign(
f[2],
f[3]))
return 1;
254 if (is_zero(
f[0]) && have_same_sign (
f[1],
f[2]) && have_opposite_sign(
f[1],
f[3]))
return 3;
255 if (is_zero(
f[0]) && have_opposite_sign(
f[1],
f[2]) && have_same_sign (
f[1],
f[3]))
return 2;
256 if (is_zero(
f[0]) && have_opposite_sign(
f[1],
f[2]) && have_opposite_sign(
f[1],
f[3]))
return 1;
257 if (is_zero(
f[1]) && have_opposite_sign(
f[0],
f[2]) && have_same_sign (
f[0],
f[3]))
return 2;
258 if (is_zero(
f[1]) && have_same_sign (
f[0],
f[2]) && have_opposite_sign(
f[0],
f[3]))
return 3;
259 if (is_zero(
f[1]) && have_opposite_sign(
f[0],
f[2]) && have_opposite_sign(
f[0],
f[3]))
return 0;
260 if (is_zero(
f[2]) && have_opposite_sign(
f[0],
f[1]) && have_same_sign (
f[0],
f[3]))
return 1;
261 if (is_zero(
f[2]) && have_same_sign (
f[0],
f[1]) && have_opposite_sign(
f[0],
f[3]))
return 3;
262 if (is_zero(
f[2]) && have_opposite_sign(
f[0],
f[1]) && have_opposite_sign(
f[0],
f[3]))
return 0;
263 if (is_zero(
f[3]) && have_opposite_sign(
f[0],
f[1]) && have_same_sign (
f[0],
f[2]))
return 1;
264 if (is_zero(
f[3]) && have_same_sign (
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]))
return 2;
265 if (is_zero(
f[3]) && have_opposite_sign(
f[0],
f[1]) && have_opposite_sign(
f[0],
f[2]))
return 0;
267 if (is_zero(
f[0]) && is_zero(
f[1]) && have_opposite_sign(
f[2],
f[3]))
return 2;
268 if (is_zero(
f[0]) && is_zero(
f[2]) && have_opposite_sign(
f[1],
f[3]))
return 1;
269 if (is_zero(
f[0]) && is_zero(
f[3]) && have_opposite_sign(
f[1],
f[2]))
return 1;
270 if (is_zero(
f[1]) && is_zero(
f[2]) && have_opposite_sign(
f[0],
f[3]))
return 0;
271 if (is_zero(
f[1]) && is_zero(
f[3]) && have_opposite_sign(
f[0],
f[2]))
return 0;
272 if (is_zero(
f[2]) && is_zero(
f[3]) && have_opposite_sign(
f[0],
f[1]))
return 0;
275 if (is_zero(
f[0]) && is_zero(
f[1]) && is_zero(
f[2]) && ! is_zero(
f[3]))
return 3;
276 if (is_zero(
f[0]) && is_zero(
f[1]) && is_zero(
f[3]) && ! is_zero(
f[2]))
return 2;
277 if (is_zero(
f[1]) && is_zero(
f[2]) && is_zero(
f[3]) && ! is_zero(
f[0]))
return 0;
283subcompute_matrix_triangle_T (
285 const std::vector<T>&
f,
286 std::vector<size_t>& j,
292 using namespace details;
294 j[0] = isolated_vertex_T (
f);
299 if (! is_zero(
f[j[0]]) && ((
f[j[0]] > 0 && j[0] % 2 == 0) || (
f[j[0]] < 0 && j[0] % 2 == 1)))
300 std::swap (j[1], j[2]);
301 T theta_1 =
f[j[1]]/(
f[j[1]]-
f[j[0]]);
302 T theta_2 =
f[j[2]]/(
f[j[2]]-
f[j[0]]);
303 T theta_3 =
f[j[3]]/(
f[j[3]]-
f[j[0]]);
306 a = theta_1*x[j[0]]+(1-theta_1)*x[j[1]];
307 b = theta_2*x[j[0]]+(1-theta_2)*x[j[2]];
308 c = theta_3*x[j[0]]+(1-theta_3)*x[j[3]];
310 if (is_zero(
f[j[1]]) && is_zero(
f[j[2]]) && is_zero(
f[j[3]])) {
317subcompute_matrix_quadrilateral_T (
319 const std::vector<T>&
f,
333 T theta_1 =
f[q[2]]/(
f[q[2]]-
f[q[0]]);
334 T theta_2 =
f[q[2]]/(
f[q[2]]-
f[q[1]]);
335 T theta_3 =
f[q[3]]/(
f[q[3]]-
f[q[0]]);
336 T theta_4 =
f[q[3]]/(
f[q[3]]-
f[q[1]]);
338 a = theta_1*x[q[0]]+(1-theta_1)*x[q[2]];
339 b = theta_2*x[q[1]]+(1-theta_2)*x[q[2]];
340 c = theta_3*x[q[0]]+(1-theta_3)*x[q[3]];
341 d = theta_4*x[q[1]]+(1-theta_4)*x[q[3]];
349template <
class T,
class M>
353 std::array<std::list<std::pair<element_type,size_t> >,
355 const communicator& comm,
366 gamma_node.resize (gamma_node_ownership);
369 iter = gamma_node_list.begin(),
370 last = gamma_node_list.end();
371 iter != last; iter++, node_iter++) {
385 gamma_side[
variant].resize (gamma_sidv_ownership, element_init);
387 for (
typename std::list<std::pair<element_type,size_type> >::const_iterator
388 iter = gamma_side_list[
variant].begin(),
389 last = gamma_side_list[
variant].end();
390 iter != last; iter++, side_iter++) {
391 *side_iter = (*iter).first;
396 const size_type undef = std::numeric_limits<size_t>::max();
397 sid_ie2bnd_ie.resize (gamma_sid_ownership, undef);
401 for (
typename std::list<std::pair<element_type,size_type> >::const_iterator
402 iter = gamma_side_list[
variant].begin(),
403 last = gamma_side_list[
variant].end();
404 iter != last; iter++, idx_iter++) {
405 *idx_iter = (*iter).second;
410 size_t variant, S_ige, k, dis_i;
411 to_solve (
size_t variant1,
size_t S_ige1=0,
size_t k1=0,
size_t dis_i1=0)
412 :
variant(variant1), S_ige(S_ige1), k(k1), dis_i(dis_i1) {}
414template <
class T,
class M>
419 std::vector<size_t>& bnd_dom_ie_list,
423 using namespace details;
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++) {
479 for (
size_type loc_idof = 0, loc_ndof =
dis_idof.size(); loc_idof < loc_ndof; 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);
500 for (
size_type loc_idof = 0, loc_ndof =
dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
513 subcompute_matrix_t (x,
f, j,
a,
b, length);
514 if (is_zero(
f[j[1]]) && is_zero(
f[j[2]])) {
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++) {
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]);
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]])) {
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++) {
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]);
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]);
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];
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));
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));
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();
819 if (!(extern_node [inod] != not_marked && marked_node [inod] == not_marked))
continue;
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;
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();
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++) {
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;
899 check_macro (S[x.k] == undef,
"external index already solved");
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;
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;
923template <
class T,
class M>
930 std::vector<size_type> bnd_dom_ie_list;
937#define _RHEOLEF_instanciation(T,M) \
938template geo_basic<T,M> level_set_internal ( \
939 const field_basic<T,M>&, \
940 const level_set_option&, \
941 std::vector<size_t>&, \
942 disarray<size_t,M>&); \
943template geo_basic<T,M> level_set ( \
944 const field_basic<T,M>&, \
945 const level_set_option&);
948#ifdef _RHEOLEF_HAVE_MPI
field::size_type size_type
see the Float page for the full documentation
see the communicator 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 T & dis_dof(size_type dis_idof) const
std::string get_approx() const
void dis_dof_update(const SetOp &=SetOp()) const
const geo_type & get_geo() const
const space_type & get_space() 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)
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
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)
This file is part of Rheolef.
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)
bool intersection_is_quadrilateral_T(const std::vector< T > &f, quadruplet &q)
std::ostream & operator<<(std::ostream &os, const catchmark &m)
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
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)
geo_basic< T, M > level_set_internal(const field_basic< T, M > &, const level_set_option &, std::vector< size_t > &, disarray< size_t, M > &)
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation