1 /*============================================================================
2 * Set of subroutines useful for intersecting entities during the joining
3 * operations:
4 * - Intersections of bounding boxes thanks to a box-tree structure,
5 * - Creation of new vertices from edge intersections,
6 * - Synchronizing intersections in parallel mode.
7 *===========================================================================*/
8
9 /*
10 This file is part of Code_Saturne, a general-purpose CFD tool.
11
12 Copyright (C) 1998-2021 EDF S.A.
13
14 This program is free software; you can redistribute it and/or modify it under
15 the terms of the GNU General Public License as published by the Free Software
16 Foundation; either version 2 of the License, or (at your option) any later
17 version.
18
19 This program is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
21 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
22 details.
23
24 You should have received a copy of the GNU General Public License along with
25 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
26 Street, Fifth Floor, Boston, MA 02110-1301, USA.
27 */
28
29 /*----------------------------------------------------------------------------*/
30
31 #include "cs_defs.h"
32
33 /*----------------------------------------------------------------------------
34 * Standard C library headers
35 *---------------------------------------------------------------------------*/
36
37 #include <stdlib.h>
38 #include <stdio.h>
39 #include <string.h>
40 #include <math.h>
41 #include <float.h>
42 #include <assert.h>
43
44 /*----------------------------------------------------------------------------
45 * Local headers
46 *---------------------------------------------------------------------------*/
47
48 #include "bft_mem.h"
49 #include "bft_printf.h"
50
51 #include "fvm_neighborhood.h"
52 #include "fvm_io_num.h"
53
54 #include "cs_all_to_all.h"
55 #include "cs_block_dist.h"
56 #include "cs_join_mesh.h"
57 #include "cs_join_set.h"
58 #include "cs_join_util.h"
59 #include "cs_log.h"
60 #include "cs_order.h"
61 #include "cs_parall.h"
62 #include "cs_search.h"
63 #include "cs_timer.h"
64
65 /*----------------------------------------------------------------------------
66 * Header for the current file
67 *---------------------------------------------------------------------------*/
68
69 #include "cs_join_intersect.h"
70
71 /*---------------------------------------------------------------------------*/
72
73 BEGIN_C_DECLS
74
75 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
76
77 /*============================================================================
78 * Macro definitions
79 *===========================================================================*/
80
81 /*============================================================================
82 * Structure and type definitions
83 *===========================================================================*/
84
85 /* -------------------------------------------------------------
86 * Definition of a structure useful for exchanging intersections
87 * between ranks (_sync_edge_inter())
88 * ------------------------------------------------------------- */
89
90 typedef struct {
91
92 cs_gnum_t vtx_gnum; /* vertex global number relative to the inter. */
93 cs_coord_t curv_abs; /* curvilinear abscissa of the intersection */
94
95 } exch_inter_t;
96
97 /*============================================================================
98 * Global variable definitions
99 *===========================================================================*/
100
101 static int _n_inter_tolerance_warnings;
102
103 static const double _cs_join_tol_eps_coef = 1.0001;
104
105 /*============================================================================
106 * Private function definitions
107 *===========================================================================*/
108
109 /*----------------------------------------------------------------------------
110 * Sort an array "a" and apply the sort to its associated array "b" (local
111 * numbering)
112 * Sort is realized thanks to a shell sort (Knuth algorithm).
113 *
114 * parameters:
115 * l <-- left bound
116 * r <-- right bound
117 * a <-> array to sort
118 * b <-> associated array
119 *---------------------------------------------------------------------------*/
120
121 inline static void
_adapted_lshellsort(cs_lnum_t l,cs_lnum_t r,cs_coord_t a[],cs_lnum_t b[])122 _adapted_lshellsort(cs_lnum_t l,
123 cs_lnum_t r,
124 cs_coord_t a[],
125 cs_lnum_t b[])
126 {
127 int i, j, h;
128 cs_lnum_t size = r - l;
129
130 if (size == 0)
131 return;
132
133 /* Compute stride */
134 for (h = 1; h <= size/9; h = 3*h+1) ;
135
136 /* Sort array */
137 for ( ; h > 0; h /= 3) {
138
139 for (i = l+h; i < r; i++) {
140
141 cs_coord_t va = a[i];
142 cs_lnum_t vb = b[i];
143
144 j = i;
145 while ( (j >= l+h) && (va < a[j-h]) ) {
146 a[j] = a[j-h];
147 b[j] = b[j-h];
148 j -= h;
149 }
150 a[j] = va;
151 b[j] = vb;
152
153 } /* Loop on array elements */
154
155 } /* End of loop on stride */
156 }
157
158 /*----------------------------------------------------------------------------
159 * Sort an array "a" and apply the sort to its associated array "b" (global
160 * numbering)
161 * Sort is realized thanks to a shell sort (Knuth algorithm).
162 *
163 * parameters:
164 * l <-- left bound
165 * r <-- right bound
166 * a <-> array to sort
167 * b <-> associated array
168 *---------------------------------------------------------------------------*/
169
170 inline static void
_adapted_gshellsort(cs_lnum_t l,cs_lnum_t r,cs_coord_t a[],cs_gnum_t b[])171 _adapted_gshellsort(cs_lnum_t l,
172 cs_lnum_t r,
173 cs_coord_t a[],
174 cs_gnum_t b[])
175 {
176 int i, j, h;
177 cs_lnum_t size = r - l;
178
179 if (size == 0)
180 return;
181
182 /* Compute stride */
183 for (h = 1; h <= size/9; h = 3*h+1) ;
184
185 /* Sort array */
186 for ( ; h > 0; h /= 3) {
187
188 for (i = l+h; i < r; i++) {
189
190 cs_coord_t va = a[i];
191 cs_gnum_t vb = b[i];
192
193 j = i;
194 while ( (j >= l+h) && (va < a[j-h]) ) {
195 a[j] = a[j-h];
196 b[j] = b[j-h];
197 j -= h;
198 }
199 a[j] = va;
200 b[j] = vb;
201
202 } /* Loop on array elements */
203
204 } /* End of loop on stride */
205 }
206
207 /*----------------------------------------------------------------------------
208 * Compute the dot product of two vectors.
209 *
210 * parameters:
211 * v1 <-- first vector
212 * v2 <-- second vector
213 *
214 * returns:
215 * the resulting dot product (v1.v2)
216 *----------------------------------------------------------------------------*/
217
218 inline static double
_dot_product(const double v1[],const double v2[])219 _dot_product(const double v1[],
220 const double v2[])
221 {
222 int i;
223 double result = 0.0;
224
225 for (i = 0; i < 3; i++)
226 result += v1[i] * v2[i];
227
228 return result;
229 }
230
231 /*----------------------------------------------------------------------------
232 * Compute the min/max coordinates for the current face in taking into account
233 * the fraction parameter.
234 *
235 * parameters:
236 * face_start <-- start index in face_vtx_lst
237 * face_end <-- end index in face_vtx_lst
238 * face_vtx_lst <-- list of vertices for "face->vertices" connectivity
239 * vertices <-- array on data associated to each selected vertex
240 * extents --> min. and max coordinates of the bounding box
241 *---------------------------------------------------------------------------*/
242
243 static void
_get_face_extents(const cs_lnum_t face_start,const cs_lnum_t face_end,const cs_lnum_t face_vtx_lst[],const cs_join_vertex_t * vertices,cs_coord_t extents[6])244 _get_face_extents(const cs_lnum_t face_start,
245 const cs_lnum_t face_end,
246 const cs_lnum_t face_vtx_lst[],
247 const cs_join_vertex_t *vertices,
248 cs_coord_t extents[6])
249 {
250 cs_lnum_t i, j;
251
252 /* Initalization */
253
254 for (i = 0; i < 3; i++) {
255 extents[i] = DBL_MAX;
256 extents[3 + i] = -DBL_MAX;
257 }
258
259 /* Loop on vertices */
260
261 for (i = face_start; i < face_end; i++) {
262
263 cs_lnum_t vtx_id = face_vtx_lst[i];
264 cs_join_vertex_t vtx = vertices[vtx_id];
265
266 for (j = 0; j < 3; j++) {
267 extents[j] = CS_MIN(extents[j], vtx.coord[j] - vtx.tolerance);
268 extents[3 + j] = CS_MAX(extents[3 + j], vtx.coord[j] + vtx.tolerance);
269 }
270
271 } /* End of loop on face->vertices connectivity */
272
273 }
274
275 /*----------------------------------------------------------------------------
276 * Compute the length of a segment between two vertices.
277 *
278 * parameters:
279 * v1 <-- cs_join_vertex_t structure for the first vertex of the segment
280 * v2 <-- cs_join_vertex_t structure for the second vertex of the segment
281 *
282 * returns:
283 * length of the segment
284 *---------------------------------------------------------------------------*/
285
286 inline static cs_real_t
_compute_length(cs_join_vertex_t v1,cs_join_vertex_t v2)287 _compute_length(cs_join_vertex_t v1,
288 cs_join_vertex_t v2)
289 {
290 cs_lnum_t k;
291 cs_real_t len = 0.0, d2 = 0.0;
292
293 for (k = 0; k < 3; k++) {
294 cs_real_t d = v1.coord[k] - v2.coord[k];
295 d2 += d * d;
296 }
297 len = sqrt(d2);
298
299 return len;
300 }
301
302 /*----------------------------------------------------------------------------
303 * Compute a new cs_join_vertex_t structure.
304 *
305 * parameters:
306 * curv_abs <-- curvilinear abscissa of the intersection
307 * gnum <-- global number associated to the new
308 * cs_join_vertex_t structure
309 * vtx_couple <-- couple of vertex numbers defining the current edge
310 * work <-- local cs_join_mesh_t structure
311 *
312 * returns:
313 * a new cs_join_vertex_t structure
314 *---------------------------------------------------------------------------*/
315
316 static cs_join_vertex_t
_get_new_vertex(cs_coord_t curv_abs,cs_gnum_t gnum,const cs_lnum_t * vtx_couple,const cs_join_mesh_t * work)317 _get_new_vertex(cs_coord_t curv_abs,
318 cs_gnum_t gnum,
319 const cs_lnum_t *vtx_couple,
320 const cs_join_mesh_t *work)
321 {
322 cs_lnum_t k;
323 cs_join_vertex_t new_vtx_data;
324
325 #if defined(DEBUG) && !defined(NDEBUG)
326 /* Avoid Valgrind warnings in byte copies due to padding */
327 memset(&new_vtx_data, 0, sizeof(cs_join_vertex_t));
328 #endif
329
330 cs_join_vertex_t v1 = work->vertices[vtx_couple[0]-1];
331 cs_join_vertex_t v2 = work->vertices[vtx_couple[1]-1];
332
333 assert(curv_abs >= 0.0);
334 assert(curv_abs <= 1.0);
335
336 /* New global number */
337
338 new_vtx_data.gnum = gnum;
339 new_vtx_data.state = CS_JOIN_STATE_NEW;
340
341 /* New tolerance */
342
343 new_vtx_data.tolerance = (1-curv_abs)*v1.tolerance + curv_abs*v2.tolerance;
344
345 /* New coordinates */
346
347 for (k = 0; k < 3; k++)
348 new_vtx_data.coord[k] = (1-curv_abs)*v1.coord[k] + curv_abs*v2.coord[k];
349
350 return new_vtx_data;
351 }
352
353 /*----------------------------------------------------------------------------
354 * Check the coherency of an equivalence. Return true if vertices are each
355 * other under their tolerance.
356 *
357 * parameters:
358 * edges <-- cs_join_edges_t structure
359 * mesh <-- cs_join_mesh_t structure associated
360 * e1_id <-- first edge implied in the equivalence
361 * curv_abs1 <-- curvilinear abscissa of the intersection on edge 1
362 * e2_id <-- second edge implied in the equivalence
363 * curv_abs2 <-- curvilinear abscissa of the intersection on edge 2
364 * verbosity <-- level of accuracy in information display
365 * logfile <-- handle to log file if verbosity > 1
366 *
367 * returns:
368 * true if the check is ok, false otherwise
369 *---------------------------------------------------------------------------*/
370
371 static bool
_check_equiv(const cs_join_edges_t * edges,const cs_join_mesh_t * mesh,cs_lnum_t e1_id,cs_coord_t curv_abs1,cs_lnum_t e2_id,cs_coord_t curv_abs2,int verbosity,FILE * logfile)372 _check_equiv(const cs_join_edges_t *edges,
373 const cs_join_mesh_t *mesh,
374 cs_lnum_t e1_id,
375 cs_coord_t curv_abs1,
376 cs_lnum_t e2_id,
377 cs_coord_t curv_abs2,
378 int verbosity,
379 FILE *logfile)
380 {
381 cs_join_vertex_t p1 = _get_new_vertex(curv_abs1, 1,
382 &(edges->def[2*e1_id]), mesh);
383 cs_join_vertex_t p2 = _get_new_vertex(curv_abs2, 2,
384 &(edges->def[2*e2_id]), mesh);
385 cs_real_t d12 = _compute_length(p1, p2);
386
387 if (p1.tolerance < d12 || p2.tolerance < d12) {
388
389 cs_lnum_t v1e1_id = edges->def[2*e1_id]-1;
390 cs_lnum_t v2e1_id = edges->def[2*e1_id+1]-1;
391 cs_lnum_t v1e2_id = edges->def[2*e2_id]-1;
392 cs_lnum_t v2e2_id = edges->def[2*e2_id+1]-1;
393
394 _n_inter_tolerance_warnings++;
395
396 if (verbosity > 3) {
397
398 fprintf(logfile,
399 "\n"
400 " Edge - Edge intersection warning between:\n"
401 " edge 1: %ld (%llu) [%ld (%llu), %ld (%llu)]\n"
402 " edge 2: %ld (%llu) [%ld (%llu), %ld (%llu)]\n"
403 " Intersection found for curv. abs. %f (e1) - %f (e2)"
404 " will be ignored.\n",
405 (long)e1_id+1, (unsigned long long)edges->gnum[e1_id],
406 (long)v1e1_id+1, (unsigned long long)mesh->vertices[v1e1_id].gnum,
407 (long)v2e1_id+1, (unsigned long long)mesh->vertices[v2e1_id].gnum,
408 (long)e2_id+1, (unsigned long long)edges->gnum[e2_id],
409 (long)v1e2_id+1, (unsigned long long)mesh->vertices[v1e2_id].gnum,
410 (long)v2e2_id+1, (unsigned long long)mesh->vertices[v2e2_id].gnum,
411 curv_abs1, curv_abs2);
412
413 if (p1.tolerance < d12 && verbosity > 4)
414 fprintf(logfile,
415 " Failure for edge 1: "
416 " Distance [v_inter1, v_inter2]: %e > v_inter1.tol: %e\n",
417 d12, p1.tolerance);
418
419 if (p2.tolerance < d12 && verbosity > 4)
420 fprintf(logfile,
421 " Failure for edge 2: "
422 " Distance [v_inter1, v_inter2]: %e > v_inter2.tol: %e\n",
423 d12, p1.tolerance);
424 }
425 return false;
426 }
427 else
428 return true;
429 }
430
431 /*----------------------------------------------------------------------------
432 * Add new equivalences between vertices in cs_join_eset_t structure.
433 *
434 * parameters:
435 * e1_id <-- data relative to edge E1
436 * e2_id <-- data relative to edge E2
437 * abs_e1 <-- curvilinear abscissa of intersection for edge E1
438 * abs_e2 <-- curvilinear abscissa of intersection for edge E2
439 * edges <-- list of edges
440 * vtx_equiv <-> pointer to a structure dealing with vertex equivalences
441 *---------------------------------------------------------------------------*/
442
443 static void
_add_trivial_equiv(cs_lnum_t e1_id,cs_lnum_t e2_id,cs_coord_t abs_e1,cs_coord_t abs_e2,const cs_join_edges_t * edges,cs_join_eset_t * vtx_equiv)444 _add_trivial_equiv(cs_lnum_t e1_id,
445 cs_lnum_t e2_id,
446 cs_coord_t abs_e1,
447 cs_coord_t abs_e2,
448 const cs_join_edges_t *edges,
449 cs_join_eset_t *vtx_equiv)
450 {
451 cs_lnum_t v1_num, v2_num;
452
453 cs_lnum_t equiv_id = vtx_equiv->n_equiv;
454
455 cs_join_eset_check_size(equiv_id, &vtx_equiv);
456
457 /* abs_ei is either 0 or 1, but we avoid exact floating-point comparisons */
458
459 if (abs_e1 < 0.5)
460 v1_num = edges->def[2*e1_id];
461 else
462 v1_num = edges->def[2*e1_id+1];
463
464 if (abs_e2 < 0.5)
465 v2_num = edges->def[2*e2_id];
466 else
467 v2_num = edges->def[2*e2_id+1];
468
469 if (v1_num < v2_num) {
470 vtx_equiv->equiv_couple[2*equiv_id] = v1_num;
471 vtx_equiv->equiv_couple[2*equiv_id+1] = v2_num;
472 }
473 else {
474 vtx_equiv->equiv_couple[2*equiv_id] = v2_num;
475 vtx_equiv->equiv_couple[2*equiv_id+1] = v1_num;
476 }
477
478 vtx_equiv->n_equiv += 1;
479 }
480
481 /*----------------------------------------------------------------------------
482 * Add a no trivial intersection in a cs_join_inter_set_t structure.
483 *
484 * parameters:
485 * e1_id <-- data relative to edge E1
486 * e2_id <-- data relative to edge E2
487 * abs_e1 <-- curvilinear abscissa of intersection(s) for edge E1
488 * abs_e2 <-- curvilinear abscissa of intersection(s) for edge E2
489 * inter_set <-> pointer to the structure maintaining data about edge-edge
490 * intersections
491 *---------------------------------------------------------------------------*/
492
493 static void
_add_inter(cs_lnum_t e1_id,cs_lnum_t e2_id,cs_coord_t abs_e1,cs_coord_t abs_e2,cs_join_inter_set_t * inter_set)494 _add_inter(cs_lnum_t e1_id,
495 cs_lnum_t e2_id,
496 cs_coord_t abs_e1,
497 cs_coord_t abs_e2,
498 cs_join_inter_set_t *inter_set)
499 {
500 cs_join_inter_t new_inter_e1, new_inter_e2;
501
502 cs_lnum_t inter_id = inter_set->n_inter;
503
504 if (inter_id + 1 > inter_set->n_max_inter) {
505
506 inter_set->n_max_inter *= 2;
507
508 BFT_REALLOC(inter_set->inter_lst,
509 2*inter_set->n_max_inter,
510 cs_join_inter_t);
511
512 }
513
514 new_inter_e1.edge_id = e1_id;
515 new_inter_e1.vtx_id = -1;
516 new_inter_e1.curv_abs = abs_e1;
517
518 new_inter_e2.edge_id = e2_id;
519 new_inter_e2.vtx_id = -1;
520 new_inter_e2.curv_abs = abs_e2;
521
522 if (e1_id < e2_id) {
523 inter_set->inter_lst[2*inter_id] = new_inter_e1;
524 inter_set->inter_lst[2*inter_id+1] = new_inter_e2;
525 }
526 else {
527 inter_set->inter_lst[2*inter_id] = new_inter_e2;
528 inter_set->inter_lst[2*inter_id+1] = new_inter_e1;
529 }
530
531 inter_set->n_inter += 1;
532 }
533
534 /*----------------------------------------------------------------------------
535 * Reduce tolerance for the weakest equivalence in order to break this
536 * equivalence
537 *
538 * parameters:
539 * n_elts <-- size of vtx_lst (n_inter on edges + the two extremities)
540 * edge_length <-- length of the edge
541 * equiv_lst <-- true if equiv between two elements else false
542 * abs_lst <-- curvilinear abscissa of each element
543 * tol_lst <-> tolerance of each element
544 *---------------------------------------------------------------------------*/
545
546 static void
_break_equivalence(cs_lnum_t n_elts,double edge_length,bool equiv_lst[],const cs_coord_t abs_lst[],const double tol_lst[])547 _break_equivalence(cs_lnum_t n_elts,
548 double edge_length,
549 bool equiv_lst[],
550 const cs_coord_t abs_lst[],
551 const double tol_lst[])
552 {
553 int i1, i2;
554 double range, _rtf, rtf12, rtf21;
555
556 int k = 0, i1_save = 0;
557 double rtf = -1.0;
558
559 for (i1 = 0; i1 < n_elts - 1; i1++) {
560
561 i2 = i1 + 1;
562
563 if (i1 == 0)
564 k = 0;
565 else
566 k += n_elts - i1;
567
568 if (equiv_lst[k] == true) {
569
570 range = (abs_lst[i2] - abs_lst[i1]) * edge_length;
571 rtf12 = range/tol_lst[i1];
572 rtf21 = range/tol_lst[i2];
573
574 assert(range >= 0.0);
575 assert(rtf12 < 1.0);
576 assert(rtf21 < 1.0);
577
578 _rtf = CS_MAX(rtf12, rtf21);
579 if (_rtf > rtf) {
580 i1_save = i1;
581 rtf = _rtf;
582 }
583
584 }
585
586 } /* End of loop on i1 */
587
588 if (rtf > 0.0) { /* We have find an equivalence to break */
589
590 /* Break equivalence between i1_save and i2_save which the weakest
591 equivalence between two consecutive vertices in the set */
592
593 for (i1 = 0; i1 < n_elts - 1; i1++) {
594
595 i2 = i1 + 1;
596
597 if (i1 == 0)
598 k = 0;
599 else
600 k += n_elts - i1;
601
602 if (i1 == i1_save) {
603 #if 0 && defined(DEBUG) && !defined(NDEBUG)
604 fprintf(cs_glob_join_log,
605 " Break equivalence between [%d, %d]\n", i1, i1_save);
606 #endif
607 for (i2 = i1 + 1; i2 < n_elts; i2++, k++)
608 equiv_lst[k] = false;
609 }
610
611 }
612
613 }
614
615 }
616
617 /*----------------------------------------------------------------------------
618 * Fill equiv_lst[] according to the tolerance related to each vertex
619 *
620 * parameters:
621 * param <-- set of user-defined parameters
622 * n_elts <-- size of vtx_lst (n_inter on edges + the two extremities)
623 * abs_lst <-- curvilinear abscissa of each element
624 * tol_lst <-> tolerance of each element
625 * equiv_lst <-> true if equiv between two elements else false
626 * tag <-- array used to mark equivalent vertices by the same tag
627 * edge_length <-- length of the edge
628 *
629 * returns:
630 * number of tolerance reductions applied
631 *---------------------------------------------------------------------------*/
632
633 static cs_lnum_t
_find_edge_equiv(cs_join_param_t param,cs_lnum_t n_elts,cs_coord_t abs_lst[],double tol_lst[],bool equiv_lst[],cs_lnum_t tag[],double edge_length)634 _find_edge_equiv(cs_join_param_t param,
635 cs_lnum_t n_elts,
636 cs_coord_t abs_lst[],
637 double tol_lst[],
638 bool equiv_lst[],
639 cs_lnum_t tag[],
640 double edge_length)
641 {
642 int i, i1, i2, k;
643 double dist;
644
645 cs_lnum_t n_loops = 0;
646 bool do_reduction = true;
647
648 /* Find equivalence between vertices sharing the same edge */
649
650 for (i1 = 0, k = 0; i1 < n_elts - 1; i1++) {
651 for (i2 = i1 + 1; i2 < n_elts; i2++, k++) {
652
653 dist = (abs_lst[i2] - abs_lst[i1]) * edge_length;
654 assert(dist >= 0.0);
655
656 /* Tag equivalence if test is true */
657
658 if (tol_lst[i1] < dist || tol_lst[i2] < dist)
659 equiv_lst[k] = false;
660 else
661 equiv_lst[k] = true;
662
663 }
664 }
665
666 while (do_reduction == true && n_loops <= param.n_max_equiv_breaks) {
667
668 /* Reduce tolerance after the first loop if necessary */
669
670 if (do_reduction == true && n_loops > 0)
671 _break_equivalence(n_elts,
672 edge_length,
673 equiv_lst,
674 abs_lst,
675 tol_lst);
676
677 /* Tag vertices in equivalence with the same even if they are
678 in equivalence thanks to transitivity */
679
680 for (i = 0; i < n_elts; i++) /* Re-initialize array */
681 tag[i] = i+1;
682
683 for (i1 = 0; i1 < n_elts-1; i1++) {
684
685 if (i1 == 0)
686 k = 0;
687 else
688 k += n_elts - i1;
689
690 if (equiv_lst[k] == true)
691 tag[i1+1] = tag[i1];
692
693 }
694
695 /* Check tolerance consistency: avoid excessive transitivity
696 on an edge */
697
698 do_reduction = false;
699 for (i1 = 0, k = 0; i1 < n_elts - 1; i1++) {
700 for (i2 = i1 + 1; i2 < n_elts; i2++, k++) {
701
702 if (equiv_lst[k] == false)
703 if (tag[i1] == tag[i2])
704 do_reduction = true;
705
706 } /* End of loop on i2 */
707 } /* End of loop on i1 */
708
709 n_loops++;
710
711 } /* End of while reduc_tol == false */
712
713 return n_loops - 1;
714 }
715
716 /*----------------------------------------------------------------------------
717 * Get 3D intersection location(s) between two edges in curvilinear abscissa.
718 *
719 * We successively try to find vertices-vertices matching,
720 * vertex-vertex matching,
721 * a "real" intersection.
722 *
723 * Each vertex owns a tolerance which drives to a sphere in which intersection
724 * is possible.
725 * When an intersection is found, we store the related curvilinear abscissa
726 * for each edge.
727 * If curvilinear abscissa : 0 => -0.01 / 1 => 1.01 in order to test an
728 * inequality rather than an equality.
729 *
730 * Let be two edges E1 (P1E1, P2E1) and E2 (P1E2, P2E2) :
731 * a point on edge E1 is defined by :
732 * P1E1 + s * (P2E1 - P1E1) with 0 <= s <= 1
733 * a point on edge B is defined by :
734 * P1E2 + t * (P2E2 - P1E2) with 0 <= t <= 1
735 *
736 * The squared distance between a point from A and a point from B is :
737 * len(s, t) = || P1E2 - P1E1 + t * (P2E2 - P1E1) - s * (P2E1 - P1E1) ||^2
738 * equivalent to :
739 * len(s, t) = a.s^2 + 2b.s.t + c.t^2 + 2d.s + 2e.t + f = d2_e1e2(s,t)
740 *
741 * We try to find (s, t) such as len(s,t) is minimal
742 * with 0 <= s <= 1 and 0 <= t <= 1
743 *
744 * -> ->
745 * a = v0 . v0 ->
746 * v0
747 * -> -> P1E1 x--------------> P2E1
748 * b = - v0 . v1 \
749 * \ P2E2
750 * -> -> \ ^
751 * c = v1 . v1 \ /
752 * -> \ /
753 * -> -> v2 \ / ->
754 * d = - v0 . v2 \ / v1
755 * \ /
756 * -> -> \ /
757 * e = v1 . v2 x
758 * P1E2
759 *
760 * ->
761 * v0 = vector (P1E1, P2E1)
762 * ->
763 * v1 = vector (P1E2, P2E2)
764 * ->
765 * v2 = vector (P1E1, P1E2)
766 *
767 * parameters:
768 * mesh <-- pointer to joining mesh
769 * edges <-- pointer to edges
770 * fraction <-- global tolerance for geometrical intersection.
771 * e1_id <-- id of edge 1
772 * abs_e1 <-- intersection location on E1 (curvilinear abscissa)
773 * e2_id <-- id of edge 2
774 * abs_e2 <-- intersection location on E2 (curvilinear abscissa)
775 * verbosity <-- level of accuracy in information display
776 * logfile <-- handle to optional log file
777 * n_inter <-> number of intersections detected.
778 *---------------------------------------------------------------------------*/
779
780 static void
_new_edge_edge_3d_inter(const cs_join_mesh_t * mesh,const cs_join_edges_t * edges,double fraction,cs_lnum_t e1_id,double abs_e1[2],cs_lnum_t e2_id,double abs_e2[2],double parall_eps2,int verbosity,FILE * logfile,cs_lnum_t * n_inter)781 _new_edge_edge_3d_inter(const cs_join_mesh_t *mesh,
782 const cs_join_edges_t *edges,
783 double fraction,
784 cs_lnum_t e1_id,
785 double abs_e1[2],
786 cs_lnum_t e2_id,
787 double abs_e2[2],
788 double parall_eps2,
789 int verbosity,
790 FILE *logfile,
791 cs_lnum_t *n_inter)
792 {
793 int k;
794 double a, b, c, d, e, f, s, t;
795 double d2_limit_e1, d_limit_e1, d2_limit_e2, d_limit_e2, d2_e1e2;
796 double inv_cross_norm2, cross_norm2;
797 double int_inter[2] = {0., 0.}, ext_inter[4], v0[3], v1[3], v2[3];
798
799 int n_int_inter = 0, n_ext_inter = 0;
800 bool int_p1e2 = false, int_p2e2 = false;
801
802 cs_lnum_t p1e1_id = edges->def[2*e1_id] - 1;
803 cs_lnum_t p2e1_id = edges->def[2*e1_id+1] - 1;
804 cs_lnum_t p1e2_id = edges->def[2*e2_id] - 1;
805 cs_lnum_t p2e2_id = edges->def[2*e2_id+1] - 1;
806
807 const cs_join_vertex_t p1e1 = mesh->vertices[p1e1_id];
808 const cs_join_vertex_t p2e1 = mesh->vertices[p2e1_id];
809 const cs_join_vertex_t p1e2 = mesh->vertices[p1e2_id];
810 const cs_join_vertex_t p2e2 = mesh->vertices[p2e2_id];
811
812 #if 0 && defined(DEBUG) && !defined(NDEBUG)
813 bool tst_dbg = (verbosity > 6 &&
814 (p1e1.gnum == 716852 || p2e1.gnum == 716852 ||
815 p1e2.gnum == 716852 || p2e2.gnum == 716852) ?
816 true : false);
817 #endif
818
819 /* Initialize parameters */
820
821 *n_inter = 0;
822
823 assert(p1e1.gnum != p2e1.gnum);
824 assert(p1e2.gnum != p2e2.gnum);
825
826 if (p1e1.gnum == p1e2.gnum || p2e1.gnum == p2e2.gnum)
827 return;
828
829 if (p1e1.gnum == p2e2.gnum || p2e1.gnum == p1e2.gnum)
830 return;
831
832 /* Compute local vectors and parameters */
833
834 for (k = 0; k < 3; k++) {
835 v0[k] = p2e1.coord[k] - p1e1.coord[k];
836 v1[k] = p2e2.coord[k] - p1e2.coord[k];
837 v2[k] = p1e2.coord[k] - p1e1.coord[k];
838 }
839
840 a = _dot_product(v0, v0);
841 b = - _dot_product(v0, v1);
842 c = _dot_product(v1, v1);
843 d = - _dot_product(v0, v2);
844 e = _dot_product(v1, v2);
845 f = _dot_product(v2, v2);
846
847 #if 0 && defined(DEBUG) && !defined(NDEBUG)
848 if (tst_dbg && logfile != NULL) {
849 fprintf(logfile,
850 "\n\np1e1 : %10llu - [%10.8e %10.8e %10.8e] - tol: %10.8g\n",
851 (unsigned long long)p1e1.gnum,
852 p1e1.coord[0], p1e1.coord[1], p1e1.coord[2], p1e1.tolerance);
853 fprintf(logfile,
854 "p2e1 : %10llu - [%10.8e %10.8e %10.8e] - tol: %10.8g\n",
855 (unsigned long long)p2e1.gnum,
856 p2e1.coord[0], p2e1.coord[1], p2e1.coord[2], p2e1.tolerance);
857 fprintf(logfile,
858 "p1e2 : %10llu - [%10.8e %10.8e %10.8e] - tol: %10.8g\n",
859 (unsigned long long)p1e2.gnum,
860 p1e2.coord[0], p1e2.coord[1], p1e2.coord[2], p1e2.tolerance);
861 fprintf(logfile,
862 "p2e2 : %10llu - [%10.8e %10.8e %10.8e] - tol: %10.8g\n\n",
863 (unsigned long long)p2e2.gnum,
864 p2e2.coord[0], p2e2.coord[1], p2e2.coord[2], p2e2.tolerance);
865 fprintf(logfile,
866 "v0 : [ %10.8e %10.8e %10.8e]\n", v0[0], v0[1], v0[2]);
867 fprintf(logfile,
868 "v1 : [ %10.8e %10.8e %10.8e]\n", v1[0], v1[1], v1[2]);
869 fprintf(logfile, "v2 : [ %10.8e %10.8e %10.8e]\n\n", v2[0], v2[1], v2[2]);
870 fprintf(logfile, "a : %10.8e, b : %10.8e, c : %10.8e\n", a, b, c);
871 fprintf(logfile, "d : %10.8e, e : %10.8e, f : %10.8e\n", d, e, f);
872 }
873 #endif
874
875 /* Check size of each edge is not equal to 0 */
876
877 assert(a > 0);
878 assert(c > 0);
879
880 /* Check computation of the tolerance */
881
882 assert(sqrt(a) * fraction * 1.00001 >= p1e1.tolerance);
883 assert(sqrt(a) * fraction * 1.00001 >= p2e1.tolerance);
884 assert(sqrt(c) * fraction * 1.00001 >= p1e2.tolerance);
885 assert(sqrt(c) * fraction * 1.00001 >= p2e2.tolerance);
886
887 /*------------------------------------------------------------------
888 * First part : Search for interior/interior intersection
889 * The only possibility is an intersection with 0 < s < 1 and
890 * 0 < t < 1.
891 *------------------------------------------------------------------*/
892
893 /* Edges are parallel if a.c - b^2 = 0.
894 * However, this term is o(length^4). So, we will prefer work with angle
895 * If two edges are not parallel => a.c - b^2 > 0.
896 *
897 * cos(theta) = v0.v1 / ||v0||.||v1||
898 * cos^2 (theta) = (v0.v1)^2 / ((v0.v0).(v1.v1))^2 = (-b)^2 / (a.c)
899 * sin^2 (theta) = 1 - b^2 / (a.c)
900 *
901 * || v0 ^ v1 ||^2 = ||v0||^2 . ||v1||^2 . sin^2(thetha)
902 *
903 */
904
905 cross_norm2 = CS_ABS(a * c - b * b);
906
907 #if 0 && defined(DEBUG) && !defined(NDEBUG)
908 if (tst_dbg && logfile != NULL)
909 fprintf(logfile,
910 " [I]: Test inter. E1 - E2: "
911 "\t cross_norm2: %12.8e - parall_limit: %12.8e\n",
912 cross_norm2, parall_eps2*a*c);
913 #endif
914
915 if (cross_norm2 > parall_eps2 * a * c) {
916
917 /* <=> 1 - b^2/(a.c) < epsilon^2 <=> sin^2(theta) < epsilon^2
918
919 Edges are not parallel.
920 Define s and t if there is an intersection on interior points.
921
922 -> -> ->
923 v2 + s.v0 = t.v1 and cross_norm2 = a.c - b^2 != 0.0
924
925 => s = (b.e - c.d)/cross_norm2
926 => t = (e.a - b.d)/cross_norm2
927 */
928
929 s = b * e - c * d;
930 t = b * d - a * e;
931 inv_cross_norm2 = 1.0 / cross_norm2;
932 s *= inv_cross_norm2;
933 t *= inv_cross_norm2;
934
935 #if 0 && defined(DEBUG) && !defined(NDEBUG)
936 if (tst_dbg && logfile != NULL)
937 fprintf(logfile,
938 " [I]-1: Test inter. E1 - E2: "
939 "\t s: %12.8e - t: %12.8e\n", s, t);
940 #endif
941
942 if (s >= 0. && s <= 1.0) {
943 if (t >= 0. && t <= 1.0) {
944
945 /* If tests are OK, we are on an interior point for E1 and E2 */
946
947 d2_e1e2 = CS_ABS(s*(a*s + 2.0*(b*t + d)) + t*(c*t + 2.0*e) + f);
948 d_limit_e1 = (1.0-s) * p1e1.tolerance + s * p2e1.tolerance;
949 d_limit_e2 = (1.0-t) * p1e2.tolerance + t * p2e2.tolerance;
950 d2_limit_e1 = d_limit_e1 * d_limit_e1;
951 d2_limit_e2 = d_limit_e2 * d_limit_e2;
952
953 #if 0 && defined(DEBUG) && !defined(NDEBUG)
954 if (tst_dbg && logfile != NULL)
955 fprintf(logfile,
956 "\t[I]-2: s = %10.8e, t = %10.8e\n"
957 "\t d2_e1e2: %10.8e, d2_limit_e1: %10.8e, "
958 "d2_limit_e2: %10.8e\n",
959 s, t, d2_e1e2, d2_limit_e1, d2_limit_e2);
960 #endif
961
962 if (d2_e1e2 <= d2_limit_e1 && d2_e1e2 <= d2_limit_e2) {
963
964 /* Under tolerance for edges E1 and E2 => intersection is possible */
965
966 assert(0.0 <= s && s <= 1.0);
967 assert(0.0 <= t && t <= 1.0);
968
969 n_int_inter = 1;
970 int_inter[0] = s;
971 int_inter[1] = t;
972
973 #if 0 && defined(DEBUG) && !defined(NDEBUG)
974 if (tst_dbg && logfile != NULL)
975 fprintf(logfile,
976 "\t[I]-3: Add int. inter. Edge-Edge (%10.8e,%10.8e)\n",
977 s, t);
978 #endif
979
980 } /* If we are under tolerance */
981
982 } /* If it is an interior point for edge E2 */
983 } /* If it is an interior point for edge E1 */
984
985 } /* End of last test : inter. E1 - E2 */
986
987 /*-----------------------------------------------------------------
988 * Second part : we try to find under the current tolerance if
989 * "parallel edges", "extremity-extremity" or "extremity-interior"
990 * are possible.
991 *-----------------------------------------------------------------*/
992
993 /* Distance between P1E1 (s = 0) and a point on E2 */
994 /* ----------------------------------------------- */
995
996 t = -e / c; /* Nearest point on E2 if s = 0 => t = -e/c */
997
998 #if 0 && defined(DEBUG) && !defined(NDEBUG)
999 if (tst_dbg && logfile != NULL)
1000 fprintf(logfile, " [II] Test inter. P1E1 - E2: t = %10.8e\n", t);
1001 #endif
1002
1003 if (t > -fraction && t < 1.0 + fraction) { /* filter */
1004
1005 d2_e1e2 = t*(c*t + 2.0*e) + f;
1006 assert(d2_e1e2 >= -1.e-8);
1007
1008 /* Tolerance for vertex p1e1 (s = 0) */
1009
1010 d2_limit_e1 = p1e1.tolerance * p1e1.tolerance;
1011
1012 if (d2_e1e2 <= d2_limit_e1) { /* Under the tolerance for vertex P1E1 */
1013
1014 /* Tolerance for a vertex on E2 located at "t" */
1015
1016 d_limit_e2 = (1-t)*p1e2.tolerance + t*p2e2.tolerance;
1017 d2_limit_e2 = d_limit_e2 * d_limit_e2;
1018
1019 if (d2_e1e2 <= d2_limit_e2) {
1020
1021 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1022 if (tst_dbg && logfile != NULL)
1023 fprintf(logfile,
1024 "\t[II]-3: Under tol. for P1E1 and a point of E2\n"
1025 "\td2_e1e2 = %10.8e, d2_lim_e1: %10.8e, d2_lim_e2: %10.8e\n",
1026 d2_e1e2, d2_limit_e1, d2_limit_e2);
1027 #endif
1028
1029 /* We are under tolerance for edge E2. We try to find if it is
1030 a vertex-vertex intersection */
1031
1032 if (t < 0.0) { /* Test intersection P1E1-P1E2 */
1033
1034 d2_e1e2 = f; /* s = t = 0.0 */
1035 d2_limit_e2 = p1e2.tolerance * p1e2.tolerance;
1036
1037 if (d2_e1e2 <= d2_limit_e1 && d2_e1e2 <= d2_limit_e2) {
1038
1039 if (_check_equiv(edges, mesh,
1040 e1_id, 0.0, e2_id, 0.0,
1041 verbosity, logfile)) {
1042
1043 /* "extremity-extremity" intersection with s = t = 0.0
1044 under the tolerance is possible */
1045
1046 int_p1e2 = true;
1047 n_ext_inter = 1;
1048 ext_inter[0] = 0.0;
1049 ext_inter[1] = 0.0;
1050
1051 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1052 if (tst_dbg && logfile != NULL)
1053 fprintf(logfile, "\t[II]-4: Add inter. Vtx-Vtx (0.00, 0.00)\n");
1054 #endif
1055 }
1056 }
1057 }
1058 else if (t > 1.0) { /* Test intersection P1E1-P2E2 */
1059
1060 d2_e1e2 = c + 2.0*e + f; /* s = 0.0 and t = 1.0 */
1061 d2_limit_e2 = p2e2.tolerance * p2e2.tolerance;
1062
1063 if (d2_e1e2 <= d2_limit_e1 && d2_e1e2 <= d2_limit_e2) {
1064
1065 if (_check_equiv(edges, mesh,
1066 e1_id, 0.0, e2_id, 1.0,
1067 verbosity, logfile)) {
1068
1069 /* "extremity-extremity" intersection with s = 0.0, t = 1.0
1070 under the tolerance is possible */
1071
1072 int_p2e2 = true;
1073 n_ext_inter = 1;
1074 ext_inter[0] = 0.0;
1075 ext_inter[1] = 1.0;
1076
1077 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1078 if (tst_dbg && logfile != NULL)
1079 fprintf(logfile,
1080 "\t[II]-5: Add inter. Vtx-Vtx (0.00, 1.00)\n");
1081 #endif
1082
1083 }
1084 }
1085 }
1086 else {
1087
1088 assert(0.0 <= t && t <= 1.0);
1089
1090 /* It's an "extremity-interior" intersection */
1091
1092 n_ext_inter = 1;
1093 ext_inter[0] = 0.0;
1094 ext_inter[1] = t;
1095
1096 if (t <= 0.0)
1097 int_p1e2 = true;
1098 if (t >= 1.0)
1099 int_p2e2 = true;
1100
1101 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1102 if (tst_dbg && logfile != NULL)
1103 fprintf(logfile,
1104 "\t[II]-6: Add inter. Vtx-Edge (0.00, %10.8e)\n", t);
1105 #endif
1106
1107 }
1108
1109 } /* End if d2_e1e2 <= d2_limit_e2 */
1110
1111 } /* End if d2_e1e2 <= d2_limit_e1 */
1112
1113 } /* End if -fraction < t < 1 + fraction */
1114
1115 /* Distance between P2E1 (s = 1) and a point on edge E2 */
1116 /* ---------------------------------------------------- */
1117
1118 t = -(b + e) / c; /* Minimum for d2_e1e2 is reached for t = -(b+e) / c
1119 when s = 1. */
1120
1121 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1122 if (tst_dbg && logfile != NULL)
1123 fprintf(logfile,
1124 " [III] Test inter. P2E1 - E2: t = %10.8e\n", t);
1125 #endif
1126
1127 if (t > -fraction && t < 1.0 + fraction) { /* filter */
1128
1129 d2_e1e2 = CS_ABS((a + 2.0*(b*t + d)) + t*(c*t + 2.0*e) + f);
1130
1131 assert(d2_e1e2 >= -1.e-8);
1132
1133 /* Tolerance for vertex p2e1 (s = 1) */
1134
1135 d2_limit_e1 = p2e1.tolerance * p2e1.tolerance;
1136
1137 if (d2_e1e2 <= d2_limit_e1) { /* Under the tolerance for vertex P2E1 */
1138
1139 d_limit_e2 = (1.0-t)*p1e2.tolerance + t*p2e2.tolerance;
1140 d2_limit_e2 = d_limit_e2 * d_limit_e2;
1141
1142 if (d2_e1e2 <= d2_limit_e2) { /* Under tolerance for edge 2 */
1143
1144 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1145 if (tst_dbg && logfile != NULL)
1146 fprintf(logfile,
1147 "\t[III]-3: Under tol. for P2E1 and a point of E2\n"
1148 "\td2_e1e2 = %10.8e, d2_lim_e1: %10.8e, d2_lim_e2: %10.8e\n",
1149 d2_e1e2, d2_limit_e1, d2_limit_e2);
1150 #endif
1151
1152 /* We are under tolerance for edge E2. We try to find if it is
1153 a vertex-vertex intersection */
1154
1155 if (t < 0.0) { /* Test intersection P2E1-P1E2 */
1156
1157 d2_e1e2 = CS_ABS(a + 2.0*d + f);
1158 d2_limit_e2 = p1e2.tolerance * p1e2.tolerance;
1159
1160 if (d2_e1e2 <= d2_limit_e1 && d2_e1e2 <= d2_limit_e2) {
1161
1162 if (_check_equiv(edges, mesh,
1163 e1_id, 1.0, e2_id, 0.0,
1164 verbosity, logfile)) {
1165
1166 /* "extremity-extremity" intersection with s = 1.0, t = 0.0
1167 under the tolerance is possible */
1168
1169 int_p1e2 = true;
1170 ext_inter[2*n_ext_inter] = 1.0;
1171 ext_inter[2*n_ext_inter+1] = 0.0;
1172 n_ext_inter += 1;
1173
1174 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1175 if (tst_dbg && logfile != NULL)
1176 fprintf(logfile,
1177 "\t[III]-4: Add inter. Vtx-Vtx (1.00, 0.00)\n");
1178 #endif
1179
1180 }
1181 }
1182 }
1183 else if (t > 1.0) { /* Test intersection P2E1-P2E2 */
1184
1185 d2_e1e2 = CS_ABS(a + 2.0*(b + d + e) + c + f);
1186 d2_limit_e2 = p2e2.tolerance * p2e2.tolerance;
1187
1188 if (d2_e1e2 <= d2_limit_e1 && d2_e1e2 <= d2_limit_e2) {
1189
1190 if (_check_equiv(edges, mesh,
1191 e1_id, 1.0, e2_id, 1.0,
1192 verbosity, logfile)) {
1193
1194 /* "extremity-extremity" intersection with s = t = 1.0
1195 under the tolerance is possible */
1196
1197 int_p2e2 = true;
1198 ext_inter[2*n_ext_inter] = 1.0;
1199 ext_inter[2*n_ext_inter+1] = 1.0;
1200 n_ext_inter += 1;
1201
1202 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1203 if (tst_dbg && logfile != NULL)
1204 fprintf(logfile,
1205 "\t[III]-5: Add inter. Vtx-Vtx (1.00, 1.00)\n");
1206 #endif
1207
1208 }
1209 }
1210 }
1211 else {
1212
1213 assert(0.0 <= t && t <= 1.0);
1214
1215 /* It's an "extremity-interior" intersection */
1216
1217 ext_inter[2*n_ext_inter] = 1.0;
1218 ext_inter[2*n_ext_inter+1] = t;
1219 n_ext_inter += 1;
1220
1221 if (t <= 0.0)
1222 int_p1e2 = true;
1223 if (t >= 1.0)
1224 int_p2e2 = true;
1225
1226 assert(n_ext_inter <= 2);
1227
1228 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1229 if (tst_dbg && logfile != NULL)
1230 fprintf(logfile,
1231 "\t[III]-6: Add inter. Vtx-Edge (1.00, %10.8e)\n", t);
1232 #endif
1233 }
1234
1235 } /* End if d2_e1e2 <= d2_limit_e2 */
1236
1237 } /* End if d2_e1e2 <= d2_limit_e1 */
1238
1239 } /* End if -fraction < t < 1 + fraction */
1240
1241 /* If vertices from edge E2 are not implied in an intersection
1242 we do a test for each vertex */
1243
1244 if (int_p1e2 == false && n_ext_inter < 2) {
1245
1246 /* Distance between P1E2 (t = 0.0) on edge E2 and a point on edge E1 */
1247 /* ----------------------------------------------------------------- */
1248
1249 s = -d / a; /* t = 0.0 */
1250
1251 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1252 if (tst_dbg && logfile != NULL)
1253 fprintf(logfile,
1254 " [IV] Test inter. P1E2 - E1: s = %10.8e\n", s);
1255 #endif
1256
1257 if (s > 0.0 && s < 1.0) { /* s = 0.0 and s = 1.0 are already tested */
1258
1259 d2_e1e2 = CS_ABS(s*(a*s + 2.0*d) + f);
1260 d2_limit_e2 = p1e2.tolerance * p1e2.tolerance;
1261
1262 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1263 if (tst_dbg && logfile != NULL)
1264 fprintf(logfile,
1265 "\t [IV-1a] d2_e1e2: %10.8e - d2_lim_e2: %10.8e\n",
1266 d2_e1e2, d2_limit_e2);
1267 #endif
1268
1269 if (d2_e1e2 <= d2_limit_e2) { /* Under the tolerance for vertex P1E2 */
1270
1271 d_limit_e1 = (1.0-s)*p1e1.tolerance + s*p2e1.tolerance;
1272 d2_limit_e1 = d_limit_e1 * d_limit_e1;
1273
1274 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1275 if (tst_dbg && logfile != NULL)
1276 fprintf(logfile,
1277 "\t [IV-1b] d2_e1e2: %10.8e - d2_lim_e1: %10.8e\n",
1278 d2_e1e2, d2_limit_e1);
1279 #endif
1280
1281 if (d2_e1e2 <= d2_limit_e1) {
1282
1283 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1284 if (tst_dbg && logfile != NULL)
1285 fprintf
1286 (logfile,
1287 "\t[IV]-2: Under tol. for P1E2 and a point of E1\n"
1288 "\t d2_e1e2 = %10.8e, d2_lim_e1: %10.8e, d2_lim_e2: %10.8e\n",
1289 d2_e1e2, d2_limit_e1, d2_limit_e2);
1290 #endif
1291
1292 /* We are under tolerance for edge E1. There is an intersection
1293 between P1E2 and a point on edge E1 */
1294
1295 ext_inter[2*n_ext_inter] = s;
1296 ext_inter[2*n_ext_inter+1] = 0.0;
1297 n_ext_inter += 1;
1298 assert(n_ext_inter <= 2);
1299
1300 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1301 if (tst_dbg && logfile != NULL)
1302 fprintf(logfile,
1303 "\t[IV]-3: Add inter. Edge-Vtx (%10.8e, 0.00)\n", s);
1304 #endif
1305
1306 } /* End if d2_e1e2 < d2_limit_e1 */
1307
1308 } /* End if d2_e1e2 < d2_limit_e2 */
1309
1310 } /* 0.0 < s < 1.0 */
1311
1312 } /* If int_p1e2 == false && n_ext_inter < 2 */
1313
1314 if (int_p2e2 == false && n_ext_inter < 2) {
1315
1316 /* Distance between P2E2 (t = 1.0) on edge E2 and a point on edge E1 */
1317 /* ----------------------------------------------------------------- */
1318
1319 s = -(b + d) / a; /* t = 1 */
1320
1321 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1322 if (tst_dbg && logfile != NULL)
1323 fprintf(logfile,
1324 " [V] Test inter. P2E2 - E1: s = %10.8e\n", s);
1325 #endif
1326
1327 if (s > 0.0 && s < 1.0) { /* s = 0.0 and s = 1.0 are already tested */
1328
1329 d2_e1e2 = CS_ABS(s*(a*s + 2.0*(b + d)) + c + 2.0*e + f);
1330 d2_limit_e2 = p2e2.tolerance * p2e2.tolerance;
1331
1332 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1333 if (tst_dbg && logfile != NULL)
1334 fprintf(logfile,
1335 "\t [V-1a] d2_e1e2: %10.8e - d2_lim_e2: %10.8e\n",
1336 d2_e1e2, d2_limit_e2);
1337 #endif
1338
1339 if (d2_e1e2 <= d2_limit_e2) { /* Under the tolerance for vertex P2E2 */
1340
1341 d_limit_e1 = (1.0-s)*p1e1.tolerance + s*p2e1.tolerance;
1342 d2_limit_e1 = d_limit_e1 * d_limit_e1;
1343
1344 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1345 if (tst_dbg && logfile != NULL)
1346 fprintf(logfile,
1347 "\t [V-1b] d2_e1e2: %10.8e - d2_lim_e1: %10.8e\n",
1348 d2_e1e2, d2_limit_e1);
1349 #endif
1350
1351 if (d2_e1e2 <= d2_limit_e1) {
1352
1353 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1354 if (tst_dbg && logfile != NULL)
1355 fprintf
1356 (logfile,
1357 "\t[V]-2: Under tol. for P2E2 and a point of E1\n"
1358 "\t d2_e1e2 = %10.8e, d2_lim_e1: %10.8e, d2_lim_e2: %10.8e\n",
1359 d2_e1e2, d2_limit_e1, d2_limit_e2);
1360 #endif
1361
1362 /* We are under tolerance for edge E1. There is an intersection
1363 between P2E2 and a point on edge E1 */
1364
1365 ext_inter[2*n_ext_inter] = s;
1366 ext_inter[2*n_ext_inter+1] = 1.0;
1367 n_ext_inter += 1;
1368 assert(n_ext_inter <= 2);
1369
1370 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1371 if (tst_dbg && logfile != NULL)
1372 fprintf(logfile,
1373 "\t[V]-3: Add inter (%10.8e, 1.01)\n", s);
1374 #endif
1375
1376 } /* End if d2_e1e2 < d2_limit_e1 */
1377
1378 } /* End if d2_e1e2 < d2_limit_e2 */
1379
1380 } /* 0.0 < s < 1.0 */
1381
1382 } /* If int_p2e2 == false && n_ext_inter < 2 */
1383
1384 /* Define intersection(s) to return */
1385
1386 if (n_int_inter == 1) {
1387
1388 if (n_ext_inter < 2) { /* Interior intersection is kept */
1389 *n_inter = 1;
1390 abs_e1[0] = int_inter[0];
1391 abs_e2[0] = int_inter[1];
1392 }
1393 else { /* A choice between interior or extremity intersection
1394 has to be done */
1395
1396 double ds = CS_ABS(ext_inter[2] - ext_inter[0]);
1397 double dt = CS_ABS(ext_inter[3] - ext_inter[1]);
1398
1399 assert(n_ext_inter == 2);
1400
1401 if (ds > fraction && dt > fraction) { /* Keep extremity inter. */
1402 *n_inter = 2;
1403 for (k = 0; k < 2; k++) {
1404 abs_e1[k] = ext_inter[2*k];
1405 abs_e2[k] = ext_inter[2*k+1];
1406 }
1407 }
1408 else { /* Keep interior inter. */
1409 *n_inter = 1;
1410 abs_e1[0] = int_inter[0];
1411 abs_e2[0] = int_inter[1];
1412 }
1413
1414 }
1415
1416 }
1417 else { /* n_int_inter == 0 */
1418
1419 *n_inter = n_ext_inter;
1420
1421 for (k = 0; k < n_ext_inter; k++) {
1422 abs_e1[k] = ext_inter[2*k];
1423 abs_e2[k] = ext_inter[2*k+1];
1424 }
1425
1426 }
1427
1428 }
1429
1430 /*----------------------------------------------------------------------------
1431 * Get 3D intersection location(s) between two edges in curvilinear abscissa.
1432 *
1433 * We successively try to find vertices-vertices matching,
1434 * vertex-vertex matching,
1435 * a "real" intersection.
1436 *
1437 * Each vertex owns a tolerance which drives to a sphere in which intersection
1438 * is possible.
1439 * When an intersection is found, we store the related curvilinear abscissa
1440 * for each edge.
1441 * If curvilinear abscissa : 0 => -0.01 / 1 => 1.01 in order to test an
1442 * inequality rather than an equality.
1443 *
1444 * Let be two edges E1 (P1E1, P2E1) and E2 (P1E2, P2E2) :
1445 * a point on edge E1 is defined by :
1446 * P1E1 + s * (P2E1 - P1E1) with 0 <= s <= 1
1447 * a point on edge B is defined by :
1448 * P1E2 + t * (P2E2 - P1E2) with 0 <= t <= 1
1449 *
1450 * The squared distance between a point from A and a point from B is :
1451 * len(s, t) = || P1E2 - P1E1 + t * (P2E2 - P1E1) - s * (P2E1 - P1E1) ||^2
1452 * equivalent to :
1453 * len(s, t) = a.s^2 + 2b.s.t + c.t^2 + 2d.s + 2e.t + f = d2_e1e2(s,t)
1454 *
1455 * We try to find (s, t) such as len(s,t) is minimal
1456 * with 0 <= s <= 1 and 0 <= t <= 1
1457 *
1458 * -> ->
1459 * a = v0 . v0 ->
1460 * v0
1461 * -> -> P1E1 x--------------> P2E1
1462 * b = - v0 . v1 \
1463 * \ P2E2
1464 * -> -> \ ^
1465 * c = v1 . v1 \ /
1466 * -> \ /
1467 * -> -> v2 \ / ->
1468 * d = - v0 . v2 \ / v1
1469 * \ /
1470 * -> -> \ /
1471 * e = v1 . v2 x
1472 * P1E2
1473 *
1474 * ->
1475 * v0 = vector (P1E1, P2E1)
1476 * ->
1477 * v1 = vector (P1E2, P2E2)
1478 * ->
1479 * v2 = vector (P1E1, P1E2)
1480 *
1481 * parameters:
1482 * mesh <-- pointer to joining mesh
1483 * edges <-- pointer to edges
1484 * fraction <-- global tolerance for geometrical intersection.
1485 * e1_id <-- id of edge 1
1486 * abs_e1 <-- intersection location on E1 (curvilinear abscissa)
1487 * e2_id <-- id of edge 2
1488 * abs_e2 <-- intersection location on E2 (curvilinear abscissa)
1489 * verbosity <-- level of detail in information display
1490 * logfile <-- handle to log file if verbosity > 1
1491 * n_inter <-> number of intersections detected.
1492 *---------------------------------------------------------------------------*/
1493
1494 static void
_edge_edge_3d_inter(const cs_join_mesh_t * mesh,const cs_join_edges_t * edges,double fraction,cs_lnum_t e1_id,double abs_e1[2],cs_lnum_t e2_id,double abs_e2[2],double parall_eps2,int verbosity,FILE * logfile,cs_lnum_t * n_inter)1495 _edge_edge_3d_inter(const cs_join_mesh_t *mesh,
1496 const cs_join_edges_t *edges,
1497 double fraction,
1498 cs_lnum_t e1_id,
1499 double abs_e1[2],
1500 cs_lnum_t e2_id,
1501 double abs_e2[2],
1502 double parall_eps2,
1503 int verbosity,
1504 FILE *logfile,
1505 cs_lnum_t *n_inter)
1506 {
1507 int k;
1508 double a, b, c, d, e, f, s, t;
1509 double d2_limit_e1, d_limit_e1, d2_limit_e2, d_limit_e2, d2_e1e2;
1510 cs_real_t v0[3], v1[3], v2[3];
1511
1512 bool int_p1e2 = false, int_p2e2 = false;
1513
1514 cs_lnum_t p1e1_id = edges->def[2*e1_id] - 1;
1515 cs_lnum_t p2e1_id = edges->def[2*e1_id+1] - 1;
1516 cs_lnum_t p1e2_id = edges->def[2*e2_id] - 1;
1517 cs_lnum_t p2e2_id = edges->def[2*e2_id+1] - 1;
1518
1519 const cs_join_vertex_t p1e1 = mesh->vertices[p1e1_id];
1520 const cs_join_vertex_t p2e1 = mesh->vertices[p2e1_id];
1521 const cs_join_vertex_t p1e2 = mesh->vertices[p1e2_id];
1522 const cs_join_vertex_t p2e2 = mesh->vertices[p2e2_id];
1523
1524 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1525 bool tst_dbg = (verbosity > 6 &&
1526 (p1e1.gnum == 716852 || p2e1.gnum == 716852 ||
1527 p1e2.gnum == 716852 || p2e2.gnum == 716852) ?
1528 true : false);
1529 #endif
1530
1531 /* Initialize parameters */
1532
1533 *n_inter = 0 ;
1534
1535 assert(p1e1.gnum != p2e1.gnum);
1536 assert(p1e2.gnum != p2e2.gnum);
1537
1538 if (p1e1.gnum == p1e2.gnum || p2e1.gnum == p2e2.gnum)
1539 return;
1540
1541 if (p1e1.gnum == p2e2.gnum || p2e1.gnum == p1e2.gnum)
1542 return;
1543
1544 /* Compute local vectors and parameters */
1545
1546 for (k = 0; k < 3; k++) {
1547 v0[k] = p2e1.coord[k] - p1e1.coord[k];
1548 v1[k] = p2e2.coord[k] - p1e2.coord[k];
1549 v2[k] = p1e2.coord[k] - p1e1.coord[k];
1550 }
1551
1552 a = _dot_product(v0, v0);
1553 b = - _dot_product(v0, v1);
1554 c = _dot_product(v1, v1);
1555 d = - _dot_product(v0, v2);
1556 e = _dot_product(v1, v2);
1557 f = _dot_product(v2, v2);
1558
1559 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1560 if (tst_dbg && logfile != NULL) {
1561 fprintf(logfile,
1562 "\n\np1e1 : %10llu - [%10.8e %10.8e %10.8e] - tol: %10.8g\n",
1563 (unsigned long long)p1e1.gnum,
1564 p1e1.coord[0], p1e1.coord[1], p1e1.coord[2], p1e1.tolerance);
1565 fprintf(logfile,
1566 "p2e1 : %10llu - [%10.8e %10.8e %10.8e] - tol: %10.8g\n",
1567 (unsigned long long)p2e1.gnum,
1568 p2e1.coord[0], p2e1.coord[1], p2e1.coord[2], p2e1.tolerance);
1569 fprintf(logfile,
1570 "p1e2 : %10llu - [%10.8e %10.8e %10.8e] - tol: %10.8g\n",
1571 (unsigned long long)p1e2.gnum,
1572 p1e2.coord[0], p1e2.coord[1], p1e2.coord[2], p1e2.tolerance);
1573 fprintf(logfile,
1574 "p2e2 : %10llu - [%10.8e %10.8e %10.8e] - tol: %10.8g\n\n",
1575 (unsigned long long)p2e2.gnum,
1576 p2e2.coord[0], p2e2.coord[1], p2e2.coord[2], p2e2.tolerance);
1577 fprintf(logfile, "v0 : [ %10.8e %10.8e %10.8e]\n", v0[0], v0[1], v0[2]);
1578 fprintf(logfile, "v1 : [ %10.8e %10.8e %10.8e]\n", v1[0], v1[1], v1[2]);
1579 fprintf(logfile, "v2 : [ %10.8e %10.8e %10.8e]\n\n", v2[0], v2[1], v2[2]);
1580 fprintf(logfile, "a : %10.8e, b : %10.8e, c : %10.8e\n", a, b, c);
1581 fprintf(logfile, "d : %10.8e, e : %10.8e, f : %10.8e\n", d, e, f);
1582 }
1583 #endif
1584
1585 /* Check size of each edge is not equal to 0 */
1586
1587 assert(a > 0);
1588 assert(c > 0);
1589
1590 /* Check computation of the tolerance */
1591
1592 assert(sqrt(a) * fraction * 1.00001 >= p1e1.tolerance);
1593 assert(sqrt(a) * fraction * 1.00001 >= p2e1.tolerance);
1594 assert(sqrt(c) * fraction * 1.00001 >= p1e2.tolerance);
1595 assert(sqrt(c) * fraction * 1.00001 >= p2e2.tolerance);
1596
1597 /*-----------------------------------------------------------------
1598 * First part : we try to find under the current tolerance if
1599 * "parallel edges", "extremity-extremity" or "extremity-interior"
1600 * are possible.
1601 *-----------------------------------------------------------------*/
1602
1603 /* Distance between P1E1 (s = 0) and a point on E2 */
1604 /* ----------------------------------------------- */
1605
1606 t = -e / c; /* Nearest point on E2 if s = 0 => t = -e/c */
1607
1608 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1609 if (tst_dbg && logfile != NULL)
1610 fprintf(logfile, " [I] Test inter. P1E1 - E2: t = %10.8e\n", t);
1611 #endif
1612
1613 if (t > -fraction && t < 1.0 + fraction) { /* filter */
1614
1615 d2_e1e2 = t*(c*t + 2.0*e) + f;
1616 assert(d2_e1e2 >= -1.e-8);
1617
1618 /* Tolerance for vertex p1e1 (s = 0) */
1619
1620 d2_limit_e1 = p1e1.tolerance * p1e1.tolerance;
1621
1622 if (d2_e1e2 <= d2_limit_e1) { /* Under the tolerance for vertex P1E1 */
1623
1624 /* Tolerance for a vertex on E2 located at "t" */
1625
1626 d_limit_e2 = (1-t)*p1e2.tolerance + t*p2e2.tolerance;
1627 d2_limit_e2 = d_limit_e2 * d_limit_e2;
1628
1629 if (d2_e1e2 <= d2_limit_e2) {
1630
1631 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1632 if (tst_dbg && logfile != NULL)
1633 fprintf(logfile,
1634 "\t[I]-3: Under tol. for P1E1 and a point of E2\n"
1635 "\td2_e1e2 = %10.8e, d2_lim_e1: %10.8e, d2_lim_e2: %10.8e\n",
1636 d2_e1e2, d2_limit_e1, d2_limit_e2);
1637 #endif
1638
1639 /* We are under tolerance for edge E2. We try to find if it is
1640 a vertex-vertex intersection */
1641
1642 if (t < 0.0) { /* Test intersection P1E1-P1E2 */
1643
1644 d2_e1e2 = f; /* s = t = 0.0 */
1645 d2_limit_e2 = p1e2.tolerance * p1e2.tolerance;
1646
1647 if (d2_e1e2 <= d2_limit_e1 && d2_e1e2 <= d2_limit_e2) {
1648
1649 if (_check_equiv(edges, mesh,
1650 e1_id, 0.0, e2_id, 0.0,
1651 verbosity, logfile)) {
1652
1653 /* "extremity-extremity" intersection with s = t = 0.0
1654 under the tolerance is possible */
1655
1656 int_p1e2 = true;
1657 abs_e1[*n_inter] = 0.0;
1658 abs_e2[*n_inter] = 0.0;
1659
1660 *n_inter += 1;
1661 assert(*n_inter == 1);
1662
1663 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1664 if (tst_dbg && logfile != NULL)
1665 fprintf(logfile,
1666 "\t[I]-4: Add inter. Vtx-Vtx (0.00, 0.00)\n");
1667 #endif
1668 }
1669 }
1670 }
1671 else if (t > 1.0) { /* Test intersection P1E1-P2E2 */
1672
1673 d2_e1e2 = c + 2.0*e + f; /* s = 0.0 and t = 1.0 */
1674 d2_limit_e2 = p2e2.tolerance * p2e2.tolerance;
1675
1676 if (d2_e1e2 <= d2_limit_e1 && d2_e1e2 <= d2_limit_e2) {
1677
1678 if (_check_equiv(edges, mesh,
1679 e1_id, 0.0, e2_id, 1.0,
1680 verbosity, logfile)) {
1681
1682 /* "extremity-extremity" intersection with s = 0.0, t = 1.0
1683 under the tolerance is possible */
1684
1685 int_p2e2 = true;
1686 abs_e1[*n_inter] = 0.0;
1687 abs_e2[*n_inter] = 1.0;
1688
1689 *n_inter += 1;
1690 assert(*n_inter == 1);
1691
1692 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1693 if (tst_dbg && logfile != NULL)
1694 fprintf(logfile,
1695 "\t[I]-5: Add inter. Vtx-Vtx (0.00, 1.00)\n");
1696 #endif
1697
1698 }
1699 }
1700 }
1701 else {
1702
1703 assert(0.0 <= t && t <= 1.0);
1704
1705 /* It's an "extremity-interior" intersection */
1706
1707 abs_e1[*n_inter] = 0.00;
1708 abs_e2[*n_inter] = t;
1709
1710 if (abs_e2[*n_inter] <= 0.0)
1711 int_p1e2 = true;
1712 if (abs_e2[*n_inter] >= 1.0)
1713 int_p2e2 = true;
1714
1715 *n_inter += 1;
1716 assert(*n_inter == 1);
1717
1718
1719 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1720 if (tst_dbg && logfile != NULL)
1721 fprintf(logfile,
1722 "\t[I]-6: Add inter. Vtx-Edge (0.00, %10.8e)\n", t);
1723 #endif
1724
1725 }
1726
1727 } /* End if d2_e1e2 <= d2_limit_e2 */
1728
1729 } /* End if d2_e1e2 <= d2_limit_e1 */
1730
1731 } /* End if -fraction < t < 1 + fraction */
1732
1733 /* Distance between P2E1 (s = 1) and a point on edge E2 */
1734 /* ---------------------------------------------------- */
1735
1736 t = -(b + e) / c; /* Minimum for d2_e1e2 is reached for t = -(b+e) / c
1737 when s = 1. */
1738
1739 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1740 if (tst_dbg && logfile != NULL)
1741 fprintf(logfile, " [II] Test inter. P2E1 - E2: t = %10.8e\n", t);
1742 #endif
1743
1744 if (t > -fraction && t < 1.0 + fraction) { /* filter */
1745
1746 d2_e1e2 = CS_ABS((a + 2.0*(b*t + d)) + t*(c*t + 2.0*e) + f);
1747
1748 assert(d2_e1e2 >= -1.e-8);
1749
1750 /* Tolerance for vertex p2e1 (s = 1) */
1751
1752 d2_limit_e1 = p2e1.tolerance * p2e1.tolerance;
1753
1754 if (d2_e1e2 <= d2_limit_e1) { /* Under the tolerance for vertex P2E1 */
1755
1756 d_limit_e2 = (1.0-t)*p1e2.tolerance + t*p2e2.tolerance;
1757 d2_limit_e2 = d_limit_e2 * d_limit_e2;
1758
1759 if (d2_e1e2 <= d2_limit_e2) { /* Under tolerance for edge 2 */
1760
1761 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1762 if (tst_dbg && logfile != NULL)
1763 fprintf(logfile,
1764 "\t[II]-3: Under tol. for P2E1 and a point of E2\n"
1765 "\td2_e1e2 = %10.8e, d2_lim_e1: %10.8e, d2_lim_e2: %10.8e\n",
1766 d2_e1e2, d2_limit_e1, d2_limit_e2);
1767 #endif
1768
1769 /* We are under tolerance for edge E2. We try to find if it is
1770 a vertex-vertex intersection */
1771
1772 if (t < 0.0) { /* Test intersection P2E1-P1E2 */
1773
1774 d2_e1e2 = CS_ABS(a + 2.0*d + f);
1775 d2_limit_e2 = p1e2.tolerance * p1e2.tolerance;
1776
1777 if (d2_e1e2 <= d2_limit_e1 && d2_e1e2 <= d2_limit_e2) {
1778
1779 if (_check_equiv(edges, mesh,
1780 e1_id, 1.0, e2_id, 0.0,
1781 verbosity, logfile)) {
1782
1783 /* "extremity-extremity" intersection with s = 1.0, t = 0.0
1784 under the tolerance is possible */
1785
1786 int_p1e2 = true;
1787 abs_e1[*n_inter] = 1.0;
1788 abs_e2[*n_inter] = 0.0;
1789
1790 *n_inter += 1;
1791
1792 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1793 if (tst_dbg && logfile != NULL)
1794 fprintf(logfile,
1795 "\t[II]-4: Add inter. Vtx-Vtx (1.00, 0.00)\n");
1796 #endif
1797
1798 }
1799 }
1800 }
1801 else if (t > 1.0) { /* Test intersection P2E1-P2E2 */
1802
1803 d2_e1e2 = CS_ABS(a + 2.0*(b + d + e) + c + f);
1804 d2_limit_e2 = p2e2.tolerance * p2e2.tolerance;
1805
1806 if (d2_e1e2 <= d2_limit_e1 && d2_e1e2 <= d2_limit_e2) {
1807
1808 if (_check_equiv(edges, mesh,
1809 e1_id, 1.0, e2_id, 1.0,
1810 verbosity, logfile)) {
1811
1812 /* "extremity-extremity" intersection with s = t = 1.0
1813 under the tolerance is possible */
1814
1815 int_p2e2 = true;
1816 abs_e1[*n_inter] = 1.00;
1817 abs_e2[*n_inter] = 1.00;
1818
1819 *n_inter += 1;
1820
1821 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1822 if (tst_dbg && logfile != NULL)
1823 fprintf(logfile,
1824 "\t[II]-5: Add inter. Vtx-Vtx (1.00, 1.00)\n");
1825 #endif
1826
1827 }
1828 }
1829 }
1830 else {
1831
1832 assert(0.0 <= t && t <= 1.0);
1833
1834 /* It's an "extremity-interior" intersection */
1835
1836 abs_e1[*n_inter] = 1.00;
1837 abs_e2[*n_inter] = t;
1838
1839 if (abs_e2[*n_inter] <= 0.0)
1840 int_p1e2 = true;
1841 if (abs_e2[*n_inter] >= 1.0)
1842 int_p2e2 = true;
1843
1844 *n_inter += 1;
1845 assert(*n_inter <= 2);
1846
1847 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1848 if (tst_dbg && logfile != NULL)
1849 fprintf(logfile,
1850 "\t[II]-6: Add inter. Vtx-Edge (1.00, %10.8e)\n", t);
1851 #endif
1852 }
1853
1854 } /* End if d2_e1e2 <= d2_limit_e2 */
1855
1856 } /* End if d2_e1e2 <= d2_limit_e1 */
1857
1858 } /* End if -fraction < t < 1 + fraction */
1859
1860 /* If two intersections are already detected, we stop here */
1861
1862 if (*n_inter == 2)
1863 return;
1864
1865 /* If vertices from edge E2 are not implied in an intersection
1866 we do a test for each vertex */
1867
1868 if (int_p1e2 == false) {
1869
1870 /* Distance between P1E2 (t = 0.0) on edge E2 and a point on edge E1 */
1871 /* ----------------------------------------------------------------- */
1872
1873 s = -d / a; /* t = 0.0 */
1874
1875 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1876 if (tst_dbg && logfile != NULL)
1877 fprintf(logfile,
1878 " [III] Test inter. P1E2 - E1: s = %10.8e\n", s);
1879 #endif
1880
1881 if (s > 0.0 && s < 1.0) { /* s = 0.0 and s = 1.0 are already tested */
1882
1883 d2_e1e2 = CS_ABS(s*(a*s + 2.0*d) + f);
1884 d2_limit_e2 = p1e2.tolerance * p1e2.tolerance;
1885
1886 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1887 if (tst_dbg && logfile != NULL)
1888 fprintf(logfile,
1889 "\t [III-1a] d2_e1e2: %10.8e - d2_lim_e2: %10.8e\n",
1890 d2_e1e2, d2_limit_e2);
1891 #endif
1892
1893 if (d2_e1e2 <= d2_limit_e2) { /* Under the tolerance for vertex P1E2 */
1894
1895 d_limit_e1 = (1.0-s)*p1e1.tolerance + s*p2e1.tolerance;
1896 d2_limit_e1 = d_limit_e1 * d_limit_e1;
1897
1898 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1899 if (tst_dbg && logfile != NULL)
1900 fprintf(logfile,
1901 "\t [IV-1b] d2_e1e2: %10.8e - d2_lim_e1: %10.8e\n",
1902 d2_e1e2, d2_limit_e1);
1903 #endif
1904
1905 if (d2_e1e2 <= d2_limit_e1) {
1906
1907 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1908 if (tst_dbg && logfile != NULL)
1909 fprintf
1910 (logfile,
1911 "\t[III]-2: Under tol. for P1E2 and a point of E1\n"
1912 "\t d2_e1e2 = %10.8e, d2_lim_e1: %10.8e, d2_lim_e2: %10.8e\n",
1913 d2_e1e2, d2_limit_e1, d2_limit_e2);
1914 #endif
1915
1916 /* We are under tolerance for edge E1. There is an intersection
1917 between P1E2 and a point on edge E1 */
1918
1919 abs_e1[*n_inter] = s;
1920 abs_e2[*n_inter] = 0.0;
1921
1922 *n_inter += 1;
1923 assert(*n_inter <= 2);
1924
1925 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1926 if (tst_dbg && logfile != NULL)
1927 fprintf(logfile,
1928 "\t[III]-3: Add inter. Edge-Vtx (%10.8e, 0.00)\n", s);
1929 #endif
1930
1931 if (*n_inter == 2)
1932 return;
1933
1934 } /* End if d2_e1e2 < d2_limit_e1 */
1935
1936 } /* End if d2_e1e2 < d2_limit_e2 */
1937
1938 } /* 0.0 < s < 1.0 */
1939
1940 } /* If int_p1e2 == false */
1941
1942 if (int_p2e2 == false) {
1943
1944 /* Distance between P2E2 (t = 1.0) on edge E2 and a point on edge E1 */
1945 /* ----------------------------------------------------------------- */
1946
1947 s = -(b + d) / a; /* t = 1 */
1948
1949 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1950 if (tst_dbg && logfile != NULL)
1951 fprintf(logfile, " [IV] Test inter. P2E2 - E1: s = %10.8e\n", s);
1952 #endif
1953
1954 if (s > 0.0 && s < 1.0) { /* s = 0.0 and s = 1.0 are already tested */
1955
1956 d2_e1e2 = CS_ABS(s*(a*s + 2.0*(b + d)) + c + 2.0*e + f);
1957 d2_limit_e2 = p2e2.tolerance * p2e2.tolerance;
1958
1959 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1960 if (tst_dbg && logfile != NULL)
1961 fprintf(logfile, "\t [IV-1a] d2_e1e2: %10.8e - d2_lim_e2: %10.8e\n",
1962 d2_e1e2, d2_limit_e2);
1963 #endif
1964
1965 if (d2_e1e2 <= d2_limit_e2) { /* Under the tolerance for vertex P2E2 */
1966
1967 d_limit_e1 = (1.0-s)*p1e1.tolerance + s*p2e1.tolerance;
1968 d2_limit_e1 = d_limit_e1 * d_limit_e1;
1969
1970 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1971 if (tst_dbg && logfile != NULL)
1972 fprintf(logfile, "\t [IV-1b] d2_e1e2: %10.8e - d2_lim_e1: %10.8e\n",
1973 d2_e1e2, d2_limit_e1);
1974 #endif
1975
1976 if (d2_e1e2 <= d2_limit_e1) {
1977
1978 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1979 if (tst_dbg && logfile != NULL)
1980 fprintf
1981 (logfile,
1982 "\t[IV]-2: Under tol. for P2E2 and a point of E1\n"
1983 "\t d2_e1e2 = %10.8e, d2_lim_e1: %10.8e, d2_lim_e2: %10.8e\n",
1984 d2_e1e2, d2_limit_e1, d2_limit_e2);
1985 #endif
1986
1987 /* We are under tolerance for edge E1. There is an intersection
1988 between P2E2 and a point on edge E1 */
1989
1990 abs_e1[*n_inter] = s;
1991 abs_e2[*n_inter] = 1.0;
1992
1993 *n_inter += 1;
1994 assert(*n_inter <= 2);
1995
1996 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1997 if (tst_dbg && logfile != NULL)
1998 fprintf(logfile, "\t[IV]-3: Add inter (%10.8e, 1.01)\n", s);
1999 #endif
2000
2001 if (*n_inter == 2)
2002 return;
2003
2004 } /* End if d2_e1e2 < d2_limit_e1 */
2005
2006 } /* End if d2_e1e2 < d2_limit_e2 */
2007
2008 } /* 0.0 < s < 1.0 */
2009
2010 } /* If int_p2e2 == false */
2011
2012 /* If there is at least one intersection, we stop here. */
2013
2014 if (*n_inter > 0)
2015 return;
2016
2017 {
2018 /*------------------------------------------------------------------
2019 * Second part : no intersection has been found ("parallel edges",
2020 * "extremity-extremity" or "extremity-interior"). The only
2021 * remaining possibility is an intersection with 0 < s < 1 and
2022 * 0 < t < 1.
2023 *------------------------------------------------------------------*/
2024
2025 /* Edges are parallel if a.c - b^2 = 0.
2026 * However, this term is o(length^4). So, we will prefer work with angle
2027 * If two edges are not parallel => a.c - b^2 > 0.
2028 *
2029 * cos(theta) = v0.v1 / ||v0||.||v1||
2030 * cos^2 (theta) = (v0.v1)^2 / ((v0.v0).(v1.v1))^2 = (-b)^2 / (a.c)
2031 * sin^2 (theta) = 1 - b^2 / (a.c)
2032 *
2033 * || v0 ^ v1 ||^2 = ||v0||^2 . ||v1||^2 . sin^2(thetha)
2034 *
2035 */
2036
2037 double inv_cross_norm2;
2038 double cross_norm2 = CS_ABS(a * c - b * b);
2039
2040 #if 0 && defined(DEBUG) && !defined(NDEBUG)
2041 if (tst_dbg && logfile != NULL)
2042 fprintf(logfile,
2043 " [V]: Test inter. E1 - E2: cross_norm2: %10.8e\n"
2044 "\ta = %10.8e, b = %10.8e, c = %10.8e, a*c = %10.8e\n",
2045 cross_norm2, a, b, c, a*c);
2046 #endif
2047
2048 if (cross_norm2 < parall_eps2 * a * c) /* <=> 1 - b^2/(a.c) < epsilon^2
2049 <=> sin^2(theta) < epsilon^2 */
2050 return;
2051
2052 /* Edges are not parallel.
2053 Define s and t if there is an intersection on interior points.
2054
2055 -> -> ->
2056 v2 + s.v0 = t.v1 and cross_norm2 = a.c - b^2 != 0.0
2057
2058 => s = (b.e - c.d)/cross_norm2
2059 => t = (e.a - b.d)/cross_norm2
2060 */
2061
2062 s = b * e - c * d;
2063 t = b * d - a * e;
2064
2065 if (s >= 0. && s <= cross_norm2) {
2066 if (t >= 0. && t <= cross_norm2) {
2067
2068 /* If tests are OK, we are on an interior point for E1 and E2 */
2069
2070 inv_cross_norm2 = 1.0 / cross_norm2;
2071 s *= inv_cross_norm2;
2072 t *= inv_cross_norm2;
2073 d2_e1e2 = CS_ABS(s*(a*s + 2.0*(b*t + d)) + t*(c*t + 2.0*e) + f);
2074 d_limit_e1 = (1.0-s) * p1e1.tolerance + s * p2e1.tolerance;
2075 d_limit_e2 = (1.0-t) * p1e2.tolerance + t * p2e2.tolerance;
2076 d2_limit_e1 = d_limit_e1 * d_limit_e1;
2077 d2_limit_e2 = d_limit_e2 * d_limit_e2;
2078
2079 #if 0 && defined(DEBUG) && !defined(NDEBUG)
2080 if (tst_dbg && logfile != NULL)
2081 fprintf(logfile,
2082 "\t[V]-2: s = %10.8e, t = %10.8e\n"
2083 "\t d2_e1e2: %10.8e, d2_limit_e1: %10.8e, "
2084 "d2_limit_e2: %10.8e\n",
2085 s, t, d2_e1e2, d2_limit_e1, d2_limit_e2);
2086 #endif
2087
2088 if (d2_e1e2 <= d2_limit_e1 && d2_e1e2 <= d2_limit_e2) {
2089
2090 /* Under tolerance for edges E1 and E2 => intersection is possible */
2091
2092 assert(0.0 <= s && s <= 1.0);
2093 assert(0.0 <= t && t <= 1.0);
2094
2095 abs_e1[*n_inter] = s;
2096 abs_e2[*n_inter] = t;
2097
2098 #if 0 && defined(DEBUG) && !defined(NDEBUG)
2099 if (tst_dbg && logfile != NULL)
2100 fprintf(logfile,
2101 "\t[V]-3: Add inter. Edge-Edge (%10.8e,%10.8e)\n", s, t);
2102 #endif
2103
2104 *n_inter += 1;
2105 assert(*n_inter == 1);
2106
2107 } /* If we are under tolerance */
2108
2109 } /* If it is an interior point for edge E2 */
2110 } /* If it is an interior point for edge E1 */
2111
2112 } /* End of last test : inter. E1 - E2 */
2113
2114 }
2115
2116 /*----------------------------------------------------------------------------
2117 * Check if we have to add the current intersection into the definition
2118 * of the reference edge description.
2119 *
2120 * parameters:
2121 * inter <-- intersection to test
2122 * ref <-- list of reference edges
2123 * edge_id <-- id of the edge to treat
2124 * cur_shift <-- shift on edges (current used positions)
2125 *
2126 * returns:
2127 * true if we should add this intersection to the cs_join_edge_inter_t
2128 * structure, false otherwise.
2129 *---------------------------------------------------------------------------*/
2130
2131 inline static bool
_need_to_add_exch_inter(exch_inter_t inter,const cs_join_inter_edges_t * ref,cs_lnum_t edge_id,const cs_lnum_t * cur_shift)2132 _need_to_add_exch_inter(exch_inter_t inter,
2133 const cs_join_inter_edges_t *ref,
2134 cs_lnum_t edge_id,
2135 const cs_lnum_t *cur_shift)
2136 {
2137 cs_lnum_t i;
2138
2139 bool ret = true;
2140
2141 for (i = ref->index[edge_id]; i < cur_shift[edge_id]; i++) {
2142
2143 if (ref->vtx_glst[i] == inter.vtx_gnum) {
2144 if (fabs(ref->abs_lst[i] - inter.curv_abs) < 1e-30) {
2145 ret = false;
2146 break;
2147 }
2148 }
2149
2150 } /* End of loop on sub-elements of the current edge */
2151
2152 return ret;
2153 }
2154
2155 /*----------------------------------------------------------------------------
2156 * Update statistics and timings for face bounding-box intersection search.
2157 *
2158 * parameters:
2159 * face_neighborhood <-- pointer to face neighborhood management structure.
2160 * box_wtime <-- bounding box construction time
2161 * box_dim --> dimension of bounding box layout
2162 * stats <-> joining statistics
2163 *---------------------------------------------------------------------------*/
2164
2165 static void
_face_bbox_search_stats(const fvm_neighborhood_t * face_neighborhood,cs_timer_counter_t box_time,int * box_dim,cs_join_stats_t * stats)2166 _face_bbox_search_stats(const fvm_neighborhood_t *face_neighborhood,
2167 cs_timer_counter_t box_time,
2168 int *box_dim,
2169 cs_join_stats_t *stats)
2170 {
2171 int i;
2172 int depth[3];
2173 cs_lnum_t _n_leaves[3], _n_boxes[3];
2174 cs_lnum_t _n_threshold_leaves[3], _n_leaf_boxes[3];
2175 size_t _mem_final[3], _mem_required[3];
2176 double build_wtime, build_cpu_time, query_wtime, query_cpu_time;
2177
2178 int dim = fvm_neighborhood_get_box_stats(face_neighborhood,
2179 depth,
2180 _n_leaves,
2181 _n_boxes,
2182 _n_threshold_leaves,
2183 _n_leaf_boxes,
2184 _mem_final,
2185 _mem_required);
2186
2187 fvm_neighborhood_get_times(face_neighborhood,
2188 &build_wtime,
2189 &build_cpu_time,
2190 &query_wtime,
2191 &query_cpu_time);
2192
2193 cs_timer_counter_t build_time, query_time;
2194
2195 build_time.nsec = build_wtime*1e9;
2196 query_time.nsec = query_wtime*1e9;
2197
2198 for (i = 0; i < 3; i++) {
2199 _mem_final[i] /= 1024;
2200 _mem_required[i] /= 1024;
2201 }
2202
2203 *box_dim = dim;
2204
2205 stats->bbox_layout = CS_MAX(stats->bbox_layout, dim);
2206
2207 if (stats->n_calls < 1) {
2208 stats->bbox_depth[1] = depth[1];
2209 stats->n_leaves[1] = _n_leaves[1];
2210 stats->n_boxes[1] = _n_boxes[1];
2211 stats->n_th_leaves[1] = _n_threshold_leaves[1];
2212 stats->n_leaf_boxes[1] = _n_leaf_boxes[1];
2213 stats->box_mem_final[1] = _mem_final[1];
2214 stats->box_mem_required[1] = _mem_required[1];
2215 }
2216
2217 stats->bbox_depth[0] += depth[0];
2218 stats->bbox_depth[1] = CS_MIN(stats->bbox_depth[1],
2219 (cs_gnum_t)depth[1]);
2220 stats->bbox_depth[2] = CS_MAX(stats->bbox_depth[2],
2221 (cs_gnum_t)depth[2]);
2222
2223 stats->n_leaves[0] += _n_leaves[0];
2224 stats->n_leaves[1] = CS_MIN(stats->n_leaves[1],
2225 (cs_gnum_t)_n_leaves[1]);
2226 stats->n_leaves[2] = CS_MAX(stats->n_leaves[2],
2227 (cs_gnum_t)_n_leaves[2]);
2228
2229 stats->n_boxes[0] += _n_boxes[0];
2230 stats->n_boxes[1] = CS_MIN(stats->n_boxes[1],
2231 (cs_gnum_t)_n_boxes[1]);
2232 stats->n_boxes[2] = CS_MAX(stats->n_boxes[2],
2233 (cs_gnum_t)_n_boxes[2]);
2234
2235 stats->n_th_leaves[0] += _n_threshold_leaves[0];
2236 stats->n_th_leaves[1] = CS_MIN(stats->n_th_leaves[1],
2237 (cs_gnum_t)_n_threshold_leaves[1]);
2238 stats->n_th_leaves[2] = CS_MAX(stats->n_th_leaves[2],
2239 (cs_gnum_t)_n_threshold_leaves[2]);
2240
2241 stats->n_leaf_boxes[0] += _n_leaf_boxes[0];
2242 stats->n_leaf_boxes[1] = CS_MIN(stats->n_leaf_boxes[1],
2243 (cs_gnum_t)_n_leaf_boxes[1]);
2244 stats->n_leaf_boxes[2] = CS_MAX(stats->n_leaf_boxes[2],
2245 (cs_gnum_t)_n_leaf_boxes[2]);
2246
2247 stats->box_mem_final[0] += _mem_final[0];
2248 stats->box_mem_final[1] = CS_MIN(stats->box_mem_final[1], _mem_final[1]);
2249 stats->box_mem_final[2] = CS_MAX(stats->box_mem_final[2], _mem_final[2]);
2250
2251 stats->box_mem_required[0] += _mem_required[0];
2252 stats->box_mem_required[1] = CS_MIN(stats->box_mem_required[1],
2253 (cs_gnum_t)_mem_required[1]);
2254 stats->box_mem_required[2] = CS_MAX(stats->box_mem_required[2],
2255 (cs_gnum_t)_mem_required[2]);
2256
2257 CS_TIMER_COUNTER_ADD(stats->t_box_build, stats->t_box_build, box_time);
2258 CS_TIMER_COUNTER_ADD(stats->t_box_build, stats->t_box_build, build_time);
2259
2260 CS_TIMER_COUNTER_ADD(stats->t_box_query, stats->t_box_query, query_time);
2261 }
2262
2263 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
2264
2265 /*============================================================================
2266 * Public function definitions
2267 *===========================================================================*/
2268
2269 /*----------------------------------------------------------------------------
2270 * Create a new cs_join_inter_set_t structure.
2271 *
2272 * parameters:
2273 * init_size <-- number of init. cs_join_inter_t structure to allocate
2274 *
2275 * returns:
2276 * a pointer to a new inter_set_t structure.
2277 *---------------------------------------------------------------------------*/
2278
2279 cs_join_inter_set_t *
cs_join_inter_set_create(cs_lnum_t init_size)2280 cs_join_inter_set_create(cs_lnum_t init_size)
2281 {
2282 cs_join_inter_set_t *new_set = NULL;
2283
2284 BFT_MALLOC(new_set, 1, cs_join_inter_set_t);
2285
2286 new_set->n_max_inter = init_size; /* default value */
2287 new_set->n_inter = 0;
2288
2289 BFT_MALLOC(new_set->inter_lst, 2*new_set->n_max_inter, cs_join_inter_t);
2290
2291 return new_set;
2292 }
2293
2294 /*----------------------------------------------------------------------------
2295 * Destroy a cs_join_inter_set_t structure.
2296 *
2297 * parameter:
2298 * inter_set <-> a pointer to the inter_set_t structure to destroy
2299 *---------------------------------------------------------------------------*/
2300
2301 void
cs_join_inter_set_destroy(cs_join_inter_set_t ** inter_set)2302 cs_join_inter_set_destroy(cs_join_inter_set_t **inter_set)
2303 {
2304 if (inter_set != NULL) {
2305 if (*inter_set != NULL) {
2306 BFT_FREE((*inter_set)->inter_lst);
2307 BFT_FREE(*inter_set);
2308 }
2309 }
2310 }
2311
2312 /*----------------------------------------------------------------------------
2313 * Dump a cs_join_inter_set_t structure.
2314 *
2315 * parameters:
2316 * f <-- handle to output file
2317 * i_set <-- cs_join_inter_set_t structure to dump
2318 * edges <-- associated cs_join_edge_t structure
2319 * mesh <-- associated cs_join_mesh_t structure
2320 *---------------------------------------------------------------------------*/
2321
2322 void
cs_join_inter_set_dump(FILE * f,const cs_join_inter_set_t * i_set,const cs_join_edges_t * edges,const cs_join_mesh_t * mesh)2323 cs_join_inter_set_dump(FILE *f,
2324 const cs_join_inter_set_t *i_set,
2325 const cs_join_edges_t *edges,
2326 const cs_join_mesh_t *mesh)
2327 {
2328 int i;
2329
2330 fprintf(f, "\n Dump an inter_set_t structure (%p)\n", (const void *)i_set);
2331
2332 fprintf(f, " n_max_inter: %10ld\n", (long)i_set->n_max_inter);
2333 fprintf(f, " n_inter : %10ld\n\n", (long)i_set->n_inter);
2334
2335 for (i = 0; i < i_set->n_inter; i++) {
2336
2337 cs_join_inter_t inter1 = i_set->inter_lst[2*i];
2338 cs_join_inter_t inter2 = i_set->inter_lst[2*i+1];
2339 cs_lnum_t v1e1_id = edges->def[2*inter1.edge_id] - 1;
2340 cs_lnum_t v2e1_id = edges->def[2*inter1.edge_id+1] - 1;
2341 cs_lnum_t v1e2_id = edges->def[2*inter2.edge_id] - 1;
2342 cs_lnum_t v2e2_id = edges->def[2*inter2.edge_id+1] - 1;
2343 cs_gnum_t v1e1 = (mesh->vertices[v1e1_id]).gnum;
2344 cs_gnum_t v2e1 = (mesh->vertices[v2e1_id]).gnum;
2345 cs_gnum_t v1e2 = (mesh->vertices[v1e2_id]).gnum;
2346 cs_gnum_t v2e2 = (mesh->vertices[v2e2_id]).gnum;
2347
2348 fprintf(f, "\n%5d - (%9llu - %9llu)\n",
2349 i, (unsigned long long)edges->gnum[inter1.edge_id],
2350 (unsigned long long)edges->gnum[inter2.edge_id]);
2351 fprintf(f, "E1 [%5llu %5llu] (%6.3f)\n",
2352 (unsigned long long)v1e1, (unsigned long long)v2e1,
2353 inter1.curv_abs);
2354 fprintf(f, "E2 [%5llu %5llu] (%6.3f)\n",
2355 (unsigned long long)v1e2, (unsigned long long)v2e2,
2356 inter2.curv_abs);
2357
2358 } /* End of loop on edge intersections */
2359
2360 fflush(f);
2361 }
2362
2363 /*----------------------------------------------------------------------------
2364 * Allocate and initialize a new cs_join_inter_edges_t structure.
2365 *
2366 * parameters:
2367 * n_edges <-- number of edges
2368 *
2369 * returns:
2370 * a pointer to the created cs_join_inter_edges_t structure.
2371 *---------------------------------------------------------------------------*/
2372
2373 cs_join_inter_edges_t *
cs_join_inter_edges_create(cs_lnum_t n_edges)2374 cs_join_inter_edges_create(cs_lnum_t n_edges)
2375 {
2376 cs_lnum_t i;
2377
2378 cs_join_inter_edges_t *inter_edges = NULL;
2379
2380
2381 /* Allocate and initialize structure */
2382
2383 BFT_MALLOC(inter_edges, 1, cs_join_inter_edges_t);
2384
2385 inter_edges->n_edges = n_edges;
2386
2387 /* Build inter_edges_idx */
2388
2389 BFT_MALLOC(inter_edges->index, n_edges + 1, cs_lnum_t);
2390
2391 for (i = 0; i < n_edges + 1; i++)
2392 inter_edges->index[i] = 0;
2393
2394 BFT_MALLOC(inter_edges->edge_gnum, n_edges, cs_gnum_t);
2395
2396 for (i = 0; i < n_edges; i++)
2397 inter_edges->edge_gnum[i] = 0;
2398
2399 inter_edges->max_sub_size = 0;
2400
2401 inter_edges->vtx_lst = NULL;
2402 inter_edges->vtx_glst = NULL;
2403 inter_edges->abs_lst = NULL;
2404
2405 return inter_edges;
2406 }
2407
2408 /*----------------------------------------------------------------------------
2409 * Build a cs_join_inter_edges_t structure (useful to find equivalence on
2410 * edges and to apply vertex merge to a cs_join_mesh_t structure).
2411 *
2412 * parameters:
2413 * edges <-- cs_join_edges_t structure
2414 * inter_set <-- structure storing data on edge intersections
2415 *
2416 * returns:
2417 * a pointer to the created cs_join_inter_edges_t structure
2418 *---------------------------------------------------------------------------*/
2419
2420 cs_join_inter_edges_t *
cs_join_inter_edges_define(const cs_join_edges_t * edges,const cs_join_inter_set_t * inter_set)2421 cs_join_inter_edges_define(const cs_join_edges_t *edges,
2422 const cs_join_inter_set_t *inter_set)
2423 {
2424 cs_lnum_t i, lst_size, n_edge_inter;
2425
2426 cs_lnum_t max_n_sub_inter = 0;
2427 cs_lnum_t *counter = NULL;
2428 cs_join_inter_edges_t *inter_edges = NULL;
2429
2430 /* Sanity checks */
2431
2432 assert(edges != NULL);
2433 assert(inter_set != NULL);
2434
2435 /* Allocate and initialize structure */
2436
2437 inter_edges = cs_join_inter_edges_create(edges->n_edges);
2438
2439 for (i = 0; i < edges->n_edges; i++)
2440 inter_edges->edge_gnum[i] = edges->gnum[i];
2441
2442 n_edge_inter = 2*inter_set->n_inter;
2443
2444 if (n_edge_inter == 0)
2445 return inter_edges;
2446
2447 for (i = 0; i < n_edge_inter; i++) {
2448
2449 cs_join_inter_t inter = inter_set->inter_lst[i];
2450 cs_lnum_t edge_id = inter.edge_id;
2451
2452 assert(edge_id < edges->n_edges);
2453
2454 if (inter.curv_abs > 0.0 && inter.curv_abs < 1.0)
2455 inter_edges->index[edge_id+1] += 1;
2456
2457 }
2458
2459 for (i = 0; i < edges->n_edges; i++) {
2460
2461 cs_lnum_t n_sub_inter = inter_edges->index[i+1];
2462
2463 max_n_sub_inter = CS_MAX(max_n_sub_inter, n_sub_inter);
2464 inter_edges->index[i+1] += inter_edges->index[i];
2465
2466 }
2467
2468 inter_edges->max_sub_size = max_n_sub_inter;
2469 lst_size = inter_edges->index[edges->n_edges];
2470
2471 /* Fill structures */
2472
2473 BFT_MALLOC(inter_edges->vtx_lst, lst_size, cs_lnum_t);
2474 BFT_MALLOC(inter_edges->abs_lst, lst_size, cs_coord_t);
2475
2476 BFT_MALLOC(counter, edges->n_edges, cs_lnum_t);
2477
2478 for (i = 0; i < edges->n_edges; i++)
2479 counter[i] = 0;
2480
2481 for (i = 0; i < n_edge_inter; i++) {
2482
2483 cs_join_inter_t inter = inter_set->inter_lst[i];
2484 cs_lnum_t edge_id = inter.edge_id;
2485
2486 if (inter.curv_abs > 0.0 && inter.curv_abs < 1.0) {
2487
2488 cs_lnum_t shift = inter_edges->index[edge_id] + counter[edge_id];
2489
2490 inter_edges->vtx_lst[shift] = inter.vtx_id+1;
2491 inter_edges->abs_lst[shift] = inter.curv_abs;
2492 counter[edge_id] += 1;
2493
2494 }
2495
2496 } /* End of loop on intersections */
2497
2498 /* Order lists */
2499
2500 for (i = 0; i < edges->n_edges; i++) {
2501
2502 cs_lnum_t start = inter_edges->index[i];
2503 cs_lnum_t end = inter_edges->index[i+1];
2504
2505 if (end - start > 1)
2506 _adapted_lshellsort(start,
2507 end,
2508 inter_edges->abs_lst,
2509 inter_edges->vtx_lst);
2510
2511 } /* End of loop on edges */
2512
2513 #if 0 && defined(DEBUG) && !defined(NDEBUG) /* Sanity check */
2514
2515 for (i = 0; i < edges->n_edges; i++) {
2516
2517 cs_lnum_t j;
2518 cs_lnum_t start = inter_edges->index[i];
2519 cs_lnum_t end = inter_edges->index[i+1];
2520
2521 if (end - start > 0) {
2522
2523 assert(inter_edges->abs_lst[start] < 1.0);
2524 assert(inter_edges->abs_lst[start] > 0.0);
2525
2526 for (j = start; j < end - 1; j++) {
2527
2528 assert(inter_edges->abs_lst[j+1] < 1.0);
2529 assert(inter_edges->abs_lst[j+1] > 0.0);
2530
2531 if (inter_edges->abs_lst[j] > inter_edges->abs_lst[j+1])
2532 bft_error(__FILE__, __LINE__, 0,
2533 _("\n Incoherency found in inter_edges structure for"
2534 " edge %d (%u):\n"
2535 " Bad ordering of curvilinear abscissa.\n"
2536 " Vertex %d (%u) of abscissa: %f is before vertex %d"
2537 " (%u) of abscissa: %f\n"),
2538 i+1, edges->gnum[i], inter_edges->vtx_lst[j],
2539 inter_edges->abs_lst[j], inter_edges->vtx_lst[j+1],
2540 inter_edges->abs_lst[j+1]);
2541
2542 if (inter_edges->vtx_lst[j] == inter_edges->vtx_lst[j+1])
2543 bft_error(__FILE__, __LINE__, 0,
2544 _("\n Incoherency found in inter_edges structure.\n"
2545 " Redundancy for edge %d (%u) :\n"
2546 " v1_num : %d\n"
2547 " v2_num : %d\n"
2548 " Vertex %d appears twice.\n"), i+1, edges->gnum[i],
2549 edges->def[2*i], edges->def[2*i+1],
2550 inter_edges->vtx_lst[j]);
2551
2552 }
2553
2554 } /* n_sub_inter > 1 */
2555
2556 } /* End of loop on edges */
2557 #endif
2558
2559 /* Free memory */
2560
2561 BFT_FREE(counter);
2562
2563 /* Return pointer */
2564
2565 return inter_edges;
2566 }
2567
2568 /*----------------------------------------------------------------------------
2569 * Destroy an cs_join_inter_edges_t structure.
2570 *
2571 * parameters:
2572 * inter_edges <-> pointer to cs_join_inter_edges_t structure to destroy
2573 *---------------------------------------------------------------------------*/
2574
2575 void
cs_join_inter_edges_destroy(cs_join_inter_edges_t ** inter_edges)2576 cs_join_inter_edges_destroy(cs_join_inter_edges_t **inter_edges)
2577 {
2578 if (inter_edges != NULL) {
2579 cs_join_inter_edges_t *ie = *inter_edges;
2580 if (ie != NULL) {
2581 BFT_FREE(ie->index);
2582 BFT_FREE(ie->edge_gnum);
2583 BFT_FREE(ie->vtx_lst);
2584 BFT_FREE(ie->vtx_glst);
2585 BFT_FREE(ie->abs_lst);
2586 BFT_FREE(*inter_edges);
2587 }
2588 }
2589 }
2590
2591 /*----------------------------------------------------------------------------
2592 * Find non-trivial equivalences between vertices sharing the same edges.
2593 *
2594 * For instance an equivalence between a vertex from the extremity of an edge
2595 * and a vertex created by an edge-edge intersection.
2596 *
2597 * parameters:
2598 * param <-- set of user-defined parameter
2599 * mesh <-- pointer to the local cs_join_mesh_t structure
2600 * which has initial vertex data
2601 * edges <-- list of edges
2602 * inter_edges <-- structure including data on edge intersections
2603 * vtx_equiv <-> structure dealing with vertex equivalences
2604 *---------------------------------------------------------------------------*/
2605
2606 void
cs_join_add_equiv_from_edges(cs_join_param_t param,cs_join_mesh_t * mesh,const cs_join_edges_t * edges,const cs_join_inter_edges_t * inter_edges,cs_join_eset_t * vtx_equiv)2607 cs_join_add_equiv_from_edges(cs_join_param_t param,
2608 cs_join_mesh_t *mesh,
2609 const cs_join_edges_t *edges,
2610 const cs_join_inter_edges_t *inter_edges,
2611 cs_join_eset_t *vtx_equiv)
2612 {
2613 cs_lnum_t i, j, k, i1, i2, size, esize, n_breaks;
2614
2615 bool *equiv_lst = NULL;
2616 cs_lnum_t *vtx_lst = NULL, *tag_lst = NULL;
2617 cs_coord_t *abs_lst = NULL;
2618 double *tol_lst = NULL;
2619 FILE *logfile = cs_glob_join_log;
2620
2621 assert(mesh != NULL);
2622 assert(edges != NULL);
2623 assert(inter_edges != NULL);
2624 assert(vtx_equiv != NULL);
2625
2626 int n_break_counter = 0, n_max_breaks = 0;
2627
2628 if (inter_edges != NULL) {
2629 if (inter_edges->index[inter_edges->n_edges] > 0) {
2630
2631 assert(inter_edges->vtx_lst != NULL);
2632 assert(inter_edges->abs_lst != NULL);
2633
2634 size = inter_edges->max_sub_size + 2;
2635 BFT_MALLOC(vtx_lst, size, cs_lnum_t);
2636 BFT_MALLOC(tag_lst, size, cs_lnum_t);
2637 BFT_MALLOC(abs_lst, size, cs_coord_t);
2638 BFT_MALLOC(tol_lst, size, double);
2639 esize = size*(size-1)/2;
2640 BFT_MALLOC(equiv_lst, esize, bool);
2641
2642 /* Main loop */
2643
2644 for (i = 0; i < inter_edges->n_edges; i++) {
2645
2646 cs_lnum_t v1_num = edges->def[2*i];
2647 cs_lnum_t v2_num = edges->def[2*i+1];
2648 cs_lnum_t v1_id = v1_num - 1;
2649 cs_lnum_t v2_id = v2_num - 1;
2650 cs_lnum_t start = inter_edges->index[i];
2651 cs_lnum_t end = inter_edges->index[i+1];
2652 cs_lnum_t n_sub_elts = 2 + end - start;
2653 double edge_length = _compute_length(mesh->vertices[v1_id],
2654 mesh->vertices[v2_id]);
2655
2656 #if 0 && defined(DEBUG) && !defined(NDEBUG)
2657 if (param.verbosity > 4) {
2658 int vid;
2659 double ds_tol;
2660 cs_gnum_t v1_gnum = (mesh->vertices[v1_num-1]).gnum;
2661 cs_gnum_t v2_gnum = (mesh->vertices[v2_num-1]).gnum;
2662
2663 fprintf(logfile,
2664 "\n%6d: [%9llu] = (%7d [%9llu] - %7d [%9llu] - len: %8.6e)\n",
2665 i, (unsigned long long)edges->gnum[i],
2666 v1_num, (unsigned long long)v1_gnum,
2667 v2_num, (unsigned long long)v2_gnum,
2668 edge_length);
2669 fflush(logfile);
2670
2671 if (inter_edges->vtx_glst == NULL) {
2672
2673 for (j = start, k = 0; j < end; j++, k++) {
2674 vid = inter_edges->vtx_lst[j] - 1;
2675 if (vid > mesh->n_vertices) {
2676 cs_join_mesh_dump(logfile, mesh);
2677 fprintf(logfile, "vid: %d - n_vertices: %d\n",
2678 vid, mesh->n_vertices);
2679 bft_error(__FILE__, __LINE__, 0,
2680 _(" Vertex number out of bounds.\n"));
2681 }
2682 ds_tol = mesh->vertices[vid].tolerance/edge_length;
2683 fprintf(logfile,
2684 " %7d (%9llu) - (%7d, s = %8.6e, ds_tol = %8.6e)\n",
2685 vid+1, (unsigned long long)mesh->vertices[vid].gnum,
2686 k+1, inter_edges->abs_lst[j], ds_tol);
2687 }
2688 }
2689 else {
2690
2691 for (j = start, k = 0; j < end; j++, k++) {
2692 vid = inter_edges->vtx_lst[j] - 1;
2693 if (vid > mesh->n_vertices) {
2694 cs_join_mesh_dump(logfile, mesh);
2695 fprintf(logfile, "vid: %d - n_vertices: %d\n",
2696 vid, mesh->n_vertices);
2697 bft_error(__FILE__, __LINE__, 0,
2698 _(" Vertex number out of bounds.\n"));
2699 }
2700 ds_tol = mesh->vertices[vid].tolerance/edge_length;
2701 fprintf(logfile,
2702 " %9llu - (%7d, s = %8.6e, ds_tol = %8.6e)\n",
2703 (unsigned long long)inter_edges->vtx_glst[j],
2704 k+1, inter_edges->abs_lst[j], ds_tol);
2705 }
2706 }
2707 } /* param.verbosity > 4 */
2708 #endif
2709
2710 /* Build temporary lists */
2711
2712 vtx_lst[0] = v1_num;
2713 abs_lst[0] = 0.0;
2714 tol_lst[0] = (mesh->vertices[v1_id]).tolerance * _cs_join_tol_eps_coef;
2715
2716 for (j = start, k = 1; j < end; j++, k++) {
2717 vtx_lst[k] = inter_edges->vtx_lst[j];
2718 abs_lst[k] = inter_edges->abs_lst[j];
2719 tol_lst[k] = (mesh->vertices[vtx_lst[k]-1]).tolerance
2720 * _cs_join_tol_eps_coef;
2721 }
2722
2723 vtx_lst[k] = v2_num;
2724 abs_lst[k] = 1.0;
2725 tol_lst[k] = (mesh->vertices[v2_id]).tolerance
2726 * _cs_join_tol_eps_coef;
2727
2728 /* Loop on couples of vertices to find if two vertices are equivalent
2729 Apply a tolerance reduction if necessary. */
2730
2731 n_breaks = _find_edge_equiv(param,
2732 n_sub_elts,
2733 abs_lst,
2734 tol_lst,
2735 equiv_lst,
2736 tag_lst,
2737 edge_length);
2738
2739 if (n_breaks > 0) {
2740 n_break_counter += 1;
2741 if (param.verbosity > 3)
2742 fprintf(logfile,
2743 " Edge %8ld: n_equiv. broken: %ld\n",
2744 (long)i+1, (long)n_breaks);
2745 }
2746
2747 n_max_breaks = CS_MAX(n_max_breaks, n_breaks);
2748
2749 /* Add new equivalences */
2750
2751 for (i1 = 0; i1 < n_sub_elts - 1; i1++) {
2752 for (i2 = i1 + 1; i2 < n_sub_elts; i2++) {
2753
2754 if (tag_lst[i1] == tag_lst[i2]) {
2755 if (vtx_lst[i1] != vtx_lst[i2]) {
2756
2757 cs_lnum_t equiv_id = vtx_equiv->n_equiv;
2758 cs_join_eset_check_size(equiv_id, &vtx_equiv);
2759
2760 #if 0 && defined(DEBUG) && !defined(NDEBUG)
2761 if (logfile != NULL)
2762 fprintf
2763 (logfile, " Add equiv %d between [%llu, %llu]\n",
2764 equiv_id+1,
2765 (unsigned long long)mesh->vertices[vtx_lst[i1]-1].gnum,
2766 (unsigned long long)mesh->vertices[vtx_lst[i2]-1].gnum);
2767 #endif
2768
2769 if (vtx_lst[i1] < vtx_lst[i2]) {
2770 vtx_equiv->equiv_couple[2*equiv_id] = vtx_lst[i1];
2771 vtx_equiv->equiv_couple[2*equiv_id+1] = vtx_lst[i2];
2772 }
2773 else {
2774 vtx_equiv->equiv_couple[2*equiv_id] = vtx_lst[i2];
2775 vtx_equiv->equiv_couple[2*equiv_id+1] = vtx_lst[i1];
2776 }
2777 vtx_equiv->n_equiv += 1;
2778
2779 }
2780 } /* Equivalence found */
2781
2782 } /* End of loop on i2 */
2783 } /* End of loop on i1 */
2784
2785 } /* End of loop on edge intersections */
2786
2787 /* Free memory */
2788
2789 BFT_FREE(vtx_lst);
2790 BFT_FREE(tag_lst);
2791 BFT_FREE(abs_lst);
2792 BFT_FREE(tol_lst);
2793 BFT_FREE(equiv_lst);
2794
2795 } /* inter_edges->index[inter_edges->n_edges] > 0 */
2796 } /* inter_edges != NULL */
2797
2798 if (param.verbosity > 0) {
2799
2800 cs_gnum_t n_g_break_counter = n_break_counter;
2801 cs_parall_counter(&n_g_break_counter, 1);
2802
2803 bft_printf(_("\n Equivalences broken for %llu edges.\n"),
2804 (unsigned long long)n_g_break_counter);
2805
2806 if (param.verbosity > 1) {
2807 cs_lnum_t g_n_max_breaks = n_max_breaks;
2808 cs_parall_counter_max(&g_n_max_breaks, 1);
2809 bft_printf(_("\n Max. number of equiv. breaks: %llu\n"),
2810 (unsigned long long)g_n_max_breaks);
2811 }
2812 }
2813
2814 }
2815
2816 #if defined(HAVE_MPI)
2817
2818 /*----------------------------------------------------------------------------
2819 * Synchronize the definition of intersections on each edge by block.
2820 *
2821 * parameters:
2822 * edges <-- cs_join_edges_t structure
2823 * mesh <-- cs_join_mesh_t structure
2824 * part <-- structure storing data on edge intersections by partition
2825 *
2826 * returns:
2827 * newly allocated cs_join_inter_edges_t, synchronized and defined on
2828 * a block
2829 *---------------------------------------------------------------------------*/
2830
2831 cs_join_inter_edges_t *
cs_join_inter_edges_part_to_block(const cs_join_mesh_t * mesh,const cs_join_edges_t * edges,const cs_join_inter_edges_t * part)2832 cs_join_inter_edges_part_to_block(const cs_join_mesh_t *mesh,
2833 const cs_join_edges_t *edges,
2834 const cs_join_inter_edges_t *part)
2835 {
2836 MPI_Comm mpi_comm = cs_glob_mpi_comm;
2837
2838 const int n_ranks = cs_glob_n_ranks;
2839 const int local_rank = cs_glob_rank_id;
2840 const cs_lnum_t n_edges = edges->n_edges;
2841
2842 /* Sanity check */
2843
2844 assert(mesh != NULL);
2845 assert(part != NULL);
2846 assert(part->n_edges == n_edges);
2847 assert(edges != NULL);
2848
2849 cs_block_dist_info_t bi
2850 = cs_block_dist_compute_sizes(local_rank,
2851 n_ranks,
2852 1,
2853 0,
2854 edges->n_g_edges);
2855
2856 cs_all_to_all_t
2857 *d = cs_all_to_all_create_from_block(n_edges,
2858 0, /* flags */
2859 part->edge_gnum,
2860 bi,
2861 mpi_comm);
2862
2863 /* Send global numbers and index */
2864
2865 cs_gnum_t *orig_gnum = cs_all_to_all_copy_array(d,
2866 CS_GNUM_TYPE,
2867 1,
2868 false, /* reverse */
2869 part->edge_gnum,
2870 NULL);
2871
2872 cs_lnum_t n_recv = cs_all_to_all_n_elts_dest(d);
2873
2874 cs_lnum_t *orig_index
2875 = cs_all_to_all_copy_index(d,
2876 false, /* reverse */
2877 part->index,
2878 NULL);
2879
2880 /* Build send_inter_list and exchange information */
2881
2882 cs_lnum_t send_inter_list_size = part->index[n_edges];
2883 exch_inter_t *send_inter_list;
2884 BFT_MALLOC(send_inter_list, send_inter_list_size, exch_inter_t);
2885
2886 for (cs_lnum_t i = 0; i < n_edges; i++) {
2887
2888 cs_lnum_t s_id = part->index[i];
2889 cs_lnum_t e_id = part->index[i+1];
2890
2891 for (cs_lnum_t j = s_id; j < e_id; j++) {
2892 cs_lnum_t vtx_id = part->vtx_lst[j] - 1;
2893 (send_inter_list[j]).vtx_gnum = (mesh->vertices[vtx_id]).gnum;
2894 (send_inter_list[j]).curv_abs = part->abs_lst[j];
2895 }
2896
2897 }
2898
2899 /* Exchange buffers;
2900 index size if adapted temporarily for the datatype (a better
2901 solution would be to allow user datatypes in cs_def.c) */
2902
2903 for (cs_lnum_t i = 0; i < n_recv+1; i++) {
2904 orig_index[i] *= sizeof(exch_inter_t);
2905 }
2906
2907 cs_lnum_t *part_index;
2908 BFT_MALLOC(part_index, n_edges+1, cs_lnum_t);
2909 for (cs_lnum_t i = 0; i < n_edges+1; i++) {
2910 part_index[i] = part->index[i] * sizeof(exch_inter_t);
2911 }
2912
2913 exch_inter_t *recv_inter_list
2914 = cs_all_to_all_copy_indexed(d,
2915 CS_CHAR,
2916 false, /* reverse */
2917 part_index,
2918 send_inter_list,
2919 orig_index,
2920 NULL);
2921
2922 BFT_FREE(part_index);
2923
2924 for (cs_lnum_t i = 0; i < n_recv+1; i++) {
2925 orig_index[i] /= sizeof(exch_inter_t);
2926 }
2927
2928 BFT_FREE(send_inter_list);
2929
2930 cs_all_to_all_destroy(&d);
2931
2932 /* Synchronize the definition of each edge.
2933 Define a new cs_join_inter_edges_t struct. */
2934
2935 cs_lnum_t block_size = 0;
2936 if (bi.gnum_range[1] > bi.gnum_range[0])
2937 block_size = bi.gnum_range[1] - bi.gnum_range[0];
2938
2939 cs_join_inter_edges_t *block
2940 = cs_join_inter_edges_create(block_size);
2941
2942 for (cs_lnum_t i = 0; i < block_size; i++) {
2943 cs_gnum_t g_id = i;
2944 block->edge_gnum[i] = bi.gnum_range[0] + g_id;
2945 }
2946
2947 /* First scan fill index_ref */
2948
2949 for (cs_lnum_t i = 0; i < n_recv; i++) {
2950
2951 cs_gnum_t num = orig_gnum[i] - bi.gnum_range[0] + 1;
2952 cs_lnum_t n_sub_elts = orig_index[i+1] - orig_index[i];
2953
2954 assert(num <= (cs_gnum_t)block_size);
2955 block->index[num] += n_sub_elts;
2956
2957 }
2958
2959 cs_lnum_t *shift_ref;
2960 BFT_MALLOC(shift_ref, block_size, cs_lnum_t);
2961
2962 for (cs_lnum_t i = 0; i < block_size; i++) {
2963 block->index[i+1] += block->index[i];
2964 shift_ref[i] = block->index[i];
2965 }
2966
2967 BFT_MALLOC(block->vtx_glst,
2968 block->index[block_size],
2969 cs_gnum_t);
2970
2971 BFT_MALLOC(block->abs_lst,
2972 block->index[block_size],
2973 cs_coord_t);
2974
2975 /* Second scan: fill buffers */
2976
2977 for (cs_lnum_t i = 0; i < n_recv; i++) {
2978
2979 cs_gnum_t block_id = orig_gnum[i] - bi.gnum_range[0];
2980
2981 assert(block_id < (cs_gnum_t)block_size);
2982
2983 cs_lnum_t s_id = orig_index[i];
2984 cs_lnum_t e_id = orig_index[i+1];
2985
2986 for (cs_lnum_t j = s_id; j < e_id; j++) {
2987
2988 exch_inter_t exch_inter = recv_inter_list[j];
2989
2990 if (_need_to_add_exch_inter(exch_inter,
2991 block,
2992 block_id,
2993 shift_ref)) {
2994
2995 cs_lnum_t _shift = shift_ref[block_id];
2996
2997 assert(_shift < block->index[block_id+1]);
2998
2999 block->vtx_glst[_shift] = exch_inter.vtx_gnum;
3000 block->abs_lst[_shift] = exch_inter.curv_abs;
3001 shift_ref[block_id] += 1;
3002
3003 } /* End of adding a new intersection in the edge definition */
3004
3005 } /* End of loop on sub_elts */
3006
3007 } /* End of loop on edge descriptions */
3008
3009 BFT_FREE(recv_inter_list);
3010
3011 /* Compact block */
3012
3013 BFT_FREE(orig_gnum);
3014 BFT_FREE(orig_index);
3015
3016 cs_lnum_t shift = 0;
3017 for (cs_lnum_t i = 0; i < block_size; i++) {
3018
3019 for (cs_lnum_t j = block->index[i]; j < shift_ref[i]; j++) {
3020
3021 block->vtx_glst[shift] = block->vtx_glst[j];
3022 block->abs_lst[shift] = block->abs_lst[j];
3023 shift++;
3024
3025 }
3026
3027 }
3028
3029 cs_lnum_t *new_index;
3030 BFT_MALLOC(new_index, block_size + 1, cs_lnum_t);
3031
3032 new_index[0] = 0;
3033 for (cs_lnum_t i = 0; i < block_size; i++)
3034 new_index[i+1] = new_index[i] + shift_ref[i] - block->index[i];
3035
3036 BFT_FREE(shift_ref);
3037 BFT_FREE(block->index);
3038
3039 block->index = new_index;
3040
3041 BFT_REALLOC(block->vtx_glst, block->index[block_size], cs_gnum_t);
3042 BFT_REALLOC(block->abs_lst, block->index[block_size], cs_coord_t);
3043
3044 /* Sort intersection by increasing curvilinear abscissa for each edge */
3045
3046 cs_lnum_t _max = 0;
3047
3048 for (cs_lnum_t i = 0; i < block_size; i++)
3049 _max = CS_MAX(_max, new_index[i+1] - new_index[i]);
3050
3051 block->max_sub_size = _max;
3052
3053 for (cs_lnum_t i = 0; i < block_size; i++) {
3054
3055 cs_lnum_t start = block->index[i];
3056 cs_lnum_t end = block->index[i+1];
3057
3058 if (end - start > 0)
3059 _adapted_gshellsort(start, end, block->abs_lst, block->vtx_glst);
3060
3061 } /* End of loop on edges of block */
3062
3063 #if 0 && defined(DEBUG) && !defined(NDEBUG)
3064 cs_join_inter_edges_dump(cs_glob_join_log, block, edges, mesh);
3065 #endif
3066
3067 return block;
3068 }
3069
3070 /*----------------------------------------------------------------------------
3071 * Synchronize the definition of intersections on each edge from a
3072 * cs_join_inter_edges_t structure defined on a block.
3073 *
3074 * parameters:
3075 * n_g_egdes <-- global number of edges
3076 * block <-- synchronized cs_join_inter_edges_t struct. by block
3077 * part <-> cs_join_inter_edges_t to synchronized on partition
3078 *---------------------------------------------------------------------------*/
3079
3080 void
cs_join_inter_edges_block_to_part(cs_gnum_t n_g_edges,const cs_join_inter_edges_t * block,cs_join_inter_edges_t * part)3081 cs_join_inter_edges_block_to_part(cs_gnum_t n_g_edges,
3082 const cs_join_inter_edges_t *block,
3083 cs_join_inter_edges_t *part)
3084 {
3085 MPI_Comm mpi_comm = cs_glob_mpi_comm;
3086
3087 const int n_ranks = cs_glob_n_ranks;
3088 const int local_rank = cs_glob_rank_id;
3089 const cs_block_dist_info_t bi = cs_block_dist_compute_sizes(local_rank,
3090 n_ranks,
3091 1,
3092 0,
3093 n_g_edges);
3094
3095 /* Sanity check */
3096
3097 assert(block != NULL);
3098 assert(part != NULL);
3099
3100 #if defined(DEBUG) && !defined(NDEBUG)
3101 {
3102 cs_lnum_t block_size = 0;
3103 if (bi.gnum_range[1] > bi.gnum_range[0])
3104 block_size = bi.gnum_range[1] - bi.gnum_range[0];
3105 assert(block->n_edges == block_size);
3106 }
3107 #endif
3108
3109 cs_all_to_all_t
3110 *d = cs_all_to_all_create_from_block(part->n_edges,
3111 CS_ALL_TO_ALL_USE_DEST_ID, /* flags */
3112 part->edge_gnum,
3113 bi,
3114 mpi_comm);
3115
3116 /* Exchange global numbers of edges which should be returned
3117 from block to partition */
3118
3119 cs_gnum_t *orig_gnum = cs_all_to_all_copy_array(d,
3120 CS_GNUM_TYPE,
3121 1,
3122 false, /* reverse */
3123 part->edge_gnum,
3124 NULL);
3125
3126 cs_lnum_t send_list_size = cs_all_to_all_n_elts_dest(d);
3127
3128 assert(send_list_size == block->n_edges);
3129
3130 /* Send the number of sub_elements for each requested edge */
3131
3132 cs_lnum_t *block_index;
3133 BFT_MALLOC(block_index, send_list_size+1, cs_lnum_t);
3134
3135 block_index[0] = 0;
3136 for (cs_lnum_t i = 0; i < send_list_size; i++) {
3137 cs_gnum_t s_id = orig_gnum[i] - bi.gnum_range[0];
3138 block_index[i+1] = block_index[i]
3139 + block->index[s_id+1] - block->index[s_id];
3140 }
3141
3142 cs_all_to_all_copy_index(d,
3143 true, /* reverse */
3144 block_index,
3145 part->index);
3146
3147 /* block struct. send the requested global edges to all the part struct
3148 on the distant ranks */
3149
3150 cs_lnum_t send_inter_list_size = block_index[send_list_size];
3151
3152 exch_inter_t *send_inter_list = NULL;
3153 BFT_MALLOC(send_inter_list, send_inter_list_size, exch_inter_t);
3154
3155 for (cs_lnum_t i = 0; i < send_list_size; i++) {
3156 cs_gnum_t b_id = orig_gnum[i] - bi.gnum_range[0];
3157
3158 cs_lnum_t s_id = block->index[b_id];
3159 cs_lnum_t e_id = block->index[b_id+1];
3160
3161 for (cs_lnum_t j = s_id; j < e_id; j++) {
3162 (send_inter_list[j]).vtx_gnum = block->vtx_glst[j];
3163 (send_inter_list[j]).curv_abs = block->abs_lst[j];
3164 }
3165
3166 } /* End of loop on elements to send */
3167
3168 BFT_FREE(orig_gnum);
3169
3170 /* Exchange buffers */
3171
3172 for (cs_lnum_t i = 0; i < send_list_size+1; i++) {
3173 block_index[i] *= sizeof(exch_inter_t);
3174 }
3175 for (cs_lnum_t i = 0; i < part->n_edges+1; i++) {
3176 part->index[i] *= sizeof(exch_inter_t);
3177 }
3178
3179 exch_inter_t *recv_inter_list
3180 = cs_all_to_all_copy_indexed(d,
3181 CS_CHAR,
3182 true, /* reverse */
3183 block_index,
3184 send_inter_list,
3185 part->index,
3186 NULL);
3187
3188 for (cs_lnum_t i = 0; i < part->n_edges+1; i++) {
3189 part->index[i] /= sizeof(exch_inter_t);
3190 }
3191
3192 /* Partial free memory */
3193
3194 BFT_FREE(send_inter_list);
3195 BFT_FREE(block_index);
3196
3197 cs_all_to_all_destroy(&d);
3198
3199 /* Update part definitions */
3200
3201 cs_lnum_t max_sub_size = 0;
3202
3203 for (cs_lnum_t i = 0; i < part->n_edges; i++) {
3204 cs_lnum_t nsub = part->index[i+1] - part->index[i];
3205 max_sub_size = CS_MAX(max_sub_size, nsub);
3206 }
3207
3208 part->max_sub_size = max_sub_size;
3209
3210 cs_lnum_t part_indexed_size = part->index[part->n_edges];
3211
3212 BFT_FREE(part->vtx_lst);
3213 part->vtx_lst = NULL;
3214 BFT_REALLOC(part->vtx_glst, part_indexed_size, cs_gnum_t);
3215 BFT_REALLOC(part->abs_lst, part_indexed_size, cs_coord_t);
3216
3217 for (cs_lnum_t i = 0; i < part_indexed_size; i++) {
3218 part->vtx_glst[i] = (recv_inter_list[i]).vtx_gnum;
3219 part->abs_lst[i] = (recv_inter_list[i]).curv_abs;
3220 }
3221
3222 /* Free memory */
3223
3224 BFT_FREE(recv_inter_list);
3225 }
3226
3227 #endif /* HAVE_MPI */
3228
3229 /*----------------------------------------------------------------------------
3230 * Redefine a cs_join_inter_edges_t structure to be consistent with the local
3231 * numbering of a given couple of cs_join_mesh_t structure and
3232 * cs_join_edges_t structure.
3233 * Add future new vertices for the face definition in cs_join_mesh_t
3234 *
3235 * parameters:
3236 * verbosity <-- verbosity level
3237 * edges <-- cs_join_edges_t structure
3238 * mesh <-> cs_join_mesh_t structure
3239 * inter_edges <-> current cs_join_inter_edges_t struct. to work with
3240 *---------------------------------------------------------------------------*/
3241
3242 void
cs_join_intersect_update_struct(int verbosity,const cs_join_edges_t * edges,cs_join_mesh_t * mesh,cs_join_inter_edges_t ** inter_edges)3243 cs_join_intersect_update_struct(int verbosity,
3244 const cs_join_edges_t *edges,
3245 cs_join_mesh_t *mesh,
3246 cs_join_inter_edges_t **inter_edges)
3247 {
3248 cs_lnum_t i, j, shift, o_id, max_size;
3249
3250 cs_lnum_t n_new_vertices = 0;
3251 cs_gnum_t *vtx_gnum = NULL, *edge_gnum = NULL;
3252 cs_lnum_t *edge_order = NULL, *vtx_order = NULL;
3253 cs_join_inter_edges_t *_inter_edges = *inter_edges;
3254 cs_join_inter_edges_t *new_inter_edges = NULL;
3255 cs_join_vertex_t *new_vertices = NULL;
3256
3257 assert(edges != NULL);
3258 assert(mesh != NULL);
3259
3260 const cs_lnum_t n_edges = edges->n_edges;
3261 const cs_lnum_t n_init_vertices = mesh->n_vertices;
3262
3263 assert(inter_edges != NULL);
3264 assert(n_edges == _inter_edges->n_edges);
3265
3266 /* Check if we have to re-order inter_edges */
3267
3268 for (i = 0; i < n_edges; i++)
3269 if (_inter_edges->edge_gnum[i] != edges->gnum[i])
3270 break;
3271
3272 if (i != n_edges) {
3273
3274 new_inter_edges = cs_join_inter_edges_create(n_edges);
3275
3276 /* Order global edge numbering */
3277
3278 BFT_MALLOC(edge_order, n_edges, cs_lnum_t);
3279 BFT_MALLOC(edge_gnum, n_edges, cs_gnum_t);
3280
3281 cs_order_gnum_allocated(NULL, edges->gnum, edge_order, n_edges);
3282
3283 for (i = 0; i < n_edges; i++)
3284 edge_gnum[i] = edges->gnum[edge_order[i]];
3285
3286 /* edge_gnum -> edge_id */
3287
3288 for (i = 0; i < n_edges; i++) {
3289
3290 cs_gnum_t e_gnum = _inter_edges->edge_gnum[i];
3291 cs_lnum_t e_id = cs_search_g_binary(n_edges, e_gnum, edge_gnum);
3292
3293 if (e_id == -1)
3294 bft_error(__FILE__, __LINE__, 0,
3295 _(" The received edge global number (%llu) is unknown"
3296 " on the current rank.\n"),
3297 (unsigned long long)_inter_edges->edge_gnum[i]);
3298
3299 o_id = edge_order[e_id];
3300 new_inter_edges->edge_gnum[o_id] = e_gnum;
3301
3302 new_inter_edges->index[o_id+1] =
3303 _inter_edges->index[i+1] - _inter_edges->index[i];
3304
3305 }
3306
3307 for (i = 0; i < n_edges; i++)
3308 new_inter_edges->index[i+1] += new_inter_edges->index[i];
3309
3310 BFT_MALLOC(new_inter_edges->vtx_glst,
3311 new_inter_edges->index[n_edges], cs_gnum_t);
3312 BFT_MALLOC(new_inter_edges->abs_lst,
3313 new_inter_edges->index[n_edges], cs_coord_t);
3314
3315 for (i = 0; i < n_edges; i++) {
3316
3317 cs_gnum_t e_gnum = _inter_edges->edge_gnum[i];
3318 cs_lnum_t e_id = cs_search_g_binary(n_edges, e_gnum, edge_gnum);
3319
3320 o_id = edge_order[e_id];
3321 shift = new_inter_edges->index[o_id];
3322
3323 for (j = _inter_edges->index[i]; j < _inter_edges->index[i+1]; j++) {
3324 new_inter_edges->vtx_glst[shift] = _inter_edges->vtx_glst[j];
3325 new_inter_edges->abs_lst[shift] = _inter_edges->abs_lst[j];
3326 shift++;
3327 }
3328
3329 }
3330
3331 /* Partial memory free */
3332
3333 BFT_FREE(edge_gnum);
3334 BFT_FREE(edge_order);
3335
3336 cs_join_inter_edges_destroy(&_inter_edges);
3337
3338 } /* End if re-order necessary */
3339
3340 else
3341 new_inter_edges = _inter_edges;
3342
3343 if (new_inter_edges->vtx_lst == NULL)
3344 BFT_MALLOC(new_inter_edges->vtx_lst,
3345 new_inter_edges->index[n_edges], cs_lnum_t);
3346
3347 /* Order global vertex numbering */
3348
3349 BFT_MALLOC(vtx_gnum, n_init_vertices, cs_gnum_t);
3350 BFT_MALLOC(vtx_order, n_init_vertices, cs_lnum_t);
3351
3352 for (i = 0; i < n_init_vertices; i++)
3353 vtx_gnum[i] = mesh->vertices[i].gnum;
3354
3355 cs_order_gnum_allocated(NULL, vtx_gnum, vtx_order, n_init_vertices);
3356
3357 for (i = 0; i < n_init_vertices; i++)
3358 vtx_gnum[i] = mesh->vertices[vtx_order[i]].gnum;
3359
3360 /* Pre-allocate a buffer to store data on possible new vertices */
3361
3362 max_size = 100;
3363 BFT_MALLOC(new_vertices, max_size, cs_join_vertex_t);
3364
3365 /* Fill vtx_lst array of the cs_join_inter_edges_t structure */
3366
3367 for (i = 0; i < n_edges; i++) {
3368
3369 cs_lnum_t start = new_inter_edges->index[i];
3370 cs_lnum_t end = new_inter_edges->index[i+1];
3371
3372 for (j = start; j < end; j++) {
3373
3374 cs_lnum_t id = cs_search_g_binary(n_init_vertices,
3375 new_inter_edges->vtx_glst[j],
3376 vtx_gnum);
3377
3378 if (id == -1) { /* New vertex to add in the mesh structure */
3379
3380 if (n_new_vertices >= max_size) {
3381 max_size *= 2;
3382 BFT_REALLOC(new_vertices, max_size, cs_join_vertex_t);
3383 }
3384
3385 new_vertices[n_new_vertices]
3386 = _get_new_vertex(new_inter_edges->abs_lst[j],
3387 new_inter_edges->vtx_glst[j],
3388 &(edges->def[2*i]),
3389 mesh);
3390
3391 /* update new_inter_edges */
3392
3393 n_new_vertices++;
3394 new_inter_edges->vtx_lst[j] = n_init_vertices + n_new_vertices;
3395
3396 }
3397 else
3398 new_inter_edges->vtx_lst[j] = vtx_order[id] + 1;
3399
3400 } /* End of loop on intersection description */
3401
3402 } /* End of loop on edges */
3403
3404 if (n_new_vertices > 0) {
3405
3406 if (verbosity > 2)
3407 fprintf(cs_glob_join_log,
3408 "\n Add %ld new vertices in the %s mesh definition"
3409 " during update of the edge definition.\n",
3410 (long)n_new_vertices, mesh->name);
3411
3412 BFT_REALLOC(mesh->vertices,
3413 n_init_vertices + n_new_vertices,
3414 cs_join_vertex_t);
3415
3416 for (i = 0; i < n_new_vertices; i++)
3417 mesh->vertices[n_init_vertices + i] = new_vertices[i];
3418
3419 mesh->n_vertices = n_init_vertices + n_new_vertices;
3420
3421 }
3422
3423 /* Free memory */
3424
3425 BFT_FREE(vtx_gnum);
3426 BFT_FREE(vtx_order);
3427 BFT_FREE(new_vertices);
3428
3429 /* Returns pointer */
3430
3431 *inter_edges = new_inter_edges;
3432 }
3433
3434 /*----------------------------------------------------------------------------
3435 * Get all real edges intersections among possible edges intersections.
3436 *
3437 * parameters:
3438 * param <-- set of user-defined parameters for the joining
3439 * edge_edge_vis <-- a pointer to a cs_join_gset_t structure
3440 * edges <-- pointer to a structure defining edges
3441 * mesh <-- pointer to the cs_join_mesh_t structure
3442 * which has the face connectivity
3443 * inter_set <-> pointer to a structure including data on edge
3444 * intersections
3445 * vtx_eset <-> pointer to a structure dealing with vertex
3446 * equivalences
3447 *
3448 * returns:
3449 * the type of joining encountered (conforming or not)
3450 *---------------------------------------------------------------------------*/
3451
3452 cs_join_type_t
cs_join_intersect_edges(cs_join_param_t param,const cs_join_gset_t * edge_edge_vis,const cs_join_edges_t * edges,const cs_join_mesh_t * mesh,cs_join_eset_t ** vtx_eset,cs_join_inter_set_t ** inter_set)3453 cs_join_intersect_edges(cs_join_param_t param,
3454 const cs_join_gset_t *edge_edge_vis,
3455 const cs_join_edges_t *edges,
3456 const cs_join_mesh_t *mesh,
3457 cs_join_eset_t **vtx_eset,
3458 cs_join_inter_set_t **inter_set)
3459 {
3460 cs_lnum_t i, j, k;
3461 double abs_e1[2], abs_e2[2];
3462
3463 cs_join_type_t join_type = CS_JOIN_TYPE_CONFORMING;
3464 cs_lnum_t n_inter = 0;
3465 cs_lnum_t n_inter_detected = 0, n_real_inter = 0, n_trivial_inter = 0;
3466 cs_gnum_t n_g_inter[3] = {0, 0, 0};
3467 cs_join_inter_set_t *_inter_set = NULL;
3468 cs_join_eset_t *_vtx_eset = NULL;
3469 FILE *logfile = cs_glob_join_log;
3470
3471 const double merge_limit = param.fraction * param.pre_merge_factor;
3472 const double parall_eps2 = 1e-6;
3473
3474 /* Sanity checks */
3475
3476 assert(mesh != NULL);
3477 assert(edges != NULL);
3478 assert(edge_edge_vis != NULL);
3479
3480 assert(vtx_eset != NULL);
3481 assert(inter_set != NULL);
3482
3483 if (param.verbosity > 3)
3484 fprintf(logfile, " Parallel intersection criterion: %8.5e\n",
3485 parall_eps2);
3486
3487 /* Initialization of structures */
3488
3489 _n_inter_tolerance_warnings = 0;
3490
3491 _inter_set = cs_join_inter_set_create(50);
3492 _vtx_eset = cs_join_eset_create(30);
3493
3494 /* Loop on edges */
3495
3496 for (i = 0; i < edge_edge_vis->n_elts; i++) {
3497
3498 int e1 = edge_edge_vis->g_elts[i]; /* This is a local number */
3499
3500 for (j = edge_edge_vis->index[i]; j < edge_edge_vis->index[i+1]; j++) {
3501
3502 int e2 = edge_edge_vis->g_list[j]; /* This is a local number */
3503 int e1_id = (e1 < e2 ? e1 - 1 : e2 - 1);
3504 int e2_id = (e1 < e2 ? e2 - 1 : e1 - 1);
3505
3506 assert(e1 != e2);
3507
3508 /* Get edge-edge intersection */
3509
3510 if (param.icm == 1)
3511 _edge_edge_3d_inter(mesh,
3512 edges,
3513 param.fraction,
3514 e1_id, abs_e1,
3515 e2_id, abs_e2,
3516 parall_eps2,
3517 param.verbosity,
3518 logfile,
3519 &n_inter);
3520
3521 else if (param.icm == 2)
3522 _new_edge_edge_3d_inter(mesh,
3523 edges,
3524 param.fraction,
3525 e1_id, abs_e1,
3526 e2_id, abs_e2,
3527 parall_eps2,
3528 param.verbosity,
3529 logfile,
3530 &n_inter);
3531
3532 n_inter_detected += n_inter;
3533
3534 #if 0 && defined(DEBUG) && !defined(NDEBUG)
3535 if (param.verbosity > 3 && n_inter > 0) {
3536
3537 cs_lnum_t v1e1 = edges->def[2*e1_id] - 1;
3538 cs_lnum_t v2e1 = edges->def[2*e1_id+1] - 1;
3539 cs_lnum_t v1e2 = edges->def[2*e2_id] - 1;
3540 cs_lnum_t v2e2 = edges->def[2*e2_id+1] - 1;
3541
3542 fprintf(logfile,
3543 "\n Intersection: "
3544 "E1 (%llu) [%llu - %llu] / E2 (%llu) [%llu - %llu]\n",
3545 (unsigned long long)edges->gnum[e1_id],
3546 (unsigned long long)mesh->vertices[v1e1].gnum,
3547 (unsigned long long)mesh->vertices[v2e1].gnum,
3548 (unsigned long long)edges->gnum[e2_id],
3549 (unsigned long long)mesh->vertices[v1e2].gnum,
3550 (unsigned long long)mesh->vertices[v2e2].gnum);
3551 fprintf(logfile, " n_inter: %d ", n_inter);
3552 for (k = 0; k < n_inter; k++)
3553 fprintf(logfile,
3554 " (%d) - s_e1 = %g, s_e2 = %g", k, abs_e1[k], abs_e2[k]);
3555 fflush(logfile);
3556 }
3557 #endif
3558
3559 for (k = 0; k < n_inter; k++) {
3560
3561 bool trivial = false;
3562
3563 if (abs_e1[k] <= merge_limit || abs_e1[k] >= 1.0 - merge_limit)
3564 if (abs_e2[k] <= merge_limit || abs_e2[k] >= 1.0 - merge_limit)
3565 trivial = true;
3566
3567 if (trivial) {
3568
3569 _add_trivial_equiv(e1_id,
3570 e2_id,
3571 abs_e1[k],
3572 abs_e2[k],
3573 edges,
3574 _vtx_eset);
3575
3576 n_trivial_inter += 1;
3577
3578 }
3579 else {
3580
3581 if (join_type == CS_JOIN_TYPE_CONFORMING)
3582 join_type = CS_JOIN_TYPE_NON_CONFORMING;
3583
3584 _add_inter(e1_id, e2_id, abs_e1[k], abs_e2[k], _inter_set);
3585
3586 }
3587
3588 } /* End of loop on detected edge_edge_vis */
3589
3590 } /* End of loop on entities intersecting elements */
3591
3592 } /* End of loop on elements in intersection list */
3593
3594 n_real_inter = n_inter_detected - n_trivial_inter;
3595
3596 if (n_inter_detected == 0)
3597 join_type = CS_JOIN_TYPE_NULL;
3598
3599 /* Order and delete redundant equivalences */
3600
3601 cs_join_eset_clean(&_vtx_eset);
3602
3603 /* Synchronize join_type and counts */
3604
3605 n_g_inter[0] = n_inter_detected;
3606 n_g_inter[1] = n_trivial_inter;
3607 n_g_inter[2] = n_real_inter;
3608
3609 #if defined(HAVE_MPI)
3610 if (cs_glob_n_ranks > 1) {
3611
3612 int tag = (int)join_type;
3613 int sync_tag = tag;
3614 cs_gnum_t tmp_inter[3];
3615
3616 MPI_Allreduce(&tag, &sync_tag, 1, MPI_INT, MPI_MAX, cs_glob_mpi_comm);
3617
3618 join_type = sync_tag;
3619
3620 MPI_Allreduce(n_g_inter, tmp_inter, 3, CS_MPI_GNUM, MPI_SUM,
3621 cs_glob_mpi_comm);
3622
3623 for (i = 0; i < 3; i++)
3624 n_g_inter[i] = tmp_inter[i];
3625 }
3626 #endif
3627
3628 if (param.verbosity > 0) {
3629
3630 bft_printf(_("\n"
3631 " Global number of intersections detected: %12llu\n"
3632 " Vertex-Vertex intersections: %12llu\n"
3633 " Other intersections: %12llu\n"),
3634 (unsigned long long)n_g_inter[0],
3635 (unsigned long long)n_g_inter[1],
3636 (unsigned long long)n_g_inter[2]);
3637
3638 if (param.verbosity > 2) {
3639 fprintf(logfile,
3640 "\n"
3641 " Local number of intersections detected: %10d\n"
3642 " Vertex-Vertex intersections: %10d\n"
3643 " Other intersections: %10d\n",
3644 (int)n_inter_detected, (int)n_trivial_inter,
3645 (int)n_real_inter);
3646 fprintf(logfile,
3647 "\n Local number of edge-edge intersection warnings: %9ld\n",
3648 (long)_n_inter_tolerance_warnings);
3649 fprintf(logfile,
3650 "\n Local number of equivalences between vertices: %9ld\n",
3651 (long)_vtx_eset->n_equiv);
3652 }
3653
3654 }
3655
3656 bft_printf_flush();
3657
3658 /* Return pointer */
3659
3660 *inter_set = _inter_set;
3661 *vtx_eset = _vtx_eset;
3662
3663 return join_type;
3664 }
3665
3666 /*----------------------------------------------------------------------------
3667 * Build a tree structure on which we associate leaves and face bounding boxes.
3668 * Create a cs_join_gset_t structure (indexed list on global numbering)
3669 * storing potential intersections between face bounding boxes.
3670 *
3671 * parameters:
3672 * param <-- set of user-defined parameters
3673 * join_mesh <-- cs_join_mesh_t structure where faces are defined
3674 * stats <-> joining statistics
3675 *
3676 * returns:
3677 * a new allocated pointer to a cs_join_gset_t structure storing the
3678 * face - face visibility.
3679 *---------------------------------------------------------------------------*/
3680
3681 cs_join_gset_t *
cs_join_intersect_faces(const cs_join_param_t param,const cs_join_mesh_t * join_mesh,cs_join_stats_t * stats)3682 cs_join_intersect_faces(const cs_join_param_t param,
3683 const cs_join_mesh_t *join_mesh,
3684 cs_join_stats_t *stats)
3685 {
3686 cs_lnum_t i;
3687
3688 int box_dim = 0;
3689 cs_coord_t *f_extents = NULL;
3690 fvm_neighborhood_t *face_neighborhood = NULL;
3691 cs_join_gset_t *face_visibility = NULL;
3692
3693 assert(join_mesh != NULL);
3694
3695 cs_timer_t t0 = cs_timer_time();
3696
3697 #if defined HAVE_MPI
3698 face_neighborhood = fvm_neighborhood_create(cs_glob_mpi_comm);
3699 #else
3700 face_neighborhood = fvm_neighborhood_create();
3701 #endif
3702
3703 fvm_neighborhood_set_options(face_neighborhood,
3704 param.tree_max_level,
3705 param.tree_n_max_boxes,
3706 param.tree_max_box_ratio,
3707 param.tree_max_box_ratio_distrib);
3708
3709 /* Allocate temporary extent arrays */
3710
3711 BFT_MALLOC(f_extents, join_mesh->n_faces*6, cs_coord_t);
3712
3713 /* Define each bounding box for the selected faces */
3714
3715 for (i = 0; i < join_mesh->n_faces; i++)
3716 _get_face_extents(join_mesh->face_vtx_idx[i],
3717 join_mesh->face_vtx_idx[i+1],
3718 join_mesh->face_vtx_lst,
3719 join_mesh->vertices,
3720 f_extents + i*6);
3721
3722 cs_timer_t t1 = cs_timer_time();
3723 cs_timer_counter_t extents_time = cs_timer_diff(&t0, &t1);
3724
3725 fvm_neighborhood_by_boxes(face_neighborhood,
3726 3, /* spatial dimension */
3727 join_mesh->n_faces,
3728 join_mesh->face_gnum,
3729 NULL,
3730 NULL,
3731 &f_extents);
3732
3733 _face_bbox_search_stats(face_neighborhood,
3734 extents_time,
3735 &box_dim,
3736 stats);
3737
3738 if (param.verbosity > 0) {
3739 bft_printf(_(" Determination of possible face intersections:\n\n"
3740 " bounding-box tree layout: %dD\n"), box_dim);
3741 bft_printf_flush();
3742 }
3743
3744 /* Retrieve face -> face visibility */
3745
3746 BFT_MALLOC(face_visibility, 1, cs_join_gset_t);
3747
3748 assert(sizeof(cs_lnum_t) == sizeof(cs_lnum_t));
3749
3750 fvm_neighborhood_transfer_data(face_neighborhood,
3751 &(face_visibility->n_elts),
3752 &(face_visibility->g_elts),
3753 &(face_visibility->index),
3754 &(face_visibility->g_list));
3755
3756 fvm_neighborhood_destroy(&face_neighborhood);
3757
3758 #if 0 && defined(DEBUG) && !defined(NDEBUG)
3759 {
3760 int len;
3761 FILE *dbg_file = NULL;
3762 char *filename = NULL;
3763
3764 len = strlen("JoinDBG_FaceVis.dat")+1+2+4;
3765 BFT_MALLOC(filename, len, char);
3766 sprintf(filename, "Join%02dDBG_FaceVis%04d.dat",
3767 param.num, CS_MAX(cs_glob_rank_id, 0));
3768 dbg_file = fopen(filename, "w");
3769
3770 cs_join_gset_dump(dbg_file, face_visibility);
3771
3772 fflush(dbg_file);
3773 BFT_FREE(filename);
3774 fclose(dbg_file);
3775 }
3776 #endif /* defined(DEBUG) && !defined(NDEBUG) */
3777
3778 return face_visibility;
3779 }
3780
3781 /*----------------------------------------------------------------------------
3782 * Transform face visibility into edge visibility (mesh->face_gnum must be
3783 * ordered).
3784 *
3785 * parameters:
3786 * mesh <-- pointer to a cs_join_mesh_t structure
3787 * edges <-- pointer to a cs_join_edges_t structure
3788 * face_visib <-- pointer to a cs_join_gset_t structure
3789 *
3790 * returns:
3791 * a new allocated cs_join_gset_t structure holding edge visibility
3792 *---------------------------------------------------------------------------*/
3793
3794 cs_join_gset_t *
cs_join_intersect_face_to_edge(const cs_join_mesh_t * mesh,const cs_join_edges_t * edges,const cs_join_gset_t * face_visib)3795 cs_join_intersect_face_to_edge(const cs_join_mesh_t *mesh,
3796 const cs_join_edges_t *edges,
3797 const cs_join_gset_t *face_visib)
3798 {
3799 cs_lnum_t i, j, k, edge_num, edge_id, shift;
3800
3801 assert(mesh != NULL);
3802 assert(edges != NULL);
3803 assert(face_visib != NULL);
3804
3805 cs_lnum_t size = 0, size_max = 0;
3806 cs_lnum_t *count = NULL, *face2edge_idx = NULL, *face2edge_lst = NULL;
3807 cs_gnum_t *tmp = NULL;
3808 cs_join_gset_t *edge_visib = NULL;
3809
3810 /* Create a local "face -> edge" connectivity
3811 First, count number of edges for each face */
3812
3813 BFT_MALLOC(face2edge_idx, mesh->n_faces + 1, cs_lnum_t);
3814
3815 face2edge_idx[0] = 0;
3816 for (i = 0; i < mesh->n_faces; i++)
3817 face2edge_idx[i+1] = mesh->face_vtx_idx[i+1] - mesh->face_vtx_idx[i];
3818
3819 for (i = 0; i < mesh->n_faces; i++)
3820 face2edge_idx[i+1] += face2edge_idx[i];
3821
3822 /* Build face2edge_lst */
3823
3824 BFT_MALLOC(face2edge_lst, face2edge_idx[mesh->n_faces], cs_lnum_t);
3825 BFT_MALLOC(count, mesh->n_faces, cs_lnum_t);
3826
3827 for (i = 0; i < mesh->n_faces; i++)
3828 count[i] = 0;
3829
3830 for (i = 0; i < mesh->n_faces; i++) {
3831
3832 cs_lnum_t start = mesh->face_vtx_idx[i];
3833 cs_lnum_t end = mesh->face_vtx_idx[i+1];
3834
3835 for (j = start; j < end - 1; j++) {
3836
3837 edge_num = cs_join_mesh_get_edge(mesh->face_vtx_lst[j] + 1,
3838 mesh->face_vtx_lst[j+1] + 1,
3839 edges);
3840
3841 shift = face2edge_idx[i] + count[i];
3842 count[i] += 1;
3843 face2edge_lst[shift] = CS_ABS(edge_num);
3844
3845 }
3846
3847 edge_num = cs_join_mesh_get_edge(mesh->face_vtx_lst[end-1] + 1,
3848 mesh->face_vtx_lst[start] + 1,
3849 edges);
3850
3851 shift = face2edge_idx[i] + count[i];
3852 count[i] += 1;
3853 face2edge_lst[shift] = CS_ABS(edge_num);
3854
3855 } /* End of loop on faces */
3856
3857 /* Transform numbering in face_visib to match numbering
3858 in the current cs_join_mesh_t structure */
3859
3860 for (i = 0; i < face_visib->n_elts; i++) {
3861
3862 cs_lnum_t face_id = cs_search_g_binary(mesh->n_faces,
3863 face_visib->g_elts[i],
3864 mesh->face_gnum);
3865
3866 face_visib->g_elts[i] = face_id;
3867
3868 for (j = face_visib->index[i]; j < face_visib->index[i+1]; j++) {
3869
3870 cs_lnum_t adj_id = cs_search_g_binary(mesh->n_faces,
3871 face_visib->g_list[j],
3872 mesh->face_gnum);
3873
3874 face_visib->g_list[j] = adj_id;
3875
3876 }
3877
3878 } /* End of loop on bounding boxes */
3879
3880 /* Create edge_visib. Unfold face_visib */
3881
3882 for (i = 0; i < face_visib->n_elts; i++) {
3883 j = face_visib->g_elts[i];
3884 size += face2edge_idx[j+1] - face2edge_idx[j];
3885 }
3886
3887 edge_visib = cs_join_gset_create(size);
3888
3889 edge_id = 0;
3890
3891 for (i = 0; i < face_visib->n_elts; i++) {
3892
3893 cs_lnum_t face_id = face_visib->g_elts[i];
3894 cs_lnum_t start = face2edge_idx[face_id];
3895 cs_lnum_t end = face2edge_idx[face_id+1];
3896
3897 size = 0;
3898 for (j = face_visib->index[i]; j < face_visib->index[i+1]; j++) {
3899
3900 cs_lnum_t adj_id = face_visib->g_list[j];
3901
3902 size += face2edge_idx[adj_id+1] - face2edge_idx[adj_id];
3903
3904 }
3905
3906 size_max = CS_MAX(size, size_max);
3907
3908 for (j = start; j < end; j++) {
3909
3910 edge_visib->g_elts[edge_id] = face2edge_lst[j];
3911 edge_visib->index[edge_id+1] = size;
3912 edge_id++;
3913
3914 }
3915
3916 } /* End of loop on bounding boxes */
3917
3918 assert(edge_id == edge_visib->n_elts);
3919
3920 /* Build index */
3921
3922 for (i = 0; i < edge_visib->n_elts; i++)
3923 edge_visib->index[i+1] += edge_visib->index[i];
3924
3925 BFT_MALLOC(edge_visib->g_list,
3926 edge_visib->index[edge_visib->n_elts],
3927 cs_gnum_t);
3928
3929 BFT_MALLOC(tmp, size_max, cs_gnum_t);
3930
3931 /* Build list */
3932
3933 edge_id = 0;
3934
3935 for (i = 0; i < face_visib->n_elts; i++) {
3936
3937 cs_lnum_t _count = 0;
3938 cs_lnum_t face_id = face_visib->g_elts[i];
3939 cs_lnum_t n_edges = face2edge_idx[face_id+1] - face2edge_idx[face_id];
3940 cs_lnum_t b_start = face_visib->index[i];
3941 cs_lnum_t b_end = face_visib->index[i+1];
3942
3943 /* Unfold face->edge connectivity for the current list of bounding boxes */
3944
3945 for (j = b_start; j < b_end; j++) {
3946
3947 cs_lnum_t adj_face_id = face_visib->g_list[j];
3948 cs_lnum_t n_adj_edges = face2edge_idx[adj_face_id+1]
3949 - face2edge_idx[adj_face_id];
3950
3951 shift = face2edge_idx[adj_face_id];
3952
3953 for (k = 0; k < n_adj_edges; k++)
3954 tmp[_count + k] = face2edge_lst[shift + k];
3955 _count += n_adj_edges;
3956
3957 }
3958
3959 for (j = 0; j < n_edges; j++, edge_id++) {
3960
3961 assert(_count == edge_visib->index[edge_id+1]-edge_visib->index[edge_id]);
3962
3963 shift = edge_visib->index[edge_id];
3964
3965 for (k = 0; k < _count; k++)
3966 edge_visib->g_list[shift + k] = tmp[k];
3967
3968 } /* End of loop on edges */
3969
3970 } /* End of loop on bounding boxes */
3971
3972 /* Free memory */
3973
3974 BFT_FREE(face2edge_idx);
3975 BFT_FREE(face2edge_lst);
3976 BFT_FREE(count);
3977 BFT_FREE(tmp);
3978
3979 /* Delete redundancies in g_elts, order g_elts and compact data */
3980
3981 cs_join_gset_merge_elts(edge_visib, 0); /* 0 = g_elts is not ordered */
3982
3983 /* Delete redundancies in g_list */
3984
3985 cs_join_gset_clean(edge_visib);
3986
3987 cs_join_gset_compress(edge_visib);
3988
3989 /* Return pointers */
3990
3991 return edge_visib;
3992 }
3993
3994 /*----------------------------------------------------------------------------
3995 * Dump a cs_join_inter_edges_t structure.
3996 *
3997 * parameters:
3998 * f <-- handle to output file
3999 * inter_edges <-- cs_join_inter_edges_t structure to dump
4000 * edges <-- list of edges
4001 * mesh <-- associated cs_join_mesh_t structure
4002 *---------------------------------------------------------------------------*/
4003
4004 void
cs_join_inter_edges_dump(FILE * f,const cs_join_inter_edges_t * inter_edges,const cs_join_edges_t * edges,const cs_join_mesh_t * mesh)4005 cs_join_inter_edges_dump(FILE *f,
4006 const cs_join_inter_edges_t *inter_edges,
4007 const cs_join_edges_t *edges,
4008 const cs_join_mesh_t *mesh)
4009 {
4010 int i, j, k;
4011
4012 fprintf(f, "\n Dump of a cs_join_inter_edges_t structure (%p)\n",
4013 (const void *)inter_edges);
4014
4015 if (inter_edges == NULL)
4016 return;
4017
4018 fprintf(f, " n_edges: %10ld\n", (long)inter_edges->n_edges);
4019 fprintf(f, " max_sub_size: %10ld\n\n", (long)inter_edges->max_sub_size);
4020
4021 for (i = 0; i < inter_edges->n_edges; i++) {
4022
4023 assert(edges != NULL);
4024 assert(mesh != NULL);
4025
4026 cs_lnum_t v1_num = edges->def[2*i];
4027 cs_lnum_t v2_num = edges->def[2*i+1];
4028 cs_gnum_t v1_gnum = (mesh->vertices[v1_num-1]).gnum;
4029 cs_gnum_t v2_gnum = (mesh->vertices[v2_num-1]).gnum;
4030 cs_lnum_t start = inter_edges->index[i];
4031 cs_lnum_t end = inter_edges->index[i+1];
4032
4033 fprintf(f, "\n%6ld: [%9llu] = (%7ld [%9llu] - %7ld [%9llu])\n",
4034 (long)i, (unsigned long long)edges->gnum[i],
4035 (long)v1_num, (unsigned long long)v1_gnum,
4036 (long)v2_num, (unsigned long long)v2_gnum);
4037
4038 fprintf(f, " n_sub_inter: %4ld - index : %7ld <-- %7ld\n",
4039 (long)end-start, (long)start, (long)end);
4040
4041 if (inter_edges->vtx_glst == NULL) {
4042
4043 for (j = start, k = 0; j < end; j++, k++)
4044 fprintf
4045 (f, " %7ld (%9ld) - (%7llu, %8.6e)\n",
4046 (long)k, (long)inter_edges->vtx_lst[j],
4047 (unsigned long long)mesh->vertices[inter_edges->vtx_lst[j]-1].gnum,
4048 inter_edges->abs_lst[j]);
4049
4050 }
4051 else {
4052
4053 for (j = start, k = 0; j < end; j++, k++)
4054 fprintf(f, " %9ld - (%7llu, %8.6e)\n",
4055 (long)k, (unsigned long long)inter_edges->vtx_glst[j],
4056 inter_edges->abs_lst[j]);
4057
4058 }
4059
4060 } /* End of loop on edge intersections */
4061
4062 fflush(f);
4063 }
4064
4065 /*---------------------------------------------------------------------------*/
4066
4067 END_C_DECLS
4068