Rheolef  7.2
an efficient C++ finite element environment
level_set.cc
Go to the documentation of this file.
1 #include "rheolef/level_set.h"
26 #include "rheolef/field.h"
27 namespace rheolef {
28 
29 // --------------------------------------------------------------------------------
30 // gestion de la numerotation locale
31 // --------------------------------------------------------------------------------
32 // TODO: move in reference_element
33 static
34 size_t
35 edge_t_iloc (size_t l, size_t m)
36 {
37  static const int edge_t_iloc_table [3][3] = {
38  {-1, 0, 2},
39  { 0,-1, 1},
40  { 2, 1,-1}};
41  return size_t(edge_t_iloc_table [l][m]);
42 }
43 static
44 size_t
45 edge_T_iloc (size_t l, size_t m)
46 {
47  static const int edge_T_iloc_table [4][4] = {
48  {-1, 0, 2, 3},
49  { 0,-1, 1, 4},
50  { 2, 1,-1, 5},
51  { 3, 4, 5,-1}};
52  return size_t(edge_T_iloc_table [l][m]);
53 }
54 static
55 size_t
56 face_T_iloc (size_t l, size_t m, size_t n)
57 {
58  static size_t face_T_iloc_table [4][4][4];
59  bool static initialized = false;
60  if (!initialized) {
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);
65  reference_element hat_K (reference_element::T);
66  for (size_t i_face = 0; i_face < hat_K.n_face(); i_face++) {
67  size_t p[3];
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;
75  }
76  }
77  return face_T_iloc_table [l][m][n];
78 }
79 // --------------------------------------------------------------------------------
80 // gestion de la precision
81 // --------------------------------------------------------------------------------
83 
84 template <class T>
85 static
86 bool
87 is_zero (const T& x) {
88  // TODO: control by level_set_option ?
89  return fabs(x) <= level_set_epsilon;
90 }
91 template <class T>
92 static
93 bool
94 have_same_sign (const T& x, const T& y) {
95  using namespace details;
96  return !is_zero(x) && !is_zero(y) && x*y > 0;
97 }
98 template <class T>
99 static
100 bool
101 have_opposite_sign (const T& x, const T& y) {
102  using namespace details;
103  return !is_zero(x) && !is_zero(y) && x*y < 0;
104 }
105 // --------------------------------------------------------------------------------
106 // 2D: fonctions locales : sur un seul element triangle
107 // --------------------------------------------------------------------------------
108 // appele lors du 1er passage qui liste les elements de la bande cas de dimension 2
109 template <class T>
110 static
111 bool
112 belongs_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;
115  // on rejette le triangle dans tous les sommets de meme signe :
116  if (is_zero(f[0]) && is_zero(f[1]) && is_zero(f[2])) return false;
117  // on rejette les triangles dont un sommet :
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;
121  return true;
122 }
123 // apellee lors du calcul des matrices M_K et A_K pour K dans la bande cas de dimension 2
124 template <class T>
125 static
126 size_t
127 isolated_vertex_t (const std::vector<T>& f)
128 {
129  using namespace details;
130  /* on retourne le sommet isole a chaque fois */
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; /* on peut retourner 2 de meme*/
135 
136  if (is_zero(f[1]) && have_opposite_sign(f[0],f[2])) return 0; /* on peut retourner 2 */
137  if (is_zero(f[2]) && have_opposite_sign(f[0],f[1])) return 0; /* on peut retourner 1 */
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;
140  return 0; /* f1 == 0 et f2 == 0 et f0 != 0 */
141 }
142 template <class T>
143 static
144 void
145 subcompute_matrix_t (
146  const std::vector<point_basic<T> >& x,
147  const std::vector<T>& f,
148  std::vector<size_t>& j,
149  point_basic<T>& a,
150  point_basic<T>& b,
151  T& S)
152 {
153  using namespace details;
154  j.resize (3);
155  j[0] = isolated_vertex_t (f);
156  j[1] = (j[0]+1) % 3;
157  j[2] = (j[1]+1) % 3;
158  // edge {j1,j2} has normal oriented as grad(f), in f>0 direction:
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]]);
162  // calcul des coordonnes d'intersection
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]];
165  S = sqrt(pow(a[0]-b[0],2)+pow(a[1]-b[1],2));
166  if (is_zero(f[j[1]]) && is_zero(f[j[2]])) {
167  S /= 2;
168  }
169 }
170 // --------------------------------------------------------------------------------
171 // 3D: fonctions locales : sur un seul element tetraedre
172 // --------------------------------------------------------------------------------
173 
174 class quadruplet {
175 public:
176  quadruplet (size_t a=0, size_t b=0, size_t c=0, size_t d=0) {
177  q[0]=a;
178  q[1]=b;
179  q[2]=c;
180  q[3]=d;
181  }
182  size_t operator[] (size_t i) const {
183  return q[i%4];
184  }
185  size_t& operator[] (size_t i) {
186  return q[i%4];
187  }
188  friend std::ostream& operator<< (std::ostream& out, const quadruplet& q) {
189  out << "((" << q[0] << "," << q[1] << "), (" << q[2] << "," << q[3] << "))";
190  return out;
191  }
192 protected:
193  size_t q[4];
194 };
195 // appele lors du 1er passage qui liste les elements de la bande cas de dimension 3
196 template <class T>
197 static
198 bool
199 belongs_to_band_T (const std::vector<T>& f)
200 {
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;
203  // cas ou 4 points s'annulent en dimension 3 est degenere
204  if (is_zero(f[0]) && is_zero(f[1]) && is_zero(f[2]) && is_zero(f[3])) return false;
205 
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;
210  // cas ou f s'annule sur 2 sommets et garde le meme signe sur les 2 autres sommets est exclu
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;
217  return true;
218 }
219 // apellee lors du calcul des matrices M_K et A_K pour T dans la bande cas de dimension
220 template <class T>
221 bool
222 intersection_is_quadrilateral_T (const std::vector<T>& f, quadruplet& q)
223 {
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);
227  return true;
228  }
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);
232  return true;
233  }
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);
237  return true;
238  }
239  return false;
240 }
241 // cas d'une intersection triangle:
242 template <class T>
243 static
244 size_t
245 isolated_vertex_T (const std::vector<T>& f)
246 {
247  using namespace details;
248  // cas ou l'intersection est un triangle
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;
252  // cas ou f s'annule sur un sommet et change de signe sur les 2 autres sommets
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;
266  // cas ou f s'annule en 2 sommets et change de signe sur les 2 autres sommets
267  if (is_zero(f[0]) && is_zero(f[1]) && have_opposite_sign(f[2],f[3])) return 2; // ou 3
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;
273 
274  // le triangle d'intersection est la face du tetradre ou f s'annule sur les 3 sommets
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;
278  return 1;
279 }
280 template <class T>
281 static
282 void
283 subcompute_matrix_triangle_T (
284  const std::vector<point_basic<T> >& x,
285  const std::vector<T>& f,
286  std::vector<size_t>& j,
287  point_basic<T>& a,
288  point_basic<T>& b,
289  point_basic<T>& c,
290  T& aire)
291 {
292  using namespace details;
293  j.resize(4);
294  j[0] = isolated_vertex_T (f);
295  j[1] = (j[0]+1) % 4;
296  j[2] = (j[1]+1) % 4;
297  j[3] = (j[2]+1) % 4;
298  // orient
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]]);
304 
305  /* calcul des coordonnees d'intersection */
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]];
309  aire = 0.5* norm (vect( b-a, c-a));
310  if (is_zero(f[j[1]]) && is_zero(f[j[2]]) && is_zero(f[j[3]])) {
311  aire /= 2;
312  }
313 }
314 template <class T>
315 static
316 void
317 subcompute_matrix_quadrilateral_T (
318  const std::vector<point_basic<T> >& x,
319  const std::vector<T>& f,
320  const quadruplet& q,
321  point_basic<T>& a,
322  point_basic<T>& b,
323  point_basic<T>& c,
324  point_basic<T>& d,
325  T& aire_Q)
326 {
327  // intersection:
328  // a = segment {x(q0) x(q2)} inter {f=0}
329  // b = segment {x(q1) x(q2)} inter {f=0}
330  // d = segment {x(q1) x(q3)} inter {f=0}
331  // c = segment {x(q0) x(q3)} inter {f=0}
332  // quadrilatere abdc = triangle(abd) union triangle(adc)
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]]);
337  /* calcul des coordonnees d'intersection */
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]];
342  aire_Q = 0.5*norm(vect(a-c,a-b)) + 0.5*norm(vect(d-c,d-b));
343 }
344 // --------------------------------------------------------------------------------
345 // level_set: compte the intersection mesh
346 // --------------------------------------------------------------------------------
348 
349 template <class T, class M>
350 void
352  const std::list<point_basic<T> >& gamma_node_list,
353  std::array<std::list<std::pair<element_type,size_t> >,
354  reference_element::max_variant> gamma_side_list,
355  const communicator& comm,
356  size_t d,
357  disarray<point_basic<T>, M>& gamma_node,
358  std::array<disarray<element_type,M>,
359  reference_element::max_variant>& gamma_side,
360  disarray<size_t,M>& sid_ie2bnd_ie)
361 {
363  // 1) nodes:
364  size_type nnod = gamma_node_list.size();
365  distributor gamma_node_ownership (distributor::decide, comm, nnod);
366  gamma_node.resize (gamma_node_ownership);
367  typename disarray<point_basic<T>, M>::iterator node_iter = gamma_node.begin();
368  for (typename std::list<point_basic<T> >::const_iterator
369  iter = gamma_node_list.begin(),
370  last = gamma_node_list.end();
371  iter != last; iter++, node_iter++) {
372  *node_iter = *iter;
373  }
374  // 2) sides:
376  size_type map_dim = d-1;
377  size_type order = 1;
378  size_type nsid = 0;
381  size_type nsidv = gamma_side_list [variant].size();
382  distributor gamma_sidv_ownership (distributor::decide, comm, nsidv);
383  nsid += nsidv;
384  element_type element_init (variant, order, alloc);
385  gamma_side[variant].resize (gamma_sidv_ownership, element_init);
386  typename disarray<element_type, M>::iterator side_iter = gamma_side[variant].begin();
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;
392  }
393  }
394  // 3) side2band correspondance
395  distributor gamma_sid_ownership (distributor::decide, comm, nsid);
396  const size_type undef = std::numeric_limits<size_t>::max();
397  sid_ie2bnd_ie.resize (gamma_sid_ownership, undef);
398  typename disarray<size_type, M>::iterator idx_iter = sid_ie2bnd_ie.begin();
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;
406  }
407  }
408 }
409 struct to_solve {
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) {}
413 };
414 template <class T, class M>
415 geo_basic<T,M>
417  const field_basic<T,M>& fh,
418  const level_set_option& opt,
419  std::vector<size_t>& bnd_dom_ie_list,
420  disarray<size_t,M>& sid_ie2bnd_ie)
421 {
422  using namespace std;
423  using namespace details;
425  level_set_epsilon = opt.epsilon; // set global variable
426  fh.dis_dof_update();
427  const geo_basic<T,M>& lambda = fh.get_geo();
428  const space_basic<T,M>& Xh = fh.get_space();
429  check_macro(lambda.order() == 1, "Only order=1 level set mesh supported");
430  check_macro(fh.get_approx() == "P1", "Only P1 level set function supported");
431  size_type order = 1;
432  std::vector<size_type> dis_idof;
433  std::vector<T> f;
434  size_type d = lambda.dimension();
435  size_type map_dim = d-1;
437  std::array<list<pair<element_type,size_type> >,
438  reference_element::max_variant> gamma_side_list;
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();
444 
445  distributor node_ownership = lambda.sizes().node_ownership;
446  size_type first_dis_inod = node_ownership.first_index();
447  disarray<size_type> marked_node (node_ownership, not_marked);
448  disarray<size_type> extern_node (node_ownership, not_marked);
449  set<size_type> ext_marked_node_idx;
450  list<to_solve> node_to_solve;
451 
452  distributor edge_ownership = lambda.sizes().ownership_by_dimension[1];
453  size_type first_dis_iedg = edge_ownership.first_index();
454  disarray<size_type> marked_edge (edge_ownership, not_marked);
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;
459 
460  distributor face_ownership = lambda.sizes().ownership_by_dimension[2];
461  size_type first_dis_ifac = face_ownership.first_index();
462  disarray<size_type> marked_face (face_ownership, not_marked);
463 
464  communicator comm = node_ownership.comm();
465  size_type my_proc = comm.rank();
466 
467  bnd_dom_ie_list.resize(0);
468 
469  // ------------------------------------------------------------
470  // 1) loop on lambda & build intersection sides
471  // ------------------------------------------------------------
472  for (size_type ie = 0, ne = lambda.size(), bnd_ie = 0; ie < ne; ie++) {
473  // ---------------------------------------------------------
474  // 1.1) fast check if there is an intersection:
475  // ---------------------------------------------------------
476  const geo_element& K = lambda [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]);
481  }
482  bool do_intersect = false;
483  switch (K.variant()) {
485  do_intersect = belongs_to_band_t (f);
486  break;
487  case reference_element::T: {
488  do_intersect = belongs_to_band_T (f);
489  break;
490  }
491  default :
492  error_macro("level set: element type `" << K.name() << "' not yet supported");
493  }
494  if (! do_intersect) continue;
495  bnd_dom_ie_list.push_back (ie);
496  // ---------------------------------------------------------
497  // 1.2) compute the intersection
498  // ---------------------------------------------------------
499  x.resize (dis_idof.size());
500  for (size_type loc_idof = 0, loc_ndof = dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
501  size_type loc_inod = loc_idof; // assume here that fh has an isoparametric approx
502  size_type dis_inod = K [loc_inod];
503  x [loc_idof] = lambda.dis_node(dis_inod);
504  }
505  element_type S (alloc);
506  switch (K.variant()) {
507  case reference_element::t: {
508  // ---------------------------------------------------------
509  // 1.2.a) triangle -> 1 edge
510  // ---------------------------------------------------------
511  point_basic<T> a, b;
512  T length;
513  subcompute_matrix_t (x, f, j, a, b, length);
514  if (is_zero(f[j[1]]) && is_zero(f[j[2]])) {
515  // ---------------------------------------------------------
516  // 1.2.a.1) edge {j1,j2} is included in the bbox mesh
517  // ---------------------------------------------------------
518  for (size_type k = 0; k < 2; k++) {
519  size_type dis_inod = K[j[k+1]];
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));
525  }
526  } else {
527  // dis_inod is owned by another partition jproc:
528  // if there is another neighbour element K' from the same partition jproc
529  // then jproc will insert dis_inod in gamma
530  // so, there is nothing to do here
531  // otherwise, dis_inod will be orphan in jproc and will be re-affected to my_proc
532  extern_node.dis_entry (dis_inod) = my_proc;
533  }
534  }
535  size_type loc_iedg = edge_t_iloc (j[1], j[2]);
536  size_type dis_iedg = K.edge (loc_iedg);
537  if (edge_ownership.is_owned (dis_iedg)) {
538  size_type iedg = dis_iedg - first_dis_iedg;
539  if (marked_edge [iedg] == not_marked) {
540  S.reset(reference_element::e, order);
541  marked_edge [iedg] = gamma_side_list [S.variant()].size();
542  size_type S_ige = marked_edge [iedg];
543  for (size_type k = 0; k < S.n_node(); k++) {
544  size_type dis_inod = K[j[k+1]];
545  if (node_ownership.is_owned (dis_inod)) {
546  size_type inod = dis_inod - first_dis_inod;
547  S[k] = marked_node [inod];
548  } else {
549  S[k] = undef;
550  node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
551  ext_marked_node_idx.insert (dis_inod);
552  }
553  }
554  gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
555  }
556  } else {
557  // intersection is an edge of lambda that is owned by another partition jproc
558  // the neighbour element K' across this e"dge belongs to the same partition jproc
559  // and will insert dis_iedg in gamma
560  // so, there is nothing to do here
561  }
562  } else {
563  // ---------------------------------------------------------
564  // 1.2.a.2) edge {j1,j2} is interior to the triangle
565  // ---------------------------------------------------------
566  S.reset(reference_element::e, order);
567  size_type S_ige = gamma_side_list [S.variant()].size ();
568  point_basic<T> xx[2] = {a,b};
569  for (size_type k = 0; k < 2; k++) {
570  if (! is_zero(f[j[k+1]]) && ! is_zero(f[j[0]])) {
571  // xk is inside edge {j0,j[k+1]} of triangle K:
572  size_type loc_iedg = edge_t_iloc (j[0], j[k+1]);
573  size_type dis_iedg = K.edge (loc_iedg);
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]);
579  }
580  S[k] = marked_edge [iedg];
581  } else {
582  S[k] = undef;
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]);
586  }
587  } else { // xk is at edge boundary: a node of the 2d mesh
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));
594  }
595  S[k] = marked_node [inod];
596  } else {
597  S[k] = undef;
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;
601  }
602  }
603  }
604  // S[0] == S[1] when is_zero(f[j[0]]) but f[j[0]] != 0, i.e. precision pbs
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));
607  }
608  break;
609  }
610  case reference_element::T: {
611  // ---------------------------------------------------------
612  // 1.2.b) tetrahedron -> 1 triangle or 1 quadrangle
613  // ---------------------------------------------------------
614  quadruplet q;
615  point_basic<T> a, b, c, d;
616  T aire;
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]])) {
620  // the full face {j1,j2,j3} is included in the surface mesh:
621  for (size_type k = 0; k < 3; k++) {
622  size_type dis_inod = K[j[k+1]];
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));
628  }
629  } else {
630  // dis_inod is owned by another partition jproc:
631  // if there is another neighbour element K' from the same partition jproc
632  // then jproc will insert dis_inod in gamma
633  // so, there is nothing to do here
634  // otherwise, dis_inod will be orphan in jproc and will be re-affected to my_proc
635  extern_node.dis_entry (dis_inod) = my_proc;
636  }
637  }
638  size_type loc_ifac = face_T_iloc (j[1], j[2], j[3]);
639  size_type dis_ifac = K.face (loc_ifac);
640  if (face_ownership.is_owned (dis_ifac)) {
641  size_type ifac = dis_ifac - first_dis_ifac;
642  if (marked_face [ifac] == not_marked) {
643  S.reset(reference_element::t, order);
644  marked_face [ifac] = gamma_side_list [S.variant()].size();
645  size_type S_ige = marked_face [ifac];
646  for (size_type k = 0; k < S.n_node(); k++) {
647  size_type dis_inod = K[j[k+1]];
648  if (node_ownership.is_owned (dis_inod)) {
649  size_type inod = dis_inod - first_dis_inod;
650  S[k] = marked_node [inod];
651  } else {
652  S[k] = undef;
653  node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
654  ext_marked_node_idx.insert (dis_inod);
655  }
656  }
657  gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
658  } else {
659  // the side will be inserted by the neighbour of K in another partition
660  // so, nothing to do
661  }
662  } else {
663  // intersection is a face of lambda that is owned by another partition jproc
664  // the neighbour element K' across this face belongs to the same partition jproc
665  // and will insert dis_ifac in gamma
666  // so, there is nothing to do here
667  }
668  } else {
669  // create the new face {j1,j2,j3} by intersections:
670  S.reset(reference_element::t, order);
671  size_type S_ige = gamma_side_list [S.variant()].size();
672  point_basic<T> xx[3] = {a,b,c};
673  for (size_type k = 0; k < 3; k++) {
674  if (! is_zero(f[j[k+1]]) && ! is_zero(f[j[0]])) {
675  // xk is inside edge {j0,j[k+1]} of triangle K:
676  size_type loc_iedg = edge_T_iloc (j[0], j[k+1]);
677  size_type dis_iedg = K.edge (loc_iedg);
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]);
683  }
684  S[k] = marked_edge [iedg];
685  } else {
686  S[k] = undef;
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]);
690  }
691  } else { // xk is at edge boundary: a node of the 2d mesh
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));
698  }
699  S[k] = marked_node [inod];
700  } else {
701  S[k] = undef;
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;
705  }
706  }
707  }
708  // S[0] == S[j] when is_zero(f[j[0]]) but f[j[0]] != 0, i.e. precision pbs
709  check_macro ((S[0] != S[1] || S[0] == undef) &&
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));
713  }
714  } else {
715  // create the new quadri face by intersections:
716  subcompute_matrix_quadrilateral_T (x, f, q, a, b, c, d, aire);
717  S.reset (reference_element::q, order);
718  size_type S_ige = gamma_side_list[reference_element::q].size();
719  size_type S1_ige = gamma_side_list[reference_element::t].size(); // or S1+S2=two 't'
720  size_type S2_ige = S1_ige + 1;
721  size_type iloc_S1[4] = {0, 1, 2, undef}; // the way to slip a q into 2 t
722  size_type iloc_S2[4] = {0, undef, 1, 2};
723  point_basic<T> xx[4] = {a,b,d,c};
724  size_type s[4] = {q[0],q[2],q[1],q[3]};
725  for (size_type k = 0; k < 4; k++) {
726  size_type k1 = (k+1) % 4;
727  if (! is_zero(f[s[k]]) && ! is_zero(f[s[k1]])) {
728  // xk is inside edge {j0,j[k+1]} of triangle K:
729  size_type loc_iedg = edge_T_iloc (s[k], s[k1]);
730  size_type dis_iedg = K.edge (loc_iedg);
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]);
736  }
737  S[k] = marked_edge [iedg];
738  } else {
739  S[k] = undef;
740  if (!opt.split_to_triangle) {
741  edge_to_solve.push_back (to_solve (reference_element::q, S_ige, k, dis_iedg));
742  } else {
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));
745  }
746  ext_marked_edge_idx.insert (dis_iedg);
747  extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
748  }
749  } else { // xk is at edge boundary: a node of the 2d mesh
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));
756  }
757  S[k] = marked_node [inod];
758  } else {
759  S[k] = undef;
760  if (!opt.split_to_triangle) {
761  node_to_solve.push_back (to_solve (reference_element::q, S_ige, k, dis_inod));
762  } else {
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));
765  }
766  ext_marked_node_idx.insert (dis_inod);
767  }
768  }
769  }
770  if (!opt.split_to_triangle) {
771  check_macro ((S[0] != S[1] || S[0] == undef) &&
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");
775  // S[0] == S[j] when is_zero(f[j[0]]) but f[j[0]] != 0, i.e. precision pbs
776  gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
777  } else {
778  // split quadri into 2 triangles
779  // one K -> two (S1,S2) faces: table element2face may return a pair of size_t
780  // but S1-> and S2->K only is required during assembly
781  element_type S1 (alloc);
782  S1.reset (reference_element::t, order);
783  for (size_type k = 0; k < 4; k++) {
784  if (iloc_S1[k] != undef) S1 [iloc_S1[k]] = S[k];
785  }
786  check_macro ((S1[0] != S1[1] || S1[0] == undef) &&
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));
790 
791  element_type S2 (alloc);
792  S2.reset (reference_element::t, order);
793  for (size_type k = 0; k < 4; k++) {
794  if (iloc_S2[k] != undef) S2 [iloc_S2[k]] = S[k];
795  }
796  check_macro ((S2[0] != S2[1] || S2[0] == undef) &&
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));
800  }
801  } // if-else
802  break;
803  }
804  default : {
805  error_macro("level set intersection: element not yet implemented: " << K.name());
806  }
807  }
808  bnd_ie++;
809  }
810  extern_node.dis_entry_assembly();
811  extern_edge.dis_entry_assembly();
812  // ------------------------------------------------------------
813  // 2) solve orphan nodes, if any
814  // ------------------------------------------------------------
815  // 2.1.a) re-affect orphan node to another process where gamma use it
816  distributor comm_ownership (comm.size(), comm, 1);
817  disarray<index_set,M> orphan_node (comm_ownership, index_set());
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;
820  // inod is orphan in this proc: not used for gamma (but used for lambda)
821  // re-affect it to a process that use it for gamma
822  size_type iproc = extern_node[inod];
823  size_type dis_inod = first_dis_inod + inod;
824  index_set dis_inod_set;
825  dis_inod_set.insert (dis_inod);
826  orphan_node.dis_entry (iproc) += dis_inod_set;
827  }
828  orphan_node.dis_entry_assembly();
829 
830  // 2.1.b) re-affect orphan edge to another process where gamma use it
831  // there could be orphan edge in 2d if lambda is a part of a regular mesh:
832  // => a boundary edge can belong to another proc. This is not yet handled
833  disarray<index_set,M> orphan_edge (comm_ownership, index_set());
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;
836  // iedg is orphan in this proc: not used for gamma (but used for lambda)
837  // re-affect it to a process that use it for gamma
838  size_type iproc = extern_edge[iedg].first;
839  size_type dis_iedg = first_dis_iedg + iedg;
840  index_set dis_iedg_set;
841  dis_iedg_set.insert (dis_iedg);
842  orphan_edge.dis_entry (iproc) += dis_iedg_set;
843  }
844  orphan_edge.dis_entry_assembly ();
845  check_macro (orphan_edge[0].size() == 0, "unexpected orphan edges");
846 
847  // 2.2) count total nodes used by gamma
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;
851  distributor gamma_node_ownership (distributor::decide, comm, gamma_nnod);
852  // 2.3) shift marked_node & marked_edge from gamma_inod local count to gamma_dis_inod one
853  size_type gamma_first_dis_inod = gamma_node_ownership.first_index();
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;
857  }
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;
861  }
862  // 2.4) append orphan node to regular one in gamma_node_list and set marked_node
863  for (index_set::const_iterator
864  iter = orphan_node[0].begin(),
865  last = orphan_node[0].end(); iter != last; ++iter) {
866  size_type dis_inod = *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));
869  }
870  marked_node.dis_entry_assembly();
871  marked_edge.dis_entry_assembly();
872  // ------------------------------------------------------------
873  // 3) convert lists to fixed size distributed arrays
874  // ------------------------------------------------------------
875  disarray<point_basic<T>, M> gamma_node;
876  std::array<disarray<element_type,M>, reference_element::max_variant> gamma_side;
877  gamma_list2disarray (gamma_node_list, gamma_side_list, comm, d, gamma_node, gamma_side, sid_ie2bnd_ie);
878  // ------------------------------------------------------------
879  // 4) replace inod to dis_inod in element lists
880  // ------------------------------------------------------------
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; // external node, will be solved later
887  S[loc_inod] += gamma_first_dis_inod;
888  }
889  }
890  }
891  // ----------------------------------------------------------------
892  // 5) solve intersection that are located on external edges & nodes
893  // ----------------------------------------------------------------
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");
900  size_type dis_inod = x.dis_i;
901  size_type iproc = node_ownership.find_owner(dis_inod);
902  size_type gamma_dis_inod = marked_node.dis_at (dis_inod);
903  S[x.k] = gamma_dis_inod;
904  }
905 
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");
912  size_type dis_iedg = x.dis_i;
913  size_type iproc = edge_ownership.find_owner(dis_iedg);
914  size_type gamma_dis_inod = marked_edge.dis_at (dis_iedg);
915  S[x.k] = gamma_dis_inod;
916  }
917  // ------------------------------------------------------------
918  // 6) convert intersection to geo
919  // ------------------------------------------------------------
920  geo_basic<T,M> gamma (lambda, gamma_node, gamma_side);
921  return gamma;
922 }
923 template <class T, class M>
924 geo_basic<T,M>
926  const field_basic<T,M>& fh,
927  const level_set_option& opt)
928 {
930  std::vector<size_type> bnd_dom_ie_list;
931  disarray<size_type,M> sid_ie2bnd_ie;
932  return level_set_internal (fh, opt, bnd_dom_ie_list, sid_ie2bnd_ie);
933 }
934 // ----------------------------------------------------------------------------
935 // instanciation in library
936 // ----------------------------------------------------------------------------
937 #define _RHEOLEF_instanciation(T,M) \
938 template 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>&); \
943 template geo_basic<T,M> level_set ( \
944  const field_basic<T,M>&, \
945  const level_set_option&);
946 
947 _RHEOLEF_instanciation(Float,sequential)
948 #ifdef _RHEOLEF_HAVE_MPI
950 #endif // _RHEOLEF_HAVE_MPI
951 
952 } // namespace
field::size_type size_type
Definition: branch.cc:430
see the Float page for the full documentation
see the disarray page for the full documentation
Definition: disarray.h:497
rep::base::iterator iterator
Definition: disarray.h:502
see the distributor page for the full documentation
Definition: distributor.h:69
size_type find_owner(size_type dis_i) const
find iproc associated to a global index dis_i: CPU=log(nproc)
Definition: distributor.cc:106
size_type size(size_type iproc) const
Definition: distributor.h:170
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
Definition: distributor.h:158
bool is_owned(size_type dis_i, size_type iproc) const
true when dis_i in [first_index(iproc):last_index(iproc)[
Definition: distributor.h:220
static const size_type decide
Definition: distributor.h:83
const communicator_type & comm() const
Definition: distributor.h:152
const space_type & get_space() const
Definition: field.h:270
const T & dis_dof(size_type dis_idof) const
Definition: field.cc:377
std::string get_approx() const
Definition: field.h:272
const geo_type & get_geo() const
Definition: field.h:271
void dis_dof_update(const SetOp &=SetOp()) const
Definition: field.h:763
see the geo_element page for the full documentation
Definition: geo_element.h:102
size_type edge(size_type i) const
Definition: geo_element.h:209
size_type face(size_type i) const
Definition: geo_element.h:210
reference_element::size_type size_type
Definition: geo_element.h:125
variant_type variant() const
Definition: geo_element.h:161
char name() const
Definition: geo_element.h:169
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
Definition: level_set.cc:347
size_t size_type
Definition: basis_get.cc:76
distributed
Definition: asr.cc:228
point_basic< T >
Definition: piola_fem.h:135
static Float level_set_epsilon
Definition: level_set.cc:82
#define error_macro(message)
Definition: dis_macros.h:49
Expr1::float_type T
Definition: field_expr.h:230
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)
Definition: level_set.cc:222
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:264
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)
Definition: level_set.cc:351
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition: catchmark.h:99
geo_basic< T, M > level_set_internal(const field_basic< T, M > &, const level_set_option &, std::vector< size_t > &, disarray< size_t, M > &)
Definition: level_set.cc:416
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
_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)
Definition: level_set.cc:925
Float gamma[][pmax+1]
Definition: cavity_dg.h:29
Definition: sphere.icc:25
Float epsilon
Expr1::memory_type M
Definition: vec_expr_v2.h:416