1 /*============================================================================
2 * Main structure associated to a mesh
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 <stdio.h>
35 #include <stdlib.h>
36 #include <string.h>
37 #include <limits.h>
38 #include <math.h>
39 #include <assert.h>
40
41 /*----------------------------------------------------------------------------
42 * Local headers
43 *----------------------------------------------------------------------------*/
44
45 #include "bft_mem.h"
46 #include "bft_printf.h"
47
48 #include "fvm_io_num.h"
49 #include "fvm_periodicity.h"
50 #include "fvm_selector.h"
51
52 #include "cs_base.h"
53 #include "cs_halo.h"
54 #include "cs_halo_perio.h"
55 #include "cs_interface.h"
56 #include "cs_log.h"
57 #include "cs_mesh_halo.h"
58 #include "cs_numbering.h"
59 #include "cs_order.h"
60 #include "cs_mesh_quantities.h"
61 #include "cs_parall.h"
62 #include "cs_ext_neighborhood.h"
63 #include "cs_timer.h"
64
65 /*----------------------------------------------------------------------------
66 * Header for the current file
67 *----------------------------------------------------------------------------*/
68
69 #include "cs_mesh.h"
70
71 /*----------------------------------------------------------------------------*/
72
73 BEGIN_C_DECLS
74
75 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
76
77 /*=============================================================================
78 * Local Macro Definitions
79 *============================================================================*/
80
81 #define CS_MESH_N_SUBS 5
82
83 /*============================================================================
84 * Local structure definitions
85 *============================================================================*/
86
87 /*=============================================================================
88 * Local Macro definitions
89 *============================================================================*/
90
91 /*============================================================================
92 * Static global variables
93 *============================================================================*/
94
95 cs_mesh_t *cs_glob_mesh = NULL; /* Pointer on the main mesh */
96
97 /*============================================================================
98 * Private function definitions
99 *============================================================================*/
100
101 /*----------------------------------------------------------------------------
102 * Update cell family array in case of parallelism.
103 *
104 * This function aims at copying main values from cells in halo (id between 1
105 * and n_cells) to ghost cells on distant ranks (id between n_cells + 1 to
106 * n_cells_with_halo).
107 *
108 * parameters:
109 * mesh <-> pointer to mesh structure
110 *----------------------------------------------------------------------------*/
111
112 static void
_sync_cell_fam(cs_mesh_t * const mesh)113 _sync_cell_fam(cs_mesh_t *const mesh)
114 {
115 cs_halo_t *halo = mesh->halo;
116
117 if (halo == NULL)
118 return;
119
120 if (mesh->verbosity > 0)
121 bft_printf(_("Synchronizing cell families\n"));
122
123 cs_halo_sync_untyped(halo, CS_HALO_EXTENDED, sizeof(int), mesh->cell_family);
124 }
125
126 /*----------------------------------------------------------------------------
127 * Compute the minimum and the maximum of a vector (locally).
128 *
129 * parameters:
130 * n_vals <-- local number of elements
131 * var <-- pointer to vector
132 * min --> minimum
133 * max --> maximum
134 *----------------------------------------------------------------------------*/
135
136 static void
_compute_local_minmax(cs_lnum_t n_vals,const cs_lnum_t var[],cs_lnum_t * min,cs_lnum_t * max)137 _compute_local_minmax(cs_lnum_t n_vals,
138 const cs_lnum_t var[],
139 cs_lnum_t *min,
140 cs_lnum_t *max)
141 {
142 cs_lnum_t i;
143 cs_lnum_t _min = var[0], _max = var[0];
144
145 for (i = 1; i < n_vals; i++) {
146 _min = CS_MIN(_min, var[i]);
147 _max = CS_MAX(_max, var[i]);
148 }
149
150 if (min != NULL) *min = _min;
151 if (max != NULL) *max = _max;
152 }
153
154 /*----------------------------------------------------------------------------
155 * Display the distribution of values of a integer vector.
156 *
157 * parameters:
158 * n_vals <-- local number of elements
159 * var <-- pointer to vector
160 *----------------------------------------------------------------------------*/
161
162 static void
_display_histograms(cs_lnum_t n_vals,const cs_lnum_t var[])163 _display_histograms(cs_lnum_t n_vals,
164 const cs_lnum_t var[])
165 {
166 cs_lnum_t i, j, k, val_max, val_min;
167 double step;
168
169 cs_lnum_t count[CS_MESH_N_SUBS];
170 int n_steps = CS_MESH_N_SUBS;
171
172 /* Compute local min and max */
173
174 if (n_vals == 0) {
175 bft_printf(_(" no value\n"));
176 return;
177 }
178
179 val_max = var[0];
180 val_min = var[0];
181 _compute_local_minmax(n_vals, var, &val_min, &val_max);
182
183 bft_printf(_(" minimum value = %10d\n"), (int)val_min);
184 bft_printf(_(" maximum value = %10d\n\n"), (int)val_max);
185
186 /* Define axis subdivisions */
187
188 for (j = 0; j < n_steps; j++)
189 count[j] = 0;
190
191 if (val_max - val_min > 0) {
192
193 if (val_max-val_min < n_steps)
194 n_steps = CS_MAX(1, floor(val_max-val_min));
195
196 step = (double)(val_max - val_min) / n_steps;
197
198 /* Loop on values */
199
200 for (i = 0; i < n_vals; i++) {
201
202 /* Associated subdivision */
203
204 for (j = 0, k = 1; k < n_steps; j++, k++) {
205 if (var[i] < val_min + k*step)
206 break;
207 }
208 count[j] += 1;
209
210 }
211
212 for (i = 0, j = 1; i < n_steps - 1; i++, j++)
213 bft_printf(" %3d : [ %10d ; %10d [ = %10d\n",
214 (int)i+1,
215 (int)(val_min + i*step),
216 (int)(val_min + j*step),
217 (int)(count[i]));
218
219 bft_printf(" %3d : [ %10d ; %10d ] = %10d\n",
220 (int)n_steps,
221 (int)(val_min + (n_steps - 1)*step),
222 (int)val_max,
223 (int)(count[n_steps - 1]));
224
225 }
226
227 else { /* if (val_max == val_min) */
228 bft_printf(" %3d : [ %10d ; %10d ] = %10d\n",
229 1, (int)(val_min), (int)val_max, (int)n_vals);
230 }
231
232 }
233
234 /*----------------------------------------------------------------------------
235 * Write a summary about halo features in log
236 *
237 * parameters:
238 * mesh <-- pointer to cs_mesh_t structure
239 * interface_time <-- time elapsed in interface build
240 * halo_time <-- time elapsed in halo build
241 * ext_neighborhood_time <-- time elapsed in extended neighborhood
242 * connectivity building
243 *----------------------------------------------------------------------------*/
244
245 static void
_print_halo_info(const cs_mesh_t * mesh,double interface_time,double halo_time,double ext_neighborhood_time)246 _print_halo_info(const cs_mesh_t *mesh,
247 double interface_time,
248 double halo_time,
249 double ext_neighborhood_time)
250 {
251 /* Summary of the computional times */
252
253 cs_log_printf(CS_LOG_PERFORMANCE,
254 _("\nHalo creation times summary\n\n"));
255
256 if (mesh->n_domains > 1 || mesh->n_init_perio > 0)
257 cs_log_printf(CS_LOG_PERFORMANCE,
258 _(" Interface creation: %.3g s\n"),
259 interface_time);
260
261 if (mesh->halo_type == CS_HALO_EXTENDED)
262 cs_log_printf(CS_LOG_PERFORMANCE,
263 _(" Extended connectivity creation: %.3g s\n"),
264 ext_neighborhood_time);
265
266 cs_log_printf(CS_LOG_PERFORMANCE,
267 _(" Halo creation: %.3g s\n\n"),
268 halo_time);
269
270 cs_log_printf(CS_LOG_PERFORMANCE,
271 _(" Total time for halo creation: %.3g s\n\n"),
272 halo_time + interface_time + ext_neighborhood_time);
273
274 cs_log_separator(CS_LOG_PERFORMANCE);
275 cs_log_printf_flush(CS_LOG_PERFORMANCE);
276 }
277
278 /*----------------------------------------------------------------------------
279 * Write a summary about cell neighbor features in log
280 *
281 * parameters:
282 * mesh <-- pointer to cs_mesh_t structure
283 *----------------------------------------------------------------------------*/
284
285 static void
_print_cell_neighbor_info(const cs_mesh_t * mesh)286 _print_cell_neighbor_info(const cs_mesh_t *mesh)
287 {
288 cs_lnum_t i, j, k;
289 float step;
290
291 int n_steps = CS_MESH_N_SUBS;
292 int n_min_neighbors = 0;
293 int n_max_neighbors = 0;
294
295 cs_gnum_t count[CS_MESH_N_SUBS];
296
297 int *n_cell_neighbors = NULL;
298
299 /* Summary of the number of cell neighbors */
300
301 BFT_MALLOC(n_cell_neighbors, mesh->n_cells_with_ghosts, int);
302
303 for (i = 0; i < mesh->n_cells_with_ghosts; i++)
304 n_cell_neighbors[i] = 0;
305
306 for (i = 0; i < mesh->n_i_faces; i++) {
307 cs_lnum_t c_id_0 = mesh->i_face_cells[i][0];
308 cs_lnum_t c_id_1 = mesh->i_face_cells[i][1];
309 n_cell_neighbors[c_id_0] += 1;
310 n_cell_neighbors[c_id_1] += 1;
311 }
312
313 if (mesh->n_cells > 0)
314 n_min_neighbors = n_cell_neighbors[0];
315
316 for (i = 0; i < mesh->n_cells; i++) {
317 if (n_cell_neighbors[i] < n_min_neighbors)
318 n_min_neighbors = n_cell_neighbors[i];
319 else if (n_cell_neighbors[i] > n_max_neighbors)
320 n_max_neighbors = n_cell_neighbors[i];
321 }
322
323 #if defined(HAVE_MPI)
324
325 if (cs_glob_n_ranks > 1) {
326
327 int n_g_min, n_g_max;
328
329 MPI_Allreduce(&n_min_neighbors, &n_g_min, 1, MPI_INT, MPI_MIN,
330 cs_glob_mpi_comm);
331 MPI_Allreduce(&n_max_neighbors, &n_g_max, 1, MPI_INT, MPI_MAX,
332 cs_glob_mpi_comm);
333
334 n_min_neighbors = n_g_min;
335 n_max_neighbors = n_g_max;
336 }
337
338 #endif /* defined(HAVE_MPI) */
339
340 bft_printf(_("\n Histogram of the number of interior faces per cell:\n\n"));
341
342 bft_printf(_(" minimum value = %10d\n"), n_min_neighbors);
343 bft_printf(_(" maximum value = %10d\n\n"), n_max_neighbors);
344
345 /* Define axis subdivisions */
346
347 for (i = 0; i < CS_MESH_N_SUBS; i++)
348 count[i] = 0;
349
350 if (n_max_neighbors - n_min_neighbors > 0) {
351
352 if (n_max_neighbors - n_min_neighbors < n_steps)
353 n_steps = n_max_neighbors - n_min_neighbors;
354
355 step = (float)(n_max_neighbors - n_min_neighbors) / n_steps;
356
357 /* Loop on values */
358
359 for (i = 0; i < mesh->n_cells; i++) {
360
361 /* Associated subdivision */
362 for (j = 0, k = 1; k < n_steps; j++, k++) {
363 if (n_cell_neighbors[i] < n_min_neighbors + k*step)
364 break;
365 }
366 count[j] += 1;
367 }
368
369 #if defined(HAVE_MPI)
370
371 if (cs_glob_n_ranks > 1) {
372
373 cs_gnum_t g_count[CS_MESH_N_SUBS];
374
375 MPI_Allreduce(count, g_count, n_steps, CS_MPI_GNUM, MPI_SUM,
376 cs_glob_mpi_comm);
377
378 for (i = 0; i < n_steps; i++)
379 count[i] = g_count[i];
380 }
381
382 #endif
383
384 for (i = 0, j = 1; i < n_steps - 1; i++, j++)
385 bft_printf(" %3d : [ %10d ; %10d [ = %10llu\n",
386 (int)i+1,
387 (int)(n_min_neighbors + i*step),
388 (int)(n_min_neighbors + j*step),
389 (unsigned long long)(count[i]));
390
391 bft_printf(" %3d : [ %10d ; %10d ] = %10llu\n",
392 n_steps,
393 (int)(n_min_neighbors + (n_steps - 1)*step),
394 n_max_neighbors,
395 (unsigned long long)(count[n_steps - 1]));
396 }
397
398 else { /* if (n_min_neighbors == n_max_neighbors) */
399 bft_printf(" %3d : [ %10d ; %10d ] = %10llu\n",
400 1, n_min_neighbors, n_max_neighbors,
401 (unsigned long long)mesh->n_g_cells);
402
403 }
404
405 bft_printf("\n ----------------------------------------------------------\n");
406
407 /* Cleanup */
408
409 BFT_FREE(n_cell_neighbors);
410 }
411
412 /*----------------------------------------------------------------------------
413 * Compare periodic couples in global numbering form (qsort function).
414 *
415 * parameters:
416 * x <-> pointer to first couple
417 * y <-> pointer to second couple
418 *
419 * returns:
420 * lexicographical
421 *----------------------------------------------------------------------------*/
422
423 static int
_compare_couples(const void * x,const void * y)424 _compare_couples(const void *x,
425 const void *y)
426 {
427 int retval = 1;
428
429 const cs_gnum_t *c0 = x;
430 const cs_gnum_t *c1 = y;
431
432 if (c0[0] < c1[0])
433 retval = -1;
434
435 else if (c0[0] == c1[0]) {
436 if (c0[1] < c1[1])
437 retval = -1;
438 else if (c0[1] == c1[1])
439 retval = 0;
440 }
441
442 return retval;
443 }
444
445 /*----------------------------------------------------------------------------
446 * Define parameters for building an vertices interface set structure on
447 * a given mesh.
448 *
449 * parameters:
450 * mesh <-- pointer to mesh structure
451 * face_ifs <-- pointer to face interface structure describing
452 * periodic face couples
453 * p_n_perio_couples --> pointer to the number of periodic couples
454 * p_perio_couples --> pointer to the periodic couple list
455 *----------------------------------------------------------------------------*/
456
457 static void
_define_perio_vtx_couples(const cs_mesh_t * mesh,cs_interface_set_t * face_ifs,cs_lnum_t * p_n_perio_couples[],cs_gnum_t ** p_perio_couples[])458 _define_perio_vtx_couples(const cs_mesh_t *mesh,
459 cs_interface_set_t *face_ifs,
460 cs_lnum_t *p_n_perio_couples[],
461 cs_gnum_t **p_perio_couples[])
462 {
463 cs_lnum_t *n_perio_couples = NULL;
464 cs_gnum_t **perio_couples = NULL;
465
466 int i, j;
467 cs_lnum_t k, l;
468 cs_lnum_t itf_start;
469
470 int perio_count = 0;
471 cs_lnum_t *send_index = NULL;
472 cs_lnum_t *recv_index = NULL;
473 cs_gnum_t *send_num = NULL, *recv_num = NULL;
474 int *tr_id = NULL;
475
476 const int n_perio = mesh->n_init_perio;
477 const int n_interfaces = cs_interface_set_size(face_ifs);
478 const cs_lnum_t n_ifs_faces = cs_interface_set_n_elts(face_ifs);
479 const cs_gnum_t *vtx_gnum = mesh->global_vtx_num;
480
481 /* List direct and reverse transforms */
482
483 BFT_MALLOC(n_perio_couples, n_perio, cs_lnum_t);
484 BFT_MALLOC(perio_couples, n_perio, cs_gnum_t *);
485
486 for (i = 0; i < n_perio; i++) {
487 n_perio_couples[i] = 0;
488 perio_couples[i] = NULL;
489 }
490
491 BFT_MALLOC(tr_id, n_perio*2, int);
492
493 for (i = 0; i < n_perio*2; i++) {
494 int rev_id = fvm_periodicity_get_reverse_id(mesh->periodicity, i);
495 if (i < rev_id) {
496 int parent_ids[2];
497 fvm_periodicity_get_parent_ids(mesh->periodicity, i, parent_ids);
498 if (parent_ids[0] < 0 && parent_ids[1] < 0) {
499 tr_id[perio_count*2] = i + 1;
500 tr_id[perio_count*2 + 1] = rev_id + 1;
501 perio_count++;
502 }
503 }
504 }
505 assert(perio_count == n_perio);
506
507 /* Count maximum vertex couples and build exchange index;
508 count only couples in direct periodicity direction, as the number
509 of faces in positive and reverse direction should match */
510
511 BFT_MALLOC(send_index, n_ifs_faces + 1, cs_lnum_t);
512 BFT_MALLOC(recv_index, n_ifs_faces + 1, cs_lnum_t);
513
514 for (k = 0; k < n_ifs_faces + 1; k++) {
515 send_index[k] = 0;
516 recv_index[k] = 0;
517 }
518
519 for (i = 0, itf_start = 0; i < n_interfaces; i++) {
520
521 const cs_interface_t *face_if = cs_interface_set_get(face_ifs, i);
522 const cs_lnum_t *tr_index = cs_interface_get_tr_index(face_if);
523 const cs_lnum_t *loc_id = cs_interface_get_elt_ids(face_if);
524
525 for (j = 0; j < n_perio; j++) {
526
527 const cs_lnum_t tr_start = tr_index[tr_id[j*2]];
528 const cs_lnum_t tr_end = tr_index[tr_id[j*2] + 1];
529 const cs_lnum_t rev_start = tr_index[tr_id[j*2+1]];
530 const cs_lnum_t rev_end = tr_index[tr_id[j*2+1] + 1];
531
532 for (k = rev_start; k < rev_end; k++) {
533 cs_lnum_t l_id = loc_id[k];
534 send_index[itf_start + k + 1]
535 = (mesh->i_face_vtx_idx[l_id+1] - mesh->i_face_vtx_idx[l_id]);
536 }
537 for (k = tr_start; k < tr_end; k++) {
538 cs_lnum_t l_id = loc_id[k];
539 cs_lnum_t n_vtx = ( mesh->i_face_vtx_idx[l_id+1]
540 - mesh->i_face_vtx_idx[l_id]);
541 recv_index[itf_start + k + 1] = n_vtx;
542 n_perio_couples[j] += n_vtx;
543 }
544 }
545
546 itf_start += cs_interface_size(face_if);
547
548 }
549
550 /* Transform count to index */
551
552 for (k = 0; k < n_ifs_faces; k++) {
553 send_index[k+1] += send_index[k];
554 recv_index[k+1] += recv_index[k];
555 }
556
557 BFT_MALLOC(send_num, send_index[n_ifs_faces], cs_gnum_t);
558 BFT_MALLOC(recv_num, recv_index[n_ifs_faces], cs_gnum_t);
559
560 /* Prepare send buffer (send reverse transformation values) */
561
562 for (i = 0, j = 0, l = 0; i < n_interfaces; i++) {
563
564 const cs_interface_t *face_if = cs_interface_set_get(face_ifs, i);
565 const cs_lnum_t *tr_index = cs_interface_get_tr_index(face_if);
566 const cs_lnum_t *loc_id = cs_interface_get_elt_ids(face_if);
567
568 for (j = 0; j < n_perio; j++) {
569
570 const cs_lnum_t start_id = tr_index[tr_id[j*2+1]];
571 const cs_lnum_t end_id = tr_index[tr_id[j*2+1] + 1];
572
573 if (vtx_gnum != NULL) {
574
575 for (k = start_id; k < end_id; k++) {
576 cs_lnum_t v_id;
577 cs_lnum_t f_id = loc_id[k];
578 for (v_id = mesh->i_face_vtx_idx[f_id];
579 v_id < mesh->i_face_vtx_idx[f_id + 1];
580 v_id++)
581 send_num[l++] = vtx_gnum[mesh->i_face_vtx_lst[v_id]];
582 }
583
584 }
585 else {
586
587 for (k = start_id; k < end_id; k++) {
588 cs_lnum_t v_id;
589 cs_lnum_t f_id = loc_id[k];
590 for (v_id = mesh->i_face_vtx_idx[f_id];
591 v_id < mesh->i_face_vtx_idx[f_id + 1];
592 v_id++)
593 send_num[l++] = mesh->i_face_vtx_lst[v_id] + 1;
594 }
595
596 }
597 }
598 }
599
600 /* Exchange global vertex numbers */
601
602 cs_interface_set_copy_indexed(face_ifs,
603 CS_GNUM_TYPE,
604 false, /* src_on_parent */
605 send_index,
606 recv_index,
607 send_num,
608 recv_num);
609
610 BFT_FREE(send_num);
611 BFT_FREE(recv_index);
612 BFT_FREE(send_index);
613
614 for (i = 0; i < n_perio; i++)
615 BFT_MALLOC(perio_couples[i], n_perio_couples[i]*2, cs_gnum_t);
616
617 /* Reset couples count */
618
619 for (i = 0; i < n_perio; i++)
620 n_perio_couples[i] = 0;
621
622 /* Now fill couples */
623
624 for (i = 0, j = 0, l = 0; i < n_interfaces; i++) {
625
626 const cs_interface_t *face_if = cs_interface_set_get(face_ifs, i);
627 const cs_lnum_t *tr_index = cs_interface_get_tr_index(face_if);
628 const cs_lnum_t *loc_ids = cs_interface_get_elt_ids(face_if);
629
630 for (j = 0; j < n_perio; j++) {
631
632 cs_lnum_t nc = n_perio_couples[j]*2;
633 const cs_lnum_t start_id = tr_index[tr_id[j*2]];
634 const cs_lnum_t end_id = tr_index[tr_id[j*2] + 1];
635
636 if (vtx_gnum != NULL) {
637
638 for (k = start_id; k < end_id; k++) {
639 cs_lnum_t v_id;
640 cs_lnum_t f_id = loc_ids[k];
641 for (v_id = mesh->i_face_vtx_idx[f_id];
642 v_id < mesh->i_face_vtx_idx[f_id + 1];
643 v_id++) {
644 perio_couples[j][nc++] = vtx_gnum[mesh->i_face_vtx_lst[v_id]];
645 perio_couples[j][nc++] = recv_num[l++];
646 }
647 }
648
649 }
650 else {
651
652 for (k = start_id; k < end_id; k++) {
653 cs_lnum_t v_id;
654 cs_lnum_t f_id = loc_ids[k];
655 for (v_id = mesh->i_face_vtx_idx[f_id];
656 v_id < mesh->i_face_vtx_idx[f_id + 1];
657 v_id++) {
658 perio_couples[j][nc++] = mesh->i_face_vtx_lst[v_id] + 1;
659 perio_couples[j][nc++] = recv_num[l++];
660 }
661 }
662 }
663
664 n_perio_couples[j] = nc/2;
665 }
666 }
667
668 BFT_FREE(recv_num);
669 BFT_FREE(tr_id);
670
671 /* Finally, remove duplicate couples */
672
673 for (i = 0; i < mesh->n_init_perio; i++) {
674
675 if (n_perio_couples[i] > 1) {
676
677 cs_lnum_t nc = n_perio_couples[i];
678 cs_gnum_t *pc = perio_couples[i];
679
680 qsort(pc, nc, sizeof(cs_gnum_t) * 2, &_compare_couples);
681
682 for (k = 1, l = 0; k < nc; k++) {
683 if ((pc[k*2] != pc[l*2]) && (pc[k*2+1] != pc[l*2+1])) {
684 l += 1;
685 pc[l*2] = pc[k*2];
686 pc[l*2+1] = pc[k*2+1];
687 }
688 }
689 n_perio_couples[i] = l+1;
690
691 BFT_REALLOC(perio_couples[i], n_perio_couples[i]*2, cs_gnum_t);
692 }
693 }
694
695 /* Set return values */
696
697 *p_n_perio_couples = n_perio_couples;
698 *p_perio_couples = perio_couples;
699 }
700
701 /*----------------------------------------------------------------------------
702 * Assign a periodicity number and optionnaly a rank id to halo elements.
703 *
704 * parameters:
705 * mesh <-- pointer to mesh structure
706 * halo_perio_num --> periodicity number associated with each halo cell,
707 * signed for direct/reverse transform.
708 * (size: mesh->n_ghost_cells)
709 * halo_rank_id --> halo local rank id associated with each halo cell,
710 * or NULL (size: mesh->n_ghost_cells)
711 *
712 * returns:
713 * the number of base periodicities defined
714 *----------------------------------------------------------------------------*/
715
716 static int
_get_halo_perio_num(const cs_mesh_t * mesh,int halo_perio_num[],int halo_rank_id[])717 _get_halo_perio_num(const cs_mesh_t *mesh,
718 int halo_perio_num[],
719 int halo_rank_id[])
720 {
721 int i, rank_id;
722 cs_lnum_t j;
723
724 int n_perio = 0;
725
726 int *tr_id = NULL;
727
728 const cs_halo_t *halo = mesh->halo;
729 const fvm_periodicity_t *periodicity = mesh->periodicity;
730
731 /* Initialization */
732
733 assert(halo != NULL);
734 assert(periodicity != NULL);
735
736 /* List direct and non-combined transforms */
737
738 BFT_MALLOC(tr_id, halo->n_transforms, int);
739 for (i = 0; i < halo->n_transforms; i++)
740 tr_id[i] = -1;
741
742 n_perio = 0;
743
744 for (i = 0; i < halo->n_transforms; i++) {
745 int rev_id = fvm_periodicity_get_reverse_id(periodicity, i);
746 if (i < rev_id) {
747 int parent_ids[2];
748 fvm_periodicity_get_parent_ids(periodicity, i, parent_ids);
749 if (parent_ids[0] < 0 && parent_ids[1] < 0) {
750 tr_id[n_perio*2] = i;
751 tr_id[n_perio*2+1] = rev_id;
752 n_perio++;
753 }
754 }
755 }
756 assert(n_perio == mesh->n_init_perio);
757
758 BFT_REALLOC(tr_id, n_perio*2, int);
759 assert(n_perio == mesh->n_init_perio);
760
761 /* Initialize halo values */
762
763 for (j = 0; j < mesh->n_ghost_cells; j++)
764 halo_perio_num[j] = 0;
765
766 if (halo_rank_id != NULL) {
767 for (j = 0; j < mesh->n_ghost_cells; j++)
768 halo_rank_id[j] = -1;
769 }
770
771 /* For each periodic face couple, we have 2 global ghost cells: one for
772 the direct periodicity, one for the reverse. To match couples only
773 once, we ignore the reverse periodicity, leaving exactly 1 ghost cell
774 per couple (which may belong to multiple periodic faces). */
775
776 /* Mark ghost cells with their periodicity number */
777
778 for (i = 0; i < n_perio; i++) {
779
780 cs_lnum_t shift, start, end;
781
782 /* Direct periodicity */
783
784 shift = 4 * halo->n_c_domains * tr_id[i*2];
785
786 for (rank_id = 0; rank_id < halo->n_c_domains; rank_id++) {
787
788 start = halo->perio_lst[shift + 4*rank_id];
789 end = start + halo->perio_lst[shift + 4*rank_id + 1];
790
791 for (j = start; j < end; j++)
792 halo_perio_num[j] = i + 1;
793
794 if (halo_rank_id != NULL) {
795 for (j = start; j < end; j++)
796 halo_rank_id[j] = rank_id;
797 }
798 }
799
800 /* Reverse periodicity */
801
802 shift = 4 * halo->n_c_domains * tr_id[i*2 + 1];
803
804 for (rank_id = 0; rank_id < halo->n_c_domains; rank_id++) {
805
806 start = halo->perio_lst[shift + 4*rank_id];
807 end = start + halo->perio_lst[shift + 4*rank_id + 1];
808
809 for (j = start; j < end; j++)
810 halo_perio_num[j] = -(i + 1);
811
812 if (halo_rank_id != NULL) {
813 for (j = start; j < end; j++)
814 halo_rank_id[j] = rank_id;
815 }
816 }
817 }
818
819 BFT_FREE(tr_id);
820
821 return n_perio;
822 }
823
824 /*----------------------------------------------------------------------------
825 * Build a cell -> face search structure for local cells with a periodic
826 * face in the direct peridicity direction.
827 *
828 * The caller is responsible for freeing the arrays allocated and returned
829 * by this function once they are no longer needed.
830 *
831 * parameters:
832 * mesh <-- pointer to mesh structure
833 * halo_perio_num <-- periodicity number of ghost cells
834 * cell_face_idx --> index (0 to n-1) of face ids for a given cell
835 * cell_face --> face ids (0 to n-1) for a given cell
836 *----------------------------------------------------------------------------*/
837
838 static void
_build_cell_face_perio_match(const cs_mesh_t * mesh,const int * halo_perio_num,cs_lnum_t ** cell_face_idx,cs_lnum_t ** cell_face)839 _build_cell_face_perio_match(const cs_mesh_t *mesh,
840 const int *halo_perio_num,
841 cs_lnum_t **cell_face_idx,
842 cs_lnum_t **cell_face)
843 {
844 cs_lnum_t i;
845
846 cs_lnum_t *_cell_face_count, *_cell_face_idx, *_cell_face;
847
848 /* Build cell -> face indexed search array containing only faces
849 adjacent to the halo. */
850
851 BFT_MALLOC(_cell_face_count, mesh->n_cells_with_ghosts, cs_lnum_t);
852 BFT_MALLOC(_cell_face_idx, mesh->n_cells_with_ghosts + 1, cs_lnum_t);
853
854 for (i = 0; i < mesh->n_cells_with_ghosts; i++)
855 _cell_face_count[i] = 0;
856
857 for (i = 0; i < mesh->n_i_faces; i++) {
858 const cs_lnum_t c_id0 = mesh->i_face_cells[i][0];
859 const cs_lnum_t c_id1 = mesh->i_face_cells[i][1];
860 if (c_id0 < mesh->n_cells && c_id1 >= mesh->n_cells) {
861 if (halo_perio_num[c_id1 - mesh->n_cells] < 0)
862 _cell_face_count[c_id0] += 1;
863 }
864 else if (c_id1 < mesh->n_cells && c_id0 >= mesh->n_cells) {
865 if (halo_perio_num[c_id0 - mesh->n_cells] < 0)
866 _cell_face_count[c_id1] += 1;
867 }
868 }
869
870 _cell_face_idx[0] = 0;
871 for (i = 0; i < mesh->n_cells_with_ghosts; i++) {
872 _cell_face_idx[i+1] = _cell_face_idx[i] + _cell_face_count[i];
873 _cell_face_count[i] = 0;
874 }
875
876 BFT_MALLOC(_cell_face, _cell_face_idx[mesh->n_cells_with_ghosts], cs_lnum_t);
877
878 for (i = 0; i < mesh->n_i_faces; i++) {
879 const cs_lnum_t c_id0 = mesh->i_face_cells[i][0];
880 const cs_lnum_t c_id1 = mesh->i_face_cells[i][1];
881 if (c_id0 < mesh->n_cells && c_id1 >= mesh->n_cells) {
882 if (halo_perio_num[c_id1 - mesh->n_cells] < 0) {
883 _cell_face[_cell_face_idx[c_id0] + _cell_face_count[c_id0]] = i;
884 _cell_face_count[c_id0] += 1;
885 }
886 }
887 else if (c_id1 < mesh->n_cells && c_id0 >= mesh->n_cells) {
888 if (halo_perio_num[c_id0 - mesh->n_cells] < 0) {
889 _cell_face[_cell_face_idx[c_id1] + _cell_face_count[c_id1]] = i;
890 _cell_face_count[c_id1] += 1;
891 }
892 }
893 }
894
895 BFT_FREE(_cell_face_count);
896
897 *cell_face_idx = _cell_face_idx;
898 *cell_face = _cell_face;
899 }
900
901 #if defined(HAVE_MPI)
902
903 /*----------------------------------------------------------------------------
904 * Get global lists of periodic faces.
905 *
906 * The caller is responsible for freeing the arrays allocated and returned
907 * by this function once they are no longer needed.
908 *
909 * parameters:
910 * mesh <-- pointer to mesh structure
911 * n_perio_face_couples --> local number of periodic couples per
912 * periodicity (size: mesh->n_init_perio)
913 * perio_face_couples --> arrays of global periodic couple face numbers,
914 * for each periodicity
915 *----------------------------------------------------------------------------*/
916
917 static void
_get_perio_faces_g(const cs_mesh_t * mesh,cs_lnum_t ** n_perio_face_couples,cs_gnum_t *** perio_face_couples)918 _get_perio_faces_g(const cs_mesh_t *mesh,
919 cs_lnum_t **n_perio_face_couples,
920 cs_gnum_t ***perio_face_couples)
921 {
922 int i, rank_id, perio_num;
923 cs_lnum_t j, k, l;
924 cs_lnum_t dest_c_id, src_c_id;
925
926 cs_lnum_t recv_size = 0;
927 cs_lnum_t *loc_cell_num = NULL;
928 cs_lnum_t *recv_count = NULL, *recv_idx = NULL;
929 cs_lnum_t *send_count = NULL, *send_idx = NULL;
930 cs_lnum_t *cell_face_idx, *cell_face;
931 cs_gnum_t *recv_vals = NULL, *send_vals = NULL;
932
933 int *halo_perio_num = NULL;
934 int *halo_rank_id = NULL;
935
936 cs_lnum_t *_n_perio_face_couples = NULL;
937 cs_gnum_t **_perio_face_couples = NULL;
938
939 int request_count = 0;
940 MPI_Request *halo_request = NULL;
941 MPI_Status *halo_status = NULL;
942
943 const cs_halo_t *halo = mesh->halo;
944
945 /* Initialization */
946
947 BFT_MALLOC(_n_perio_face_couples, mesh->n_init_perio, cs_lnum_t);
948 BFT_MALLOC(_perio_face_couples, mesh->n_init_perio, cs_gnum_t *);
949
950 for (i = 0; i < mesh->n_init_perio; i++) {
951 _n_perio_face_couples[i] = 0;
952 _perio_face_couples[i] = NULL;
953 }
954
955 assert(halo != NULL);
956
957 /* Mark ghost cells with their periodicity number and halo rank id */
958
959 BFT_MALLOC(halo_perio_num, mesh->n_ghost_cells, int);
960 BFT_MALLOC(halo_rank_id, mesh->n_ghost_cells, int);
961
962 _get_halo_perio_num(mesh, halo_perio_num, halo_rank_id);
963
964 /* Exchange local cell numbers, so that cell-face reverse matching
965 is made easier later */
966
967 BFT_MALLOC(loc_cell_num, mesh->n_cells_with_ghosts, cs_lnum_t);
968 for (j = 0; j < mesh->n_cells; j++)
969 loc_cell_num[j] = j+1;
970 for (j = mesh->n_cells; j < mesh->n_cells_with_ghosts; j++)
971 loc_cell_num[j] = 0;
972
973 cs_halo_sync_num(halo, CS_HALO_STANDARD, loc_cell_num);
974
975 /* Now compute send and receive counts and build indexes */
976
977 BFT_MALLOC(send_count, halo->n_c_domains, cs_lnum_t);
978 BFT_MALLOC(recv_count, halo->n_c_domains, cs_lnum_t);
979 BFT_MALLOC(send_idx, halo->n_c_domains + 1, cs_lnum_t);
980 BFT_MALLOC(recv_idx, halo->n_c_domains + 1, cs_lnum_t);
981
982 for (i = 0; i < halo->n_c_domains; i++) {
983 send_count[i] = 0;
984 recv_count[i] = 0;
985 }
986
987 for (j = 0; j < mesh->n_i_faces; j++) {
988 const cs_lnum_t h_id0 = mesh->i_face_cells[j][0] - mesh->n_cells;
989 const cs_lnum_t h_id1 = mesh->i_face_cells[j][1] - mesh->n_cells;
990 if (h_id0 >= 0) {
991 perio_num = halo_perio_num[h_id0];
992 if (perio_num > 0)
993 send_count[halo_rank_id[h_id0]] += 1;
994 else if (perio_num < 0)
995 recv_count[halo_rank_id[h_id0]] += 1;
996 }
997 else if (h_id1 >= 0) {
998 perio_num = halo_perio_num[h_id1];
999 if (perio_num > 0)
1000 send_count[halo_rank_id[h_id1]] += 1;
1001 else if (perio_num < 0)
1002 recv_count[halo_rank_id[h_id1] ]+= 1;
1003 }
1004 }
1005
1006 send_idx[0] = 0;
1007 recv_idx[0] = 0;
1008 for (i = 0; i < halo->n_c_domains; i++) {
1009 send_idx[i+1] = send_idx[i] + send_count[i];
1010 recv_idx[i+1] = recv_idx[i] + recv_count[i];
1011 send_count[i] = 0;
1012 }
1013
1014 /* Assemble the tuples to send; a tuple contains:
1015 destination cell id (cell matched by direct periodicity),
1016 source cell (cell matched by reverse periodicity),
1017 periodicity number,
1018 global face number (2nd face in a couple). */
1019
1020 BFT_MALLOC(send_vals, send_idx[halo->n_c_domains]*4, cs_gnum_t);
1021
1022 for (j = 0; j < mesh->n_i_faces; j++) {
1023
1024 const cs_lnum_t c_id0 = mesh->i_face_cells[j][0];
1025 const cs_lnum_t c_id1 = mesh->i_face_cells[j][1];
1026
1027 dest_c_id = -1;
1028 if (c_id0 >= mesh->n_cells) {
1029 perio_num = halo_perio_num[c_id0 - mesh->n_cells];
1030 if (perio_num > 0) {
1031 dest_c_id = c_id1;
1032 src_c_id = c_id0;
1033 rank_id = halo_rank_id[c_id0 - mesh->n_cells];
1034 }
1035 }
1036 else if (c_id1 >= mesh->n_cells) {
1037 perio_num = halo_perio_num[c_id1 - mesh->n_cells];
1038 if (perio_num > 0) {
1039 dest_c_id = c_id0;
1040 src_c_id = c_id1;
1041 rank_id = halo_rank_id[c_id1 - mesh->n_cells];
1042 }
1043 }
1044 if (dest_c_id != -1) {
1045 k = send_idx[rank_id] + send_count[rank_id];
1046 send_vals[k*4] = loc_cell_num[dest_c_id];
1047 send_vals[k*4+1] = loc_cell_num[src_c_id];
1048 send_vals[k*4+2] = perio_num;
1049 send_vals[k*4+3] = mesh->global_i_face_num[j];
1050 send_count[rank_id] += 1;
1051 }
1052 }
1053
1054 BFT_FREE(halo_rank_id);
1055
1056 /* Now exchange face matching data */
1057
1058 BFT_MALLOC(halo_request, halo->n_c_domains*2, MPI_Request);
1059 BFT_MALLOC(halo_status, halo->n_c_domains*2, MPI_Status);
1060
1061 BFT_MALLOC(recv_vals, recv_idx[halo->n_c_domains]*4, cs_gnum_t);
1062
1063 /* Receive data from distant ranks */
1064
1065 for (rank_id = 0; rank_id < halo->n_c_domains; rank_id++) {
1066 if (recv_count[rank_id] > 0)
1067 MPI_Irecv(recv_vals + (recv_idx[rank_id]*4),
1068 recv_count[rank_id]*4,
1069 CS_MPI_GNUM,
1070 halo->c_domain_rank[rank_id],
1071 halo->c_domain_rank[rank_id],
1072 cs_glob_mpi_comm,
1073 &(halo_request[request_count++]));
1074 }
1075
1076 /* Send data to distant ranks */
1077
1078 for (rank_id = 0; rank_id < halo->n_c_domains; rank_id++) {
1079 if (send_count[rank_id] > 0)
1080 MPI_Isend(send_vals + (send_idx[rank_id]*4),
1081 send_count[rank_id]*4,
1082 CS_MPI_GNUM,
1083 halo->c_domain_rank[rank_id],
1084 cs_glob_rank_id,
1085 cs_glob_mpi_comm,
1086 &(halo_request[request_count++]));
1087 }
1088
1089 /* Wait for all exchanges */
1090
1091 MPI_Waitall(request_count, halo_request, halo_status);
1092
1093 BFT_FREE(halo_request);
1094 BFT_FREE(halo_status);
1095
1096 BFT_FREE(send_vals);
1097 BFT_FREE(send_count);
1098 BFT_FREE(send_idx);
1099
1100 recv_size = recv_idx[halo->n_c_domains];
1101
1102 BFT_FREE(recv_count);
1103 BFT_FREE(recv_idx);
1104
1105 /* Count periodic face couples */
1106
1107 for (j = 0; j < recv_size; j++) {
1108 perio_num = recv_vals[j*4+2];
1109 assert(perio_num > 0);
1110 _n_perio_face_couples[perio_num - 1] += 1;
1111 }
1112
1113 /* Allocate periodic couple lists */
1114
1115 for (i = 0; i < mesh->n_init_perio; i++) {
1116 BFT_MALLOC(_perio_face_couples[i],
1117 _n_perio_face_couples[i]*2,
1118 cs_gnum_t);
1119 _n_perio_face_couples[i] = 0;
1120 }
1121
1122 /* Build cell -> face indexed search array containing only faces
1123 adjacent to the halo. */
1124
1125 _build_cell_face_perio_match(mesh,
1126 halo_perio_num,
1127 &cell_face_idx,
1128 &cell_face);
1129
1130 /* Now match faces; faces already matched are marked so as not to be
1131 re-matched, in case faces have been split and multiple faces thus share
1132 the same cells (in which case we assume that subfaces are locally ordered
1133 identically). */
1134
1135 for (j = 0; j < recv_size; j++) {
1136
1137 cs_lnum_t dest_c_num = recv_vals[j*4];
1138 cs_lnum_t src_c_num = recv_vals[j*4+1];
1139 int perio_id = recv_vals[j*4+2] - 1;
1140 cs_gnum_t perio_face_gnum = recv_vals[j*4+3];
1141
1142 for (k = cell_face_idx[src_c_num-1]; k < cell_face_idx[src_c_num]; k++) {
1143
1144 l = cell_face[k]; /* Face _id */
1145
1146 if (l > -1) {
1147 cs_lnum_t loc_c_num, dist_c_id;
1148 if (mesh->i_face_cells[l][0] < mesh->i_face_cells[l][1]) {
1149 loc_c_num = mesh->i_face_cells[l][0] + 1;
1150 dist_c_id = mesh->i_face_cells[l][1];
1151 }
1152 else {
1153 loc_c_num = mesh->i_face_cells[l][1] + 1;
1154 dist_c_id = mesh->i_face_cells[l][0];
1155 }
1156 assert(loc_c_num <= mesh->n_cells && dist_c_id >= mesh->n_cells);
1157
1158 /* If we have a match, update periodic couples */
1159
1160 if ( src_c_num == loc_c_num
1161 && dest_c_num == loc_cell_num[dist_c_id]
1162 && halo_perio_num[dist_c_id-mesh->n_cells] == - (perio_id + 1)) {
1163 cs_lnum_t couple_id = _n_perio_face_couples[perio_id];
1164 cell_face[k] = -1; /* Mark as used */
1165 _perio_face_couples[perio_id][couple_id*2]
1166 = mesh->global_i_face_num[l];
1167 _perio_face_couples[perio_id][couple_id*2 + 1] = perio_face_gnum;
1168 _n_perio_face_couples[perio_id] += 1;
1169 }
1170 }
1171 }
1172 }
1173
1174 BFT_FREE(halo_perio_num);
1175
1176 BFT_FREE(cell_face_idx);
1177 BFT_FREE(cell_face);
1178 BFT_FREE(recv_vals);
1179 BFT_FREE(loc_cell_num);
1180
1181 /* Set return values */
1182
1183 *n_perio_face_couples = _n_perio_face_couples;
1184 *perio_face_couples = _perio_face_couples;
1185 }
1186
1187 #endif /* defined(HAVE_MPI) */
1188
1189 /*----------------------------------------------------------------------------
1190 * Get global lists of periodic faces.
1191 *
1192 * This version is used for serial mode.
1193 *
1194 * The caller is responsible for freeing the arrays allocated and returned
1195 * by this function once they are no longer needed.
1196 *
1197 * parameters:
1198 * mesh <-- pointer to mesh structure
1199 * n_perio_face_couples --> global number of periodic couples per
1200 * periodicity (size: mesh->n_init_perio)
1201 * perio_face_couples --> arrays of global periodic couple face numbers,
1202 * for each periodicity
1203 *----------------------------------------------------------------------------*/
1204
1205 static void
_get_perio_faces_l(const cs_mesh_t * mesh,cs_lnum_t ** n_perio_face_couples,cs_gnum_t *** perio_face_couples)1206 _get_perio_faces_l(const cs_mesh_t *mesh,
1207 cs_lnum_t **n_perio_face_couples,
1208 cs_gnum_t ***perio_face_couples)
1209 {
1210 int i, perio_num;
1211 cs_lnum_t j, k, l;
1212 cs_lnum_t dest_c_id, src_c_id;
1213
1214 cs_lnum_t *loc_cell_num = NULL;
1215 cs_lnum_t *cell_face_idx, *cell_face;
1216
1217 int *halo_perio_num = NULL;
1218
1219 cs_lnum_t *_n_perio_face_couples = NULL;
1220 cs_gnum_t **_perio_face_couples = NULL;
1221
1222 const cs_halo_t *halo = mesh->halo;
1223
1224 /* Initialization */
1225
1226 BFT_MALLOC(_n_perio_face_couples, mesh->n_init_perio, cs_lnum_t);
1227 BFT_MALLOC(_perio_face_couples, mesh->n_init_perio, cs_gnum_t *);
1228
1229 for (i = 0; i < mesh->n_init_perio; i++) {
1230 _n_perio_face_couples[i] = 0;
1231 _perio_face_couples[i] = NULL;
1232 }
1233
1234 assert(halo != NULL);
1235 assert(halo->n_c_domains == 1);
1236
1237 /* Mark ghost cells with their periodicity number */
1238
1239 BFT_MALLOC(halo_perio_num, mesh->n_ghost_cells, int);
1240
1241 _get_halo_perio_num(mesh, halo_perio_num, NULL);
1242
1243 /* Exchange local cell numbers, so that cell-face reverse matching
1244 is made easier later */
1245
1246 BFT_MALLOC(loc_cell_num, mesh->n_cells_with_ghosts, cs_lnum_t);
1247 for (j = 0; j < mesh->n_cells; j++)
1248 loc_cell_num[j] = j+1;
1249 for (j = mesh->n_cells; j < mesh->n_cells_with_ghosts; j++)
1250 loc_cell_num[j] = 0;
1251
1252 cs_halo_sync_num(halo, CS_HALO_STANDARD, loc_cell_num);
1253
1254 /* Now compute number of periodic couples */
1255
1256 for (j = 0; j < mesh->n_i_faces; j++) {
1257 const cs_lnum_t h_id0 = mesh->i_face_cells[j][0] - mesh->n_cells;
1258 const cs_lnum_t h_id1 = mesh->i_face_cells[j][1] - mesh->n_cells;
1259 if (h_id0 >= 0) {
1260 if (halo_perio_num[h_id0] > 0)
1261 _n_perio_face_couples[halo_perio_num[h_id0] - 1] += 1;
1262 }
1263 else if (h_id1 >= 0) {
1264 if (halo_perio_num[h_id1] > 0)
1265 _n_perio_face_couples[halo_perio_num[h_id1] - 1] += 1;
1266 }
1267 }
1268
1269 /* Allocate periodic couple lists */
1270
1271 for (i = 0; i < mesh->n_init_perio; i++) {
1272 BFT_MALLOC(_perio_face_couples[i],
1273 _n_perio_face_couples[i]*2,
1274 cs_gnum_t);
1275 _n_perio_face_couples[i] = 0;
1276 }
1277
1278 /* Build cell -> face indexed search array containing only faces
1279 adjacent to the halo. */
1280
1281 _build_cell_face_perio_match(mesh,
1282 halo_perio_num,
1283 &cell_face_idx,
1284 &cell_face);
1285
1286 /* Now match faces; faces already matched are marked so as not to be
1287 re-matched, in case faces have been split and multiple faces thus share
1288 the same cells (in which case we assume that subfaces are locally ordered
1289 identically). */
1290
1291 for (j = 0; j < mesh->n_i_faces; j++) {
1292
1293 const cs_lnum_t c_id0 = mesh->i_face_cells[j][0];
1294 const cs_lnum_t c_id1 = mesh->i_face_cells[j][1];
1295
1296 dest_c_id = -1;
1297 if (c_id0 >= mesh->n_cells) {
1298 perio_num = halo_perio_num[c_id0 - mesh->n_cells];
1299 if (perio_num > 0) {
1300 dest_c_id = c_id1;
1301 src_c_id = c_id0;
1302 }
1303 }
1304 else if (c_id1 >= mesh->n_cells) {
1305 perio_num = halo_perio_num[c_id1 - mesh->n_cells];
1306 if (perio_num > 0) {
1307 dest_c_id = c_id0;
1308 src_c_id = c_id1;
1309 }
1310 }
1311 if (dest_c_id != -1) {
1312
1313 cs_lnum_t dest_c_num = loc_cell_num[dest_c_id];
1314 cs_lnum_t src_c_num = loc_cell_num[src_c_id];
1315
1316 for (k = cell_face_idx[src_c_num-1]; k < cell_face_idx[src_c_num]; k++) {
1317
1318 l = cell_face[k]; /* Face _id */
1319
1320 if (l > -1) {
1321 cs_lnum_t loc_c_num, dist_c_id;
1322 if (mesh->i_face_cells[l][0] < mesh->i_face_cells[l][1]) {
1323 loc_c_num = mesh->i_face_cells[l][0] + 1;
1324 dist_c_id = mesh->i_face_cells[l][1];
1325 }
1326 else {
1327 loc_c_num = mesh->i_face_cells[l][1] + 1;
1328 dist_c_id = mesh->i_face_cells[l][0];
1329 }
1330 assert(loc_c_num <= mesh->n_cells && dist_c_id >= mesh->n_cells);
1331
1332 /* If we have a match, update periodic couples */
1333
1334 if ( src_c_num == loc_c_num
1335 && dest_c_num == loc_cell_num[dist_c_id]
1336 && halo_perio_num[dist_c_id-mesh->n_cells] == -perio_num) {
1337 int perio_id = perio_num - 1;
1338 cs_lnum_t couple_id = _n_perio_face_couples[perio_id];
1339 cell_face[k] = -1; /* Mark as used */
1340 _perio_face_couples[perio_id][couple_id*2] = l+1;
1341 _perio_face_couples[perio_id][couple_id*2 + 1] = j+1;
1342 _n_perio_face_couples[perio_id] += 1;
1343 }
1344 }
1345 }
1346 }
1347 }
1348
1349 if (mesh->global_i_face_num != NULL) {
1350 for (i = 0; i < mesh->n_init_perio; i++) {
1351 for (j = 0; j < _n_perio_face_couples[i]*2; j++)
1352 _perio_face_couples[i][j]
1353 = mesh->global_i_face_num[_perio_face_couples[i][j] - 1];
1354 }
1355 }
1356
1357 BFT_FREE(halo_perio_num);
1358
1359 BFT_FREE(cell_face_idx);
1360 BFT_FREE(cell_face);
1361 BFT_FREE(loc_cell_num);
1362
1363 /* Set return values */
1364
1365 *n_perio_face_couples = _n_perio_face_couples;
1366 *perio_face_couples = _perio_face_couples;
1367 }
1368
1369 /*----------------------------------------------------------------------------
1370 * Check for free (unreferenced) vertices from a mesh.
1371 *
1372 * parameters:
1373 * mesh <-> pointer to mesh structure
1374 *
1375 * returns:
1376 * number of free vertices
1377 *----------------------------------------------------------------------------*/
1378
1379 static cs_gnum_t
_check_free_vertices(cs_mesh_t * mesh)1380 _check_free_vertices(cs_mesh_t *mesh)
1381 {
1382 char *ref = NULL;
1383
1384 cs_gnum_t retval = 0;
1385
1386 /* Mark vertices */
1387
1388 BFT_MALLOC(ref, mesh->n_vertices, char);
1389
1390 for (cs_lnum_t i = 0; i < mesh->n_vertices; i++)
1391 ref[i] = 0;
1392
1393 for (cs_lnum_t i = 0; i < mesh->i_face_vtx_connect_size; i++)
1394 ref[mesh->i_face_vtx_lst[i]] = 1;
1395
1396 for (cs_lnum_t i = 0; i < mesh->b_face_vtx_connect_size; i++)
1397 ref[mesh->b_face_vtx_lst[i]] = 1;
1398
1399 /* Transform marker to mapping */
1400
1401 for (cs_lnum_t i = 0; i < mesh->n_vertices; i++) {
1402 if (ref[i] == 0)
1403 retval += 1;
1404 }
1405
1406 #if defined(HAVE_MPI)
1407 if (cs_glob_n_ranks > 1) {
1408 cs_gnum_t n_g_count = retval;
1409 MPI_Allreduce(&n_g_count, &retval, 1, CS_MPI_GNUM, MPI_SUM,
1410 cs_glob_mpi_comm);
1411 }
1412 #endif
1413
1414 BFT_FREE(ref);
1415
1416 return retval;
1417 }
1418
1419 /*----------------------------------------------------------------------------
1420 * Discard free (unreferenced) vertices from a mesh.
1421 *
1422 * parameters:
1423 * mesh <-> pointer to mesh structure
1424 *----------------------------------------------------------------------------*/
1425
1426 static void
_discard_free_vertices(cs_mesh_t * mesh)1427 _discard_free_vertices(cs_mesh_t *mesh)
1428 {
1429 cs_lnum_t i;
1430 cs_lnum_t n_vertices = 0;
1431 cs_lnum_t *new_vertex_id = NULL;
1432
1433 /* Mark vertices */
1434
1435 BFT_MALLOC(new_vertex_id, mesh->n_vertices, cs_lnum_t);
1436
1437 for (i = 0; i < mesh->n_vertices; i++)
1438 new_vertex_id[i] = -1;
1439
1440 for (i = 0; i < mesh->i_face_vtx_connect_size; i++)
1441 new_vertex_id[mesh->i_face_vtx_lst[i]] = 0;
1442
1443 for (i = 0; i < mesh->b_face_vtx_connect_size; i++)
1444 new_vertex_id[mesh->b_face_vtx_lst[i]] = 0;
1445
1446 /* Transform marker to mapping */
1447
1448 for (i = 0; i < mesh->n_vertices; i++) {
1449 if (new_vertex_id[i] != -1)
1450 new_vertex_id[i] = n_vertices++;
1451 }
1452
1453 /* Update local mesh structure if necessary */
1454
1455 if (n_vertices < mesh->n_vertices) {
1456
1457 /* Update face connectivity */
1458
1459 for (i = 0; i < mesh->i_face_vtx_connect_size; i++)
1460 mesh->i_face_vtx_lst[i] = new_vertex_id[mesh->i_face_vtx_lst[i]];
1461
1462 for (i = 0; i < mesh->b_face_vtx_connect_size; i++)
1463 mesh->b_face_vtx_lst[i] = new_vertex_id[mesh->b_face_vtx_lst[i]];
1464
1465 /* Update vertex connectivity and global numbering */
1466
1467 for (i = 0; i < mesh->n_vertices; i++) {
1468 int l;
1469 const cs_lnum_t j = new_vertex_id[i];
1470 if (j != -1) {
1471 for (l = 0; l < 3; l++)
1472 mesh->vtx_coord[j*3 + l] = mesh->vtx_coord[i*3 + l];
1473 if (mesh->global_vtx_num != NULL)
1474 mesh->global_vtx_num[j] = mesh->global_vtx_num[i];
1475 }
1476 }
1477
1478 /* Update extended connectivity */
1479
1480 if (mesh->gcell_vtx_lst != NULL) {
1481 const cs_lnum_t n = mesh->gcell_vtx_idx[mesh->n_ghost_cells];
1482 for (i = 0; i < n; i++)
1483 mesh->gcell_vtx_lst[i] = new_vertex_id[mesh->gcell_vtx_lst[i]];
1484 }
1485
1486 /* Update mesh structure */
1487
1488 mesh->n_vertices = n_vertices;
1489
1490 BFT_REALLOC(mesh->vtx_coord, n_vertices*3, cs_real_t);
1491
1492 if (mesh->global_vtx_num != NULL)
1493 BFT_REALLOC(mesh->global_vtx_num, n_vertices, cs_gnum_t);
1494 }
1495
1496 if (mesh->vtx_interfaces != NULL)
1497 cs_interface_set_renumber(mesh->vtx_interfaces, new_vertex_id);
1498
1499 BFT_FREE(new_vertex_id);
1500
1501 /* Build an I/O numbering to compact the global numbering */
1502
1503 if (cs_glob_n_ranks > 1) {
1504
1505 fvm_io_num_t *tmp_num = fvm_io_num_create(NULL,
1506 mesh->global_vtx_num,
1507 mesh->n_vertices,
1508 0);
1509
1510 if (mesh->n_vertices > 0)
1511 memcpy(mesh->global_vtx_num,
1512 fvm_io_num_get_global_num(tmp_num),
1513 mesh->n_vertices*sizeof(cs_gnum_t));
1514
1515 mesh->n_g_vertices = fvm_io_num_get_global_count(tmp_num);
1516
1517 assert(fvm_io_num_get_local_count(tmp_num) == (cs_lnum_t)mesh->n_vertices);
1518
1519 tmp_num = fvm_io_num_destroy(tmp_num);
1520
1521 }
1522 else { /* if (cs_glob_ranks == 1) */
1523
1524 assert(mesh->global_vtx_num == NULL);
1525 mesh->n_g_vertices = mesh->n_vertices;
1526
1527 }
1528 }
1529
1530 /*----------------------------------------------------------------------------
1531 * Remove selector and associated structures.
1532 *
1533 * mesh <-> pointer to a mesh structure
1534 *----------------------------------------------------------------------------*/
1535
1536 static void
_free_selectors(cs_mesh_t * mesh)1537 _free_selectors(cs_mesh_t *mesh)
1538 {
1539 if (mesh->select_cells != NULL)
1540 mesh->select_cells = fvm_selector_destroy(mesh->select_cells);
1541 if (mesh->select_i_faces != NULL)
1542 mesh->select_i_faces = fvm_selector_destroy(mesh->select_i_faces);
1543 if (mesh->select_b_faces != NULL)
1544 mesh->select_b_faces = fvm_selector_destroy(mesh->select_b_faces);
1545
1546 /* Destroy group class set after selectors, who reference it */
1547
1548 if (mesh->class_defs != NULL)
1549 mesh->class_defs
1550 = fvm_group_class_set_destroy(mesh->class_defs);
1551 }
1552
1553 /*----------------------------------------------------------------------------
1554 * Prepare local processor count for mesh stats.
1555 *
1556 * parameters:
1557 * mesh <-- pointer to mesh structure
1558 * n_group_elts --> number of cells, interior faces, boundary faces,
1559 * and isolated faces belonging to each group
1560 * (size: mesh->n_groups * 4, interleaved)
1561 * n_perio_faces --> number of faces belonging to each group
1562 * (size: mesh->n_init_perio)
1563 *----------------------------------------------------------------------------*/
1564
1565 static void
_prepare_mesh_group_stats(const cs_mesh_t * mesh,cs_gnum_t n_group_elts[],cs_gnum_t n_perio_faces[])1566 _prepare_mesh_group_stats(const cs_mesh_t *mesh,
1567 cs_gnum_t n_group_elts[],
1568 cs_gnum_t n_perio_faces[])
1569 {
1570 cs_lnum_t i, j;
1571
1572 int *i_face_flag = NULL;
1573 cs_lnum_t *f_count = NULL;
1574
1575 /* In case of parallel and/or periodic faces, build a flag so as to count
1576 purely parallel faces only once, but periodic faces twice.
1577 The flag is defined to the periodicity number for periodic faces,
1578 0 for regular or purely parallel faces with an "outside" ghost cell,
1579 -1 for purely parallel faces with an "inside" ghost cell. */
1580
1581 if (mesh->halo != NULL) {
1582
1583 BFT_MALLOC(i_face_flag, mesh->n_i_faces, int);
1584
1585 if (mesh->n_init_perio > 0) {
1586 cs_mesh_get_face_perio_num(mesh, i_face_flag);
1587 for (i = 0; i < mesh->n_i_faces; i++)
1588 if (i_face_flag[i] < 0)
1589 i_face_flag[i] = - i_face_flag[i];
1590 }
1591 else {
1592 for (i = 0; i < mesh->n_i_faces; i++)
1593 i_face_flag[i] = 0;
1594 }
1595
1596 for (i = 0; i < mesh->n_i_faces; i++) {
1597 cs_lnum_t c_num_0 = mesh->i_face_cells[i][0] + 1;
1598 if (i_face_flag[i] == 0 && c_num_0 > mesh->n_cells)
1599 i_face_flag[i] = -1;
1600 }
1601
1602 }
1603
1604 /* Set counts to zero */
1605
1606 BFT_MALLOC(f_count, mesh->n_families*4, cs_lnum_t);
1607 for (i = 0; i < mesh->n_families*4; i++)
1608 f_count[i] = 0;
1609
1610 for (i = 0; i < mesh->n_init_perio; i++)
1611 n_perio_faces[i] = 0;
1612
1613 /* Cell counters */
1614
1615 for (i = 0; i < mesh->n_cells; i++)
1616 f_count[(mesh->cell_family[i] - 1) * 4] += 1;
1617
1618 /* Interior face counters */
1619
1620 if (i_face_flag != NULL) {
1621 for (i = 0; i < mesh->n_i_faces; i++) {
1622 const int flag = i_face_flag[i];
1623 if (flag > 0)
1624 n_perio_faces[flag - 1] += 1;
1625 if (flag > -1)
1626 f_count[(mesh->i_face_family[i] - 1) * 4 + 1] += 1;
1627 }
1628 BFT_FREE(i_face_flag);
1629 }
1630 else {
1631 for (i = 0; i < mesh->n_i_faces; i++)
1632 f_count[(mesh->i_face_family[i] - 1) * 4 + 1] += 1;
1633 }
1634
1635 /* Boundary face counters */
1636
1637 if (mesh->n_b_faces > 0) {
1638 for (i = 0; i < mesh->n_b_faces; i++) {
1639 int col = (mesh->b_face_cells[i] > -1) ? 2 : 3;
1640 f_count[(mesh->b_face_family[i] - 1) * 4 + col] += 1;
1641 }
1642 }
1643
1644 /* Now transfer count from group classes to groups */
1645
1646 for (i = 0; i < mesh->n_groups*4; i++)
1647 n_group_elts[i] = 0;
1648
1649 for (i = 0; i < mesh->n_families; i++) {
1650 for (j = 0; j < mesh->n_max_family_items; j++) {
1651 int group_id = -1 - mesh->family_item[mesh->n_families*j + i];
1652 if (group_id >= 0) {
1653 int k;
1654 for (k = 0; k < 4; k++)
1655 n_group_elts[group_id*4 + k] += f_count[i*4 + k];
1656 }
1657 }
1658 }
1659
1660 BFT_FREE(f_count);
1661 }
1662
1663 /*----------------------------------------------------------------------------
1664 * Print mesh group statistics.
1665 *
1666 * parameters:
1667 * mesh <-- pointer to mesh structure
1668 *----------------------------------------------------------------------------*/
1669
1670 static void
_print_mesh_group_stats(const cs_mesh_t * mesh)1671 _print_mesh_group_stats(const cs_mesh_t *mesh)
1672 {
1673 int i;
1674 size_t count_size = mesh->n_groups*4 + mesh->n_init_perio;
1675 cs_gnum_t *count = NULL, *n_elt_groups = NULL, *n_perio_faces = NULL;
1676
1677 if (count_size == 0)
1678 return;
1679
1680 BFT_MALLOC(count, count_size, cs_gnum_t);
1681
1682 n_elt_groups = count;
1683 n_perio_faces = count + mesh->n_groups*4;
1684
1685 _prepare_mesh_group_stats(mesh, n_elt_groups, n_perio_faces);
1686
1687 #if defined(HAVE_MPI)
1688 if (cs_glob_n_ranks > 1) {
1689 cs_gnum_t *_count = NULL;
1690 BFT_MALLOC(_count, count_size, cs_gnum_t);
1691 MPI_Allreduce(count, _count, count_size, CS_MPI_GNUM, MPI_SUM,
1692 cs_glob_mpi_comm);
1693 memcpy(count, _count, sizeof(cs_gnum_t)*count_size);
1694 BFT_FREE(_count);
1695 }
1696 #endif
1697
1698 /* Print perodic faces information */
1699
1700 if (mesh->n_init_perio > 0) {
1701
1702
1703 bft_printf(_("\n Periodic faces (which are also interior faces):\n"));
1704
1705 for (i = 0; i < mesh->n_init_perio; i++)
1706 bft_printf(_(" Periodicity %2d: %llu face couples\n"),
1707 i+1, (unsigned long long)(n_perio_faces[i]/2));
1708
1709 }
1710
1711 /* Print group information */
1712
1713 if (mesh->n_groups > 0) {
1714
1715 bft_printf(_("\n Groups:\n"));
1716
1717 for (i = 0; i < mesh->n_groups; i++) {
1718 bft_printf(" \"%s\"\n", mesh->group + mesh->group_idx[i]);
1719 if (n_elt_groups[i*4] > 0)
1720 bft_printf(_(" cells: %12llu\n"),
1721 (unsigned long long)(n_elt_groups[i*4]));
1722 if (n_elt_groups[i*4 + 1] > 0)
1723 bft_printf(_(" interior faces: %12llu\n"),
1724 (unsigned long long)(n_elt_groups[i*4 + 1]));
1725 if (n_elt_groups[i*4 + 2] > 0)
1726 bft_printf(_(" boundary faces: %12llu\n"),
1727 (unsigned long long)(n_elt_groups[i*4 + 2]));
1728 if (n_elt_groups[i*4 + 3] > 0)
1729 bft_printf(_(" isolated faces: %12llu\n"),
1730 (unsigned long long)(n_elt_groups[i*4 + 2]));
1731 }
1732
1733 }
1734
1735 BFT_FREE(count);
1736
1737 if (mesh->n_init_perio > 0 || mesh->n_groups > 0)
1738 bft_printf("\n");
1739 }
1740
1741 /*----------------------------------------------------------------------------
1742 * Write a summary about mesh features in log file
1743 *
1744 * parameters:
1745 * mesh <-- pointer to cs_mesh_t structure
1746 *----------------------------------------------------------------------------*/
1747
1748 static void
_print_mesh_info(cs_mesh_t * mesh)1749 _print_mesh_info(cs_mesh_t *mesh)
1750 {
1751 cs_halo_t *halo = mesh->halo;
1752
1753 cs_lnum_t *rank_buffer = NULL;
1754
1755 /* Summary of cell and ghost cell distribution */
1756
1757 if (mesh->n_domains > 1) {
1758
1759 BFT_MALLOC(rank_buffer, mesh->n_domains, cs_lnum_t);
1760
1761 #if defined(HAVE_MPI)
1762 MPI_Allgather(&(mesh->n_cells), 1, CS_MPI_LNUM,
1763 rank_buffer , 1, CS_MPI_LNUM, cs_glob_mpi_comm);
1764 #endif
1765
1766 bft_printf(_("\n Histogram of the number of cells per rank:\n\n"));
1767
1768 _display_histograms(mesh->n_domains, rank_buffer);
1769
1770 } /* End if n_domains > 1 */
1771
1772 else
1773 bft_printf
1774 (_(" Number of cells: %ld\n"),
1775 (long)mesh->n_cells);
1776
1777 bft_printf("\n ----------------------------------------------------------\n");
1778
1779 if (mesh->n_domains > 1) {
1780
1781 #if defined(HAVE_MPI)
1782 MPI_Allgather(&(mesh->n_cells_with_ghosts), 1, CS_MPI_LNUM,
1783 rank_buffer, 1, CS_MPI_LNUM, cs_glob_mpi_comm);
1784 #endif
1785
1786 bft_printf(_("\n Histogram of the number of standard + halo cells "
1787 "per rank:\n\n"));
1788
1789 _display_histograms(mesh->n_domains, rank_buffer);
1790
1791 } /* End if n_domains > 1 */
1792
1793 else
1794 bft_printf
1795 (_(" Number of cells + halo cells: %ld\n\n"),
1796 (long)mesh->n_cells_with_ghosts);
1797
1798
1799 bft_printf("\n ----------------------------------------------------------\n");
1800
1801 if (halo != NULL) {
1802
1803 cs_lnum_t n_std_ghost_cells = halo->n_elts[CS_HALO_STANDARD];
1804
1805 if (mesh->n_domains > 1) {
1806
1807 cs_lnum_t n_gcells = mesh->n_ghost_cells;
1808
1809 #if defined(HAVE_MPI)
1810 MPI_Allgather(&n_gcells , 1, CS_MPI_LNUM,
1811 rank_buffer, 1, CS_MPI_LNUM, cs_glob_mpi_comm);
1812 #endif
1813
1814 bft_printf
1815 (_("\n Histogram of the number of ghost cells per rank:\n\n"));
1816
1817 _display_histograms(mesh->n_domains, rank_buffer);
1818
1819 if (mesh->halo_type == CS_HALO_EXTENDED) {
1820
1821 #if defined(HAVE_MPI)
1822 MPI_Allgather(&n_std_ghost_cells, 1, CS_MPI_LNUM,
1823 rank_buffer, 1, CS_MPI_LNUM, cs_glob_mpi_comm);
1824 #endif
1825
1826 bft_printf
1827 (_("\n"
1828 " Histogram of the number of ghost cells\n"
1829 " in the standard neighborhood per rank:\n\n"));
1830
1831 _display_histograms(mesh->n_domains, rank_buffer);
1832 }
1833
1834 } /* If n_ranks > 1 */
1835
1836 else {
1837 bft_printf(_("\n Number of ghost cells: %10ld\n"),
1838 (long)mesh->n_ghost_cells);
1839 if (mesh->halo_type == CS_HALO_EXTENDED)
1840 bft_printf(_(" in the standard neighborhood: %10ld\n"),
1841 (long)n_std_ghost_cells);
1842 }
1843
1844 bft_printf("\n"
1845 " ----------------------------------------------------------\n");
1846
1847 } /* End if halo != NULL */
1848
1849 /* Summary of faces distribution */
1850
1851 if (mesh->n_domains > 1) {
1852
1853 #if defined(HAVE_MPI)
1854 MPI_Allgather(&(mesh->n_i_faces), 1, CS_MPI_LNUM,
1855 rank_buffer, 1, CS_MPI_LNUM, cs_glob_mpi_comm);
1856 #endif
1857
1858 bft_printf
1859 (_("\n Histogram of the number of interior faces per rank:\n\n"));
1860
1861 _display_histograms(mesh->n_domains, rank_buffer);
1862
1863 } /* End if n_domains > 1 */
1864
1865 else
1866 bft_printf
1867 (_(" Number of interior faces: %ld\n"),
1868 (long)mesh->n_i_faces);
1869
1870 bft_printf("\n ----------------------------------------------------------\n");
1871
1872 if (mesh->n_domains > 1) {
1873
1874 #if defined(HAVE_MPI)
1875 MPI_Allgather(&(mesh->n_b_faces), 1, CS_MPI_LNUM,
1876 rank_buffer, 1, CS_MPI_LNUM, cs_glob_mpi_comm);
1877 #endif
1878
1879 bft_printf
1880 (_("\n Histogram of the number of boundary faces per rank:\n\n"));
1881
1882 _display_histograms(mesh->n_domains, rank_buffer);
1883
1884 } /* End if n_domains > 1 */
1885
1886 else
1887 bft_printf
1888 (_(" Number of boundary faces: %ld\n"),
1889 (long)mesh->n_b_faces);
1890
1891 bft_printf("\n ----------------------------------------------------------\n");
1892
1893 /* Add cell neighbor info */
1894
1895 _print_cell_neighbor_info(mesh);
1896
1897 #if defined(HAVE_MPI)
1898
1899 /* Summary of the number of neighbors */
1900
1901 if (mesh->n_domains > 1) {
1902
1903 cs_lnum_t n_c_domains = halo->n_c_domains;
1904
1905 MPI_Allgather(&n_c_domains, 1, CS_MPI_LNUM,
1906 rank_buffer , 1, CS_MPI_LNUM, cs_glob_mpi_comm);
1907
1908 bft_printf(_("\n Histogram of the number of neighboring domains "
1909 "per rank:\n\n"));
1910
1911 _display_histograms(mesh->n_domains, rank_buffer);
1912
1913 bft_printf("\n"
1914 " ----------------------------------------------------------\n");
1915
1916 } /* If n_domains > 1 */
1917
1918 #endif
1919
1920 /* Cleanup */
1921
1922 if (mesh->n_domains > 1)
1923 BFT_FREE(rank_buffer);
1924
1925 bft_printf_flush();
1926 }
1927
1928 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1929
1930 /*============================================================================
1931 * Public function definitions for Fortran API
1932 *============================================================================*/
1933
1934 /*----------------------------------------------------------------------------
1935 * Update a scalar array in case of parallelism and/or periodicity.
1936 *
1937 * Fortran interface:
1938 *
1939 * subroutine synsca(var)
1940 * *****************
1941 *
1942 * var : <-> : scalar array
1943 *----------------------------------------------------------------------------*/
1944
CS_PROCF(synsca,SYNSCA)1945 void CS_PROCF(synsca, SYNSCA)
1946 (
1947 cs_real_t var[]
1948 )
1949 {
1950 cs_mesh_sync_var_scal(var);
1951 }
1952
1953 /*----------------------------------------------------------------------------
1954 * Update a scalar array in case of parallelism and/or periodicity,
1955 * using an extended halo.
1956 *
1957 * Fortran interface:
1958 *
1959 * subroutine synsce(var)
1960 * *****************
1961 *
1962 * var : <-> : scalar array
1963 *----------------------------------------------------------------------------*/
1964
CS_PROCF(synsce,SYNSCE)1965 void CS_PROCF(synsce, SYNSCE)
1966 (
1967 cs_real_t var[]
1968 )
1969 {
1970 cs_mesh_sync_var_scal_ext(var);
1971 }
1972
1973 /*----------------------------------------------------------------------------
1974 * Update a vector array in case of parallelism and/or periodicity.
1975 *
1976 * Fortran interface:
1977 *
1978 * subroutine synvin(var)
1979 * *****************
1980 *
1981 * var : <-> : interleaved vector (of dimension 3)
1982 *----------------------------------------------------------------------------*/
1983
CS_PROCF(synvin,SYNVIN)1984 void CS_PROCF(synvin, SYNVIN)
1985 (
1986 cs_real_t var[]
1987 )
1988 {
1989 cs_mesh_sync_var_vect(var);
1990 }
1991
1992 /*----------------------------------------------------------------------------
1993 * Update a vector array in case of parallelism and/or periodicity,
1994 * using an extended halo.
1995 *
1996 * Fortran interface:
1997 *
1998 * subroutine synvin(var)
1999 * *****************
2000 *
2001 * var : <-> : interleaved vector (of dimension 3)
2002 *----------------------------------------------------------------------------*/
2003
CS_PROCF(synvie,SYNVIE)2004 void CS_PROCF(synvie, SYNVIE)
2005 (
2006 cs_real_t var[]
2007 )
2008 {
2009 cs_mesh_sync_var_vect_ext(var);
2010 }
2011
2012 /*----------------------------------------------------------------------------
2013 * Update a tensor array in case of parallelism and/or periodicity.
2014 *
2015 * Fortran interface:
2016 *
2017 * subroutine syntin(var)
2018 * *****************
2019 *
2020 * var : <-> : interleaved tensor (of dimension 3x3)
2021 *----------------------------------------------------------------------------*/
2022
CS_PROCF(syntin,SYNTIN)2023 void CS_PROCF(syntin, SYNTIN)
2024 (
2025 cs_real_t var[]
2026 )
2027 {
2028 cs_mesh_sync_var_tens(var);
2029 }
2030
2031 /*----------------------------------------------------------------------------
2032 * Update a symmetric tensor array in case of parallelism and/or periodicity.
2033 *
2034 * Fortran interface:
2035 *
2036 * subroutine syntis(var)
2037 * *****************
2038 *
2039 * var : <-> : interleaved symmetric tensor (of dimension 6)
2040 *----------------------------------------------------------------------------*/
2041
CS_PROCF(syntis,SYNTIS)2042 void CS_PROCF(syntis, SYNTIS)
2043 (
2044 cs_real_t var[]
2045 )
2046 {
2047 cs_mesh_sync_var_sym_tens((cs_real_6_t *)var);
2048 }
2049
2050 /*=============================================================================
2051 * Public function definitions
2052 *============================================================================*/
2053
2054 /*----------------------------------------------------------------------------
2055 * Create an empty mesh structure
2056 *
2057 * returns:
2058 * pointer to created mesh structure
2059 *----------------------------------------------------------------------------*/
2060
2061 cs_mesh_t *
cs_mesh_create(void)2062 cs_mesh_create(void)
2063 {
2064 cs_mesh_t *mesh = NULL;
2065
2066 BFT_MALLOC (mesh, 1, cs_mesh_t);
2067
2068 /* General features */
2069
2070 mesh->dim = 3;
2071 mesh->domain_num = cs_glob_rank_id + 1;
2072 mesh->n_domains = 0;
2073
2074 mesh->time_dep = CS_MESH_FIXED;
2075
2076 /* Local dimensions */
2077
2078 mesh->n_cells = 0;
2079 mesh->n_i_faces = 0;
2080 mesh->n_b_faces = 0;
2081 mesh->n_vertices = 0;
2082 mesh->i_face_vtx_connect_size = 0;
2083 mesh->b_face_vtx_connect_size = 0;
2084
2085 /* Global dimensions */
2086
2087 mesh->n_g_cells = 0;
2088 mesh->n_g_i_faces = 0;
2089 mesh->n_g_b_faces = 0;
2090 mesh->n_g_vertices = 0;
2091
2092 mesh->n_g_i_c_faces = 0;
2093
2094 /* Local structures */
2095
2096 mesh->vtx_coord = NULL;
2097 mesh->i_face_cells = NULL;
2098 mesh->b_face_cells = NULL;
2099 mesh->i_face_vtx_idx = NULL;
2100 mesh->b_face_vtx_idx = NULL;
2101 mesh->i_face_vtx_lst = NULL;
2102 mesh->b_face_vtx_lst = NULL;
2103
2104 /* Global numbering */
2105
2106 mesh->global_cell_num = NULL;
2107 mesh->global_i_face_num = NULL;
2108 mesh->global_b_face_num = NULL;
2109 mesh->global_vtx_num = NULL;
2110
2111 /* Periodic features */
2112
2113 mesh->n_init_perio = 0;
2114 mesh->n_transforms = 0;
2115 mesh->have_rotation_perio = 0;
2116 mesh->periodicity = NULL;
2117
2118 /* Halo features */
2119
2120 mesh->vtx_interfaces = NULL;
2121 mesh->halo_type = CS_HALO_N_TYPES;
2122 mesh->vtx_range_set = NULL;
2123 mesh->n_ghost_cells = 0;
2124 mesh->n_cells_with_ghosts = 0;
2125 mesh->halo = NULL;
2126
2127 /* Extended connectivity and neighborhood */
2128
2129 mesh->n_b_cells = 0;
2130 mesh->b_cells = NULL;
2131
2132 mesh->cell_cells_idx = NULL;
2133 mesh->cell_cells_lst = NULL;
2134 mesh->gcell_vtx_idx = NULL;
2135 mesh->gcell_vtx_lst = NULL;
2136
2137 /* Numbering info */
2138
2139 mesh->cell_numbering = NULL;
2140 mesh->i_face_numbering = NULL;
2141 mesh->b_face_numbering = NULL;
2142 mesh->vtx_numbering = NULL;
2143
2144 /* Group and family features */
2145
2146 mesh->n_groups = 0;
2147 mesh->group_idx = NULL;
2148 mesh->group = NULL;
2149
2150 mesh->n_max_family_items = 0;
2151 mesh->n_families = 0;
2152
2153 mesh->family_item = NULL;
2154 mesh->cell_family = NULL;
2155 mesh->i_face_family = NULL;
2156 mesh->b_face_family = NULL;
2157
2158 /* Refinement */
2159
2160 mesh->i_face_r_gen = NULL;
2161
2162 /* Selector features */
2163
2164 mesh->class_defs = NULL;
2165
2166 mesh->select_cells = NULL;
2167 mesh->select_i_faces = NULL;
2168 mesh->select_b_faces = NULL;
2169
2170 /* Status flags */
2171
2172 mesh->n_g_free_faces = 0;
2173
2174 mesh->n_g_b_faces_all = 0;
2175 mesh->n_b_faces_all = 0;
2176
2177 mesh->verbosity = 1;
2178 mesh->modified = 0;
2179 mesh->save_if_modified = 1;
2180
2181 return (mesh);
2182 }
2183
2184 /*----------------------------------------------------------------------------
2185 * Destroy a mesh structure.
2186 *
2187 * mesh <-> pointer to a mesh structure
2188 *
2189 * returns:
2190 * NULL pointer
2191 *----------------------------------------------------------------------------*/
2192
2193 cs_mesh_t *
cs_mesh_destroy(cs_mesh_t * mesh)2194 cs_mesh_destroy(cs_mesh_t *mesh)
2195 {
2196 BFT_FREE(mesh->vtx_coord);
2197 BFT_FREE(mesh->i_face_cells);
2198 BFT_FREE(mesh->b_face_cells);
2199 BFT_FREE(mesh->i_face_vtx_idx);
2200 BFT_FREE(mesh->b_face_vtx_idx);
2201 BFT_FREE(mesh->i_face_vtx_lst);
2202 BFT_FREE(mesh->b_face_vtx_lst);
2203
2204 BFT_FREE(mesh->global_cell_num);
2205 BFT_FREE(mesh->global_i_face_num);
2206 BFT_FREE(mesh->global_b_face_num);
2207 BFT_FREE(mesh->global_vtx_num);
2208
2209 BFT_FREE(mesh->group_idx);
2210 BFT_FREE(mesh->group);
2211
2212 BFT_FREE(mesh->family_item);
2213 BFT_FREE(mesh->cell_family);
2214 BFT_FREE(mesh->i_face_family);
2215 BFT_FREE(mesh->b_face_family);
2216
2217 BFT_FREE(mesh->i_face_r_gen);
2218
2219 /* Free periodic structures */
2220
2221 if (mesh->n_init_perio > 0)
2222 mesh->periodicity = fvm_periodicity_destroy(mesh->periodicity);
2223
2224 cs_mesh_free_rebuildable(mesh, true);
2225
2226 BFT_FREE(mesh);
2227
2228 return mesh;
2229 }
2230
2231 /*----------------------------------------------------------------------------
2232 * Update (compactify) an array of global numbers.
2233 *
2234 * parameters:
2235 * n_elts <-> number of local elements
2236 * elt_gnum <-> global element numbers
2237 *
2238 * return:
2239 * associated global number of elements
2240 *----------------------------------------------------------------------------*/
2241
2242 cs_gnum_t
cs_mesh_compact_gnum(cs_lnum_t n_elts,cs_gnum_t * elt_gnum)2243 cs_mesh_compact_gnum(cs_lnum_t n_elts,
2244 cs_gnum_t *elt_gnum)
2245 {
2246 cs_gnum_t n_g_elts = n_elts;
2247
2248 if (cs_glob_n_ranks > 1 || elt_gnum != NULL) {
2249
2250 fvm_io_num_t *tmp_num = fvm_io_num_create(NULL, elt_gnum, n_elts, 0);
2251
2252 if (n_elts > 0)
2253 memcpy(elt_gnum,
2254 fvm_io_num_get_global_num(tmp_num),
2255 n_elts*sizeof(cs_gnum_t));
2256
2257 n_g_elts = fvm_io_num_get_global_count(tmp_num);
2258
2259 assert(fvm_io_num_get_local_count(tmp_num) == n_elts);
2260
2261 tmp_num = fvm_io_num_destroy(tmp_num);
2262
2263 }
2264
2265 return n_g_elts;
2266 }
2267
2268 /*----------------------------------------------------------------------------
2269 * Remove arrays and structures that may be rebuilt.
2270 *
2271 * mesh <-> pointer to a mesh structure
2272 * free_halos <-- if true, free halos and parallel/periodic interface
2273 * structures
2274 *----------------------------------------------------------------------------*/
2275
2276 void
cs_mesh_free_rebuildable(cs_mesh_t * mesh,bool free_halos)2277 cs_mesh_free_rebuildable(cs_mesh_t *mesh,
2278 bool free_halos)
2279 {
2280 /* Free structures that may be rebuilt */
2281
2282 BFT_FREE(mesh->b_cells);
2283
2284 if (mesh->cell_cells_idx != NULL) {
2285 BFT_FREE(mesh->cell_cells_idx);
2286 BFT_FREE(mesh->cell_cells_lst);
2287 }
2288
2289 if (mesh->gcell_vtx_idx != NULL) {
2290 BFT_FREE(mesh->gcell_vtx_idx);
2291 BFT_FREE(mesh->gcell_vtx_lst);
2292 }
2293
2294 /* Free halo, interface and range set structures */
2295
2296 if (free_halos) {
2297
2298 if (mesh->vtx_interfaces != NULL)
2299 cs_interface_set_destroy(&(mesh->vtx_interfaces));
2300 if (mesh->halo != NULL)
2301 cs_halo_destroy(&(mesh->halo));
2302 if (mesh->vtx_range_set != NULL)
2303 cs_range_set_destroy(&(mesh->vtx_range_set));
2304
2305 }
2306
2307 /* Free numbering info */
2308
2309 if (mesh->cell_numbering != NULL)
2310 cs_numbering_destroy(&(mesh->cell_numbering));
2311 if (mesh->i_face_numbering != NULL)
2312 cs_numbering_destroy(&(mesh->i_face_numbering));
2313 if (mesh->b_face_numbering != NULL)
2314 cs_numbering_destroy(&(mesh->b_face_numbering));
2315 if (mesh->vtx_numbering != NULL)
2316 cs_numbering_destroy(&(mesh->vtx_numbering));
2317
2318 /* Free selection structures */
2319
2320 _free_selectors(mesh);
2321 }
2322
2323 /*----------------------------------------------------------------------------
2324 * Discard free (isolated) faces from a mesh.
2325 *
2326 * This should always be done before using the mesh for computation.
2327 *
2328 * parameters:
2329 * mesh <-> pointer to mesh structure
2330 *----------------------------------------------------------------------------*/
2331
2332 void
cs_mesh_discard_free_faces(cs_mesh_t * mesh)2333 cs_mesh_discard_free_faces(cs_mesh_t *mesh)
2334 {
2335 cs_lnum_t i;
2336 cs_gnum_t n_g_b_faces_old = 0, n_g_vertices_old = 0;
2337 cs_lnum_t j = 0, k = 0, l = 0;
2338
2339 if (mesh->n_g_free_faces == 0)
2340 return;
2341
2342 n_g_b_faces_old = mesh->n_g_b_faces;
2343 n_g_vertices_old = mesh->n_g_vertices;
2344
2345 for (i = 0; i < mesh->n_b_faces; i++) {
2346
2347 if (mesh->b_face_cells[i] > -1) {
2348
2349 mesh->b_face_cells[j] = mesh->b_face_cells[i];
2350 mesh->b_face_family[j] = mesh->b_face_family[i];
2351
2352 mesh->b_face_vtx_idx[j] = l;
2353 for (k = mesh->b_face_vtx_idx[i]; k < mesh->b_face_vtx_idx[i+1]; k++)
2354 mesh->b_face_vtx_lst[l++] = mesh->b_face_vtx_lst[k];
2355
2356 if (mesh->global_b_face_num != NULL)
2357 mesh->global_b_face_num[j] = mesh->global_b_face_num[i];
2358
2359 j += 1;
2360 }
2361
2362 }
2363
2364 mesh->b_face_vtx_idx[j] = l;
2365 mesh->b_face_vtx_connect_size = l;
2366
2367 /* Resize arrays if necessary */
2368
2369 if (j < i) {
2370
2371 BFT_REALLOC(mesh->b_face_cells, j, cs_lnum_t);
2372 BFT_REALLOC(mesh->b_face_family, j, int);
2373 BFT_REALLOC(mesh->b_face_vtx_idx, j+1, cs_lnum_t);
2374 BFT_REALLOC(mesh->b_face_vtx_lst, k, cs_lnum_t);
2375 if (mesh->global_b_face_num != NULL) {
2376 BFT_REALLOC(mesh->global_b_face_num, j, cs_gnum_t);
2377 }
2378
2379 mesh->n_b_faces = j;
2380 }
2381
2382 /* Build an I/O numbering on boundary faces to compact the global numbering */
2383
2384 mesh->n_g_b_faces = cs_mesh_compact_gnum(mesh->n_b_faces,
2385 mesh->global_b_face_num);
2386
2387 /* Now also clean-up unreferenced vertices */
2388
2389 _discard_free_vertices(mesh);
2390
2391 bft_printf(_("\n"
2392 " Removed %llu isolated faces\n"
2393 " Number of initial vertices: %llu\n"
2394 " Number of vertices: %llu\n\n"),
2395 (unsigned long long)(n_g_b_faces_old - mesh->n_g_b_faces),
2396 (unsigned long long)(n_g_vertices_old),
2397 (unsigned long long)(mesh->n_g_vertices));
2398
2399 mesh->n_g_free_faces = 0;
2400
2401 mesh->modified |= CS_MESH_MODIFIED;
2402 }
2403
2404 /*----------------------------------------------------------------------------
2405 * Discard free (isolated) vertices from a mesh.
2406 *
2407 * This is recommended before using the mesh for computation.
2408 *
2409 * parameters:
2410 * mesh <-> pointer to mesh structure
2411 *----------------------------------------------------------------------------*/
2412
2413 void
cs_mesh_discard_free_vertices(cs_mesh_t * mesh)2414 cs_mesh_discard_free_vertices(cs_mesh_t *mesh)
2415 {
2416 cs_gnum_t n_g_f_vertices = _check_free_vertices(mesh);
2417
2418 if (n_g_f_vertices > 0) {
2419 cs_gnum_t n_g_vertices_old = mesh->n_g_vertices;
2420 _discard_free_vertices(mesh);
2421 bft_printf(_("\n"
2422 " Removed isolated vertices\n"
2423 " Number of initial vertices: %llu\n"
2424 " Number of vertices: %llu\n\n"),
2425 (unsigned long long)(n_g_vertices_old),
2426 (unsigned long long)(mesh->n_g_vertices));
2427
2428 mesh->modified |= CS_MESH_MODIFIED;
2429 }
2430 }
2431
2432 /*----------------------------------------------------------------------------
2433 * Compute global face connectivity size.
2434 *
2435 * Faces on simple parallel boundaries are counted only once, but periodic
2436 * faces are counted twice.
2437 *
2438 * parameters:
2439 * mesh <-- pointer to a cs_mesh_t structure
2440 * g_i_face_vertices_size --> global interior face connectivity size, or NULL
2441 * g_b_face_vertices_size --> global boundary face connectivity size, or NULL
2442 *----------------------------------------------------------------------------*/
2443
2444 void
cs_mesh_g_face_vertices_sizes(const cs_mesh_t * mesh,cs_gnum_t * g_i_face_vertices_size,cs_gnum_t * g_b_face_vertices_size)2445 cs_mesh_g_face_vertices_sizes(const cs_mesh_t *mesh,
2446 cs_gnum_t *g_i_face_vertices_size,
2447 cs_gnum_t *g_b_face_vertices_size)
2448 {
2449 cs_gnum_t _g_face_vertices_size[2] = {0, 0};
2450
2451 #if defined(HAVE_MPI)
2452
2453 if (cs_glob_n_ranks > 1) {
2454
2455 cs_lnum_t i;
2456 cs_gnum_t _l_face_vertices_size[2] = {0, 0};
2457
2458 /* Without periodicity, count faces bordering halo once */
2459
2460 if (mesh->n_init_perio == 0) {
2461 for (i = 0; i < mesh->n_i_faces; i++) {
2462 if (mesh->i_face_cells[i][0] < mesh->n_cells)
2463 _l_face_vertices_size[0] += ( mesh->i_face_vtx_idx[i+1]
2464 - mesh->i_face_vtx_idx[i]);
2465 }
2466 }
2467
2468 /* With periodicity, count faces bordering halo once for purely parallel
2469 halo */
2470
2471 else {
2472
2473 cs_lnum_t rank_id, t_id, shift;
2474 cs_lnum_t start = 0, end = 0;
2475 int *perio_flag = NULL;
2476
2477 const cs_halo_t *halo = mesh->halo;
2478 const cs_lnum_t n_transforms = halo->n_transforms;
2479
2480 BFT_MALLOC(perio_flag, mesh->n_ghost_cells, int);
2481 for (i = 0; i < mesh->n_ghost_cells; i++)
2482 perio_flag[i] = 0;
2483
2484 for (t_id = 0; t_id < n_transforms; t_id++) {
2485 shift = 4 * halo->n_c_domains * t_id;
2486 for (rank_id = 0; rank_id < halo->n_c_domains; rank_id++) {
2487 start = halo->perio_lst[shift + 4*rank_id];
2488 end = start + halo->perio_lst[shift + 4*rank_id + 1];
2489 for (i = start; i < end; i++)
2490 perio_flag[i] = 1;
2491 }
2492 }
2493
2494 for (i = 0; i < mesh->n_i_faces; i++) {
2495 int count_face = 1;
2496 if (mesh->i_face_cells[i][0] >= mesh->n_cells) {
2497 if (perio_flag[mesh->i_face_cells[i][0] - mesh->n_cells] == 0)
2498 count_face = 0;
2499 }
2500 if (count_face)
2501 _l_face_vertices_size[0] += ( mesh->i_face_vtx_idx[i+1]
2502 - mesh->i_face_vtx_idx[i]);
2503 }
2504
2505 BFT_FREE(perio_flag);
2506 }
2507
2508 _l_face_vertices_size[1] = mesh->b_face_vtx_connect_size;
2509
2510 MPI_Allreduce(&_l_face_vertices_size, &_g_face_vertices_size, 2,
2511 CS_MPI_GNUM, MPI_SUM, cs_glob_mpi_comm);
2512 }
2513
2514 #endif /* defined(HAVE_MPI) */
2515
2516 if (cs_glob_n_ranks == 1) {
2517 _g_face_vertices_size[0] = mesh->i_face_vtx_connect_size;
2518 _g_face_vertices_size[1] = mesh->b_face_vtx_connect_size;
2519 }
2520
2521 if (g_i_face_vertices_size != NULL)
2522 *g_i_face_vertices_size = _g_face_vertices_size[0];
2523 if (g_b_face_vertices_size != NULL)
2524 *g_b_face_vertices_size = _g_face_vertices_size[1];
2525 }
2526
2527 /*----------------------------------------------------------------------------
2528 * Generate or update list of mesh boundary cells.
2529 *
2530 * parameters:
2531 * mesh <-> pointer to a cs_mesh_t structure
2532 *----------------------------------------------------------------------------*/
2533
2534 void
cs_mesh_update_b_cells(cs_mesh_t * mesh)2535 cs_mesh_update_b_cells(cs_mesh_t *mesh)
2536 {
2537 cs_lnum_t i;
2538
2539 cs_lnum_t n_b_cells = 0;
2540 bool *flag = NULL;
2541
2542 BFT_MALLOC(flag, mesh->n_cells, bool);
2543
2544 for (i = 0; i < mesh->n_cells; i++)
2545 flag[i] = false;
2546
2547 for (i = 0; i < mesh->n_b_faces; i++) {
2548 if (mesh->b_face_cells[i] > -1)
2549 flag[mesh->b_face_cells[i]] = true;
2550 }
2551
2552 for (i = 0, n_b_cells = 0; i < mesh->n_cells; i++) {
2553 if (flag[i] == true)
2554 n_b_cells++;
2555 }
2556
2557 mesh->n_b_cells = n_b_cells;
2558 BFT_REALLOC(mesh->b_cells, mesh->n_b_cells, cs_lnum_t);
2559
2560 for (i = 0, n_b_cells = 0; i < mesh->n_cells; i++) {
2561 if (flag[i] == true)
2562 mesh->b_cells[n_b_cells++] = i;
2563 }
2564
2565 BFT_FREE(flag);
2566 }
2567
2568 /*----------------------------------------------------------------------------
2569 * Compute or update mesh structure members that depend on other members,
2570 * but whose results may be reused, such as global number of elements
2571 * (cells, vertices, internal and border faces) and sync cell family.
2572 *
2573 * parameters:
2574 * mesh <-> pointer to a cs_mesh_t structure
2575 *----------------------------------------------------------------------------*/
2576
2577 void
cs_mesh_update_auxiliary(cs_mesh_t * mesh)2578 cs_mesh_update_auxiliary(cs_mesh_t *mesh)
2579 {
2580 cs_lnum_t i;
2581
2582 #if defined(HAVE_MPI)
2583
2584 if (cs_glob_n_ranks > 1) {
2585
2586 cs_gnum_t n_g_elts[4], max_elt_num[4];
2587
2588 if (mesh->verbosity > 0)
2589 bft_printf(_("\n Global definition of the number of elements "
2590 "(cells, vertices, faces...)\n"));
2591
2592 /* Global dimensions of the mesh */
2593
2594 max_elt_num[0] = mesh->n_cells;
2595 MPI_Allreduce(max_elt_num, n_g_elts, 1, CS_MPI_GNUM, MPI_SUM,
2596 cs_glob_mpi_comm);
2597
2598 max_elt_num[1] = 0;
2599 for (i = 0; i < mesh->n_i_faces; i++) {
2600 if (mesh->global_i_face_num[i] > max_elt_num[1])
2601 max_elt_num[1] = mesh->global_i_face_num[i];
2602 }
2603
2604 max_elt_num[2] = 0;
2605 for (i = 0; i < mesh->n_b_faces; i++) {
2606 if (mesh->global_b_face_num[i] > max_elt_num[2])
2607 max_elt_num[2] = mesh->global_b_face_num[i];
2608 }
2609
2610 max_elt_num[3] = 0;
2611 for (i = 0; i < mesh->n_vertices; i++) {
2612 if (mesh->global_vtx_num[i] > max_elt_num[3])
2613 max_elt_num[3] = mesh->global_vtx_num[i];
2614 }
2615
2616 MPI_Allreduce(max_elt_num + 1, n_g_elts + 1, 3, CS_MPI_GNUM, MPI_MAX,
2617 cs_glob_mpi_comm);
2618
2619 mesh->n_g_cells = n_g_elts[0];
2620 mesh->n_g_i_faces = n_g_elts[1];
2621 mesh->n_g_b_faces = n_g_elts[2];
2622 mesh->n_g_vertices = n_g_elts[3];
2623
2624 }
2625
2626 #endif
2627
2628 if (cs_glob_n_ranks == 1) {
2629 mesh->n_g_cells = mesh->n_cells;
2630 mesh->n_g_i_faces = mesh->n_i_faces;
2631 mesh->n_g_b_faces = mesh->n_b_faces;
2632 mesh->n_g_vertices = mesh->n_vertices;
2633 }
2634
2635 mesh->n_g_i_c_faces = mesh->n_g_i_faces;
2636
2637 if (mesh->n_init_perio > 0) {
2638
2639 const cs_lnum_t n_cells = mesh->n_cells;
2640 cs_gnum_t n_g_i_c_faces = 0;
2641 for (i = 0; i < mesh->n_i_faces; i++) {
2642 if (mesh->i_face_cells[i][0] < n_cells)
2643 n_g_i_c_faces++;
2644 }
2645
2646 if (cs_glob_n_ranks == 1)
2647 mesh->n_g_i_c_faces = n_g_i_c_faces;
2648 #if defined(HAVE_MPI)
2649 else if (cs_glob_n_ranks > 1)
2650 MPI_Allreduce(&n_g_i_c_faces, &(mesh->n_g_i_c_faces), 1,
2651 CS_MPI_GNUM, MPI_SUM, cs_glob_mpi_comm);
2652 #endif
2653 }
2654
2655 /* Sync cell family array (also in case of periodicity) */
2656
2657 _sync_cell_fam(mesh);
2658
2659 /* Number of boundary cells */
2660
2661 cs_mesh_update_b_cells(mesh);
2662 }
2663
2664 /*----------------------------------------------------------------------------*/
2665 /*!
2666 * Creation and initialization of halo structures.
2667 *
2668 * Treatment of parallel and/or periodic halos for standard and extended
2669 * ghost cells according to halo type requested by global options.
2670 *
2671 * parameters:
2672 * \param[in, out] mesh pointer to mesh structure
2673 * \param[in, out] mb pointer to mesh builder
2674 * (for periodicity)
2675 * \param[in] halo_type type of halo (standard or extended)
2676 * \param[in] verbosity verbosity
2677 * \param[in] rebuild_vtx_interface also rebuild vertex interfaces ?
2678 *----------------------------------------------------------------------------*/
2679
2680 void
cs_mesh_init_halo(cs_mesh_t * mesh,cs_mesh_builder_t * mb,cs_halo_type_t halo_type,int verbosity,bool rebuild_vtx_interface)2681 cs_mesh_init_halo(cs_mesh_t *mesh,
2682 cs_mesh_builder_t *mb,
2683 cs_halo_type_t halo_type,
2684 int verbosity,
2685 bool rebuild_vtx_interface)
2686 {
2687 cs_lnum_t i;
2688 double t1, t2;
2689 double halo_time = 0, interface_time = 0, ext_neighborhood_time = 0;
2690
2691 int *perio_num = NULL;
2692 cs_lnum_t *gcell_vtx_idx = NULL, *gcell_vtx_lst = NULL;
2693 cs_lnum_t *n_periodic_couples = NULL;
2694 cs_gnum_t *g_i_face_num = NULL;
2695 cs_gnum_t *g_vertex_num = NULL;
2696 cs_gnum_t **periodic_couples = NULL;
2697
2698 cs_interface_set_t *face_interfaces = NULL;
2699
2700 const cs_lnum_t n_i_faces = mesh->n_i_faces;
2701 const cs_lnum_t n_vertices = mesh->n_vertices;
2702
2703 int s_verbosity = mesh->verbosity;
2704 s_verbosity = verbosity;
2705
2706 /* Build halo */
2707
2708 #if 0
2709 /* TODO: activating this option needs a better upstream setting
2710 of extended neighborhood type */
2711 if (cs_ext_neighborhood_get_type() > CS_EXT_NEIGHBORHOOD_NONE)
2712 halo_type = CS_HALO_EXTENDED;
2713 #endif
2714
2715 mesh->halo_type = halo_type;
2716
2717 mesh->n_ghost_cells = 0;
2718 mesh->n_cells_with_ghosts = mesh->n_cells + mesh->n_ghost_cells;
2719
2720 if (mesh->n_domains > 1 || mesh->n_init_perio > 0) {
2721
2722 t1 = cs_timer_wtime();
2723
2724 if (verbosity > 0)
2725 bft_printf("\n"
2726 " ----------------------------------------------------------"
2727 "\n");
2728
2729 if (mesh->periodicity != NULL) {
2730 if (fvm_periodicity_get_n_levels(mesh->periodicity) == 1) {
2731 if (verbosity > 0)
2732 bft_printf(_(" Composing periodicities\n"));
2733 fvm_periodicity_combine(mesh->periodicity, 0);
2734 }
2735 }
2736
2737 if (verbosity > 0) {
2738 if (halo_type == CS_HALO_EXTENDED)
2739 bft_printf(_("\n Halo construction with extended neighborhood\n"
2740 " ============================================\n\n"));
2741
2742 else
2743 bft_printf(_("\n Halo construction with standard neighborhood\n"
2744 " ============================================\n\n"));
2745 }
2746
2747 /* Build periodic numbering */
2748
2749 if (mesh->n_init_perio > 0) {
2750 BFT_MALLOC(perio_num, mesh->n_init_perio, int);
2751 for (i = 0; i < mesh->n_init_perio; i++)
2752 perio_num[i] = i+1;
2753 }
2754
2755 if (verbosity > 0)
2756 bft_printf(_(" Face interfaces creation\n"));
2757
2758 /* Build purely parallel cs_interface_set structure */
2759
2760 g_i_face_num = mesh->global_i_face_num;
2761 g_vertex_num = mesh->global_vtx_num;
2762
2763 if (g_i_face_num == NULL) {
2764 BFT_MALLOC(g_i_face_num, n_i_faces, cs_gnum_t);
2765 for (i = 0; i < n_i_faces; i++)
2766 g_i_face_num[i] = (cs_gnum_t)i+1;
2767 }
2768
2769 if (mb != NULL) {
2770 face_interfaces = cs_interface_set_create(n_i_faces,
2771 NULL,
2772 g_i_face_num,
2773 mesh->periodicity,
2774 mesh->n_init_perio,
2775 perio_num,
2776 mb->n_per_face_couples,
2777 (const cs_gnum_t *const *)(mb->per_face_couples));
2778 for (i = 0; i < mb->n_perio; i++)
2779 BFT_FREE(mb->per_face_couples[i]);
2780 BFT_FREE(mb->per_face_couples);
2781 BFT_FREE(mb->n_per_face_couples);
2782 mb->n_perio = 0;
2783 }
2784 else
2785 face_interfaces = cs_interface_set_create(n_i_faces,
2786 NULL,
2787 g_i_face_num,
2788 mesh->periodicity,
2789 0,
2790 NULL,
2791 NULL,
2792 NULL);
2793
2794 if (mesh->global_i_face_num != g_i_face_num)
2795 BFT_FREE(g_i_face_num);
2796
2797 if (g_vertex_num == NULL) {
2798 BFT_MALLOC(g_vertex_num, n_vertices, cs_gnum_t);
2799 for (i = 0; i < n_vertices; i++)
2800 g_vertex_num[i] = (cs_gnum_t)i+1;
2801 }
2802
2803 if (mesh->n_init_perio > 0) {
2804
2805 mesh->n_transforms = fvm_periodicity_get_n_transforms(mesh->periodicity);
2806
2807 bft_printf(_(" Definition of periodic vertices\n"));
2808
2809 _define_perio_vtx_couples(mesh,
2810 face_interfaces,
2811 &n_periodic_couples,
2812 &periodic_couples);
2813
2814 #if 0 /* For debugging purposes */
2815 for (i = 0; i < mesh->n_init_perio; i++) {
2816 cs_lnum_t j;
2817 bft_printf("\n\n Periodicity number: %d\n", perio_num[i]);
2818 bft_printf(" Number of couples : %d\n", n_periodic_couples[i]);
2819 for (j = 0; j < n_periodic_couples[i]; j++)
2820 bft_printf("%12d --> %12d\n",
2821 periodic_couples[i][2*j], periodic_couples[i][2*j + 1]);
2822 }
2823 fvm_periodicity_dump(mesh->periodicity);
2824 #endif
2825
2826 }
2827
2828 if (rebuild_vtx_interface) {
2829
2830 if (verbosity > 0)
2831 bft_printf(_(" Vertex interfaces creation\n"));
2832
2833 if (mesh->vtx_interfaces != NULL)
2834 cs_interface_set_destroy(&mesh->vtx_interfaces);
2835
2836 mesh->vtx_interfaces = cs_interface_set_create(n_vertices,
2837 NULL,
2838 g_vertex_num,
2839 mesh->periodicity,
2840 mesh->n_init_perio,
2841 perio_num,
2842 n_periodic_couples,
2843 (const cs_gnum_t *const *)periodic_couples);
2844
2845 }
2846
2847 if (mesh->global_vtx_num != g_vertex_num)
2848 BFT_FREE(g_vertex_num);
2849
2850 BFT_FREE(perio_num);
2851 BFT_FREE(n_periodic_couples);
2852
2853 for (i = 0; i < mesh->n_init_perio; i++)
2854 BFT_FREE(periodic_couples[i]);
2855 BFT_FREE(periodic_couples);
2856
2857 #if 0 /* For debugging purposes */
2858 bft_printf("Dump final vertices interface:\n");
2859 cs_interface_set_add_match_ids(mesh->vtx_interfaces);
2860 cs_interface_set_dump(mesh->vtx_interfaces);
2861 cs_interface_set_free_match_ids(mesh->vtx_interfaces);
2862 #endif
2863
2864 t2 = cs_timer_wtime();
2865 interface_time = t2-t1;
2866
2867 t1 = cs_timer_wtime();
2868
2869 /* Creation of the cs_halo_t structure. */
2870
2871 if (verbosity > 0) {
2872 bft_printf(_(" Halo creation\n"));
2873 bft_printf_flush();
2874 }
2875
2876 mesh->halo = cs_halo_create(mesh->vtx_interfaces);
2877
2878 if (verbosity > 0) {
2879 bft_printf(_(" Halo definition\n"));
2880 bft_printf_flush();
2881 }
2882
2883 cs_mesh_halo_define(mesh,
2884 face_interfaces,
2885 mesh->vtx_interfaces,
2886 &gcell_vtx_idx,
2887 &gcell_vtx_lst);
2888
2889 cs_interface_set_destroy(&face_interfaces);
2890
2891 cs_halo_create_complete(mesh->halo);
2892
2893 t2 = cs_timer_wtime();
2894 halo_time = t2-t1;
2895
2896 } /* end if (mesh->n_domains > 1 || mesh->n_init_perio > 0) */
2897
2898 /* Define a cell -> cells connectivity for the extended neighborhood
2899 if necessary */
2900
2901 if (halo_type == CS_HALO_EXTENDED) {
2902
2903 t1 = cs_timer_wtime();
2904 if (verbosity > 0) {
2905 bft_printf(_(" Extended neighborhood structures definition\n"));
2906 bft_printf_flush();
2907 }
2908
2909 mesh->gcell_vtx_idx = gcell_vtx_idx;
2910 mesh->gcell_vtx_lst = gcell_vtx_lst;
2911
2912 cs_ext_neighborhood_define(mesh);
2913
2914 bft_printf_flush();
2915 t2 = cs_timer_wtime();
2916 ext_neighborhood_time = t2-t1;
2917
2918 }
2919 else {
2920 BFT_FREE(gcell_vtx_idx);
2921 BFT_FREE(gcell_vtx_lst);
2922 }
2923
2924 /* Output for log */
2925
2926 if (verbosity > 0) {
2927
2928 if (mesh->halo_type != CS_HALO_N_TYPES)
2929 _print_halo_info(mesh,
2930 interface_time,
2931 halo_time,
2932 ext_neighborhood_time);
2933
2934 else if (halo_type == CS_HALO_EXTENDED) {
2935 cs_log_printf(CS_LOG_PERFORMANCE,
2936 _("\nExtended connectivity creation (%.3g s)\n"),
2937 ext_neighborhood_time);
2938 cs_log_printf(CS_LOG_PERFORMANCE, "\n");
2939 cs_log_separator(CS_LOG_PERFORMANCE);
2940 }
2941
2942 _print_mesh_info(mesh);
2943
2944 }
2945
2946 mesh->verbosity = s_verbosity;
2947 }
2948
2949 /*----------------------------------------------------------------------------
2950 * Get the global number of ghost cells.
2951 *
2952 * parameters:
2953 * mesh <-- pointer to a mesh structure
2954 *
2955 * returns:
2956 * global number of ghost cells
2957 *---------------------------------------------------------------------------*/
2958
2959 cs_gnum_t
cs_mesh_n_g_ghost_cells(cs_mesh_t * mesh)2960 cs_mesh_n_g_ghost_cells(cs_mesh_t *mesh)
2961 {
2962 cs_gnum_t n_g_ghost_cells = mesh->n_ghost_cells;
2963
2964 #if defined(HAVE_MPI)
2965 if (cs_glob_n_ranks > 1) {
2966 cs_gnum_t _n_g_ghost_cells = n_g_ghost_cells;
2967 MPI_Allreduce(&_n_g_ghost_cells, &n_g_ghost_cells, 1, CS_MPI_GNUM,
2968 MPI_SUM, cs_glob_mpi_comm);
2969 }
2970 #endif
2971
2972 return n_g_ghost_cells;
2973 }
2974
2975 /*----------------------------------------------------------------------------
2976 * Order family numbers and remove duplicates
2977 *
2978 * parameters
2979 * mesh <-> pointer to mesh structure
2980 *----------------------------------------------------------------------------*/
2981
2982 void
cs_mesh_clean_families(cs_mesh_t * mesh)2983 cs_mesh_clean_families(cs_mesh_t *mesh)
2984 {
2985 size_t i, j, gc_id, gc_id_prev;
2986 int max_val = 0;
2987 size_t gc_count = 0;
2988 size_t n_gc = mesh->n_families;
2989 size_t n_gc_vals = mesh->n_max_family_items;
2990 size_t size_tot = n_gc * n_gc_vals;
2991 cs_gnum_t *interlaced = NULL;
2992 cs_lnum_t *order = NULL;
2993 int *renum = NULL;
2994
2995 if (mesh->n_families < 2)
2996 return;
2997
2998 /* Build and order interlaced copy with only positive values */
2999
3000 BFT_MALLOC(interlaced, size_tot, cs_gnum_t);
3001
3002 for (i = 0; i < size_tot; i++) {
3003 int val = mesh->family_item[i];
3004 if (val > max_val)
3005 max_val = val;
3006 }
3007
3008 for (i = 0; i < n_gc; i++) {
3009 for (j = 0; j < n_gc_vals; j++) {
3010 int val = mesh->family_item[j*n_gc + i];
3011 if (val < 0)
3012 val = -val + max_val;
3013 interlaced[i*n_gc_vals + j] = val;
3014 }
3015 }
3016
3017 order = cs_order_gnum_s(NULL, interlaced, n_gc_vals, n_gc);
3018
3019 /* Prepare removal of duplicates and renumbering */
3020
3021 BFT_MALLOC(renum, n_gc, int);
3022
3023 gc_id = order[0];
3024 gc_id_prev = gc_id;
3025 gc_count = 1;
3026 renum[gc_id] = 0;
3027
3028 for (i = 1; i < n_gc; i++) {
3029 char is_same = '\1';
3030 gc_id = order[i];
3031 for (j = 0; j < n_gc_vals; j++) {
3032 if ( interlaced[gc_id_prev*n_gc_vals + j]
3033 != interlaced[gc_id*n_gc_vals + j])
3034 is_same = '\0';
3035 }
3036 if (is_same != '\1') {
3037 gc_id_prev = gc_id;
3038 gc_count += 1;
3039 }
3040 renum[gc_id] = gc_count - 1;
3041 }
3042
3043 /* Update definitions */
3044
3045 mesh->n_families = gc_count;
3046 BFT_REALLOC(mesh->family_item, gc_count*n_gc_vals, int);
3047
3048 for (i = 0; i < n_gc; i++) {
3049 gc_id = renum[i];
3050 for (j = 0; j < n_gc_vals; j++)
3051 mesh->family_item[j*gc_count + gc_id] = interlaced[i*n_gc_vals + j];
3052 }
3053
3054 size_tot = gc_count * n_gc_vals;
3055 for (i = 0; i < size_tot; i++) {
3056 int val = mesh->family_item[i];
3057 if (val > max_val)
3058 val = -(val - max_val);
3059 mesh->family_item[i] = val;
3060 }
3061
3062 BFT_FREE(interlaced);
3063 BFT_FREE(order);
3064
3065 /* Update references */
3066
3067 if (mesh->cell_family != NULL) {
3068 for (i = 0; i < (size_t)(mesh->n_cells); i++) {
3069 int val = mesh->cell_family[i];
3070 if (val != 0)
3071 mesh->cell_family[i] = renum[val -1] + 1;
3072 }
3073 }
3074
3075 if (mesh->i_face_family != NULL) {
3076 for (i = 0; i < (size_t)(mesh->n_i_faces); i++) {
3077 int val = mesh->i_face_family[i];
3078 if (val != 0)
3079 mesh->i_face_family[i] = renum[val - 1] + 1;
3080 }
3081 }
3082
3083 if (mesh->b_face_family != NULL) {
3084 for (i = 0; i < (size_t)(mesh->n_b_faces); i++) {
3085 int val = mesh->b_face_family[i];
3086 if (val != 0)
3087 mesh->b_face_family[i] = renum[val - 1] + 1;
3088 }
3089 }
3090
3091 BFT_FREE(renum);
3092 }
3093
3094 /*----------------------------------------------------------------------------
3095 * Create group classes based on a mesh's family definitions.
3096 *
3097 * parameters:
3098 * mesh <-- pointer to mesh structure
3099 *
3100 * returns:
3101 * pointer to group classes structure based on mesh's family definitions
3102 *----------------------------------------------------------------------------*/
3103
3104 fvm_group_class_set_t *
cs_mesh_create_group_classes(cs_mesh_t * mesh)3105 cs_mesh_create_group_classes(cs_mesh_t *mesh)
3106 {
3107 int i, j;
3108 int grp_nbr, grp_num;
3109
3110 char **group = NULL;
3111
3112 fvm_group_class_set_t *class_defs = fvm_group_class_set_create();
3113
3114 /* Construction of the fvm_group_class structure */
3115
3116 BFT_MALLOC(group, mesh->n_max_family_items, char*);
3117
3118 for (i = 0; i < mesh->n_families; i++) {
3119
3120 grp_nbr = 0;
3121
3122 for (j = 0; j < mesh->n_max_family_items; j++) {
3123
3124 if (mesh->family_item[j * mesh->n_families + i] < 0) {
3125 /* Fortran formulation */
3126 grp_num = -mesh->family_item[j*mesh->n_families + i] -1;
3127 group[grp_nbr++] = mesh->group + mesh->group_idx[grp_num];
3128 }
3129
3130 }
3131
3132 fvm_group_class_set_add(class_defs,
3133 grp_nbr,
3134 (const char **)group);
3135
3136 } /* End of loop on families */
3137
3138 BFT_FREE(group);
3139
3140 return class_defs;
3141 }
3142
3143 /*----------------------------------------------------------------------------
3144 * Define group classes for a mesh based on its family definitions.
3145 *
3146 * parameters:
3147 * mesh <-> pointer to mesh structure
3148 *----------------------------------------------------------------------------*/
3149
3150 void
cs_mesh_init_group_classes(cs_mesh_t * mesh)3151 cs_mesh_init_group_classes(cs_mesh_t *mesh)
3152 {
3153 if (mesh->class_defs != NULL)
3154 mesh->class_defs = fvm_group_class_set_destroy(mesh->class_defs);
3155
3156 mesh->class_defs = cs_mesh_create_group_classes(mesh);
3157 }
3158
3159 /*----------------------------------------------------------------------------
3160 * Assign selectors to global mesh.
3161 *
3162 * Should be called once the mesh is fully built.
3163 *----------------------------------------------------------------------------*/
3164
3165 void
cs_mesh_init_selectors(void)3166 cs_mesh_init_selectors(void)
3167 {
3168 if (cs_glob_mesh->class_defs == NULL)
3169 cs_mesh_init_group_classes(cs_glob_mesh);
3170
3171 /* Construction of the selectors */
3172
3173 cs_glob_mesh->select_cells
3174 = fvm_selector_create(cs_glob_mesh->dim,
3175 cs_glob_mesh->n_cells,
3176 cs_glob_mesh->class_defs,
3177 cs_glob_mesh->cell_family,
3178 1,
3179 cs_glob_mesh_quantities->cell_cen,
3180 NULL);
3181
3182 cs_glob_mesh->select_b_faces
3183 = fvm_selector_create(cs_glob_mesh->dim,
3184 cs_glob_mesh->n_b_faces,
3185 cs_glob_mesh->class_defs,
3186 cs_glob_mesh->b_face_family,
3187 1,
3188 cs_glob_mesh_quantities->b_face_cog,
3189 cs_glob_mesh_quantities->b_face_normal);
3190
3191 cs_glob_mesh->select_i_faces
3192 = fvm_selector_create(cs_glob_mesh->dim,
3193 cs_glob_mesh->n_i_faces,
3194 cs_glob_mesh->class_defs,
3195 cs_glob_mesh->i_face_family,
3196 1,
3197 cs_glob_mesh_quantities->i_face_cog,
3198 cs_glob_mesh_quantities->i_face_normal);
3199
3200 }
3201
3202 /*----------------------------------------------------------------------------
3203 * Update selector and associated structures.
3204 *
3205 * parameters:
3206 * mesh <-> pointer to a mesh structure
3207 *----------------------------------------------------------------------------*/
3208
3209 void
cs_mesh_update_selectors(cs_mesh_t * mesh)3210 cs_mesh_update_selectors(cs_mesh_t *mesh)
3211 {
3212 _free_selectors(mesh);
3213 cs_mesh_init_selectors();
3214 }
3215
3216 /*----------------------------------------------------------------------------
3217 * Update a scalar array in case of parallelism and/or periodicity.
3218 *
3219 * Note: this function is only present so that a C equivalent to the
3220 * Fortran wrappers is available. In C code, directly using the
3221 * cs_halo_sync_var() is preferred.
3222 *
3223 * parameters:
3224 * var <-> scalar array
3225 *----------------------------------------------------------------------------*/
3226
3227 void
cs_mesh_sync_var_scal(cs_real_t * var)3228 cs_mesh_sync_var_scal(cs_real_t *var)
3229 {
3230 const cs_halo_t *halo = cs_glob_mesh->halo;
3231
3232 if (halo != NULL)
3233 cs_halo_sync_var(halo, CS_HALO_STANDARD, var);
3234 }
3235
3236 /*----------------------------------------------------------------------------
3237 * Update a scalar array in case of parallelism and/or periodicity,
3238 * using an extended halo.
3239 *
3240 * Note: this function is only present so that a C equivalent to the
3241 * Fortran wrappers is available. In C code, directly using the
3242 * cs_halo_sync_var() is preferred.
3243 *
3244 * parameters:
3245 * var <-> scalar array
3246 *----------------------------------------------------------------------------*/
3247
3248 void
cs_mesh_sync_var_scal_ext(cs_real_t * var)3249 cs_mesh_sync_var_scal_ext(cs_real_t *var)
3250 {
3251 const cs_halo_t *halo = cs_glob_mesh->halo;
3252
3253 if (halo != NULL)
3254 cs_halo_sync_var(halo, CS_HALO_EXTENDED, var);
3255 }
3256
3257 /*----------------------------------------------------------------------------
3258 * Update a vector array in case of parallelism and/or periodicity.
3259 *
3260 * parameters:
3261 * var <-> interleaved vector (of dimension 3)
3262 *----------------------------------------------------------------------------*/
3263
3264 void
cs_mesh_sync_var_vect(cs_real_t * var)3265 cs_mesh_sync_var_vect(cs_real_t *var)
3266 {
3267 const cs_halo_t *halo = cs_glob_mesh->halo;
3268
3269 if (halo == NULL) return;
3270
3271 cs_halo_sync_var_strided(halo, CS_HALO_STANDARD, var, 3);
3272
3273 if (cs_glob_mesh->n_init_perio > 0)
3274 cs_halo_perio_sync_var_vect(halo,
3275 CS_HALO_STANDARD,
3276 var,
3277 3);
3278 }
3279
3280 /*----------------------------------------------------------------------------
3281 * Update a vector array in case of parallelism and/or periodicity,
3282 * using an extended halo.
3283 *
3284 * parameters:
3285 * var <-> interleaved vector (of dimension 3)
3286 *----------------------------------------------------------------------------*/
3287
3288 void
cs_mesh_sync_var_vect_ext(cs_real_t * var)3289 cs_mesh_sync_var_vect_ext(cs_real_t *var)
3290 {
3291 const cs_halo_t *halo = cs_glob_mesh->halo;
3292
3293 if (halo == NULL) return;
3294
3295 cs_halo_sync_var_strided(halo, CS_HALO_EXTENDED, var, 3);
3296
3297 if (cs_glob_mesh->n_init_perio > 0)
3298 cs_halo_perio_sync_var_vect(halo,
3299 CS_HALO_EXTENDED,
3300 var,
3301 3);
3302 }
3303
3304 /*----------------------------------------------------------------------------
3305 * Update a tensor array in case of parallelism and/or periodicity.
3306 *
3307 * parameters:
3308 * var <-> interleaved tensor (of dimension 3x3)
3309 *----------------------------------------------------------------------------*/
3310
3311 void
cs_mesh_sync_var_tens(cs_real_t * var)3312 cs_mesh_sync_var_tens(cs_real_t *var)
3313 {
3314 const cs_halo_t *halo = cs_glob_mesh->halo;
3315
3316 if (halo == NULL) return;
3317
3318 cs_halo_sync_var_strided(halo, CS_HALO_STANDARD, var, 9);
3319
3320 if (cs_glob_mesh->n_init_perio > 0)
3321 cs_halo_perio_sync_var_tens(halo,
3322 CS_HALO_STANDARD,
3323 var);
3324 }
3325
3326 /*----------------------------------------------------------------------------
3327 * Update a symmetric tensor array in case of parallelism and/or periodicity.
3328 *
3329 * parameters:
3330 * var <-> symmetric interleaved tensor (of dimension 6)
3331 *----------------------------------------------------------------------------*/
3332
3333 void
cs_mesh_sync_var_sym_tens(cs_real_6_t * var)3334 cs_mesh_sync_var_sym_tens(cs_real_6_t *var)
3335 {
3336 const cs_halo_t *halo = cs_glob_mesh->halo;
3337
3338 if (halo == NULL) return;
3339
3340 cs_halo_sync_var_strided(halo, CS_HALO_STANDARD, (cs_real_t *)var, 6);
3341
3342 if (cs_glob_mesh->n_init_perio > 0)
3343 cs_halo_perio_sync_var_sym_tens(halo,
3344 CS_HALO_STANDARD,
3345 (cs_real_t *)var);
3346 }
3347
3348 /*----------------------------------------------------------------------------
3349 * Get global lists of periodic face couples.
3350 *
3351 * In parallel, each face couple may appear on only one rank.
3352 *
3353 * The caller is responsible for freeing the arrays allocated and returned
3354 * by this function once they are no onger needed.
3355 *
3356 * parameters:
3357 * mesh <-- pointer to mesh structure
3358 * n_perio_face_couples --> global number of periodic couples per
3359 * periodicity (size: mesh->n_init_perio)
3360 * perio_face_couples --> arrays of global periodic couple face numbers,
3361 * for each periodicity
3362 *----------------------------------------------------------------------------*/
3363
3364 void
cs_mesh_get_perio_faces(const cs_mesh_t * mesh,cs_lnum_t ** n_perio_face_couples,cs_gnum_t *** perio_face_couples)3365 cs_mesh_get_perio_faces(const cs_mesh_t *mesh,
3366 cs_lnum_t **n_perio_face_couples,
3367 cs_gnum_t ***perio_face_couples)
3368 {
3369 assert (mesh != NULL);
3370
3371 if (mesh->n_init_perio == 0)
3372 return;
3373
3374 #if defined(HAVE_MPI)
3375
3376 if (cs_glob_n_ranks > 1)
3377 _get_perio_faces_g(mesh,
3378 n_perio_face_couples,
3379 perio_face_couples);
3380
3381 #endif /* defined(HAVE_MPI) */
3382
3383 if (cs_glob_n_ranks == 1)
3384 _get_perio_faces_l(mesh,
3385 n_perio_face_couples,
3386 perio_face_couples);
3387
3388 /* Ensure couples are sorted (not normally necessary, but safe) */
3389
3390 if (n_perio_face_couples != NULL) {
3391 int i;
3392 for (i = 0; i < mesh->n_init_perio; i++) {
3393 if ((*n_perio_face_couples)[i] > 1)
3394 qsort((*perio_face_couples)[i],
3395 (*n_perio_face_couples)[i],
3396 sizeof(cs_gnum_t) * 2,
3397 &_compare_couples);
3398 }
3399 }
3400 }
3401
3402 /*----------------------------------------------------------------------------
3403 * Build global cell numbering array extended to ghost cell values.
3404 *
3405 * If the blank_perio flag is nonzero, periodic ghost cell numbers
3406 * are set to zero instead of the value of the matching cell.
3407 *
3408 * The caller is responsible for freeing the returned array when it
3409 * is no longer useful.
3410 *
3411 * parameters:
3412 * mesh <-- pointer to mesh structure
3413 * blank_perio <-- flag to zeroe periodic cell values
3414 *----------------------------------------------------------------------------*/
3415
3416 cs_gnum_t *
cs_mesh_get_cell_gnum(const cs_mesh_t * mesh,int blank_perio)3417 cs_mesh_get_cell_gnum(const cs_mesh_t *mesh,
3418 int blank_perio)
3419 {
3420 cs_lnum_t i;
3421 cs_gnum_t *cell_gnum = NULL;
3422
3423 assert (mesh != NULL);
3424
3425 /* Allocate array */
3426
3427 BFT_MALLOC(cell_gnum, mesh->n_cells_with_ghosts, cs_gnum_t);
3428
3429 /* Build global cell numbering including parallel halos */
3430
3431 for (i = 0; i < mesh->n_cells; i++)
3432 cell_gnum[i] = mesh->global_cell_num[i];
3433 for (i = mesh->n_cells; i < mesh->n_cells_with_ghosts; i++)
3434 cell_gnum[i] = 0;
3435
3436 if (mesh->halo != NULL) {
3437
3438 cs_halo_sync_untyped(mesh->halo,
3439 CS_HALO_EXTENDED,
3440 sizeof(cs_gnum_t),
3441 cell_gnum);
3442
3443 if (blank_perio) {
3444
3445 const cs_halo_t *halo = mesh->halo;
3446
3447 cs_lnum_t rank_id, t_id, shift;
3448 cs_lnum_t start = 0, end = 0;
3449
3450 const cs_lnum_t n_transforms = halo->n_transforms;
3451 const cs_lnum_t n_elts = halo->n_local_elts;
3452
3453 for (t_id = 0; t_id < n_transforms; t_id++) {
3454
3455 shift = 4 * halo->n_c_domains * t_id;
3456
3457 for (rank_id = 0; rank_id < halo->n_c_domains; rank_id++) {
3458
3459 start = halo->perio_lst[shift + 4*rank_id];
3460 end = start + halo->perio_lst[shift + 4*rank_id + 1];
3461 for (i = start; i < end; i++)
3462 cell_gnum[n_elts+i] = 0;
3463
3464 start = halo->perio_lst[shift + 4*rank_id + 2];
3465 end = start + halo->perio_lst[shift + 4*rank_id + 3];
3466 for (i = start; i < end; i++)
3467 cell_gnum[n_elts+i] = 0;
3468
3469 } /* End of loop on ranks */
3470
3471 } /* End of loop on transformation */
3472
3473 } /* End for blank_perio */
3474
3475 }
3476
3477 /* Return global cell number */
3478
3479 return cell_gnum;
3480 }
3481
3482 /*----------------------------------------------------------------------------
3483 * Mark interior faces with the number of their associated periodic
3484 * transform id.
3485 *
3486 * parameters:
3487 * mesh <-- pointer to mesh structure
3488 * perio_num --> periodicity number associated with each face, signed for
3489 * direct/reverse transform, 0 for non-periodic faces
3490 * (size: mesh->n_i_faces)
3491 *----------------------------------------------------------------------------*/
3492
3493 void
cs_mesh_get_face_perio_num(const cs_mesh_t * mesh,int perio_num[])3494 cs_mesh_get_face_perio_num(const cs_mesh_t *mesh,
3495 int perio_num[])
3496 {
3497 cs_lnum_t i;
3498
3499 for (i = 0; i < mesh->n_i_faces; i++)
3500 perio_num[i] = 0;
3501
3502 if (mesh->n_init_perio > 0) {
3503
3504 int *halo_perio_num = NULL;
3505
3506 /* Mark ghost cells with their periodicity number */
3507
3508 BFT_MALLOC(halo_perio_num, mesh->n_ghost_cells, int);
3509
3510 _get_halo_perio_num(mesh, halo_perio_num, NULL);
3511
3512 for (i = 0; i < mesh->n_i_faces; i++) {
3513 const cs_lnum_t h_id0 = mesh->i_face_cells[i][0] - mesh->n_cells;
3514 const cs_lnum_t h_id1 = mesh->i_face_cells[i][1] - mesh->n_cells;
3515 if (h_id0 >= 0) {
3516 if (halo_perio_num[h_id0] != 0)
3517 perio_num[i] = halo_perio_num[h_id0];
3518 }
3519 else if (h_id1 >= 0) {
3520 if (halo_perio_num[h_id1] != 0)
3521 perio_num[i] = halo_perio_num[h_id1];
3522 }
3523 }
3524
3525 BFT_FREE(halo_perio_num);
3526 }
3527 }
3528
3529 /*----------------------------------------------------------------------------
3530 * Print information on a mesh structure.
3531 *
3532 * parameters:
3533 * mesh <-- pointer to mesh structure.
3534 * name <-- associated name.
3535 *----------------------------------------------------------------------------*/
3536
3537 void
cs_mesh_print_info(const cs_mesh_t * mesh,const char * name)3538 cs_mesh_print_info(const cs_mesh_t *mesh,
3539 const char *name)
3540 {
3541 if (mesh->n_g_vertices > 0) {
3542
3543 cs_lnum_t i, vtx_id;
3544 cs_lnum_t dim = mesh->dim;
3545 cs_real_t min_xyz[3] = { 1.e127, 1.e127, 1.e127};
3546 cs_real_t max_xyz[3] = {-1.e127, -1.e127, -1.e127};
3547
3548 for (vtx_id = 0 ; vtx_id < mesh->n_vertices ; vtx_id++) {
3549
3550 for (i = 0 ; i < dim ; i++) {
3551
3552 if (mesh->vtx_coord[vtx_id*dim + i] < min_xyz[i])
3553 min_xyz[i] = mesh->vtx_coord[vtx_id*dim + i];
3554
3555 if (mesh->vtx_coord[vtx_id*dim + i] > max_xyz[i])
3556 max_xyz[i] = mesh->vtx_coord[vtx_id*dim + i];
3557
3558 }
3559
3560 } /* End of loop on vertices */
3561
3562 #if defined(HAVE_MPI)
3563
3564 if (cs_glob_n_ranks > 1) {
3565 cs_real_t g_min_xyz[3];
3566 cs_real_t g_max_xyz[3];
3567 MPI_Allreduce(min_xyz, g_min_xyz, dim, CS_MPI_REAL, MPI_MIN,
3568 cs_glob_mpi_comm);
3569 MPI_Allreduce(max_xyz, g_max_xyz, dim, CS_MPI_REAL, MPI_MAX,
3570 cs_glob_mpi_comm);
3571 for (i = 0 ; i < dim ; i++) {
3572 min_xyz[i] = g_min_xyz[i];
3573 max_xyz[i] = g_max_xyz[i];
3574 }
3575 }
3576
3577 #endif
3578
3579 bft_printf(_("\n"
3580 " Mesh coordinates: minimum and maximum\n"
3581 " X : %14.7e %14.7e\n"
3582 " Y : %14.7e %14.7e\n"
3583 " Z : %14.7e %14.7e\n"),
3584 min_xyz[0], max_xyz[0], min_xyz[1], max_xyz[1],
3585 min_xyz[2], max_xyz[2]);
3586
3587 }
3588
3589 bft_printf(_(" %s\n"
3590 " Number of cells: %llu\n"
3591 " Number of interior faces: %llu\n"
3592 " Number of boundary faces: %llu\n"
3593 " Number of vertices: %llu\n"),
3594 name,
3595 (unsigned long long)(mesh->n_g_cells),
3596 (unsigned long long)(mesh->n_g_i_faces),
3597 (unsigned long long)(mesh->n_g_b_faces - mesh->n_g_free_faces),
3598 (unsigned long long)(mesh->n_g_vertices));
3599
3600 if (mesh->n_g_free_faces > 0)
3601 bft_printf(_("\n"
3602 " Number of isolated faces: %llu\n"),
3603 (unsigned long long)(mesh->n_g_free_faces));
3604
3605 _print_mesh_group_stats(mesh);
3606 }
3607
3608 /*----------------------------------------------------------------------------
3609 * Print statistics about mesh selectors usage to log.
3610 *
3611 * parameters:
3612 * mesh <-- pointer to a mesh structure
3613 *----------------------------------------------------------------------------*/
3614
3615 void
cs_mesh_selector_stats(cs_mesh_t * mesh)3616 cs_mesh_selector_stats(cs_mesh_t *mesh)
3617 {
3618 int n_calls[3] = {0, 0, 0};
3619 double wtimes[3] = {0., 0., 0.};
3620
3621 if (mesh->select_cells != NULL)
3622 fvm_selector_get_stats(mesh->select_cells, n_calls, wtimes);
3623 if (mesh->select_i_faces != NULL)
3624 fvm_selector_get_stats(mesh->select_i_faces, n_calls + 1, wtimes + 1);
3625 if (mesh->select_b_faces != NULL)
3626 fvm_selector_get_stats(mesh->select_b_faces, n_calls + 2, wtimes + 2);
3627
3628 #if defined(HAVE_MPI)
3629 if (cs_glob_n_ranks > 1) {
3630 int i;
3631 double wtimes_glob[3];
3632 MPI_Allreduce(wtimes, wtimes_glob, 3, MPI_DOUBLE, MPI_MAX,
3633 cs_glob_mpi_comm);
3634 for (i = 0; i < 3; i++)
3635 wtimes[i] = wtimes_glob[i];
3636 }
3637 #endif
3638
3639 cs_log_printf(CS_LOG_PERFORMANCE,
3640 _("\n"
3641 "Mesh entity selections by criteria statistics:\n\n"
3642 " entity type evaluations time\n"
3643 " -----------------------------------------\n"
3644 " cells %10d %12.5f\n"
3645 " interior faces %10d %12.5f\n"
3646 " boundary faces %10d %12.5f\n"),
3647 n_calls[0], wtimes[0],
3648 n_calls[1], wtimes[1],
3649 n_calls[2], wtimes[2]);
3650
3651 cs_log_printf(CS_LOG_PERFORMANCE, "\n");
3652 cs_log_separator(CS_LOG_PERFORMANCE);
3653 }
3654
3655 /*----------------------------------------------------------------------------
3656 * Dump of a mesh structure.
3657 *
3658 * parameters:
3659 * mesh <-- pointer to mesh structure.
3660 *----------------------------------------------------------------------------*/
3661
3662 void
cs_mesh_dump(const cs_mesh_t * mesh)3663 cs_mesh_dump(const cs_mesh_t *mesh)
3664 {
3665 bft_printf("\n\nDUMP OF THE MESH STRUCTURE: %p\n\n", (const void *)mesh);
3666
3667 bft_printf("space dim : %d\n"
3668 "n_domains : %d\n"
3669 "domain_num: %d\n",
3670 (int)mesh->dim, mesh->n_domains, mesh->domain_num);
3671
3672 bft_printf("\nNumber of families: %3d\n",mesh->n_families);
3673
3674 bft_printf("\nLocal dimensions:\n"
3675 "n_cells: %ld\n"
3676 "n_cells_with_ghosts: %ld\n"
3677 "n_vertices: %ld\n"
3678 "n_i_faces: %ld\n"
3679 "n_b_faces: %ld\n",
3680 (long)mesh->n_cells, (long)mesh->n_cells_with_ghosts,
3681 (long)mesh->n_vertices,
3682 (long)mesh->n_i_faces, (long)mesh->n_b_faces);
3683
3684 bft_printf("\nGlobal dimensions:\n"
3685 "n_g_cells: %llu\n"
3686 "n_g_vertices: %llu\n"
3687 "n_g_i_faces: %llu\n"
3688 "n_g_b_faces: %llu\n",
3689 (unsigned long long)mesh->n_g_cells,
3690 (unsigned long long)mesh->n_g_vertices,
3691 (unsigned long long)mesh->n_g_i_faces,
3692 (unsigned long long)mesh->n_g_b_faces);
3693
3694 bft_printf("\n\n --------"
3695 " Vertices"
3696 " --------\n\n");
3697
3698 bft_printf("\nVertex coordinates:\n");
3699 for (cs_lnum_t i = 0; i < mesh->n_vertices; i++)
3700 bft_printf(" <%3ld > %10.3f %10.3f %10.3f\n",
3701 (long)i, mesh->vtx_coord[3*i], mesh->vtx_coord[3*i+1],
3702 mesh->vtx_coord[3*i+2]);
3703
3704 if (mesh->global_vtx_num != NULL) {
3705 bft_printf("\nGlobal vertex numbering:\n");
3706 for (cs_lnum_t i = 0; i < mesh->n_vertices; i++)
3707 bft_printf(" <%7ld > %10llu\n",
3708 (long)i, (unsigned long long)(mesh->global_vtx_num[i]));
3709 }
3710
3711 bft_printf("\n\n ---------------------------"
3712 " Internal faces connectivity"
3713 " ---------------------------\n\n");
3714
3715 bft_printf("\nInternal faces -> Cells connectivity:\n");
3716 for (cs_lnum_t i = 0; i < mesh->n_i_faces; i++)
3717 bft_printf(" < %7ld > %7ld <----> %7ld\n", (long)i,
3718 (long)mesh->i_face_cells[i][0],
3719 (long)mesh->i_face_cells[i][1]);
3720
3721 bft_printf("\nInternal faces -> vertices connectivity:\n");
3722 for (cs_lnum_t i = 0; i < mesh->n_i_faces; i++) {
3723 bft_printf(" < %7ld >", (long)i);
3724 for (cs_lnum_t j = mesh->i_face_vtx_idx[i];
3725 j < mesh->i_face_vtx_idx[i+1];
3726 j++)
3727 bft_printf(" %7ld ", (long)mesh->i_face_vtx_lst[j]);
3728 bft_printf("\n");
3729 }
3730
3731 bft_printf("\nFamily of each internal face:\n");
3732 for (cs_lnum_t i = 0; i < mesh->n_i_faces; i++)
3733 bft_printf(" < %3ld > %5d\n", (long)i, mesh->i_face_family[i]);
3734
3735 if (mesh->i_face_r_gen != NULL) {
3736 bft_printf("Refinement generation of each internal face:\n");
3737 for (cs_lnum_t i = 0; i < mesh->n_i_faces; i++)
3738 bft_printf(" < %3ld > %5d\n", (long)i, (int)(mesh->i_face_r_gen[i]));
3739 }
3740
3741 if (mesh->global_i_face_num != NULL) {
3742
3743 bft_printf("\nInternal faces global numbering:\n");
3744 for (cs_lnum_t i = 0; i < mesh->n_i_faces; i++)
3745 bft_printf(" < %7ld > %12llu\n",
3746 (long)i, (unsigned long long)(mesh->global_i_face_num[i]));
3747 bft_printf("\n");
3748
3749 }
3750
3751 bft_printf("\n\n -------------------------"
3752 " Border faces connectivity"
3753 " -------------------------\n\n");
3754
3755 bft_printf("\nBorder faces -> Cells connectivity:\n");
3756 for (cs_lnum_t i = 0; i < mesh->n_b_faces; i++)
3757 bft_printf(" < %7ld > %7ld\n", (long)i, (long)mesh->b_face_cells[i]);
3758
3759 bft_printf("\nBorder faces -> vertices connectivity:\n");
3760 for (cs_lnum_t i = 0; i < mesh->n_b_faces; i++) {
3761 bft_printf(" < %7ld >", (long)i);
3762 for (cs_lnum_t j = mesh->b_face_vtx_idx[i];
3763 j < mesh->b_face_vtx_idx[i+1];
3764 j++)
3765 bft_printf(" %7ld ", (long)mesh->b_face_vtx_lst[j]);
3766 bft_printf("\n");
3767 }
3768
3769 bft_printf("\nFamily of each boundary face:\n");
3770 for (cs_lnum_t i = 0; i < mesh->n_b_faces; i++)
3771 bft_printf(" < %3ld > %5d\n", (long)i, mesh->b_face_family[i]);
3772
3773 if (mesh->global_b_face_num != NULL) {
3774
3775 bft_printf("\nBoundary faces global numbering:\n");
3776 for (cs_lnum_t i = 0; i < mesh->n_b_faces; i++)
3777 bft_printf(" < %7ld > %12llu\n",
3778 (long)i, (unsigned long long)(mesh->global_b_face_num[i]));
3779 bft_printf("\n");
3780
3781 }
3782
3783 bft_printf("\n\n -------------------------"
3784 " Cells"
3785 " -------------------------\n\n");
3786
3787 if (mesh->global_cell_num != NULL) {
3788
3789 bft_printf("\nCell global numbering:\n");
3790 for (cs_lnum_t i = 0; i < mesh->n_cells; i++)
3791 bft_printf(" < %7ld > %12llu\n", (long)i,
3792 (unsigned long long)(mesh->global_cell_num[i]));
3793 bft_printf("\n");
3794
3795 }
3796
3797 bft_printf("Family of each cell:\n");
3798 for (cs_lnum_t i = 0; i < mesh->n_cells_with_ghosts; i++)
3799 bft_printf(" < %3ld > %5d\n", (long)i, mesh->cell_family[i]);
3800
3801 if (mesh->halo != NULL) {
3802
3803 cs_halo_t *halo = mesh->halo;
3804
3805 bft_printf("\nHalo information: %p\n", (const void *)halo);
3806
3807 bft_printf("n_c_domains: %d\n", halo->n_c_domains);
3808 bft_printf("n_ghost_cells: %ld\n", (long)mesh->n_ghost_cells);
3809 bft_printf("n_std_ghost_cells: %ld\n",
3810 (long)halo->n_elts[CS_HALO_STANDARD]);
3811 bft_printf("n_ext_ghost_cells: %ld\n",
3812 (long)( halo->n_elts[CS_HALO_EXTENDED]
3813 - halo->n_elts[CS_HALO_STANDARD]));
3814
3815 for (int i = 0; i < halo->n_c_domains; i++) {
3816
3817 bft_printf("\n\nRank id: %d\n"
3818 "Halo index start: %ld end: %ld\n"
3819 "Send index start: %ld end: %ld\n"
3820 "Send cell ids:\n",
3821 halo->c_domain_rank[i],
3822 (long)halo->index[2*i], (long)halo->index[2*i+2],
3823 (long)halo->send_index[2*i], (long)halo->send_index[2*i+2]);
3824 for (cs_lnum_t j = halo->send_index[2*i]; j < halo->send_index[2*i+2]; j++)
3825 bft_printf(" %10ld : %10ld\n", (long)j, (long)halo->send_list[j]);
3826
3827 } /* End of loop on the frontiers of halo */
3828
3829 if (mesh->n_init_perio > 0 && halo->perio_lst != NULL) {
3830
3831 const cs_lnum_t n_c_domains = halo->n_c_domains;
3832 const cs_lnum_t n_transforms = mesh->n_transforms;
3833
3834 bft_printf("\n\nHalo's data in case of periodicity:\n");
3835 bft_printf("n_transforms: %d\n",mesh->n_transforms);
3836
3837 bft_printf("\nData in the standard halo\n");
3838 for (int i = 0; i < n_transforms; i++)
3839 for (int j = 0; j < n_c_domains; j++)
3840 bft_printf("< rank:%3d >< transform:%2d > start_idx: %5ld"
3841 " n_elts: %5ld\n",
3842 halo->c_domain_rank[j], i,
3843 (long)halo->perio_lst[4*n_c_domains*i + 4*j],
3844 (long)halo->perio_lst[4*n_c_domains*i + 4*j+1]);
3845
3846 bft_printf("\nData in the extended halo\n");
3847 for (int i = 0; i < n_transforms; i++)
3848 for (int j = 0; j < n_c_domains; j++)
3849 bft_printf("< rank:%3d >< transform:%2d > "
3850 "start_idx: %5ld, n_elts: %5ld\n",
3851 halo->c_domain_rank[j], i,
3852 (long)halo->perio_lst[4*n_c_domains*i + 4*j+2],
3853 (long)halo->perio_lst[4*n_c_domains*i + 4*j+3]);
3854
3855 } /* End if n_perio > 0 */
3856
3857 } /* End if mesh->halo != NULL */
3858
3859 if (mesh->cell_cells_idx != NULL) {
3860
3861 bft_printf("\n\nCell -> cells connectivity for extended neighborhood\n\n");
3862 for (cs_lnum_t i = 0; i < mesh->n_cells; i++) {
3863 bft_printf("< cell id:%3ld> ", (long)i);
3864 for (cs_lnum_t j = mesh->cell_cells_idx[i];
3865 j < mesh->cell_cells_idx[i+1];
3866 j++)
3867 bft_printf("%ld ", (long)mesh->cell_cells_lst[j]);
3868 bft_printf("\n");
3869 }
3870
3871 }
3872
3873 if (mesh->gcell_vtx_idx != NULL) {
3874
3875 bft_printf("\n\nGhost cell -> vertices connectivity "
3876 "for extended neighborhood\n\n");
3877 for (cs_lnum_t i = 0; i < mesh->n_ghost_cells; i++) {
3878 bft_printf("< gcell id:%3ld> ", (long)i + mesh->n_cells);
3879 for (cs_lnum_t j = mesh->gcell_vtx_idx[i];
3880 j < mesh->gcell_vtx_idx[i+1];
3881 j++)
3882 bft_printf("%ld ", (long)mesh->gcell_vtx_lst[j]);
3883 bft_printf("\n");
3884 }
3885
3886 }
3887
3888 /* Dump numbering info */
3889
3890 cs_numbering_dump(mesh->cell_numbering);
3891 cs_numbering_dump(mesh->i_face_numbering);
3892 cs_numbering_dump(mesh->b_face_numbering);
3893 cs_numbering_dump(mesh->vtx_numbering);
3894
3895 /* Modification flag */
3896
3897 bft_printf("\nModification flag:\n");
3898 bft_printf("modified: %d\n",mesh->modified);
3899
3900 bft_printf("\n\nEND OF DUMP OF MESH STRUCTURE\n\n");
3901 bft_printf_flush();
3902 }
3903
3904 /*----------------------------------------------------------------------------*/
3905
3906 END_C_DECLS
3907