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