2
/// This file is part of Rheolef.
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
6
/// Rheolef is free software; you can redistribute it and/or modify
7
/// it under the terms of the GNU General Public License as published by
8
/// the Free Software Foundation; either version 2 of the License, or
9
/// (at your option) any later version.
11
/// Rheolef is distributed in the hope that it will be useful,
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14
/// GNU General Public License for more details.
16
/// You should have received a copy of the GNU General Public License
17
/// along with Rheolef; if not, write to the Free Software
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20
/// =========================================================================
21
/// Banded level set routines
23
/// Authors: Lara Aborm, Jocelyn Etienne, Pierre Saramito
25
#include "rheolef/level_set.h"
28
// --------------------------------------------------------------------------------
29
// gestion de la numerotation locale
30
// --------------------------------------------------------------------------------
31
// TODO: move in reference_element
34
edge_t_iloc (size_t l, size_t m)
36
static const int edge_t_iloc_table [3][3] = {
40
return size_t(edge_t_iloc_table [l][m]);
44
edge_T_iloc (size_t l, size_t m)
46
static const int edge_T_iloc_table [4][4] = {
51
return size_t(edge_T_iloc_table [l][m]);
55
face_T_iloc (size_t l, size_t m, size_t n)
57
static size_t face_T_iloc_table [4][4][4];
58
bool static initialized = false;
60
for (size_t i = 0; i < 4; i++)
61
for (size_t j = 0; j < 4; j++)
62
for (size_t k = 0; k < 4; k++)
63
face_T_iloc_table [i][j][k] = size_t(-1);
64
reference_element hat_K (reference_element::T);
65
for (size_t i_face = 0; i_face < hat_K.n_face(); i_face++) {
67
for (size_t k = 0; k < 3; k++) p[k] = hat_K.subgeo_local_vertex(2,i_face,k);
68
face_T_iloc_table [p[0]][p[1]][p[2]] = i_face;
69
face_T_iloc_table [p[0]][p[2]][p[1]] = i_face;
70
face_T_iloc_table [p[1]][p[0]][p[2]] = i_face;
71
face_T_iloc_table [p[1]][p[2]][p[0]] = i_face;
72
face_T_iloc_table [p[2]][p[0]][p[1]] = i_face;
73
face_T_iloc_table [p[2]][p[1]][p[0]] = i_face;
76
return face_T_iloc_table [l][m][n];
78
// --------------------------------------------------------------------------------
79
// gestion de la precision
80
// --------------------------------------------------------------------------------
81
static Float level_set_epsilon = 100*std::numeric_limits<Float>::epsilon();
86
is_zero (const T& x) {
87
// TODO: control by level_set_option_type ?
88
return fabs(x) <= level_set_epsilon;
93
have_same_sign (const T& x, const T& y) {
94
using namespace details;
95
return !is_zero(x) && !is_zero(y) && x*y > 0;
100
have_opposite_sign (const T& x, const T& y) {
101
using namespace details;
102
return !is_zero(x) && !is_zero(y) && x*y < 0;
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
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;
122
// apellee lors du calcul des matrices M_K et A_K pour K dans la bande cas de dimension 2
126
isolated_vertex_t (const std::vector<T>& f)
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*/
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 */
144
subcompute_matrix_t (
145
const std::vector<point_basic<T> >& x,
146
const std::vector<T>& f,
147
std::vector<size_t>& j,
152
using namespace details;
154
j[0] = isolated_vertex_t (f);
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íntersection
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]])) {
169
// --------------------------------------------------------------------------------
170
// 3D: fonctions locales : sur un seul element tetraedre
171
// --------------------------------------------------------------------------------
175
quadruplet (size_t a=0, size_t b=0, size_t c=0, size_t d=0) {
181
size_t operator[] (size_t i) const {
184
size_t& operator[] (size_t i) {
187
friend std::ostream& operator<< (std::ostream& out, const quadruplet& q) {
188
out << "((" << q[0] << "," << q[1] << "), (" << q[2] << "," << q[3] << "))";
194
// appele lors du 1er passage qui liste les elements de la bande cas de dimension 3
198
belongs_to_band_T (const std::vector<T>& f)
200
using namespace details;
201
if (have_same_sign(f[0],f[1]) && have_same_sign(f[0],f[2]) && have_same_sign(f[2],f[3])) return false;
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;
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;
218
// apellee lors du calcul des matrices M_K et A_K pour T dans la bande cas de dimension
221
intersection_is_quadrilateral_T (const std::vector<T>& f, quadruplet& q)
223
if (have_same_sign(f[0],f[1]) && have_opposite_sign(f[0],f[2]) && have_same_sign(f[2],f[3])) {
224
if (f[0] > 0) q = quadruplet(0,1, 2,3);
225
else q = quadruplet(2,3, 0,1);
228
if (have_opposite_sign(f[0],f[1]) && have_same_sign(f[0],f[2]) && have_opposite_sign(f[2],f[3])) {
229
if (f[0] < 0) q = quadruplet(0,2, 1,3);
230
else q = quadruplet(1,3, 0,2);
233
if (have_opposite_sign(f[0],f[1]) && have_opposite_sign(f[0],f[2]) && have_opposite_sign(f[2],f[3])) {
234
if (f[0] > 0) q = quadruplet(0,3, 1,2);
235
else q = quadruplet(1,2, 0,3);
240
// cas d'une intersection triangle:
244
isolated_vertex_T (const std::vector<T>& f)
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;
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;
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,
291
using namespace details;
293
j[0] = isolated_vertex_T (f);
298
if (! is_zero(f[j[0]]) && ((f[j[0]] > 0 && j[0] % 2 == 0) || (f[j[0]] < 0 && j[0] % 2 == 1)))
299
std::swap (j[1], j[2]);
300
T theta_1 = f[j[1]]/(f[j[1]]-f[j[0]]);
301
T theta_2 = f[j[2]]/(f[j[2]]-f[j[0]]);
302
T theta_3 = f[j[3]]/(f[j[3]]-f[j[0]]);
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]])) {
316
subcompute_matrix_quadrilateral_T (
317
const std::vector<point_basic<T> >& x,
318
const std::vector<T>& f,
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));
343
// --------------------------------------------------------------------------------
344
// level_set: compte the intersection mesh
345
// --------------------------------------------------------------------------------
346
typedef geo_element_auto<heap_allocator<geo_element::size_type> > element_type;
348
template <class T, class M>
351
const std::list<point_basic<T> >& gamma_node_list,
352
boost::array<std::list<std::pair<element_type,size_t> >,
353
reference_element::max_variant> gamma_side_list,
354
const communicator& comm,
356
array<point_basic<T>, M>& gamma_node,
357
boost::array<array<element_type,M>,
358
reference_element::max_variant>& gamma_side,
359
array<size_t,M>& sid_ie2bnd_ie)
361
typedef geo_element::size_type size_type;
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 array<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++) {
374
heap_allocator<size_type> alloc;
375
size_type map_dim = d-1;
378
for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
379
variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
380
size_type nsidv = gamma_side_list [variant].size();
381
distributor gamma_sidv_ownership (distributor::decide, comm, nsidv);
383
element_type element_init (variant, order, alloc);
384
gamma_side[variant].resize (gamma_sidv_ownership, element_init);
385
typename array<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;
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 array<size_type, M>::iterator idx_iter = sid_ie2bnd_ie.begin();
398
for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
399
variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
400
for (typename std::list<std::pair<element_type,size_type> >::const_iterator
401
iter = gamma_side_list[variant].begin(),
402
last = gamma_side_list[variant].end();
403
iter != last; iter++, idx_iter++) {
404
*idx_iter = (*iter).second;
409
template <class T, class M>
412
const array<point_basic<T>, M>& gamma_node,
413
const boost::array<array<element_type,M>, reference_element::max_variant>& gamma_side,
417
typedef geo_element::size_type size_type;
418
dout << "mesh" << endl
421
<< " dimension " << d << endl
422
<< " nodes " << gamma_node.dis_size() << endl;
424
if (gamma_side[reference_element::t].dis_size()!= 0) {
425
dout << " triangles " << gamma_side[reference_element::t].dis_size() << endl;
427
if (gamma_side[reference_element::q].dis_size()!= 0) {
428
dout << " quadrangles " << gamma_side[reference_element::q].dis_size() << endl;
431
dout << " edges " << gamma_side[reference_element::e].dis_size() << endl;
433
dout << "end header" << endl
435
<< setprecision(7); // avoid non-regr test to depend upon mach prec
436
gamma_node.put_values (dout, _point_put<T>(d));
438
size_type map_dim = d-1;
439
for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
440
variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
441
dout << gamma_side [variant];
446
size_t variant, S_ige, k, dis_i;
447
to_solve (size_t variant1, size_t S_ige1=0, size_t k1=0, size_t dis_i1=0)
448
: variant(variant1), S_ige(S_ige1), k(k1), dis_i(dis_i1) {}
450
template <class T, class M>
453
const field_basic<T,M>& fh,
454
const level_set_option_type& opt,
455
std::vector<size_t>& bnd_dom_ie_list,
456
array<size_t,M>& sid_ie2bnd_ie)
459
using namespace details;
460
typedef geo_element::size_type size_type;
461
level_set_epsilon = opt.epsilon; // set global variable
463
const geo_basic<T,M>& lambda = fh.get_geo();
464
const space_basic<T,M>& Xh = fh.get_space();
465
check_macro(lambda.order() == 1, "Only order=1 level set mesh supported");
466
check_macro(fh.get_approx() == "P1", "Only P1 level set function supported");
468
std::vector<size_type> dis_idof;
470
size_type d = lambda.dimension();
471
size_type map_dim = d-1;
472
heap_allocator<size_type> alloc;
473
boost::array<list<pair<element_type,size_type> >,
474
reference_element::max_variant> gamma_side_list;
475
list<point_basic<T> > gamma_node_list;
476
std::vector<size_type> j(d+1);
477
std::vector<point_basic<T> > x(d+1);
478
const size_type not_marked = numeric_limits<size_t>::max();
479
const size_type undef = numeric_limits<size_t>::max();
481
distributor node_ownership = lambda.sizes().node_ownership;
482
size_type first_dis_inod = node_ownership.first_index();
483
array<size_type> marked_node (node_ownership, not_marked);
484
array<size_type> extern_node (node_ownership, not_marked);
485
set<size_type> ext_marked_node_idx;
486
list<to_solve> node_to_solve;
488
distributor edge_ownership = lambda.sizes().ownership_by_dimension[1];
489
size_type first_dis_iedg = edge_ownership.first_index();
490
array<size_type> marked_edge (edge_ownership, not_marked);
491
array<std::pair<size_type,point_basic<T> > >
492
extern_edge (edge_ownership, std::make_pair(not_marked,point_basic<T>()));
493
set<size_type> ext_marked_edge_idx;
494
list<to_solve> edge_to_solve;
496
distributor face_ownership = lambda.sizes().ownership_by_dimension[2];
497
size_type first_dis_ifac = face_ownership.first_index();
498
array<size_type> marked_face (face_ownership, not_marked);
500
communicator comm = node_ownership.comm();
501
size_type my_proc = comm.rank();
503
bnd_dom_ie_list.resize(0);
505
// ------------------------------------------------------------
506
// 1) loop on lambda & build intersection sides
507
// ------------------------------------------------------------
508
for (size_type ie = 0, ne = lambda.size(), bnd_ie = 0; ie < ne; ie++) {
509
// ---------------------------------------------------------
510
// 1.1) fast check if there is an intersection:
511
// ---------------------------------------------------------
512
const geo_element& K = lambda [ie];
513
Xh.dis_idof (K, dis_idof);
514
f.resize (dis_idof.size());
515
for (size_type loc_idof = 0, loc_ndof = dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
516
f [loc_idof] = fh.dis_dof (dis_idof[loc_idof]);
518
bool do_intersect = false;
519
switch (K.variant()) {
520
case reference_element::t:
521
do_intersect = belongs_to_band_t (f);
523
case reference_element::T: {
524
do_intersect = belongs_to_band_T (f);
528
error_macro("level set: element type `" << K.name() << "' not yet supported");
530
if (! do_intersect) continue;
531
bnd_dom_ie_list.push_back (ie);
532
// ---------------------------------------------------------
533
// 1.2) compute the intersection
534
// ---------------------------------------------------------
535
x.resize (dis_idof.size());
536
for (size_type loc_idof = 0, loc_ndof = dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
537
size_type loc_inod = loc_idof; // assume here that fh has an isoparametric approx
538
size_type dis_inod = K [loc_inod];
539
x [loc_idof] = lambda.dis_node(dis_inod);
541
element_type S (alloc);
542
switch (K.variant()) {
543
case reference_element::t: {
544
// ---------------------------------------------------------
545
// 1.2.a) triangle -> 1 edge
546
// ---------------------------------------------------------
549
subcompute_matrix_t (x, f, j, a, b, length);
550
if (is_zero(f[j[1]]) && is_zero(f[j[2]])) {
551
// ---------------------------------------------------------
552
// 1.2.a.1) edge {j1,j2} is included in the bbox mesh
553
// ---------------------------------------------------------
554
for (size_type k = 0; k < 2; k++) {
555
size_type dis_inod = K[j[k+1]];
556
if (node_ownership.is_owned (dis_inod)) {
557
size_type inod = dis_inod - first_dis_inod;
558
if (marked_node [inod] == not_marked) {
559
marked_node [inod] = gamma_node_list.size();
560
gamma_node_list.push_back (lambda.node(inod));
563
// dis_inod is owned by another partition jproc:
564
// if there is another neighbour element K' from the same partition jproc
565
// then jproc will insert dis_inod in gamma
566
// so, there is nothing to do here
567
// otherwise, dis_inod will be orphan in jproc and will be re-affected to my_proc
568
extern_node.dis_entry (dis_inod) = my_proc;
571
size_type loc_iedg = edge_t_iloc (j[1], j[2]);
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
S.reset(reference_element::e, order);
577
marked_edge [iedg] = gamma_side_list [S.variant()].size();
578
size_type S_ige = marked_edge [iedg];
579
for (size_type k = 0; k < S.n_node(); k++) {
580
size_type dis_inod = K[j[k+1]];
581
if (node_ownership.is_owned (dis_inod)) {
582
size_type inod = dis_inod - first_dis_inod;
583
S[k] = marked_node [inod];
586
node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
587
ext_marked_node_idx.insert (dis_inod);
590
gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
593
// intersection is an edge of lambda that is owned by another partition jproc
594
// the neighbour element K' accross this e"dge belongs to the same partition jproc
595
// and will insert dis_iedg in gamma
596
// so, there is nothing to do here
599
// ---------------------------------------------------------
600
// 1.2.a.2) edge {j1,j2} is interior to the triangle
601
// ---------------------------------------------------------
602
S.reset(reference_element::e, order);
603
size_type S_ige = gamma_side_list [S.variant()].size ();
604
point_basic<T> xx[2] = {a,b};
605
for (size_type k = 0; k < 2; k++) {
606
if (! is_zero(f[j[k+1]]) && ! is_zero(f[j[0]])) {
607
// xk is inside edge {j0,j[k+1]} of triangle K:
608
size_type loc_iedg = edge_t_iloc (j[0], j[k+1]);
609
size_type dis_iedg = K.edge (loc_iedg);
610
if (edge_ownership.is_owned (dis_iedg)) {
611
size_type iedg = dis_iedg - first_dis_iedg;
612
if (marked_edge [iedg] == not_marked) {
613
marked_edge [iedg] = gamma_node_list.size();
614
gamma_node_list.push_back (xx[k]);
616
S[k] = marked_edge [iedg];
619
edge_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_iedg));
620
ext_marked_edge_idx.insert (dis_iedg);
621
extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
623
} else { // xk is at edge boundary: a node of the 2d mesh
624
size_type dis_inod = (!is_zero(f[j[0]])) ? K[j[k+1]] : K[j[0]];
625
if (node_ownership.is_owned (dis_inod)) {
626
size_type inod = dis_inod - first_dis_inod;
627
if (marked_node [inod] == not_marked) {
628
marked_node [inod] = gamma_node_list.size();
629
gamma_node_list.push_back (lambda.node(inod));
631
S[k] = marked_node [inod];
634
node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
635
ext_marked_node_idx.insert (dis_inod);
636
extern_node.dis_entry (dis_inod) = my_proc;
640
// S[0] == S[1] when is_zero(f[j[0]]) but f[j[0]] != 0, i.e. precision pbs
641
check_macro (S[0] != S[1] || S[0] == undef, "degenerate 2d intersection");
642
gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
646
case reference_element::T: {
647
// ---------------------------------------------------------
648
// 1.2.b) tetrahedron -> 1 triangle or 1 quadrangle
649
// ---------------------------------------------------------
651
point_basic<T> a, b, c, d;
653
if (!intersection_is_quadrilateral_T (f, q)) {
654
subcompute_matrix_triangle_T (x, f, j, a, b, c, aire);
655
if (is_zero(f[j[1]]) && is_zero(f[j[2]]) && is_zero(f[j[3]])) {
656
// the full face {j1,j2,j3} is included in the surface mesh:
657
for (size_type k = 0; k < 3; k++) {
658
size_type dis_inod = K[j[k+1]];
659
if (node_ownership.is_owned (dis_inod)) {
660
size_type inod = dis_inod - first_dis_inod;
661
if (marked_node [inod] == not_marked) {
662
marked_node [inod] = gamma_node_list.size();
663
gamma_node_list.push_back (lambda.node(inod));
666
// dis_inod is owned by another partition jproc:
667
// if there is another neighbour element K' from the same partition jproc
668
// then jproc will insert dis_inod in gamma
669
// so, there is nothing to do here
670
// otherwise, dis_inod will be orphan in jproc and will be re-affected to my_proc
671
extern_node.dis_entry (dis_inod) = my_proc;
674
size_type loc_ifac = face_T_iloc (j[1], j[2], j[3]);
675
size_type dis_ifac = K.face (loc_ifac);
676
if (face_ownership.is_owned (dis_ifac)) {
677
size_type ifac = dis_ifac - first_dis_ifac;
678
if (marked_face [ifac] == not_marked) {
679
S.reset(reference_element::t, order);
680
marked_face [ifac] = gamma_side_list [S.variant()].size();
681
size_type S_ige = marked_face [ifac];
682
for (size_type k = 0; k < S.n_node(); k++) {
683
size_type dis_inod = K[j[k+1]];
684
if (node_ownership.is_owned (dis_inod)) {
685
size_type inod = dis_inod - first_dis_inod;
686
S[k] = marked_node [inod];
689
node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
690
ext_marked_node_idx.insert (dis_inod);
693
gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
695
// the side will be inserted by the neighbour of K in another partition
699
// intersection is a face of lambda that is owned by another partition jproc
700
// the neighbour element K' accross this face belongs to the same partition jproc
701
// and will insert dis_ifac in gamma
702
// so, there is nothing to do here
705
// create the new face {j1,j2,j3} by intersections:
706
S.reset(reference_element::t, order);
707
size_type S_ige = gamma_side_list [S.variant()].size();
708
point_basic<T> xx[3] = {a,b,c};
709
for (size_type k = 0; k < 3; k++) {
710
if (! is_zero(f[j[k+1]]) && ! is_zero(f[j[0]])) {
711
// xk is inside edge {j0,j[k+1]} of triangle K:
712
size_type loc_iedg = edge_T_iloc (j[0], j[k+1]);
713
size_type dis_iedg = K.edge (loc_iedg);
714
if (edge_ownership.is_owned (dis_iedg)) {
715
size_type iedg = dis_iedg - first_dis_iedg;
716
if (marked_edge [iedg] == not_marked) {
717
marked_edge [iedg] = gamma_node_list.size();
718
gamma_node_list.push_back (xx[k]);
720
S[k] = marked_edge [iedg];
723
edge_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_iedg));
724
ext_marked_edge_idx.insert (dis_iedg);
725
extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
727
} else { // xk is at edge boundary: a node of the 2d mesh
728
size_type dis_inod = (!is_zero(f[j[0]])) ? K[j[k+1]] : K[j[0]];
729
if (node_ownership.is_owned (dis_inod)) {
730
size_type inod = dis_inod - first_dis_inod;
731
if (marked_node [inod] == not_marked) {
732
marked_node [inod] = gamma_node_list.size();
733
gamma_node_list.push_back (lambda.node(inod));
735
S[k] = marked_node [inod];
738
node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
739
ext_marked_node_idx.insert (dis_inod);
740
extern_node.dis_entry (dis_inod) = my_proc;
744
// S[0] == S[j] when is_zero(f[j[0]]) but f[j[0]] != 0, i.e. precision pbs
745
check_macro ((S[0] != S[1] || S[0] == undef) &&
746
(S[1] != S[2] || S[1] == undef) &&
747
(S[2] != S[0] || S[2] == undef), "degenerate 3d intersection");
748
gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
751
// create the new quadri face by intersections:
752
subcompute_matrix_quadrilateral_T (x, f, q, a, b, c, d, aire);
753
S.reset (reference_element::q, order);
754
size_type S_ige = gamma_side_list[reference_element::q].size();
755
size_type S1_ige = gamma_side_list[reference_element::t].size(); // or S1+S2=two 't'
756
size_type S2_ige = S1_ige + 1;
757
size_type iloc_S1[4] = {0, 1, 2, undef}; // the way to slip a q into 2 t
758
size_type iloc_S2[4] = {0, undef, 1, 2};
759
point_basic<T> xx[4] = {a,b,d,c};
760
size_type s[4] = {q[0],q[2],q[1],q[3]};
761
for (size_type k = 0; k < 4; k++) {
762
size_type k1 = (k+1) % 4;
763
if (! is_zero(f[s[k]]) && ! is_zero(f[s[k1]])) {
764
// xk is inside edge {j0,j[k+1]} of triangle K:
765
size_type loc_iedg = edge_T_iloc (s[k], s[k1]);
766
size_type dis_iedg = K.edge (loc_iedg);
767
if (edge_ownership.is_owned (dis_iedg)) {
768
size_type iedg = dis_iedg - first_dis_iedg;
769
if (marked_edge [iedg] == not_marked) {
770
marked_edge [iedg] = gamma_node_list.size();
771
gamma_node_list.push_back (xx[k]);
773
S[k] = marked_edge [iedg];
776
if (!opt.split_to_triangle) {
777
edge_to_solve.push_back (to_solve (reference_element::q, S_ige, k, dis_iedg));
779
if (iloc_S1[k] != undef) edge_to_solve.push_back (to_solve (reference_element::t, S1_ige, iloc_S1[k], dis_iedg));
780
if (iloc_S2[k] != undef) edge_to_solve.push_back (to_solve (reference_element::t, S2_ige, iloc_S2[k], dis_iedg));
782
ext_marked_edge_idx.insert (dis_iedg);
783
extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
785
} else { // xk is at edge boundary: a node of the 2d mesh
786
size_type dis_inod = is_zero(f[s[k]]) ? K[s[k]] : K[s[k1]];
787
if (node_ownership.is_owned (dis_inod)) {
788
size_type inod = dis_inod - first_dis_inod;
789
if (marked_node [inod] == not_marked) {
790
marked_node [inod] = gamma_node_list.size();
791
gamma_node_list.push_back (lambda.node(inod));
793
S[k] = marked_node [inod];
796
if (!opt.split_to_triangle) {
797
node_to_solve.push_back (to_solve (reference_element::q, S_ige, k, dis_inod));
799
if (iloc_S1[k] != undef) node_to_solve.push_back (to_solve (reference_element::t, S1_ige, iloc_S1[k], dis_inod));
800
if (iloc_S2[k] != undef) node_to_solve.push_back (to_solve (reference_element::t, S2_ige, iloc_S2[k], dis_inod));
802
ext_marked_node_idx.insert (dis_inod);
806
if (!opt.split_to_triangle) {
807
check_macro ((S[0] != S[1] || S[0] == undef) &&
808
(S[1] != S[2] || S[1] == undef) &&
809
(S[2] != S[3] || S[2] == undef) &&
810
(S[3] != S[0] || S[3] == undef), "degenerate 3d intersection");
811
// S[0] == S[j] when is_zero(f[j[0]]) but f[j[0]] != 0, i.e. precision pbs
812
gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
814
// split quadri into 2 triangles
815
// one K -> two (S1,S2) faces: table element2face may return a pair of size_t
816
// but S1-> and S2->K only is required in form_element
817
element_type S1 (alloc);
818
S1.reset (reference_element::t, order);
819
for (size_type k = 0; k < 4; k++) {
820
if (iloc_S1[k] != undef) S1 [iloc_S1[k]] = S[k];
822
check_macro ((S1[0] != S1[1] || S1[0] == undef) &&
823
(S1[1] != S1[2] || S1[1] == undef) &&
824
(S1[2] != S1[0] || S1[2] == undef), "degenerate 3d intersection");
825
gamma_side_list [S1.variant()].push_back (make_pair(S1,bnd_ie));
827
element_type S2 (alloc);
828
S2.reset (reference_element::t, order);
829
for (size_type k = 0; k < 4; k++) {
830
if (iloc_S2[k] != undef) S2 [iloc_S2[k]] = S[k];
832
check_macro ((S2[0] != S2[1] || S2[0] == undef) &&
833
(S2[1] != S2[2] || S2[1] == undef) &&
834
(S2[2] != S2[0] || S2[2] == undef), "degenerate 3d intersection");
835
gamma_side_list [S2.variant()].push_back (make_pair(S2,bnd_ie));
841
error_macro("level set intersection: element not yet implemented: " << K.name());
846
extern_node.dis_entry_assembly();
847
extern_edge.dis_entry_assembly();
848
// ------------------------------------------------------------
849
// 2) solve orphan nodes, if any
850
// ------------------------------------------------------------
851
// 2.1.a) re-affect orphan node to another process where gamma use it
852
distributor comm_ownership (comm.size(), comm, 1);
853
array<index_set,M> orphan_node (comm_ownership, index_set());
854
for (size_type inod = 0, nnod = node_ownership.size(); inod < nnod; inod++) {
855
if (!(extern_node [inod] != not_marked && marked_node [inod] == not_marked)) continue;
856
// inod is orphan in this proc: not used for gamma (but used for lambda)
857
// re-affect it to a process that use it for gamma
858
size_type iproc = extern_node[inod];
859
size_type dis_inod = first_dis_inod + inod;
860
index_set dis_inod_set;
861
dis_inod_set.insert (dis_inod);
862
orphan_node.dis_entry (iproc) += dis_inod_set;
864
orphan_node.dis_entry_assembly ();
866
// 2.1.b) re-affect orphan edge to another process where gamma use it
867
// there could be orphan edge in 2d if lambda is a part of a regular mesh:
868
// => a boundary edge can belong to another proc. This is not yet handled
869
array<index_set,M> orphan_edge (comm_ownership, index_set());
870
for (size_type iedg = 0, nedg = edge_ownership.size(); iedg < nedg; iedg++) {
871
if (!(extern_edge [iedg].first != not_marked && marked_edge [iedg] == not_marked)) continue;
872
// iedg is orphan in this proc: not used for gamma (but used for lambda)
873
// re-affect it to a process that use it for gamma
874
size_type iproc = extern_edge[iedg].first;
875
size_type dis_iedg = first_dis_iedg + iedg;
876
index_set dis_iedg_set;
877
dis_iedg_set.insert (dis_iedg);
878
orphan_edge.dis_entry (iproc) += dis_iedg_set;
880
orphan_edge.dis_entry_assembly ();
881
check_macro (orphan_edge[0].size() == 0, "unexpected orphan edges");
883
// 2.2) count total nodes used by gamma
884
size_type orphan_gamma_nnod = orphan_node[0].size();
885
size_type regular_gamma_nnod = gamma_node_list.size();
886
size_type gamma_nnod = regular_gamma_nnod + orphan_gamma_nnod;
887
distributor gamma_node_ownership (distributor::decide, comm, gamma_nnod);
888
// 2.3) shift marked_node & marked_edge from gamma_inod local count to gamma_dis_inod one
889
size_type gamma_first_dis_inod = gamma_node_ownership.first_index();
890
for (size_type inod = 0, nnod = node_ownership.size(); inod < nnod; inod++) {
891
if (marked_node [inod] == not_marked) continue;
892
marked_node [inod] += gamma_first_dis_inod;
894
for (size_type iedg = 0, nedg = edge_ownership.size(); iedg < nedg; iedg++) {
895
if (marked_edge [iedg] == not_marked) continue;
896
marked_edge [iedg] += gamma_first_dis_inod;
898
// 2.4) append orphan node to regular one in gamma_node_list and set marked_node
899
for (index_set::const_iterator
900
iter = orphan_node[0].begin(),
901
last = orphan_node[0].end(); iter != last; ++iter) {
902
size_type dis_inod = *iter;
903
marked_node.dis_entry(dis_inod) = gamma_first_dis_inod + gamma_node_list.size();
904
gamma_node_list.push_back (lambda.dis_node(dis_inod));
906
marked_node.dis_entry_assembly();
907
marked_edge.dis_entry_assembly();
908
// ------------------------------------------------------------
909
// 3) convert lists to fixed size distributed arrays
910
// ------------------------------------------------------------
911
array<point_basic<T>, M> gamma_node;
912
boost::array<array<element_type,M>, reference_element::max_variant> gamma_side;
913
gamma_list2array (gamma_node_list, gamma_side_list, comm, d, gamma_node, gamma_side, sid_ie2bnd_ie);
914
// ------------------------------------------------------------
915
// 4) replace inod to dis_inod in element lists
916
// ------------------------------------------------------------
917
for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
918
variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
919
for (size_type ige = 0, nge = gamma_side[variant].size(); ige < nge; ige++) {
920
element_type& S = gamma_side[variant][ige];
921
for (size_type loc_inod = 0, loc_nnod = S.n_node(); loc_inod < loc_nnod; loc_inod++) {
922
if (S[loc_inod] == undef) continue; // external node, will be solved later
923
S[loc_inod] += gamma_first_dis_inod;
927
// ----------------------------------------------------------------
928
// 5) solve intersection that are located on external edges & nodes
929
// ----------------------------------------------------------------
930
marked_node.set_dis_indexes (ext_marked_node_idx);
931
for (list<to_solve>::const_iterator iter = node_to_solve.begin(),
932
last = node_to_solve.end(); iter != last; iter++) {
933
const to_solve& x = *iter;
934
element_type& S = gamma_side[x.variant][x.S_ige];
935
check_macro (S[x.k] == undef, "external index already solved");
936
size_type dis_inod = x.dis_i;
937
size_type iproc = node_ownership.find_owner(dis_inod);
938
size_type gamma_dis_inod = marked_node.dis_at (dis_inod);
939
S[x.k] = gamma_dis_inod;
942
marked_edge.set_dis_indexes (ext_marked_edge_idx);
943
for (list<to_solve>::const_iterator iter = edge_to_solve.begin(),
944
last = edge_to_solve.end(); iter != last; iter++) {
945
const to_solve& x = *iter;
946
element_type& S = gamma_side[x.variant][x.S_ige];
947
check_macro (S[x.k] == undef, "external index already solved");
948
size_type dis_iedg = x.dis_i;
949
size_type iproc = edge_ownership.find_owner(dis_iedg);
950
size_type gamma_dis_inod = marked_edge.dis_at (dis_iedg);
951
S[x.k] = gamma_dis_inod;
953
// ------------------------------------------------------------
954
// 6) convert intersection to geo
955
// ------------------------------------------------------------
956
geo_basic<T,M> gamma (lambda, gamma_node, gamma_side);
959
template <class T, class M>
962
const field_basic<T,M>& fh,
963
const level_set_option_type& opt)
965
typedef geo_element::size_type size_type;
966
std::vector<size_type> bnd_dom_ie_list;
967
array<size_type,M> sid_ie2bnd_ie;
968
return level_set_internal (fh, opt, bnd_dom_ie_list, sid_ie2bnd_ie);
970
// ----------------------------------------------------------------------------
971
// instanciation in library
972
// ----------------------------------------------------------------------------
973
#define _RHEOLEF_instanciation(T,M) \
974
template geo_basic<T,M> level_set_internal ( \
975
const field_basic<T,M>&, \
976
const level_set_option_type&, \
977
std::vector<size_t>&, \
979
template geo_basic<T,M> level_set ( \
980
const field_basic<T,M>&, \
981
const level_set_option_type&);
983
_RHEOLEF_instanciation(Float,sequential)
984
#ifdef _RHEOLEF_HAVE_MPI
985
_RHEOLEF_instanciation(Float,distributed)
986
#endif // _RHEOLEF_HAVE_MPI
added
removed
nfem/plib/level_set.cc
