1 /*============================================================================
2 * Manipulation of low-level structures for the joining operations
3 *===========================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *---------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <errno.h>
35 #include <stdio.h>
36 #include <string.h>
37 #include <math.h>
38
39 /*----------------------------------------------------------------------------
40 * Local headers
41 *---------------------------------------------------------------------------*/
42
43 #include "bft_mem.h"
44 #include "bft_printf.h"
45
46 #include "fvm_defs.h"
47 #include "fvm_io_num.h"
48
49 #include "cs_join_util.h"
50 #include "cs_file.h"
51 #include "cs_mesh.h"
52 #include "cs_order.h"
53 #include "cs_search.h"
54 #include "cs_sort.h"
55
56 /*----------------------------------------------------------------------------
57 * Header for the current file
58 *---------------------------------------------------------------------------*/
59
60 #include "cs_join_set.h"
61
62 /*---------------------------------------------------------------------------*/
63
64 BEGIN_C_DECLS
65
66 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
67
68 /*============================================================================
69 * Static global variables
70 *===========================================================================*/
71
72 int cs_glob_join_count = 0;
73 int cs_glob_n_joinings = 0;
74 cs_join_t **cs_glob_join_array = NULL;
75
76 FILE *cs_glob_join_log = NULL;
77
78 /*============================================================================
79 * Macro and type definitions
80 *===========================================================================*/
81
82 /* Directory name separator
83 (historically, '/' for Unix/Linux, '\' for Windows, ':' for Mac
84 but '/' should work for all on modern systems) */
85
86 #define DIR_SEPARATOR '/'
87
88 /*============================================================================
89 * Private function definitions
90 *===========================================================================*/
91
92 /*----------------------------------------------------------------------------
93 * Initialize a cs_join_param_t structure.
94 *
95 * parameters:
96 * join_num <-- number of the current joining operation
97 * fraction <-- value of the fraction parameter
98 * plane <-- value of the plane parameter
99 * perio_type <-- periodicity type (FVM_PERIODICITY_NULL if none)
100 * perio_matrix <-- periodicity transformation matrix
101 * verbosity <-- level of verbosity required
102 * visualization <-- level of visualization required
103 * preprocessing <-- is joining part of the preprocessing stage ?
104 *
105 * returns:
106 * a pointer to a cs_join_param_t structure
107 *---------------------------------------------------------------------------*/
108
109 static cs_join_param_t
_join_param_define(int join_num,float fraction,float plane,fvm_periodicity_type_t perio_type,double perio_matrix[3][4],int verbosity,int visualization,bool preprocessing)110 _join_param_define(int join_num,
111 float fraction,
112 float plane,
113 fvm_periodicity_type_t perio_type,
114 double perio_matrix[3][4],
115 int verbosity,
116 int visualization,
117 bool preprocessing)
118 {
119 double cplane;
120 cs_join_param_t param;
121
122 param.num = join_num;
123
124 /* Possible periodicity information */
125
126 param.perio_type = perio_type;
127
128 if (param.perio_type == FVM_PERIODICITY_NULL) {
129 int i, j;
130 for (i = 0; i < 3; i++) {
131 for (j = 0; j < 4; j++)
132 param.perio_matrix[i][j] = 0.;
133 }
134 }
135 else
136 memcpy(param.perio_matrix, perio_matrix, sizeof(double)*12);
137
138 /* geometric parameters */
139
140 /* parameter used to compute the tolerance associated to each vertex.
141 Also used for finding equivalent vertices during edge intersections */
142
143 param.fraction = fraction;
144
145 /* parameter used to judge if two faces are in the same plane (during
146 the face splitting) */
147
148 param.plane = plane;
149 cplane = cos(param.plane *acos(-1.0)/180.);
150 param.plane_criteria = cplane * cplane;
151
152 /* Coef. used to modify the tolerance associated to each vertex BEFORE the
153 merge operation.
154 If coef = 0.0 => no vertex merge
155 If coef < 1.0 => reduce vertex merge
156 If coef = 1.0 => no change
157 If coef > 1.0 => increase vertex merge */
158
159 param.merge_tol_coef = 1.0;
160
161 /* Coef. used to compute a limit under which two vertices are merged
162 before the merge step.
163 Default value: 1e-3; Should be small. [1e-4, 1e-2] */
164
165 param.pre_merge_factor = 0.05;
166
167 /* Maximum number of equivalence breaks */
168
169 param.n_max_equiv_breaks = 500;
170
171 /* Tolerance computation mode: tcm
172 1: (default) tol = min. edge length related to a vertex * fraction
173 2: tolerance is computed like in mode 1 with in addition, the
174 multiplication by a coef. which is equal to the max sin(e1, e2)
175 where e1 and e2 are two edges sharing the same vertex V for which
176 we want to compute the tolerance
177 11: like 1 but only in taking into account only the selected faces
178 12: like 2 but only in taking into account only the selected faces
179 */
180
181 param.tcm = 1;
182
183 /* Intersection computation mode: icm
184 1: (default) Original algorithm. Try to clip intersection on extremity
185 2: New intersection algorithm. Avoid to clip intersection on extremity
186 */
187
188 param.icm = 1;
189
190 /* Maximum number of sub-faces for an initial selected face
191 Default value: 200 */
192
193 param.max_sub_faces = 200;
194
195 /* Deepest level reachable during tree building
196 Default value: 30 */
197
198 param.tree_max_level = 30;
199
200 /* Max. number of boxes which can be related to a leaf of the tree
201 if level != tree_max_level
202 Default value: 25 */
203
204 param.tree_n_max_boxes = 25;
205
206 /* Stop tree building if: n_linked_boxes > tree_max_box_ratio*n_init_boxes
207 Default value: 5.0 */
208
209 param.tree_max_box_ratio = 5.0;
210 param.tree_max_box_ratio_distrib = 2.0;
211
212 /* Level of display */
213
214 param.verbosity = verbosity;
215 param.visualization = visualization;
216
217 /* Is joining part of preprocessing ? */
218
219 param.preprocessing = preprocessing;
220
221 return param;
222 }
223
224 /*----------------------------------------------------------------------------
225 * Initialize a cs_join_stats_t structure.
226 *
227 * parameters:
228 * join_num <-- number of the current joining operation
229 * fraction <-- value of the fraction parameter
230 * plane <-- value of the plane parameter
231 * perio_type <-- periodicity type (FVM_PERIODICITY_NULL if none)
232 * perio_matrix <-- periodicity transformation matrix
233 * verbosity <-- level of verbosity required
234 * visualization <-- level of visualization required
235 * preprocessing <-- is joining part of the preprocessing stage ?
236 *
237 * returns:
238 * a pointer to a cs_join_param_t structure
239 *---------------------------------------------------------------------------*/
240
241 static cs_join_stats_t
_join_stats_init(void)242 _join_stats_init(void)
243 {
244 cs_join_stats_t stats;
245
246 memset(&stats, 0, sizeof(cs_join_stats_t));
247
248 CS_TIMER_COUNTER_INIT(stats.t_box_build);
249 CS_TIMER_COUNTER_INIT(stats.t_box_query);
250 CS_TIMER_COUNTER_INIT(stats.t_inter_sort);
251
252 CS_TIMER_COUNTER_INIT(stats.t_l_join_mesh);
253 CS_TIMER_COUNTER_INIT(stats.t_edge_inter);
254 CS_TIMER_COUNTER_INIT(stats.t_new_vtx);
255 CS_TIMER_COUNTER_INIT(stats.t_merge_vtx);
256 CS_TIMER_COUNTER_INIT(stats.t_u_merge_vtx);
257 CS_TIMER_COUNTER_INIT(stats.t_split_faces);
258
259 CS_TIMER_COUNTER_INIT(stats.t_total);
260
261 return stats;
262 }
263
264 /*----------------------------------------------------------------------------
265 * Initialize a structure for the synchronization of single
266 * elements
267 *
268 * returns:
269 * a pointer to a new structure used for synchronizing single elements
270 *----------------------------------------------------------------------------*/
271
272 static cs_join_sync_t *
_create_join_sync(void)273 _create_join_sync(void)
274 {
275 cs_join_sync_t *sync = NULL;
276
277 BFT_MALLOC(sync, 1, cs_join_sync_t);
278
279 sync->n_elts = 0;
280 sync->n_ranks = 0;
281 sync->ranks = NULL;
282 sync->index = NULL;
283 sync->array = NULL;
284
285 return sync;
286 }
287
288 /*----------------------------------------------------------------------------
289 * Destroy a structure for the synchronization of single elements.
290 *
291 * parameters:
292 * sync <-> pointer to a structure used for synchronizing single elements
293 *----------------------------------------------------------------------------*/
294
295 static void
_destroy_join_sync(cs_join_sync_t ** sync)296 _destroy_join_sync(cs_join_sync_t **sync)
297 {
298 cs_join_sync_t *_sync = *sync;
299
300 if (_sync == NULL)
301 return;
302
303 if (_sync->array != NULL)
304 BFT_FREE(_sync->array);
305 if (_sync->ranks != NULL)
306 BFT_FREE(_sync->ranks);
307 BFT_FREE(_sync->index);
308
309 BFT_FREE(_sync);
310
311 *sync = _sync;
312 }
313
314 /*----------------------------------------------------------------------------
315 * Reduce numbering for the selected boundary faces.
316 * After this function, we have a compact global face numbering for the
317 * selected faces.
318 *
319 * parameters:
320 * n_select_faces <-- number of selected faces
321 * p_reduce_gnum <-> pointer to the reduced numbering
322 * p_reduce_gnum_index <-> pointer to an index on ranks for the reduce
323 * numbering
324 *---------------------------------------------------------------------------*/
325
326 static void
_compact_face_gnum_selection(cs_lnum_t n_select_faces,cs_gnum_t * reduce_gnum[],cs_gnum_t * reduce_gnum_index[])327 _compact_face_gnum_selection(cs_lnum_t n_select_faces,
328 cs_gnum_t *reduce_gnum[],
329 cs_gnum_t *reduce_gnum_index[])
330 {
331 cs_lnum_t i;
332
333 cs_gnum_t shift = 0;
334 cs_gnum_t *_reduce_gnum = *reduce_gnum;
335 cs_gnum_t *_reduce_gnum_index = *reduce_gnum_index;
336
337 const int n_ranks = cs_glob_n_ranks;
338 const int local_rank = CS_MAX(cs_glob_rank_id, 0);
339
340 assert(_reduce_gnum_index == NULL);
341
342 BFT_MALLOC(_reduce_gnum_index, n_ranks + 1, cs_gnum_t);
343
344 for (i = 0; i < n_ranks; i++)
345 _reduce_gnum_index[i] = 0;
346
347 if (n_ranks > 1) {
348 #if defined(HAVE_MPI)
349 MPI_Comm mpi_comm = cs_glob_mpi_comm;
350 cs_gnum_t _n_faces = n_select_faces;
351
352 MPI_Allgather(&_n_faces, 1, CS_MPI_GNUM,
353 &(_reduce_gnum_index[1]), 1, CS_MPI_GNUM, mpi_comm);
354 #endif
355
356 for (i = 0; i < n_ranks; i++)
357 _reduce_gnum_index[i+1] += _reduce_gnum_index[i];
358
359 shift = _reduce_gnum_index[local_rank];
360
361 }
362 else {
363
364 assert(n_ranks == 1);
365 _reduce_gnum_index[n_ranks] = (cs_gnum_t)n_select_faces;
366
367 }
368
369 BFT_MALLOC(_reduce_gnum, n_select_faces, cs_gnum_t);
370
371 for (i = 0; i < n_select_faces; i++)
372 _reduce_gnum[i] = shift + i + 1;
373
374 /* Returns pointer */
375
376 *reduce_gnum = _reduce_gnum;
377 *reduce_gnum_index = _reduce_gnum_index;
378 }
379
380 /*----------------------------------------------------------------------------
381 * Extract faces implied in the current joining operation.
382 * These are faces which share at least one vertex which is in the
383 * select_vertices array.
384 *
385 * parameters:
386 * n_vertices <-- number of vertices in the whole mesh
387 * selection <-- pointer to a selection structure
388 * n_faces <-- number of faces in the whole mesh
389 * f2v_idx <-- "face -> vertex" connect. index
390 * f2v_lst <-- "face -> vertex" connect. list
391 * n_contig_faces <-> pointer to the number of contiguous faces
392 * contig_faces <-> pointer to the list of contiguous faces
393 *---------------------------------------------------------------------------*/
394
395 static void
_extract_contig_faces(cs_lnum_t n_vertices,cs_join_select_t * selection,cs_lnum_t n_faces,const cs_lnum_t f2v_idx[],const cs_lnum_t f2v_lst[],cs_lnum_t * n_contig_faces,cs_lnum_t * contig_faces[])396 _extract_contig_faces(cs_lnum_t n_vertices,
397 cs_join_select_t *selection,
398 cs_lnum_t n_faces,
399 const cs_lnum_t f2v_idx[],
400 const cs_lnum_t f2v_lst[],
401 cs_lnum_t *n_contig_faces,
402 cs_lnum_t *contig_faces[])
403 {
404 cs_lnum_t i, j, vtx_id, shift;
405
406 cs_lnum_t _n_contig_faces = 0;
407 cs_lnum_t *_contig_faces = NULL, *counter = NULL;
408 cs_lnum_t *v2f_idx = NULL, *v2f_lst = NULL;
409
410 const cs_lnum_t n_select_vertices = selection->n_vertices;
411 const cs_lnum_t n_single_vertices = selection->s_vertices->n_elts;
412 const cs_lnum_t *select_vertices = selection->vertices;
413 const cs_lnum_t *single_vertices = selection->s_vertices->array;
414
415 if (n_select_vertices + n_single_vertices == 0)
416 return;
417
418 /* Reverse face -> vertex connectivity */
419
420 BFT_MALLOC(counter, n_vertices, cs_lnum_t);
421
422 for (i = 0; i < n_vertices; i++)
423 counter[i] = 0;
424
425 for (i = 0; i < n_faces; i++) {
426 for (j = f2v_idx[i]; j < f2v_idx[i+1]; j++) {
427 vtx_id = f2v_lst[j];
428 counter[vtx_id] += 1;
429 }
430 } /* End of loop on faces */
431
432 /* Define v2f_idx */
433
434 BFT_MALLOC(v2f_idx, n_vertices + 1, cs_lnum_t);
435
436 v2f_idx[0] = 0;
437 for (i = 0; i < n_vertices; i++)
438 v2f_idx[i+1] = v2f_idx[i] + counter[i];
439
440 for (i = 0; i < n_vertices; i++)
441 counter[i] = 0;
442
443 /* Define v2f_lst */
444
445 BFT_MALLOC(v2f_lst, v2f_idx[n_vertices], cs_lnum_t);
446
447 for (i = 0; i < n_faces; i++) {
448
449 for (j = f2v_idx[i]; j < f2v_idx[i+1]; j++) {
450
451 vtx_id = f2v_lst[j];
452 shift = v2f_idx[vtx_id] + counter[vtx_id];
453 v2f_lst[shift] = i+1;
454 counter[vtx_id] += 1;
455
456 }
457
458 } /* End of loop on faces */
459
460 BFT_REALLOC(counter, n_faces, cs_lnum_t);
461
462 for (i = 0; i < n_faces; i++)
463 counter[i] = 0;
464
465 /* Count the number of contiguous faces */
466
467 for (i = 0; i < n_select_vertices; i++) {
468
469 vtx_id = select_vertices[i] - 1;
470
471 for (j = v2f_idx[vtx_id]; j < v2f_idx[vtx_id+1]; j++)
472 counter[v2f_lst[j]-1] = 1;
473
474 }
475
476 for (i = 0; i < n_single_vertices; i++) {
477
478 vtx_id = single_vertices[i] - 1;
479
480 for (j = v2f_idx[vtx_id]; j < v2f_idx[vtx_id+1]; j++)
481 counter[v2f_lst[j]-1] = 1;
482
483 }
484
485 for (i = 0; i < n_faces; i++)
486 _n_contig_faces += counter[i];
487
488 /* Define contig_faces */
489
490 BFT_MALLOC(_contig_faces, _n_contig_faces, cs_lnum_t);
491
492 _n_contig_faces = 0;
493 for (i = 0; i < n_faces; i++) {
494 if (counter[i] == 1) {
495 _contig_faces[_n_contig_faces] = i+1;
496 _n_contig_faces += 1;
497 }
498 }
499
500 /* Free memory */
501
502 BFT_FREE(v2f_idx);
503 BFT_FREE(v2f_lst);
504 BFT_FREE(counter);
505
506 /* Return pointers */
507
508 *n_contig_faces = _n_contig_faces;
509 *contig_faces = _contig_faces;
510 }
511
512 #if defined(HAVE_MPI)
513
514 /*----------------------------------------------------------------------------
515 * Define a structure used for synchronizing "single" vertices.
516 * Use a cs_interface_t structure to help the build.
517 *
518 * parameters:
519 * interfaces <-- pointer to a cs_interface_set_t structure
520 * var_size <-- number of elements in var buffer
521 * count <-> counter buffer (0: not selected, 1 otherwise)
522 * related_ranks <-> rank associated to each single vertex (size: var_size)
523 * single <-> data about the distribution of single vertices
524 * verbosity <-- verbosity level
525 *----------------------------------------------------------------------------*/
526
527 static void
_add_single_vertices(cs_interface_set_t * interfaces,cs_lnum_t var_size,cs_lnum_t * count,int * related_ranks,cs_join_sync_t * single,int verbosity)528 _add_single_vertices(cs_interface_set_t *interfaces,
529 cs_lnum_t var_size,
530 cs_lnum_t *count,
531 int *related_ranks,
532 cs_join_sync_t *single,
533 int verbosity)
534 {
535 int distant_rank, n_interfaces, last_found_rank;
536 cs_lnum_t id, ii;
537
538 int n_max_ranks = 0;
539 cs_lnum_t count_size = 0, n_entities = 0, total_size = 0, n_max_elts = 0;
540 cs_lnum_t *buf = NULL, *recv_buf = NULL;
541
542 const cs_lnum_t *local_id = NULL;
543 const cs_interface_t *interface = NULL;
544
545 assert(count != NULL);
546 assert(single != NULL);
547
548 /* Initialize and allocate */
549
550 n_interfaces = cs_interface_set_size(interfaces);
551
552 total_size = cs_interface_set_n_elts(interfaces);
553
554 BFT_MALLOC(buf, total_size, cs_lnum_t);
555
556 /* Exchange with distant ranks */
557
558 cs_interface_set_copy_array(interfaces,
559 CS_LNUM_TYPE,
560 1,
561 true,
562 count,
563 buf);
564
565 /* Now we estimate the max. number of ranks and elements involved */
566
567 count_size = 0;
568 last_found_rank = -1;
569
570 for (id = 0; id < n_interfaces; id++) {
571
572 /* Scan data */
573
574 interface = cs_interface_set_get(interfaces, id);
575 distant_rank = cs_interface_rank(interface);
576 n_entities = cs_interface_size(interface);
577 local_id = cs_interface_get_elt_ids(interface);
578
579 recv_buf = buf + count_size;
580
581 for (ii = 0; ii < n_entities; ii++) {
582
583 int vtx_id = local_id[ii];
584
585 assert(vtx_id < var_size);
586
587 if (count[vtx_id] == 0 && recv_buf[ii] > 0) { /* Find a single vertex */
588
589 if (last_found_rank != distant_rank) {
590 last_found_rank = distant_rank;
591 n_max_ranks++;
592 }
593 n_max_elts++;
594
595 }
596
597 }
598
599 count_size += n_entities;
600
601 }
602
603 if (n_max_elts > 0) { /* We have found single vertices */
604
605 BFT_MALLOC(single->ranks, n_max_ranks, int);
606 BFT_MALLOC(single->index, n_max_ranks + 1, cs_lnum_t);
607 BFT_MALLOC(single->array, n_max_elts, cs_lnum_t);
608
609 count_size = 0;
610 last_found_rank = -1;
611 single->index[0] = 0;
612 single->n_elts = 0;
613 single->n_ranks = 0;
614
615 for (id = 0; id < n_interfaces; id++) {
616
617 /* Scan data */
618
619 interface = cs_interface_set_get(interfaces, id);
620 distant_rank = cs_interface_rank(interface);
621 n_entities = cs_interface_size(interface);
622 local_id = cs_interface_get_elt_ids(interface);
623
624 recv_buf = buf + count_size;
625
626 for (ii = 0; ii < n_entities; ii++) {
627
628 int vtx_id = local_id[ii];
629
630 if (count[vtx_id] == 0 && recv_buf[ii] > 0) {
631
632 if (last_found_rank != distant_rank) {
633 last_found_rank = distant_rank;
634 single->ranks[single->n_ranks++] = distant_rank;
635 }
636
637 single->array[single->n_elts++] = local_id[ii] + 1;
638 single->index[single->n_ranks] = single->n_elts;
639
640 related_ranks[vtx_id] = distant_rank;
641 count[vtx_id] = recv_buf[ii];
642
643 }
644
645 }
646
647 count_size += n_entities;
648
649 } /* End of loop on interfaces */
650
651 BFT_REALLOC(single->ranks, single->n_ranks, int);
652 BFT_REALLOC(single->index, single->n_ranks + 1, cs_lnum_t);
653 BFT_REALLOC(single->array, single->n_elts, cs_lnum_t);
654
655 } /* End if n_max_elts > 0 */
656
657 BFT_FREE(buf);
658
659 if (cs_glob_join_log != NULL && verbosity > 3) {
660 int i, j;
661 FILE *logfile = cs_glob_join_log;
662
663 fprintf(logfile,
664 "\n Single vertices for the joining operation: (%p)\n",
665 (void *)single);
666 fprintf(logfile,
667 " Single vertices: n_elts : %8ld\n"
668 " Single vertices: n_ranks: %8d\n",
669 (long)single->n_elts, single->n_ranks);
670
671 if (single->n_elts > 0) {
672 for (i = 0; i < single->n_ranks; i++)
673 for (j = single->index[i]; j < single->index[i+1]; j++)
674 fprintf(logfile, " %9ld | %6d | %9ld\n",
675 (long)j, single->ranks[i], (long)single->array[j]);
676 fprintf(logfile, "\n");
677 }
678 fflush(logfile);
679 }
680
681 }
682
683 /*-----------------------------------------------------------------------------
684 * Define a structure used for synchronizing "single" vertices.
685 * Use a cs_interface_t structure to help the build.
686 *
687 * parameters:
688 * interfaces --> pointer to a cs_interface_set_t structure
689 * var_size --> number of elements in var buffer
690 * related_ranks <-> rank buffer for synchronization
691 * coupled <-> pointer to a structure to build on coupled vertices
692 * verbosity <-- verbosity level
693 *----------------------------------------------------------------------------*/
694
695 static void
_add_coupled_vertices(cs_interface_set_t * interfaces,cs_lnum_t var_size,int * related_ranks,cs_join_sync_t * coupled,int verbosity)696 _add_coupled_vertices(cs_interface_set_t *interfaces,
697 cs_lnum_t var_size,
698 int *related_ranks,
699 cs_join_sync_t *coupled,
700 int verbosity)
701 {
702 int distant_rank, n_interfaces, last_found_rank;
703 cs_lnum_t id, ii;
704
705 cs_lnum_t n_entities = 0, count_size = 0, total_size = 0;
706 int *buf = NULL, *recv_buf = NULL;
707
708 cs_datatype_t int_type = (sizeof(int) == 8) ? CS_INT64 : CS_INT32;
709
710 const int local_rank = CS_MAX(cs_glob_rank_id, 0);
711 const cs_lnum_t *local_id = NULL;
712 const cs_interface_t *interface = NULL;
713
714 assert(related_ranks != NULL);
715 assert(coupled != NULL);
716
717 /* Initialize and allocate */
718
719 n_interfaces = cs_interface_set_size(interfaces);
720
721 n_interfaces = cs_interface_set_size(interfaces);
722
723 total_size = cs_interface_set_n_elts(interfaces);
724
725 BFT_MALLOC(buf, total_size, int);
726
727 /* Exchange with distant ranks */
728
729 cs_interface_set_copy_array(interfaces,
730 int_type,
731 1,
732 true,
733 related_ranks,
734 buf);
735
736 /* Now we scan each part to build coupled */
737
738 count_size = 0;
739 coupled->n_elts = 0;
740 coupled->n_ranks = 0;
741 last_found_rank = -1;
742
743 for (id = 0; id < n_interfaces; id++) {
744
745 /* Scan data */
746
747 interface = cs_interface_set_get(interfaces, id);
748 distant_rank = cs_interface_rank(interface);
749 n_entities = cs_interface_size(interface);
750
751 recv_buf = buf + count_size;
752
753 for (ii = 0; ii < n_entities; ii++) {
754
755 if (recv_buf[ii] == local_rank) { /* Coupled vertex found */
756 coupled->n_elts++;
757 if (last_found_rank != distant_rank) {
758 last_found_rank = distant_rank;
759 coupled->n_ranks++;
760 }
761 }
762
763 }
764 count_size += n_entities;
765
766 }
767
768 if (coupled->n_elts > 0) {
769
770 int rank_shift = 0, vtx_shift = 0;
771
772 BFT_MALLOC(coupled->array, coupled->n_elts, cs_lnum_t);
773 BFT_MALLOC(coupled->index, coupled->n_ranks + 1, cs_lnum_t);
774 BFT_MALLOC(coupled->ranks, coupled->n_ranks, int);
775
776 coupled->index[0] = 0;
777
778 count_size = 0;
779 last_found_rank = -1;
780
781 for (id = 0; id < n_interfaces; id++) {
782
783 /* Retrieve data */
784
785 interface = cs_interface_set_get(interfaces, id);
786 distant_rank = cs_interface_rank(interface);
787 n_entities = cs_interface_size(interface);
788 local_id = cs_interface_get_elt_ids(interface);
789
790 recv_buf = buf + count_size;
791
792 for (ii = 0; ii < n_entities; ii++) {
793
794 if (recv_buf[ii] == local_rank) {
795
796 if (last_found_rank != distant_rank) {
797 last_found_rank = distant_rank;
798 coupled->ranks[rank_shift++] = distant_rank;
799 }
800
801 assert(local_id[ii] < var_size);
802 coupled->array[vtx_shift++] = local_id[ii] + 1;
803 coupled->index[rank_shift] = vtx_shift;
804 }
805
806 }
807 count_size += n_entities;
808
809 }
810
811 } /* End if coupled->n_elts > 0 */
812
813 BFT_FREE(buf);
814
815 if (cs_glob_join_log != NULL && verbosity > 3) {
816
817 int i, j;
818 FILE *logfile = cs_glob_join_log;
819
820 fprintf(logfile, "\n Coupled vertices for the joining operation: (%p)\n",
821 (void *)coupled);
822 fprintf(logfile,
823 " Coupled vertices: n_elts : %8ld\n"
824 " Coupled vertices: n_ranks: %8d\n",
825 (long)coupled->n_elts, coupled->n_ranks);
826
827 if (coupled->n_elts > 0) {
828 for (i = 0; i < coupled->n_ranks; i++)
829 for (j = coupled->index[i]; j < coupled->index[i+1]; j++)
830 fprintf(logfile, " %9ld | %6d | %9ld\n",
831 (long)j, coupled->ranks[i], (long)coupled->array[j]);
832 fprintf(logfile, "\n");
833 }
834 fflush(logfile);
835
836 }
837
838 }
839
840 /*----------------------------------------------------------------------------
841 * Get the full selection of vertices to extract. Only in parallel. Somme
842 * vertices may have been selected on another rank and not on the local rank
843 * but you have to take them into account to have a good update of the mesh.
844 *
845 * parameters:
846 * n_vertices <-- number of vertices in the parent mesh
847 * ifs <-- pointer to a cs_interface_set_t struct.
848 * p_vtx_tag <-> pointer to an array on vertices. tag=1 if selected
849 * join_select <-> pointer to a fvm_join_selection_t structure
850 * verbosity <-- verbosity level
851 *---------------------------------------------------------------------------*/
852
853 static void
_get_missing_vertices(cs_lnum_t n_vertices,cs_interface_set_t * ifs,cs_lnum_t * p_vtx_tag[],cs_join_select_t * selection,int verbosity)854 _get_missing_vertices(cs_lnum_t n_vertices,
855 cs_interface_set_t *ifs,
856 cs_lnum_t *p_vtx_tag[],
857 cs_join_select_t *selection,
858 int verbosity)
859 {
860 cs_lnum_t i;
861 cs_gnum_t n_l_elts, n_g_elts;
862
863 cs_lnum_t *vtx_tag = NULL;
864 int *related_ranks = NULL;
865
866 /* Define a counter on vertices. 1 if selected, 0 otherwise */
867
868 BFT_MALLOC(vtx_tag, n_vertices, cs_lnum_t);
869 BFT_MALLOC(related_ranks, n_vertices, int);
870
871 for (i = 0; i < n_vertices; i++) {
872 vtx_tag[i] = 0;
873 related_ranks[i] = -1;
874 }
875
876 for (i = 0; i < selection->n_vertices; i++)
877 vtx_tag[selection->vertices[i] - 1] = 1;
878
879 /* Has to be done before the search of single edges because
880 we need a synchronized vtx_tag */
881
882 _add_single_vertices(ifs,
883 n_vertices,
884 vtx_tag,
885 related_ranks,
886 selection->s_vertices,
887 verbosity);
888
889 n_l_elts = selection->s_vertices->n_elts;
890
891 MPI_Allreduce(&n_l_elts, &n_g_elts, 1, CS_MPI_GNUM,
892 MPI_SUM, cs_glob_mpi_comm);
893
894 if (n_g_elts > 0) {
895
896 bft_printf(_("\n Global number of single vertices found: %6llu\n"),
897 (unsigned long long)n_g_elts);
898 bft_printf_flush();
899 selection->do_single_sync = true;
900
901 _add_coupled_vertices(ifs,
902 n_vertices,
903 related_ranks,
904 selection->c_vertices,
905 verbosity);
906
907 } /* End if n_g_elts > 0 */
908
909 /* Return pointers */
910
911 *p_vtx_tag = vtx_tag;
912
913 /* Free memory */
914
915 BFT_FREE(related_ranks);
916 }
917
918 /*----------------------------------------------------------------------------
919 * Define a vertex -> vertex connectivity for vertices belonging to the
920 * selected boundary faces.
921 *
922 * parameters:
923 * n_vertices <-- number of vertices in the parent mesh
924 * selection <-- pointer to a fvm_join_selection_t structure
925 * b_f2v_idx <-- boundary "face -> vertex" connect. index
926 * b_f2v_lst <-- boundary "face -> vertex" connect. list
927 * p_v2v_idx <-> vertex -> vertex connect. index
928 * p_v2v_lst <-> vertex -> vertex connect. list
929 *---------------------------------------------------------------------------*/
930
931 static void
_get_select_v2v_connect(cs_lnum_t n_vertices,cs_join_select_t * selection,cs_lnum_t b_f2v_idx[],cs_lnum_t b_f2v_lst[],cs_lnum_t * p_v2v_idx[],cs_lnum_t * p_v2v_lst[])932 _get_select_v2v_connect(cs_lnum_t n_vertices,
933 cs_join_select_t *selection,
934 cs_lnum_t b_f2v_idx[],
935 cs_lnum_t b_f2v_lst[],
936 cs_lnum_t *p_v2v_idx[],
937 cs_lnum_t *p_v2v_lst[])
938 {
939 cs_lnum_t i, j, save, s, e, n_sel_edges, shift;
940
941 cs_lnum_t *count = NULL, *sel_v2v_idx = NULL, *sel_v2v_lst = NULL;
942
943 /* Build a vertex -> vertex connectivity for the selected boundary faces */
944
945 BFT_MALLOC(sel_v2v_idx, n_vertices + 1, cs_lnum_t);
946
947 for (i = 0; i < n_vertices + 1; i++)
948 sel_v2v_idx[i] = 0;
949
950 cs_join_build_edges_idx(selection->n_faces,
951 selection->faces,
952 b_f2v_idx,
953 b_f2v_lst,
954 sel_v2v_idx);
955
956 BFT_MALLOC(count, n_vertices, cs_lnum_t);
957
958 for (i = 0; i < n_vertices; i++) {
959 sel_v2v_idx[i+1] += sel_v2v_idx[i];
960 count[i] = 0;
961 }
962
963 BFT_MALLOC(sel_v2v_lst, sel_v2v_idx[n_vertices], cs_lnum_t);
964
965 cs_join_build_edges_lst(selection->n_faces,
966 selection->faces,
967 b_f2v_idx,
968 b_f2v_lst,
969 count,
970 sel_v2v_idx,
971 sel_v2v_lst);
972
973 BFT_FREE(count);
974
975 /* Order sub-lists and then clean repeated elements */
976
977 for (i = 0; i < n_vertices; i++)
978 cs_sort_shell(sel_v2v_idx[i], sel_v2v_idx[i+1], sel_v2v_lst);
979
980 save = sel_v2v_idx[0];
981 shift = 0;
982
983 for (i = 0; i < n_vertices; i++) {
984
985 s = save;
986 e = sel_v2v_idx[i+1];
987
988 if (e - s > 0) {
989 sel_v2v_lst[shift++] = sel_v2v_lst[s];
990 for (j = s + 1; j < e; j++)
991 if (sel_v2v_lst[j-1] != sel_v2v_lst[j])
992 sel_v2v_lst[shift++] = sel_v2v_lst[j];
993 }
994
995 save = e;
996 sel_v2v_idx[i+1] = shift;
997
998 }
999
1000 n_sel_edges = sel_v2v_idx[n_vertices];
1001 BFT_REALLOC(sel_v2v_lst, n_sel_edges, cs_lnum_t);
1002
1003 /* Return pointers */
1004
1005 *p_v2v_idx = sel_v2v_idx;
1006 *p_v2v_lst = sel_v2v_lst;
1007
1008 }
1009
1010 /*----------------------------------------------------------------------------
1011 * Get the related edge id from a couple of vertex ids.
1012 *
1013 * parameters:
1014 * v1_id <-- first vertex id
1015 * v2_id <-- second vertex id
1016 * v2v_idx <-- vertex -> vertex connect. index
1017 * v2v_lst <-- vertex -> vertex connect. list
1018 *
1019 * returns:
1020 * related edge_id in cs_join_edges_t structure
1021 *---------------------------------------------------------------------------*/
1022
1023 inline static cs_lnum_t
_get_edge_id(cs_lnum_t v1_id,cs_lnum_t v2_id,const cs_lnum_t v2v_idx[],const cs_lnum_t v2v_lst[])1024 _get_edge_id(cs_lnum_t v1_id,
1025 cs_lnum_t v2_id,
1026 const cs_lnum_t v2v_idx[],
1027 const cs_lnum_t v2v_lst[])
1028 {
1029 int i, va, vb;
1030
1031 cs_lnum_t edge_id = -1;
1032
1033 if (v1_id < v2_id)
1034 va = v1_id, vb = v2_id;
1035 else
1036 vb = v1_id, va = v2_id;
1037
1038 for (i = v2v_idx[va]; i < v2v_idx[va+1]; i++)
1039 if (v2v_lst[i] == vb + 1)
1040 break;
1041
1042 if (i != v2v_idx[va+1])
1043 edge_id = i;
1044
1045 return edge_id;
1046 }
1047
1048 /*----------------------------------------------------------------------------
1049 * Get the full selection of faces related to "single" vertices.
1050 * Done only if the run is parallel.
1051 *
1052 * parameters:
1053 * vertex_tag <-- tag to know if a vertices is in selection
1054 * v1_id <-- first vertex id
1055 * v2_id <-- second vertex id
1056 * ref_v2v_idx <-- vertex -> vertex connect. index
1057 * ref_v2v_lst <-- vertex -> vertex connect. list
1058 * p_tmp_size <-> pointer to the current number of single edges
1059 * p_max_size <-> pointer to the max allocated size of new_s_vertices
1060 * p_tmp_size <-> pointer to the single edges definition
1061 *---------------------------------------------------------------------------*/
1062
1063 static void
_add_s_edge(cs_lnum_t vertex_tag[],cs_lnum_t v1_id,cs_lnum_t v2_id,const cs_lnum_t sel_v2v_idx[],const cs_lnum_t sel_v2v_lst[],cs_lnum_t * p_tmp_size,cs_lnum_t * p_max_size,cs_lnum_t * p_tmp_edges[])1064 _add_s_edge(cs_lnum_t vertex_tag[],
1065 cs_lnum_t v1_id,
1066 cs_lnum_t v2_id,
1067 const cs_lnum_t sel_v2v_idx[],
1068 const cs_lnum_t sel_v2v_lst[],
1069 cs_lnum_t *p_tmp_size,
1070 cs_lnum_t *p_max_size,
1071 cs_lnum_t *p_tmp_edges[])
1072 {
1073 cs_lnum_t i, a, b, edge_id;
1074
1075 if (vertex_tag[v1_id] > 0 && vertex_tag[v2_id] > 0) {
1076
1077 bool is_found = false;
1078
1079 cs_lnum_t tmp_size = *p_tmp_size;
1080 cs_lnum_t max_size = *p_max_size;
1081 cs_lnum_t *tmp_edges = *p_tmp_edges;
1082
1083 edge_id = _get_edge_id(v1_id, v2_id, sel_v2v_idx, sel_v2v_lst);
1084
1085 if (edge_id == -1) { /* Edge not found among the selected edges */
1086
1087 assert(v1_id != v2_id);
1088 if (v1_id < v2_id)
1089 a = v1_id, b = v2_id;
1090 else
1091 a = v2_id, b = v1_id;
1092
1093 /* n_edges is assumed to be small */
1094 for (i = 0; i < tmp_size; i++) {
1095 if (tmp_edges[2*i] == a) {
1096 if (tmp_edges[2*i+1] == b) {
1097 is_found = true;
1098 break;
1099 }
1100 }
1101 }
1102
1103 if (!is_found) {
1104
1105 tmp_edges[2*tmp_size] = a;
1106 tmp_edges[2*tmp_size+1] = b;
1107 tmp_size += 1;
1108
1109 if (tmp_size >= max_size) {
1110 max_size *= 2;
1111 BFT_REALLOC(tmp_edges, 2*max_size, cs_lnum_t);
1112 }
1113
1114 }
1115
1116 } /* this edge should be among the selected edges */
1117
1118 /* Return pointers */
1119
1120 *p_max_size = max_size;
1121 *p_tmp_size = tmp_size;
1122 *p_tmp_edges = tmp_edges;
1123
1124 }
1125 }
1126
1127 /*----------------------------------------------------------------------------
1128 * Copy selected edge definitions at parallel interfaces (only in parallel).
1129 *
1130 * The returned arrays (sel_v2v_recv_idx and sel_v2v_recv_lst) both use
1131 * 0 to n-1 numbering.
1132
1133 * The caller is responsible for freeing those returned arrays.
1134 *
1135 * parameters:
1136 * ifs <-- pointer to the interface set on vertices
1137 * sel_v2v_idx <-- vertex -> vertex connect. index
1138 * sel_v2v_lst <-- vertex -> vertex connect. list
1139 * sel_v2v_recv_idx --> pointer to received edges index
1140 * sel_v2v_recv_lst --> pointer to received edges connectivity
1141 *---------------------------------------------------------------------------*/
1142
1143 static void
_copy_interface_edges(cs_interface_set_t * ifs,cs_lnum_t sel_v2v_idx[],cs_lnum_t sel_v2v_lst[],cs_lnum_t ** sel_v2v_recv_idx,cs_lnum_t ** sel_v2v_recv_lst)1144 _copy_interface_edges(cs_interface_set_t *ifs,
1145 cs_lnum_t sel_v2v_idx[],
1146 cs_lnum_t sel_v2v_lst[],
1147 cs_lnum_t **sel_v2v_recv_idx,
1148 cs_lnum_t **sel_v2v_recv_lst)
1149 {
1150 int i;
1151 cs_lnum_t j, if_shift;
1152
1153 /* Exchange v2v info;
1154 only info relative to edges where both vertices are available
1155 on the distant side is sent. */
1156
1157 int n_interfaces = cs_interface_set_size(ifs);
1158
1159 cs_lnum_t ifs_tot_size = cs_interface_set_n_elts(ifs);
1160
1161 cs_lnum_t *_sel_v2v_send_idx = NULL, *_sel_v2v_recv_idx = NULL;
1162 cs_lnum_t *_sel_v2v_send_lst = NULL, *_sel_v2v_recv_lst = NULL;
1163
1164 BFT_MALLOC(_sel_v2v_send_idx, ifs_tot_size + 1, cs_lnum_t);
1165 BFT_MALLOC(_sel_v2v_recv_idx, ifs_tot_size + 1, cs_lnum_t);
1166
1167 /* Counting pass for v2v_info exchange
1168 (send/receive indexes are prepared as counts) */
1169
1170 if_shift = 0;
1171 _sel_v2v_send_idx[0] = 0;
1172 _sel_v2v_recv_idx[0] = 0;
1173
1174 for (i = 0; i < n_interfaces; i++) {
1175
1176 const cs_interface_t *interface = cs_interface_set_get(ifs, i);
1177 const cs_lnum_t *local_id = cs_interface_get_elt_ids(interface);
1178 const cs_lnum_t if_size = cs_interface_size(interface);
1179
1180 for (j = 0; j < if_size; j++) {
1181
1182 cs_lnum_t k = local_id[j];
1183 cs_lnum_t s_id = sel_v2v_idx[k];
1184 cs_lnum_t e_id = sel_v2v_idx[k+1];
1185
1186 _sel_v2v_send_idx[if_shift + j + 1] = 0;
1187
1188 for (cs_lnum_t l = s_id; l < e_id; l++) {
1189 if (cs_search_binary(if_size, sel_v2v_lst[l] - 1, local_id) > -1)
1190 _sel_v2v_send_idx[if_shift + j + 1] += 1;
1191 }
1192
1193 }
1194
1195 if_shift += if_size;
1196
1197 }
1198
1199 cs_interface_set_copy_array(ifs,
1200 CS_LNUM_TYPE,
1201 1,
1202 false,
1203 _sel_v2v_send_idx + 1,
1204 _sel_v2v_recv_idx + 1);
1205
1206 for (j = 0; j < ifs_tot_size; j++) {
1207 _sel_v2v_send_idx[j+1] += _sel_v2v_send_idx[j];
1208 _sel_v2v_recv_idx[j+1] += _sel_v2v_recv_idx[j];
1209 }
1210
1211 /* Data pass for v2v_info exchange */
1212
1213 cs_interface_set_add_match_ids(ifs);
1214
1215 BFT_MALLOC(_sel_v2v_send_lst, _sel_v2v_send_idx[ifs_tot_size], cs_lnum_t);
1216 BFT_MALLOC(_sel_v2v_recv_lst, _sel_v2v_recv_idx[ifs_tot_size], cs_lnum_t);
1217
1218 if_shift = 0;
1219
1220 for (i = 0; i < n_interfaces; i++) {
1221
1222 const cs_interface_t *interface = cs_interface_set_get(ifs, i);
1223 const cs_lnum_t if_size = cs_interface_size(interface);
1224
1225 const cs_lnum_t *local_id = cs_interface_get_elt_ids(interface);
1226 const cs_lnum_t *distant_id = cs_interface_get_match_ids(interface);
1227
1228 for (j = 0; j < if_size; j++) {
1229
1230 cs_lnum_t k = local_id[j];
1231 cs_lnum_t s_id = sel_v2v_idx[k];
1232 cs_lnum_t e_id = sel_v2v_idx[k+1];
1233 cs_lnum_t v_send_size = 0;
1234
1235 for (cs_lnum_t l = s_id; l < e_id; l++) {
1236 cs_lnum_t m = cs_search_binary(if_size, sel_v2v_lst[l] - 1, local_id);
1237 if (m > -1) {
1238 _sel_v2v_send_lst[_sel_v2v_send_idx[if_shift + j] + v_send_size]
1239 = distant_id[m];
1240 v_send_size += 1;
1241 }
1242 }
1243
1244 }
1245
1246 if_shift += if_size;
1247
1248 }
1249
1250 cs_interface_set_copy_indexed(ifs,
1251 CS_LNUM_TYPE,
1252 false,
1253 _sel_v2v_send_idx,
1254 _sel_v2v_recv_idx,
1255 _sel_v2v_send_lst,
1256 _sel_v2v_recv_lst);
1257
1258 /* Free temporary data and set return pointers */
1259
1260 BFT_FREE(_sel_v2v_send_idx);
1261 BFT_FREE(_sel_v2v_send_lst);
1262
1263 cs_interface_set_free_match_ids(ifs);
1264
1265 *sel_v2v_recv_idx = _sel_v2v_recv_idx;
1266 *sel_v2v_recv_lst = _sel_v2v_recv_lst;
1267 }
1268
1269 /*----------------------------------------------------------------------------
1270 * Get the full selection of single edges. Done only if the run is parallel.
1271 *
1272 * parameters:
1273 * ifs <-- pointer to the interface set on vertices
1274 * vertex_tag <-- tag to know if a vertices is in selection
1275 * selection <-- pointer to a fvm_join_selection_t structure
1276 * sel_v2v_idx <-- vertex -> vertex connect. index
1277 * sel_v2v_lst <-- vertex -> vertex connect. list
1278 * b_f2v_idx <-- boundary "face -> vertex" connect. index
1279 * b_f2v_lst <-- boundary "face -> vertex" connect. list
1280 * i_f2v_idx <-- interior "face -> vertex" connect. index
1281 * i_f2v_lst <-- interior "face -> vertex" connect. list
1282 * i_face_cells <-- interior face -> cells connect.
1283 * s_edges <-> pointer to the single edges structure to define
1284 * verbosity <-- verbosity level
1285 *---------------------------------------------------------------------------*/
1286
1287 static void
_add_single_edges(cs_interface_set_t * ifs,cs_lnum_t vertex_tag[],cs_join_select_t * selection,cs_lnum_t sel_v2v_idx[],cs_lnum_t sel_v2v_lst[],cs_lnum_t b_f2v_idx[],cs_lnum_t b_f2v_lst[],cs_lnum_t i_f2v_idx[],cs_lnum_t i_f2v_lst[],cs_lnum_2_t i_face_cells[],cs_join_sync_t * s_edges,int verbosity)1288 _add_single_edges(cs_interface_set_t *ifs,
1289 cs_lnum_t vertex_tag[],
1290 cs_join_select_t *selection,
1291 cs_lnum_t sel_v2v_idx[],
1292 cs_lnum_t sel_v2v_lst[],
1293 cs_lnum_t b_f2v_idx[],
1294 cs_lnum_t b_f2v_lst[],
1295 cs_lnum_t i_f2v_idx[],
1296 cs_lnum_t i_f2v_lst[],
1297 cs_lnum_2_t i_face_cells[],
1298 cs_join_sync_t *s_edges,
1299 int verbosity)
1300 {
1301 int i;
1302 cs_lnum_t j, fid, s, e;
1303
1304 cs_lnum_t tmp_size = 0, max_size = 10;
1305 cs_lnum_t *sel_v2v_recv_idx = NULL, *sel_v2v_recv_lst = NULL;
1306 cs_lnum_t *tmp_edges = NULL;
1307
1308 assert(s_edges != NULL);
1309
1310 /* Exchange v2v info;
1311 only info relative to edges where both vertices are available
1312 on the distant side is sent. */
1313
1314 _copy_interface_edges(ifs,
1315 sel_v2v_idx,
1316 sel_v2v_lst,
1317 &sel_v2v_recv_idx,
1318 &sel_v2v_recv_lst);
1319
1320 /* Scan adjacent faces to find new "single" edges */
1321
1322 BFT_MALLOC(tmp_edges, 2*max_size, cs_lnum_t);
1323
1324 for (i = 0; i < selection->n_b_adj_faces; i++) {
1325
1326 fid = selection->b_adj_faces[i] - 1;
1327 s = b_f2v_idx[fid];
1328 e = b_f2v_idx[fid+1];
1329
1330 for (j = s; j < e - 1; j++)
1331 _add_s_edge(vertex_tag,
1332 b_f2v_lst[j],
1333 b_f2v_lst[j+1],
1334 sel_v2v_idx,
1335 sel_v2v_lst,
1336 &tmp_size,
1337 &max_size,
1338 &tmp_edges);
1339
1340 _add_s_edge(vertex_tag,
1341 b_f2v_lst[e-1],
1342 b_f2v_lst[s],
1343 sel_v2v_idx,
1344 sel_v2v_lst,
1345 &tmp_size,
1346 &max_size,
1347 &tmp_edges);
1348
1349 }
1350
1351 for (i = 0; i < selection->n_i_adj_faces; i++) {
1352
1353 fid = selection->i_adj_faces[i] - 1;
1354
1355 if (i_face_cells[fid][0] == -1 || i_face_cells[fid][1] == -1) {
1356
1357 s = i_f2v_idx[fid];
1358 e = i_f2v_idx[fid+1];
1359
1360 for (j = s; j < e - 1; j++)
1361 _add_s_edge(vertex_tag,
1362 i_f2v_lst[j],
1363 i_f2v_lst[j+1],
1364 sel_v2v_idx,
1365 sel_v2v_lst,
1366 &tmp_size,
1367 &max_size,
1368 &tmp_edges);
1369
1370 _add_s_edge(vertex_tag,
1371 i_f2v_lst[e-1],
1372 i_f2v_lst[s],
1373 sel_v2v_idx,
1374 sel_v2v_lst,
1375 &tmp_size,
1376 &max_size,
1377 &tmp_edges);
1378
1379 } /* Face on a parallel boundary */
1380
1381 }
1382
1383 /* Find the related ranks for synchronization */
1384
1385 if (tmp_size > 0) {
1386
1387 int distant_rank;
1388
1389 cs_lnum_t if_shift;
1390 int last_found_rank = -1;
1391 cs_lnum_t n_s_edges = 0;
1392
1393 bool *edge_tag = NULL;
1394
1395 const int n_interfaces = cs_interface_set_size(ifs);
1396
1397 s_edges->n_elts = tmp_size;
1398 BFT_REALLOC(tmp_edges, 2*tmp_size, cs_lnum_t);
1399
1400 BFT_MALLOC(edge_tag, tmp_size, bool);
1401 BFT_MALLOC(s_edges->array, 2*tmp_size, cs_lnum_t);
1402 BFT_MALLOC(s_edges->index, n_interfaces + 1, cs_lnum_t);
1403 BFT_MALLOC(s_edges->ranks, n_interfaces, int);
1404
1405 for (i = 0; i < tmp_size; i++)
1406 edge_tag[i] = false; /* Not matched */
1407
1408 /* Loop on interfaces (naturally ordered by rank) */
1409
1410 if_shift = 0;
1411
1412 for (i = 0; i < n_interfaces; i++) {
1413
1414 const cs_interface_t *interface = cs_interface_set_get(ifs, i);
1415 const cs_lnum_t if_size = cs_interface_size(interface);
1416
1417 const cs_lnum_t *local_id = cs_interface_get_elt_ids(interface);
1418
1419 distant_rank = cs_interface_rank(interface);
1420
1421 /* Loop on single edges */
1422
1423 for (j = 0; j < tmp_size; j++) {
1424
1425 if (edge_tag[j] == false) { /* Not already treated */
1426
1427 cs_lnum_t i0 = cs_search_binary(if_size,
1428 tmp_edges[2*j],
1429 local_id);
1430
1431 if (i0 > -1) {
1432
1433 cs_lnum_t s_id = sel_v2v_recv_idx[if_shift + i0];
1434 cs_lnum_t e_id = sel_v2v_recv_idx[if_shift + i0 + 1];
1435
1436 /* Search on short, unsorted array */
1437
1438 for (cs_lnum_t l = s_id; l < e_id; l++) {
1439
1440 if (sel_v2v_recv_lst[l] == tmp_edges[2*j+1]) {
1441
1442 if (last_found_rank != distant_rank) {
1443 s_edges->ranks[s_edges->n_ranks] = distant_rank;
1444 s_edges->index[s_edges->n_ranks] = n_s_edges;
1445 s_edges->n_ranks++;
1446 last_found_rank = distant_rank;
1447 }
1448
1449 edge_tag[j] = true; /* Tag as done */
1450 s_edges->array[2*n_s_edges] = tmp_edges[2*j] + 1;
1451 s_edges->array[2*n_s_edges+1] = tmp_edges[2*j+1] + 1;
1452 n_s_edges++;
1453
1454 break;
1455 }
1456 }
1457
1458 }
1459
1460 } /* Not matched yet */
1461
1462 } /* End of loop on single edges */
1463
1464 if_shift += if_size;
1465
1466 } /* End of loop on interfaces */
1467
1468 s_edges->index[s_edges->n_ranks] = n_s_edges;
1469 s_edges->n_elts = n_s_edges;
1470
1471 /* Memory management */
1472
1473 if (s_edges->n_elts != tmp_size)
1474 BFT_REALLOC(s_edges->array, 2*s_edges->n_elts, cs_lnum_t);
1475
1476 BFT_REALLOC(s_edges->ranks, s_edges->n_ranks, int);
1477 BFT_REALLOC(s_edges->index, s_edges->n_ranks + 1, cs_lnum_t);
1478 BFT_FREE(edge_tag);
1479
1480 } /* End if tmp_size > 0 */
1481
1482 /* Free memory */
1483
1484 BFT_FREE(tmp_edges);
1485 BFT_FREE(sel_v2v_recv_idx);
1486 BFT_FREE(sel_v2v_recv_lst);
1487
1488 /* Logging */
1489
1490 if (cs_glob_join_log != NULL && verbosity > 3) {
1491
1492 FILE *logfile = cs_glob_join_log;
1493
1494 fprintf(logfile, "\n Single edges for the joining operation: (%p)\n",
1495 (void *)s_edges);
1496 fprintf(logfile,
1497 " Single edges: n_elts : %8ld\n"
1498 " Single edges: n_ranks: %8d\n",
1499 (long)s_edges->n_elts, s_edges->n_ranks);
1500
1501 if (s_edges->n_elts > 0) {
1502 for (i = 0; i < s_edges->n_ranks; i++)
1503 for (j = s_edges->index[i]; j < s_edges->index[i+1]; j++)
1504 fprintf(logfile, " %9ld | %6d | (%9ld, %9ld)\n",
1505 (long)j, s_edges->ranks[i], (long)s_edges->array[2*j],
1506 (long)s_edges->array[2*j+1]);
1507 fprintf(logfile, "\n");
1508 }
1509 fflush(logfile);
1510 }
1511 }
1512
1513 /*----------------------------------------------------------------------------
1514 * Define a structure for the coupled edges. Only done if single edges have
1515 * been detected and only in parallel.
1516 *
1517 * parameters:
1518 * ifs <-- pointer to the interface set on vertices
1519 * s_edges <-- single edges structure used to build coupled_edges
1520 * c_edges <-> pointer to the coupled edges structure to define
1521 * verbosity <-- verbosity level
1522 *---------------------------------------------------------------------------*/
1523
1524 static void
_add_coupled_edges(cs_interface_set_t * ifs,cs_join_sync_t * s_edges,cs_join_sync_t * c_edges,int verbosity)1525 _add_coupled_edges(cs_interface_set_t *ifs,
1526 cs_join_sync_t *s_edges,
1527 cs_join_sync_t *c_edges,
1528 int verbosity)
1529 {
1530 cs_lnum_t i, j, id, n_entities;
1531 int request_count, rank_shift, distant_rank;
1532
1533 cs_lnum_t *buf = NULL, *recv_buf = NULL, *send_buf = NULL;
1534 MPI_Request *request = NULL;
1535 MPI_Status *status = NULL;
1536 MPI_Comm mpi_comm = cs_glob_mpi_comm;
1537
1538 const int local_rank = CS_MAX(cs_glob_rank_id, 0);
1539 const cs_lnum_t *local_id = NULL, *distant_id = NULL;
1540 const int n_interfaces = cs_interface_set_size(ifs);
1541 const cs_interface_t *interface = NULL;
1542
1543 assert(s_edges != NULL);
1544 assert(c_edges != NULL);
1545
1546 /* Exchange number of single edges */
1547
1548 BFT_MALLOC(request, n_interfaces * 2, MPI_Request);
1549 BFT_MALLOC(status, n_interfaces * 2, MPI_Status);
1550 BFT_MALLOC(buf, 2*n_interfaces, cs_lnum_t);
1551
1552 for (i = 0; i < 2*n_interfaces; i++)
1553 buf[i] = 0;
1554
1555 request_count = 0;
1556
1557 for (id = 0; id < n_interfaces; id++) {
1558
1559 interface = cs_interface_set_get(ifs, id);
1560 distant_rank = cs_interface_rank(interface);
1561
1562 MPI_Irecv(&(buf[id]),
1563 1,
1564 CS_MPI_LNUM,
1565 distant_rank,
1566 distant_rank,
1567 mpi_comm,
1568 &(request[request_count++]));
1569
1570 }
1571
1572 /* Send */
1573
1574 rank_shift = 0;
1575
1576 for (id = 0; id < n_interfaces; id++) {
1577
1578 /* Preparation of data to send */
1579
1580 interface = cs_interface_set_get(ifs, id);
1581 distant_rank = cs_interface_rank(interface);
1582
1583 if (rank_shift < s_edges->n_ranks) {
1584 if (s_edges->ranks[rank_shift] == distant_rank) {
1585 buf[n_interfaces + id] =
1586 s_edges->index[rank_shift+1] - s_edges->index[rank_shift];
1587 rank_shift++;
1588 }
1589 }
1590
1591 MPI_Isend(&(buf[n_interfaces + id]),
1592 1,
1593 CS_MPI_LNUM,
1594 distant_rank,
1595 local_rank,
1596 mpi_comm,
1597 &(request[request_count++]));
1598
1599 } /* End of loop on interfaces */
1600
1601 /* Sync after each rank had received all the messages */
1602
1603 MPI_Waitall(request_count, request, status);
1604
1605 /* Define c_edges */
1606
1607 for (i = 0; i < n_interfaces; i++) {
1608 if (buf[i] > 0) {
1609 c_edges->n_elts += buf[i];
1610 c_edges->n_ranks += 1;
1611 }
1612 }
1613
1614 BFT_MALLOC(c_edges->ranks, c_edges->n_ranks, int);
1615 BFT_MALLOC(c_edges->index, c_edges->n_ranks + 1, cs_lnum_t);
1616 BFT_MALLOC(c_edges->array, 2*c_edges->n_elts, cs_lnum_t);
1617
1618 for (i = 0; i < c_edges->n_ranks + 1; i++)
1619 c_edges->index[i] = 0;
1620
1621 rank_shift = 0;
1622
1623 for (i = 0; i < n_interfaces; i++) {
1624 if (buf[i] > 0) {
1625
1626 interface = cs_interface_set_get(ifs, i);
1627 distant_rank = cs_interface_rank(interface);
1628 c_edges->ranks[rank_shift++] = distant_rank;
1629 c_edges->index[rank_shift] = buf[i];
1630
1631 }
1632 }
1633
1634 assert(rank_shift == c_edges->n_ranks);
1635
1636 for (i = 0; i < c_edges->n_ranks; i++)
1637 c_edges->index[i+1] += c_edges->index[i];
1638
1639 assert(c_edges->index[c_edges->n_ranks] == c_edges->n_elts);
1640
1641 BFT_FREE(buf);
1642
1643 /* Exchange couple of vertices defining coupled edges */
1644
1645 request_count = 0;
1646 rank_shift = 0;
1647
1648 for (id = 0; id < n_interfaces; id++) {
1649
1650 interface = cs_interface_set_get(ifs, id);
1651 distant_rank = cs_interface_rank(interface);
1652
1653 if (rank_shift < c_edges->n_ranks) {
1654 if (c_edges->ranks[rank_shift] == distant_rank) {
1655
1656 n_entities = c_edges->index[rank_shift+1] - c_edges->index[rank_shift];
1657 n_entities *= 2; /* couple of vertices */
1658 recv_buf = c_edges->array + 2*c_edges->index[rank_shift];
1659 rank_shift++;
1660
1661 MPI_Irecv(recv_buf, /* receive distant num */
1662 n_entities,
1663 CS_MPI_LNUM,
1664 distant_rank,
1665 distant_rank,
1666 mpi_comm,
1667 &(request[request_count++]));
1668
1669 }
1670 }
1671
1672 } /* End of loop on interfaces */
1673
1674 /* Send */
1675
1676 rank_shift = 0;
1677
1678 for (id = 0; id < n_interfaces; id++) {
1679
1680 interface = cs_interface_set_get(ifs, id);
1681 distant_rank = cs_interface_rank(interface);
1682
1683 if (rank_shift < s_edges->n_ranks) {
1684 if (s_edges->ranks[rank_shift] == distant_rank) {
1685
1686 n_entities = s_edges->index[rank_shift+1] - s_edges->index[rank_shift];
1687 n_entities *= 2; /* couple of vertices */
1688 send_buf = s_edges->array + 2*s_edges->index[rank_shift];
1689 rank_shift++;
1690
1691 MPI_Isend(send_buf,
1692 n_entities,
1693 CS_MPI_LNUM,
1694 distant_rank,
1695 local_rank,
1696 mpi_comm,
1697 &(request[request_count++]));
1698
1699 }
1700 }
1701
1702 } /* End of loop on interfaces */
1703
1704 /* Sync after each rank had received all the messages */
1705
1706 MPI_Waitall(request_count, request, status);
1707
1708 BFT_FREE(request);
1709 BFT_FREE(status);
1710
1711 rank_shift = 0;
1712
1713 /* Switch received couple vertices from distant id. to local id. */
1714
1715 cs_interface_set_add_match_ids(ifs);
1716
1717 for (i = 0; i < n_interfaces; i++) {
1718
1719 interface = cs_interface_set_get(ifs, i);
1720 distant_rank = cs_interface_rank(interface);
1721 distant_id = cs_interface_get_match_ids(interface);
1722
1723 n_entities = cs_interface_size(interface);
1724
1725 if (rank_shift < c_edges->n_ranks) {
1726 if (c_edges->ranks[rank_shift] == distant_rank) {
1727
1728 local_id = cs_interface_get_elt_ids(interface);
1729
1730 for (j = c_edges->index[rank_shift];
1731 j < c_edges->index[rank_shift+1];
1732 j++) {
1733
1734 id = cs_search_binary(n_entities, c_edges->array[2*j] - 1, distant_id);
1735 assert(id != -1);
1736 c_edges->array[2*j] = local_id[id] + 1;
1737
1738 id = cs_search_binary(n_entities,
1739 c_edges->array[2*j+1] - 1,
1740 distant_id);
1741 assert(id != -1);
1742 c_edges->array[2*j+1] = local_id[id] + 1;
1743
1744 } /* Loop on couple edges */
1745
1746 rank_shift++;
1747
1748 } /* This rank is in the list of related ranks */
1749 }
1750
1751 } /* Loop on interfaces */
1752
1753 cs_interface_set_free_match_ids(ifs);
1754
1755 if (cs_glob_join_log != NULL && verbosity > 3) {
1756
1757 FILE *logfile = cs_glob_join_log;
1758
1759 fprintf(logfile, "\n Coupled edges for the joining operation: (%p)\n",
1760 (void *)c_edges);
1761 fprintf(logfile,
1762 " Coupled edges: n_elts : %8ld\n"
1763 " Coupled edges: n_ranks: %8d\n",
1764 (long)c_edges->n_elts, c_edges->n_ranks);
1765
1766 if (c_edges->n_elts > 0) {
1767 int shift;
1768 for (i = 0, shift = 0; i < c_edges->n_ranks; i++) {
1769 for (j = c_edges->index[i]; j < c_edges->index[i+1]; j++) {
1770 fprintf(logfile, " %9ld | %6d | (%9ld, %9ld)\n",
1771 (long)shift, c_edges->ranks[i],
1772 (long)c_edges->array[2*j], (long)c_edges->array[2*j+1]);
1773 shift++;
1774 }
1775 }
1776 fprintf(logfile, "\n");
1777 fflush(logfile);
1778 }
1779 }
1780
1781 }
1782
1783 /*----------------------------------------------------------------------------
1784 * Get the full selection of single edges. Done only if the run is parallel.
1785 *
1786 * parameters:
1787 * selection <-> pointer to a fvm_join_selection_t structure
1788 * sel_v2v_idx <-- vertex -> vertex connect. index
1789 * sel_v2v_lst <-- vertex -> vertex connect. list
1790 *---------------------------------------------------------------------------*/
1791
1792 static void
_filter_edge_element(cs_join_select_t * selection,const cs_lnum_t sel_v2v_idx[],const cs_lnum_t sel_v2v_lst[])1793 _filter_edge_element(cs_join_select_t *selection,
1794 const cs_lnum_t sel_v2v_idx[],
1795 const cs_lnum_t sel_v2v_lst[])
1796 {
1797 cs_lnum_t i, j, vid1, vid2, edge_id, request_count, shift, save;
1798
1799 int *c_edge_tag = NULL, *s_edge_tag = NULL;
1800 cs_join_sync_t *s_edges = selection->s_edges;
1801 cs_join_sync_t *c_edges = selection->c_edges;
1802
1803 MPI_Request *request = NULL;
1804 MPI_Status *status = NULL;
1805 MPI_Comm mpi_comm = cs_glob_mpi_comm;
1806
1807 const int loc_rank = CS_MAX(cs_glob_rank_id, 0);
1808
1809 assert(cs_glob_n_ranks > 1);
1810 assert(c_edges != NULL);
1811 assert(s_edges != NULL);
1812
1813 /* Allocate MPI buffers used for exchanging data */
1814
1815 BFT_MALLOC(request, c_edges->n_ranks + s_edges->n_ranks, MPI_Request);
1816 BFT_MALLOC(status, c_edges->n_ranks + s_edges->n_ranks, MPI_Status);
1817
1818 BFT_MALLOC(c_edge_tag, c_edges->n_elts, int);
1819 BFT_MALLOC(s_edge_tag, s_edges->n_elts, int);
1820
1821 for (i = 0; i < c_edges->n_elts; i++)
1822 c_edge_tag[i] = 1; /* define as selected */
1823
1824 for (i = 0; i < s_edges->n_elts; i++)
1825 s_edge_tag[i] = 1; /* define as selected */
1826
1827 for (i = 0; i < c_edges->n_elts; i++) {
1828
1829 vid1 = c_edges->array[2*i] - 1;
1830 vid2 = c_edges->array[2*i+1] - 1;
1831 edge_id = _get_edge_id(vid1, vid2, sel_v2v_idx, sel_v2v_lst);
1832
1833 if (edge_id == -1)
1834 c_edge_tag[i] = 0; /* unselect this coupled edge */
1835
1836 } /* End of loop on c_edges */
1837
1838 /* Exchange between ranks the status of the coupled edges.
1839 If one receive a tag equal to 0 => delete the related single edge */
1840
1841 request_count = 0;
1842
1843 for (i = 0; i < s_edges->n_ranks; i++) {
1844
1845 int distant_rank = s_edges->ranks[i];
1846 int length = s_edges->index[i+1] - s_edges->index[i];
1847 int *recv_buf = s_edge_tag + s_edges->index[i];
1848
1849 MPI_Irecv(recv_buf,
1850 length,
1851 MPI_INT,
1852 distant_rank,
1853 distant_rank,
1854 mpi_comm,
1855 &(request[request_count++]));
1856
1857 }
1858
1859 /* Send data to distant ranks */
1860
1861 for (i = 0; i < c_edges->n_ranks; i++) {
1862
1863 int distant_rank = c_edges->ranks[i];
1864 int length = c_edges->index[i+1] - c_edges->index[i];
1865 int *send_buf = c_edge_tag + c_edges->index[i];
1866
1867 MPI_Isend(send_buf,
1868 length,
1869 MPI_INT,
1870 distant_rank,
1871 loc_rank,
1872 mpi_comm,
1873 &(request[request_count++]));
1874
1875 }
1876
1877 /* Wait for all exchanges */
1878
1879 MPI_Waitall(request_count, request, status);
1880
1881 /* Delete unselected elements */
1882
1883 shift = 0;
1884 if (c_edges-> n_elts > 0) {
1885
1886 save = c_edges->index[0];
1887
1888 for (i = 0; i < c_edges->n_ranks; i++) {
1889
1890 for (j = save; j < c_edges->index[i+1]; j++) {
1891 if (c_edge_tag[j] == 1) {
1892 c_edges->array[2*shift] = c_edges->array[2*j];
1893 c_edges->array[2*shift+1] = c_edges->array[2*j+1];
1894 shift++;
1895 }
1896 }
1897 save = c_edges->index[i+1];
1898 c_edges->index[i+1] = shift;
1899
1900 }
1901 c_edges->n_elts = shift;
1902
1903 } /* c_edges->n_elts > 0 */
1904
1905 shift = 0;
1906 if (s_edges->n_elts > 0) {
1907
1908 save = s_edges->index[0];
1909
1910 for (i = 0; i < s_edges->n_ranks; i++) {
1911
1912 for (j = save; j < s_edges->index[i+1]; j++) {
1913 if (s_edge_tag[j] == 1) {
1914 s_edges->array[2*shift] = s_edges->array[2*j];
1915 s_edges->array[2*shift+1] = s_edges->array[2*j+1];
1916 shift++;
1917 }
1918 }
1919 save = s_edges->index[i+1];
1920 s_edges->index[i+1] = shift;
1921
1922 }
1923 s_edges->n_elts = shift;
1924
1925 } /* s_edges->n_elts > 0 */
1926
1927 /* Free memory */
1928
1929 BFT_FREE(c_edge_tag);
1930 BFT_FREE(s_edge_tag);
1931 BFT_FREE(request);
1932 BFT_FREE(status);
1933 }
1934
1935 /*----------------------------------------------------------------------------
1936 * Get the full selection of vertices to extract. Only in parallel. Somme
1937 * vertices may have been selected on another rank and not on the local rank
1938 * but you have to take them into account to have a good update of the mesh
1939 *
1940 * parameters:
1941 * b_f2v_idx <-- boundary "face -> vertex" connect. index
1942 * b_f2v_lst <-- boundary "face -> vertex" connect. list
1943 * i_f2v_idx <-- interior "face -> vertex" connect. index
1944 * i_f2v_lst <-- interior "face -> vertex" connect. list
1945 * n_vertices <-- number of vertices in the parent mesh
1946 * vtx_tag <-- tag on vertices. 1 if selected, 0 otherwise
1947 * ifs <-- pointer to a cs_interface_set_t struct.
1948 * i_face_cells <-- interior face -> cells connect.
1949 * join_select <-> pointer to a fvm_join_selection_t structure
1950 * verbosity <-- verbosity level
1951 *---------------------------------------------------------------------------*/
1952
1953 static void
_get_missing_edges(cs_lnum_t b_f2v_idx[],cs_lnum_t b_f2v_lst[],cs_lnum_t i_f2v_idx[],cs_lnum_t i_f2v_lst[],cs_lnum_t n_vertices,cs_lnum_t vtx_tag[],cs_interface_set_t * ifs,cs_lnum_2_t i_face_cells[],cs_join_select_t * selection,int verbosity)1954 _get_missing_edges(cs_lnum_t b_f2v_idx[],
1955 cs_lnum_t b_f2v_lst[],
1956 cs_lnum_t i_f2v_idx[],
1957 cs_lnum_t i_f2v_lst[],
1958 cs_lnum_t n_vertices,
1959 cs_lnum_t vtx_tag[],
1960 cs_interface_set_t *ifs,
1961 cs_lnum_2_t i_face_cells[],
1962 cs_join_select_t *selection,
1963 int verbosity)
1964 {
1965 cs_gnum_t n_l_elts, n_g_elts;
1966
1967 cs_lnum_t *sel_v2v_idx = NULL, *sel_v2v_lst = NULL;
1968
1969 assert(ifs != NULL);
1970 assert(selection != NULL);
1971
1972 /* Define single edge element */
1973
1974 _get_select_v2v_connect(n_vertices,
1975 selection,
1976 b_f2v_idx,
1977 b_f2v_lst,
1978 &sel_v2v_idx,
1979 &sel_v2v_lst);
1980
1981 _add_single_edges(ifs,
1982 vtx_tag,
1983 selection,
1984 sel_v2v_idx,
1985 sel_v2v_lst,
1986 b_f2v_idx,
1987 b_f2v_lst,
1988 i_f2v_idx,
1989 i_f2v_lst,
1990 i_face_cells,
1991 selection->s_edges,
1992 verbosity);
1993
1994 n_l_elts = selection->s_edges->n_elts;
1995
1996 MPI_Allreduce(&n_l_elts, &n_g_elts, 1, CS_MPI_GNUM,
1997 MPI_SUM, cs_glob_mpi_comm);
1998
1999 if (n_g_elts > 0) {
2000
2001 bft_printf(_(" Global number of single edges found: %6llu\n"),
2002 (unsigned long long)n_g_elts);
2003 bft_printf_flush();
2004 selection->do_single_sync = true;
2005
2006 _add_coupled_edges(ifs,
2007 selection->s_edges,
2008 selection->c_edges,
2009 verbosity);
2010
2011 _filter_edge_element(selection,
2012 sel_v2v_idx,
2013 sel_v2v_lst);
2014
2015 } /* End if n_g_s_elts > 0 */
2016
2017 if (cs_glob_join_log != NULL && verbosity > 3) {
2018 int i, j;
2019 cs_join_sync_t *s_edges = selection->s_edges;
2020 cs_join_sync_t *c_edges = selection->c_edges;
2021 FILE *logfile = cs_glob_join_log;
2022
2023 fprintf(logfile, "\n After possible filtering\n");
2024
2025 fprintf(logfile, "\n Coupled edges for the joining operation: (%p)\n",
2026 (void *)c_edges);
2027 fprintf(logfile,
2028 " Coupled edges: n_elts : %8ld\n"
2029 " Coupled edges: n_ranks: %8d\n",
2030 (long)c_edges->n_elts, c_edges->n_ranks);
2031
2032 if (c_edges->n_elts > 0) {
2033 int shift;
2034 for (i = 0, shift = 0; i < c_edges->n_ranks; i++) {
2035 for (j = c_edges->index[i]; j < c_edges->index[i+1]; j++) {
2036 fprintf(logfile, " %9ld | %6d | (%9ld, %9ld)\n",
2037 (long)shift, c_edges->ranks[i],
2038 (long)c_edges->array[2*j], (long)c_edges->array[2*j+1]);
2039 shift++;
2040 }
2041 }
2042 fprintf(logfile, "\n");
2043 }
2044
2045 fprintf(logfile, "\n Single edges for the joining operation: (%p)\n",
2046 (void *)s_edges);
2047 fprintf(logfile,
2048 " Single edges: n_elts : %8ld\n"
2049 " Single edges: n_ranks: %8d\n",
2050 (long)s_edges->n_elts, s_edges->n_ranks);
2051
2052 if (s_edges->n_elts > 0) {
2053 for (i = 0; i < s_edges->n_ranks; i++)
2054 for (j = s_edges->index[i]; j < s_edges->index[i+1]; j++)
2055 fprintf(logfile, " %9ld | %6d | (%9ld, %9ld)\n",
2056 (long)j, s_edges->ranks[i], (long)s_edges->array[2*j],
2057 (long)s_edges->array[2*j+1]);
2058 fprintf(logfile, "\n");
2059 }
2060
2061 fflush(logfile);
2062 bft_printf_flush();
2063 }
2064
2065 /* Free memory */
2066
2067 BFT_FREE(sel_v2v_idx);
2068 BFT_FREE(sel_v2v_lst);
2069
2070 }
2071 #endif /* defined(HAVE_MPI) */
2072
2073 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
2074
2075 /*============================================================================
2076 * Public function definitions
2077 *===========================================================================*/
2078
2079 /*----------------------------------------------------------------------------
2080 * Create and initialize a cs_join_t structure.
2081 *
2082 * parameters:
2083 * join_number <-- number related to the joining operation
2084 * sel_criteria <-- boundary face selection criteria
2085 * fraction <-- value of the fraction parameter
2086 * plane <-- value of the plane parameter
2087 * perio_type <-- periodicity type (FVM_PERIODICITY_NULL if none)
2088 * perio_matrix <-- periodicity transformation matrix
2089 * verbosity <-- level of verbosity required
2090 * visualization <-- level of visualization required
2091 * preprocessing <-- is joining part of the preprocessing stage ?
2092 *
2093 * returns:
2094 * a pointer to a new allocated cs_join_t structure
2095 *---------------------------------------------------------------------------*/
2096
2097 cs_join_t *
cs_join_create(int join_number,const char * sel_criteria,float fraction,float plane,fvm_periodicity_type_t perio_type,double perio_matrix[3][4],int verbosity,int visualization,bool preprocessing)2098 cs_join_create(int join_number,
2099 const char *sel_criteria,
2100 float fraction,
2101 float plane,
2102 fvm_periodicity_type_t perio_type,
2103 double perio_matrix[3][4],
2104 int verbosity,
2105 int visualization,
2106 bool preprocessing)
2107 {
2108 size_t l;
2109
2110 cs_join_t *join = NULL;
2111
2112 /* Check main parameter values */
2113
2114 if (fraction < 0.0 || fraction >= 1.0)
2115 bft_error(__FILE__, __LINE__, 0,
2116 _("Mesh joining:"
2117 " Forbidden value for the fraction parameter.\n"
2118 " It must be between [0.0, 1.0[ and is here: %f\n"),
2119 fraction);
2120
2121 if (plane < 0.0 || plane >= 90.0)
2122 bft_error(__FILE__, __LINE__, 0,
2123 _("Mesh joining:"
2124 " Forbidden value for the plane parameter.\n"
2125 " It must be between [0, 90] and is here: %f\n"),
2126 plane);
2127
2128 /* Initialize structure */
2129
2130 BFT_MALLOC(join, 1, cs_join_t);
2131
2132 join->selection = NULL;
2133
2134 join->param = _join_param_define(join_number,
2135 fraction,
2136 plane,
2137 perio_type,
2138 perio_matrix,
2139 verbosity,
2140 visualization,
2141 preprocessing);
2142
2143 join->stats = _join_stats_init();
2144
2145 join->log_name = NULL;
2146
2147 /* Copy the selection criteria for future use */
2148
2149 l = strlen(sel_criteria);
2150 BFT_MALLOC(join->criteria, l + 1, char);
2151 strcpy(join->criteria, sel_criteria);
2152
2153 /* Initialize log file if necessary */
2154
2155 if (verbosity > 2) {
2156 char logname[80];
2157 char dir[] = "log";
2158 char rank_add[16] = "";
2159 char perio_add[16] = "";
2160 if (cs_file_isdir(dir) == 0) {
2161 if (cs_glob_rank_id < 1)
2162 if (cs_file_mkdir_default(dir) != 0)
2163 bft_error(__FILE__, __LINE__, 0,
2164 _("The log directory cannot be created"));
2165 #if defined(HAVE_MPI)
2166 if (cs_glob_n_ranks > 1)
2167 MPI_Barrier(cs_glob_mpi_comm); /* to avoid race conditions */
2168 #endif
2169 }
2170 if (perio_type != FVM_PERIODICITY_NULL)
2171 strcpy(perio_add, "_perio");
2172 if (cs_glob_n_ranks > 1)
2173 sprintf(rank_add, "_r%04d", cs_glob_rank_id);
2174 sprintf(logname, "log%cjoin_%02d%s%s.log", DIR_SEPARATOR,
2175 join_number, perio_add, rank_add);
2176 BFT_MALLOC(join->log_name, strlen(logname) + 1, char);
2177 strcpy(join->log_name, logname);
2178 }
2179
2180 return join;
2181 }
2182
2183 /*----------------------------------------------------------------------------
2184 * Destroy a cs_join_t structure.
2185 *
2186 * parameters:
2187 * join <-> pointer to the cs_join_t structure to destroy
2188 *---------------------------------------------------------------------------*/
2189
2190 void
cs_join_destroy(cs_join_t ** join)2191 cs_join_destroy(cs_join_t **join)
2192 {
2193 if (*join != NULL) {
2194
2195 cs_join_t *_join = *join;
2196
2197 BFT_FREE(_join->log_name);
2198 BFT_FREE(_join->criteria);
2199
2200 BFT_FREE(_join);
2201 *join = NULL;
2202
2203 }
2204 }
2205
2206 /*----------------------------------------------------------------------------
2207 * Create and initialize a cs_join_select_t structure.
2208 *
2209 * parameters:
2210 * selection_criteria <-- pointer to a cs_mesh_select_t structure
2211 * perio_type <-- periodicity type (FVM_PERIODICITY_NULL if none)
2212 * verbosity <-- level of verbosity required
2213 *
2214 * returns:
2215 * pointer to a newly created cs_join_select_t structure
2216 *---------------------------------------------------------------------------*/
2217
2218 cs_join_select_t *
cs_join_select_create(const char * selection_criteria,fvm_periodicity_type_t perio_type,int verbosity)2219 cs_join_select_create(const char *selection_criteria,
2220 fvm_periodicity_type_t perio_type,
2221 int verbosity)
2222 {
2223 cs_lnum_t *vtx_tag = NULL;
2224 cs_join_select_t *selection = NULL;
2225 cs_lnum_t *order = NULL, *ordered_faces = NULL;
2226 cs_interface_set_t *ifs = NULL;
2227 cs_mesh_t *mesh = cs_glob_mesh;
2228 FILE *logfile = cs_glob_join_log;
2229
2230 const int n_ranks = cs_glob_n_ranks;
2231
2232 assert(mesh != NULL);
2233
2234 /* Initialize cs_join_select_t struct. */
2235
2236 BFT_MALLOC(selection, 1, cs_join_select_t);
2237
2238 selection->n_init_b_faces = mesh->n_b_faces;
2239 selection->n_init_i_faces = mesh->n_i_faces;
2240 selection->n_init_vertices = mesh->n_vertices;
2241
2242 selection->n_faces = 0;
2243 selection->n_g_faces = 0;
2244
2245 selection->faces = NULL;
2246 selection->compact_face_gnum = NULL;
2247 selection->compact_rank_index = NULL;
2248
2249 selection->n_vertices = 0;
2250 selection->n_g_vertices = 0;
2251 selection->vertices = NULL;
2252
2253 selection->n_b_adj_faces = 0;
2254 selection->b_adj_faces = NULL;
2255 selection->n_i_adj_faces = 0;
2256 selection->i_adj_faces = NULL;
2257
2258 selection->b_face_state = NULL;
2259 selection->i_face_state = NULL;
2260
2261 selection->n_couples = 0;
2262 selection->per_v_couples = NULL;
2263
2264 /*
2265 Single elements (Only possible in parallel. It appears
2266 when the domain splitting has a poor quality and elements
2267 on the joining interface are prisms or tetraedrals)
2268 s = single / c = coupled
2269 */
2270
2271 selection->do_single_sync = false;
2272 selection->s_vertices = _create_join_sync();
2273 selection->c_vertices = _create_join_sync();
2274 selection->s_edges = _create_join_sync();
2275 selection->c_edges = _create_join_sync();
2276
2277 /* Extract selected boundary faces */
2278
2279 BFT_MALLOC(selection->faces, mesh->n_b_faces, cs_lnum_t);
2280
2281 cs_selector_get_b_face_num_list(selection_criteria,
2282 &(selection->n_faces),
2283 selection->faces);
2284
2285 /* In case of periodicity, ensure no isolated faces are
2286 selected */
2287
2288 if (perio_type != FVM_PERIODICITY_NULL) {
2289 cs_lnum_t j = 0;
2290 for (cs_lnum_t i = 0; i < selection->n_faces; i++) {
2291 cs_lnum_t f_id = selection->faces[i]-1;
2292 if (mesh->b_face_cells[f_id] > -1)
2293 selection->faces[j++] = f_id+1;
2294 }
2295 selection->n_faces = j;
2296 }
2297
2298 BFT_MALLOC(order, selection->n_faces, cs_lnum_t);
2299 BFT_MALLOC(ordered_faces, selection->n_faces, cs_lnum_t);
2300
2301 cs_order_gnum_allocated(selection->faces, NULL, order, selection->n_faces);
2302
2303 for (cs_lnum_t i = 0; i < selection->n_faces; i++)
2304 ordered_faces[i] = selection->faces[order[i]];
2305
2306 BFT_FREE(order);
2307 BFT_FREE(selection->faces);
2308 selection->faces = ordered_faces;
2309
2310 if (n_ranks == 1)
2311 selection->n_g_faces = selection->n_faces;
2312
2313 #if defined(HAVE_MPI)
2314 if (n_ranks > 1) { /* Parallel treatment */
2315
2316 cs_gnum_t n_l_faces = selection->n_faces;
2317
2318 MPI_Allreduce(&n_l_faces, &(selection->n_g_faces),
2319 1, CS_MPI_GNUM, MPI_SUM, cs_glob_mpi_comm);
2320
2321 }
2322 #endif
2323
2324 if (verbosity > 0)
2325 bft_printf
2326 (_(" Global number of boundary faces selected for joining: %10llu\n"),
2327 (unsigned long long)selection->n_g_faces);
2328
2329 /* Define a compact global numbering on selected boundary faces and
2330 build an index on ranks on this compact numbering */
2331
2332 _compact_face_gnum_selection(selection->n_faces,
2333 &(selection->compact_face_gnum),
2334 &(selection->compact_rank_index));
2335
2336 assert(selection->n_g_faces == selection->compact_rank_index[n_ranks]);
2337
2338 /* Extract selected vertices from the selected boundary faces */
2339
2340 cs_join_extract_vertices(selection->n_faces,
2341 selection->faces,
2342 mesh->b_face_vtx_idx,
2343 mesh->b_face_vtx_lst,
2344 mesh->n_vertices,
2345 &(selection->n_vertices),
2346 &(selection->vertices));
2347
2348 #if defined(HAVE_MPI)
2349 if (n_ranks > 1) { /* Search for missing vertices */
2350
2351 assert(mesh->global_vtx_num != NULL);
2352
2353 ifs = cs_interface_set_create(mesh->n_vertices,
2354 NULL,
2355 mesh->global_vtx_num,
2356 NULL,
2357 0,
2358 NULL,
2359 NULL,
2360 NULL);
2361
2362 assert(ifs != NULL);
2363
2364 _get_missing_vertices(mesh->n_vertices,
2365 ifs,
2366 &vtx_tag,
2367 selection,
2368 verbosity);
2369
2370 }
2371 #endif
2372
2373 /* Extract list of boundary faces contiguous to the selected vertices */
2374
2375 _extract_contig_faces(mesh->n_vertices,
2376 selection,
2377 mesh->n_b_faces,
2378 mesh->b_face_vtx_idx,
2379 mesh->b_face_vtx_lst,
2380 &(selection->n_b_adj_faces),
2381 &(selection->b_adj_faces));
2382
2383 /* Remove boundary faces already defined in selection->faces */
2384
2385 cs_join_clean_selection(&(selection->n_b_adj_faces),
2386 &(selection->b_adj_faces),
2387 selection->n_faces,
2388 selection->faces);
2389
2390 /* Extract list of interior faces contiguous to the selected vertices */
2391
2392 _extract_contig_faces(mesh->n_vertices,
2393 selection,
2394 mesh->n_i_faces,
2395 mesh->i_face_vtx_idx,
2396 mesh->i_face_vtx_lst,
2397 &(selection->n_i_adj_faces),
2398 &(selection->i_adj_faces));
2399
2400 /* Check if there is no forgotten vertex in the selection.
2401 Otherwise define structures to enable future synchronization.
2402 Only possible in parallel. */
2403
2404 #if defined(HAVE_MPI)
2405 if (n_ranks > 1) {
2406
2407 _get_missing_edges(mesh->b_face_vtx_idx,
2408 mesh->b_face_vtx_lst,
2409 mesh->i_face_vtx_idx,
2410 mesh->i_face_vtx_lst,
2411 mesh->n_vertices,
2412 vtx_tag,
2413 ifs,
2414 mesh->i_face_cells,
2415 selection,
2416 verbosity);
2417
2418 BFT_FREE(vtx_tag);
2419 cs_interface_set_destroy(&ifs);
2420
2421 }
2422 #endif
2423
2424 /* Face state setting */
2425
2426 BFT_MALLOC(selection->b_face_state, mesh->n_b_faces, cs_join_state_t);
2427 BFT_MALLOC(selection->i_face_state, mesh->n_i_faces, cs_join_state_t);
2428
2429 for (cs_lnum_t i = 0; i < mesh->n_b_faces; i++)
2430 selection->b_face_state[i] = CS_JOIN_STATE_UNDEF;
2431
2432 for (cs_lnum_t i = 0; i < mesh->n_i_faces; i++)
2433 selection->i_face_state[i] = CS_JOIN_STATE_UNDEF;
2434
2435 for (cs_lnum_t i = 0; i < selection->n_faces; i++)
2436 selection->b_face_state[selection->faces[i]-1] = CS_JOIN_STATE_ORIGIN;
2437
2438 for (cs_lnum_t i = 0; i < selection->n_b_adj_faces; i++)
2439 selection->b_face_state[selection->b_adj_faces[i]-1] = CS_JOIN_STATE_ORIGIN;
2440
2441 for (cs_lnum_t i = 0; i < selection->n_i_adj_faces; i++)
2442 selection->i_face_state[selection->i_adj_faces[i]-1] = CS_JOIN_STATE_ORIGIN;
2443
2444 /* Display information according to the level of verbosity */
2445
2446 if (verbosity > 2) {
2447
2448 assert(logfile != NULL);
2449
2450 fprintf(logfile,
2451 "\n Local information about selection structure:\n");
2452 fprintf(logfile,
2453 " number of faces: %8ld\n",
2454 (long)selection->n_faces);
2455 fprintf(logfile,
2456 " number of vertices: %8ld\n",
2457 (long)selection->n_vertices);
2458 fprintf(logfile,
2459 " number of adj. boundary faces: %8ld\n",
2460 (long)selection->n_b_adj_faces);
2461 fprintf(logfile,
2462 " number of adj. interior faces: %8ld\n",
2463 (long)selection->n_i_adj_faces);
2464
2465 if (selection->do_single_sync == true) {
2466 fprintf(logfile,
2467 "\n Information on single/coupled elements:\n");
2468 fprintf(logfile,
2469 " Number of single vertices : %6ld with %3d related ranks\n",
2470 (long)selection->s_vertices->n_elts,
2471 selection->s_vertices->n_ranks);
2472 fprintf(logfile,
2473 " Number of coupled vertices: %6ld with %3d related ranks\n",
2474 (long)selection->c_vertices->n_elts,
2475 selection->c_vertices->n_ranks);
2476 fprintf(logfile,
2477 " Number of single edges : %6ld with %3d related ranks\n",
2478 (long)selection->s_edges->n_elts,
2479 selection->s_edges->n_ranks);
2480 fprintf(logfile,
2481 " Number of coupled edges : %6ld with %3d related ranks\n",
2482 (long)selection->c_edges->n_elts,
2483 selection->c_edges->n_ranks);
2484 }
2485
2486 if (verbosity > 3) {
2487 fprintf(logfile,
2488 "\n Compact index on ranks for the selected faces:\n");
2489 for (int i = 0; i < n_ranks + 1; i++)
2490 fprintf(logfile,
2491 " %5d | %11llu\n", i,
2492 (unsigned long long)selection->compact_rank_index[i]);
2493 fprintf(logfile, "\n");
2494 }
2495
2496 /* Debug level logging */
2497
2498 if (verbosity > 4) {
2499
2500 fprintf(logfile,
2501 "\n Selected faces for the joining operation:\n");
2502 for (cs_lnum_t i = 0; i < selection->n_faces; i++)
2503 fprintf(logfile,
2504 " %9ld | %9ld | %10llu\n",
2505 (long)i, (long)selection->faces[i],
2506 (unsigned long long)selection->compact_face_gnum[i]);
2507 fprintf(logfile, "\n");
2508 }
2509
2510 if (verbosity > 4) {
2511 fprintf(logfile,
2512 "\n Selected vertices for the joining operation:\n");
2513 for (cs_lnum_t i = 0; i < selection->n_vertices; i++)
2514 fprintf(logfile,
2515 " %9ld | %9ld\n", (long)i, (long)selection->vertices[i]);
2516 fprintf(logfile, "\n");
2517 }
2518
2519 if (verbosity > 4) {
2520 fprintf(logfile,
2521 "\n Contiguous boundary faces for the joining operation:\n");
2522 for (cs_lnum_t i = 0; i < selection->n_b_adj_faces; i++)
2523 fprintf(logfile,
2524 " %9ld | %9ld\n", (long)i, (long)selection->b_adj_faces[i]);
2525 fprintf(logfile, "\n");
2526 }
2527
2528 if (verbosity > 4) {
2529 fprintf(logfile,
2530 "\n Contiguous interior faces for the joining operation:\n");
2531 for (cs_lnum_t i = 0; i < selection->n_i_adj_faces; i++)
2532 fprintf(logfile, " %9ld | %9ld\n",
2533 (long)i, (long)selection->i_adj_faces[i]);
2534 fprintf(logfile, "\n");
2535 }
2536
2537 fflush(logfile);
2538
2539 bft_printf_flush();
2540
2541 } /* End if verbosity > 2 */
2542
2543 return selection;
2544 }
2545
2546 /*----------------------------------------------------------------------------
2547 * Destroy a cs_join_select_t structure.
2548 *
2549 * parameters:
2550 * param <-- user-defined joining parameters
2551 * join_select <-- pointer to pointer to structure to destroy
2552 *---------------------------------------------------------------------------*/
2553
2554 void
cs_join_select_destroy(cs_join_param_t param,cs_join_select_t ** join_select)2555 cs_join_select_destroy(cs_join_param_t param,
2556 cs_join_select_t **join_select)
2557 {
2558 if (*join_select != NULL) {
2559
2560 cs_join_select_t *_js = *join_select;
2561
2562 BFT_FREE(_js->faces);
2563 BFT_FREE(_js->compact_face_gnum);
2564 BFT_FREE(_js->compact_rank_index);
2565 BFT_FREE(_js->vertices);
2566 BFT_FREE(_js->b_adj_faces);
2567 BFT_FREE(_js->i_adj_faces);
2568
2569 BFT_FREE(_js->b_face_state);
2570 BFT_FREE(_js->i_face_state);
2571
2572 if (param.perio_type != FVM_PERIODICITY_NULL)
2573 BFT_FREE(_js->per_v_couples);
2574
2575 _destroy_join_sync(&(_js->s_vertices));
2576 _destroy_join_sync(&(_js->c_vertices));
2577 _destroy_join_sync(&(_js->s_edges));
2578 _destroy_join_sync(&(_js->c_edges));
2579
2580 BFT_FREE(*join_select);
2581 *join_select = NULL;
2582
2583 }
2584 }
2585
2586 /*----------------------------------------------------------------------------
2587 * Extract vertices from a selection of faces.
2588 *
2589 * parameters:
2590 * n_select_faces <-- number of selected faces
2591 * select_faces <-- list of faces selected
2592 * f2v_idx <-- "face -> vertex" connect. index
2593 * f2v_lst <-- "face -> vertex" connect. list
2594 * n_vertices <-- number of vertices
2595 * n_select_vertices <-> pointer to the number of selected vertices
2596 * select_vertices <-> pointer to the list of selected vertices
2597 *---------------------------------------------------------------------------*/
2598
2599 void
cs_join_extract_vertices(cs_lnum_t n_select_faces,const cs_lnum_t * select_faces,const cs_lnum_t * f2v_idx,const cs_lnum_t * f2v_lst,cs_lnum_t n_vertices,cs_lnum_t * n_select_vertices,cs_lnum_t * select_vertices[])2600 cs_join_extract_vertices(cs_lnum_t n_select_faces,
2601 const cs_lnum_t *select_faces,
2602 const cs_lnum_t *f2v_idx,
2603 const cs_lnum_t *f2v_lst,
2604 cs_lnum_t n_vertices,
2605 cs_lnum_t *n_select_vertices,
2606 cs_lnum_t *select_vertices[])
2607 {
2608 int i, j, face_id;
2609
2610 cs_lnum_t _n_select_vertices = 0;
2611 cs_lnum_t *counter = NULL, *_select_vertices = NULL;
2612
2613 if (n_select_faces > 0) {
2614
2615 BFT_MALLOC(counter, n_vertices, cs_lnum_t);
2616
2617 for (i = 0; i < n_vertices; i++)
2618 counter[i] = 0;
2619
2620 for (i = 0; i < n_select_faces; i++) {
2621
2622 face_id = select_faces[i] - 1;
2623
2624 for (j = f2v_idx[face_id]; j < f2v_idx[face_id+1]; j++)
2625 counter[f2v_lst[j]] = 1;
2626
2627 }
2628
2629 for (i = 0; i < n_vertices; i++)
2630 _n_select_vertices += counter[i];
2631
2632 BFT_MALLOC(_select_vertices, _n_select_vertices, cs_lnum_t);
2633
2634 _n_select_vertices = 0;
2635 for (i = 0; i < n_vertices; i++)
2636 if (counter[i] == 1)
2637 _select_vertices[_n_select_vertices++] = i + 1;
2638
2639 assert(_n_select_vertices > 0);
2640
2641 BFT_FREE(counter);
2642
2643 } /* End if n_select_faces > 0 */
2644
2645 /* Return pointers */
2646
2647 *n_select_vertices = _n_select_vertices;
2648 *select_vertices = _select_vertices;
2649 }
2650
2651 /*----------------------------------------------------------------------------
2652 * Eliminate redundancies found between two lists of elements.
2653 * Delete elements in elts[] and keep elements in the reference list.
2654 *
2655 * parameters:
2656 * n_elts <-> number of elements in the list to clean
2657 * elts <-> list of elements in the list to clean
2658 * n_ref_elts <-- number of elements in the reference list
2659 * ref_elts <-- list of reference elements
2660 *---------------------------------------------------------------------------*/
2661
2662 void
cs_join_clean_selection(cs_lnum_t * n_elts,cs_lnum_t * elts[],cs_lnum_t n_ref_elts,cs_lnum_t ref_elts[])2663 cs_join_clean_selection(cs_lnum_t *n_elts,
2664 cs_lnum_t *elts[],
2665 cs_lnum_t n_ref_elts,
2666 cs_lnum_t ref_elts[])
2667 {
2668 cs_lnum_t i = 0, j = 0;
2669 cs_lnum_t _n_elts = 0;
2670 cs_lnum_t *_elts = *elts;
2671
2672 while (i < *n_elts && j < n_ref_elts) {
2673
2674 if (_elts[i] < ref_elts[j])
2675 _elts[_n_elts++] = _elts[i++];
2676 else if (_elts[i] > ref_elts[j])
2677 j++;
2678 else
2679 i++, j++;
2680
2681 }
2682
2683 for (;i < *n_elts; i++, _n_elts++)
2684 _elts[_n_elts] = _elts[i];
2685
2686 BFT_REALLOC(_elts, _n_elts, cs_lnum_t);
2687
2688 *n_elts = _n_elts;
2689 *elts = _elts;
2690 }
2691
2692 /*----------------------------------------------------------------------------
2693 * Build vertex -> vertex index for a selection of faces.
2694 *
2695 * "v2v_idx" is already allocated to the number of vertices in the mesh.
2696 * At this stage, it is just a counter.
2697 *
2698 * parameters:
2699 * n_faces <-- number of selected faces
2700 * faces <-- list of selected faces
2701 * f2v_idx <-- face -> vertex connectivity index
2702 * f2v_lst <-- face -> vertex connectivity list
2703 * v2v_idx <-> index to build (already allocated and may be used again)
2704 *---------------------------------------------------------------------------*/
2705
2706 void
cs_join_build_edges_idx(cs_lnum_t n_faces,const cs_lnum_t faces[],const cs_lnum_t f2v_idx[],const cs_lnum_t f2v_lst[],cs_lnum_t v2v_idx[])2707 cs_join_build_edges_idx(cs_lnum_t n_faces,
2708 const cs_lnum_t faces[],
2709 const cs_lnum_t f2v_idx[],
2710 const cs_lnum_t f2v_lst[],
2711 cs_lnum_t v2v_idx[])
2712 {
2713 /* Loop on all selected faces. No need to loop on other faces because
2714 the selected vertices are all found with this only step. */
2715
2716 for (cs_lnum_t i = 0; i < n_faces; i++) {
2717
2718 cs_lnum_t v1, v2;
2719
2720 cs_lnum_t fid = faces[i] - 1;
2721 cs_lnum_t s = f2v_idx[fid];
2722 cs_lnum_t e = f2v_idx[fid+1];
2723
2724 for (cs_lnum_t j = s; j < e - 1; j++) { /* scan edges */
2725
2726 v1 = f2v_lst[j] + 1;
2727 v2 = f2v_lst[j+1] + 1;
2728
2729 if (v1 < v2)
2730 v2v_idx[v1] += 1;
2731 else if (v2 < v1)
2732 v2v_idx[v2] += 1;
2733 else
2734 bft_error(__FILE__, __LINE__, 0,
2735 _(" Inconsistent mesh definition. Cannot build edges.\n"
2736 " Face %ld has the same vertex %ld twice.\n"),
2737 (long)fid+1, (long)v1);
2738
2739 }
2740
2741 /* Last edge */
2742
2743 v1 = f2v_lst[e-1] + 1;
2744 v2 = f2v_lst[s] + 1;
2745
2746 if (v1 < v2)
2747 v2v_idx[v1] += 1;
2748 else if (v2 < v1)
2749 v2v_idx[v2] += 1;
2750 else
2751 bft_error(__FILE__, __LINE__, 0,
2752 _(" Inconsistent mesh definition. Cannot build edges.\n"
2753 " Face %ld has the same vertex %ld twice.\n"),
2754 (long)fid+1, (long)v1);
2755
2756 } /* End of loop on selected faces */
2757 }
2758
2759 /*----------------------------------------------------------------------------
2760 * Build vertex -> vertex list for a selection of faces.
2761 * "count" and "v2v_lst" are already allocated to the number of vertices in
2762 * the mesh.
2763 *
2764 * parameters:
2765 * n_faces <-- number of selected faces
2766 * faces <-- list of selected faces
2767 * f2v_idx <-- face -> vertex connectivity index
2768 * f2v_lst <-- face -> vertex connectivity list
2769 * count <-> array used to count the number of values already added
2770 * v2v_idx <-- vertex -> vertex connect. index
2771 * v2v_lst <-> vertex -> vertex connect. list to build (can be used again)
2772 *---------------------------------------------------------------------------*/
2773
2774 void
cs_join_build_edges_lst(cs_lnum_t n_faces,const cs_lnum_t faces[],const cs_lnum_t f2v_idx[],const cs_lnum_t f2v_lst[],cs_lnum_t count[],const cs_lnum_t v2v_idx[],cs_lnum_t v2v_lst[])2775 cs_join_build_edges_lst(cs_lnum_t n_faces,
2776 const cs_lnum_t faces[],
2777 const cs_lnum_t f2v_idx[],
2778 const cs_lnum_t f2v_lst[],
2779 cs_lnum_t count[],
2780 const cs_lnum_t v2v_idx[],
2781 cs_lnum_t v2v_lst[])
2782 {
2783 cs_lnum_t i, j, v1_id, v2_id, fid, s, e, shift;
2784
2785 for (i = 0; i < n_faces; i++) {
2786
2787 fid = faces[i] - 1;
2788 s = f2v_idx[fid];
2789 e = f2v_idx[fid+1];
2790
2791 for (j = s; j < e - 1; j++) { /* Scan edges */
2792
2793 v1_id = f2v_lst[j];
2794 v2_id = f2v_lst[j+1];
2795
2796 if (v1_id < v2_id) {
2797 shift = v2v_idx[v1_id] + count[v1_id];
2798 v2v_lst[shift] = v2_id + 1;
2799 count[v1_id] += 1;
2800 }
2801 else if (v2_id < v1_id) {
2802 shift = v2v_idx[v2_id] + count[v2_id];
2803 v2v_lst[shift] = v1_id + 1;
2804 count[v2_id] += 1;
2805 }
2806
2807 }
2808
2809 /* Last edge */
2810
2811 v1_id = f2v_lst[e-1];
2812 v2_id = f2v_lst[s];
2813
2814 if (v1_id < v2_id) {
2815 shift = v2v_idx[v1_id] + count[v1_id];
2816 v2v_lst[shift] = v2_id + 1;
2817 count[v1_id] += 1;
2818 }
2819 else if (v2_id < v1_id) {
2820 shift = v2v_idx[v2_id] + count[v2_id];
2821 v2v_lst[shift] = v1_id + 1;
2822 count[v2_id] += 1;
2823 }
2824
2825 } /* End of loop on selected faces */
2826
2827 }
2828
2829 /*----------------------------------------------------------------------------*/
2830
2831 END_C_DECLS
2832