Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
level_set.cc
Go to the documentation of this file.
1
25#include "rheolef/level_set.h"
26#include "rheolef/field.h"
27namespace rheolef {
28
29// --------------------------------------------------------------------------------
30// gestion de la numerotation locale
31// --------------------------------------------------------------------------------
32// TODO: move in reference_element
33static
34size_t
35edge_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}
43static
44size_t
45edge_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}
54static
55size_t
56face_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// --------------------------------------------------------------------------------
82static Float level_set_epsilon = 100*std::numeric_limits<Float>::epsilon();
83
84template <class T>
85static
86bool
87is_zero (const T& x) {
88 // TODO: control by level_set_option ?
89 return fabs(x) <= level_set_epsilon;
90}
91template <class T>
92static
93bool
94have_same_sign (const T& x, const T& y) {
95 using namespace details;
96 return !is_zero(x) && !is_zero(y) && x*y > 0;
97}
98template <class T>
99static
100bool
101have_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
109template <class T>
110static
111bool
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;
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
124template <class T>
125static
126size_t
127isolated_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}
142template <class T>
143static
144void
145subcompute_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
174class quadruplet {
175public:
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 }
192protected:
193 size_t q[4];
194};
195// appele lors du 1er passage qui liste les elements de la bande cas de dimension 3
196template <class T>
197static
198bool
199belongs_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
220template <class T>
221bool
222intersection_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:
242template <class T>
243static
244size_t
245isolated_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}
280template <class T>
281static
282void
283subcompute_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}
314template <class T>
315static
316void
317subcompute_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
349template <class T, class M>
350void
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>,
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;
380 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
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();
400 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
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}
409struct 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};
414template <class T, class M>
415geo_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;
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()) {
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 }
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 // ------------------------------------------------------------
882 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
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}
923template <class T, class M>
924geo_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) \
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&);
946
948#ifdef _RHEOLEF_HAVE_MPI
949_RHEOLEF_instanciation(Float,distributed)
950#endif // _RHEOLEF_HAVE_MPI
951
952} // namespace
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
see the communicator 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)
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
Definition distributor.h:83
const communicator_type & comm() const
const T & dis_dof(size_type dis_idof) const
Definition field.cc:377
std::string get_approx() const
Definition field.h:272
void dis_dof_update(const SetOp &=SetOp()) const
Definition field.h:763
const geo_type & get_geo() const
Definition field.h:271
const space_type & get_space() const
Definition field.h:270
generic mesh with rerefence counting
Definition geo.h:1089
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
the finite element space
Definition space.h:382
geo_element_auto< heap_allocator< geo_element::size_type > > element_type
Definition level_set.cc:347
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)")
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)
Definition level_set.cc:351
bool intersection_is_quadrilateral_T(const std::vector< T > &f, quadruplet &q)
Definition level_set.cc:222
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
Definition point.h:264
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120
geo_basic< T, M > level_set(const field_basic< T, M > &fh, const level_set_option &opt)
Definition level_set.cc:925
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
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
STL namespace.
Definition cavity_dg.h:29
Definition sphere.icc:25
Expr1::memory_type M