1 /*============================================================================
2 * Main structure for an I/O numbering scheme associated with mesh entities
3 * (such as cells, faces, and vertices);
4 *
5 * In parallel mode, such a scheme is important so as to redistribute
6 * locally numbered entities on n processes to files written by p
7 * processes, with p <= n.
8 *
9 * Only the case where p = 1 is presently implemented, so the numbering
10 * scheme is simply based on entity's global labels.
11 *
12 * For p > 1, it would probably be necessary to extend the numbering
13 * schemes so as to account for the fact that a given entity may have
14 * a main index on its main associated domain, but may be present
15 * as a ghost entity with another index on neighboring domains.
16 *============================================================================*/
17
18 /*
19 This file is part of Code_Saturne, a general-purpose CFD tool.
20
21 Copyright (C) 1998-2021 EDF S.A.
22
23 This program is free software; you can redistribute it and/or modify it under
24 the terms of the GNU General Public License as published by the Free Software
25 Foundation; either version 2 of the License, or (at your option) any later
26 version.
27
28 This program is distributed in the hope that it will be useful, but WITHOUT
29 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
30 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
31 details.
32
33 You should have received a copy of the GNU General Public License along with
34 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
35 Street, Fifth Floor, Boston, MA 02110-1301, USA.
36 */
37
38 /*----------------------------------------------------------------------------*/
39
40 #include "cs_defs.h"
41
42 /*----------------------------------------------------------------------------
43 * Standard C library headers
44 *----------------------------------------------------------------------------*/
45
46 #include <assert.h>
47 #include <float.h>
48 #include <math.h>
49 #include <stdio.h>
50 #include <string.h>
51
52 /*----------------------------------------------------------------------------
53 * Local headers
54 *----------------------------------------------------------------------------*/
55
56 #include "bft_mem.h"
57 #include "bft_printf.h"
58
59 #include "fvm_hilbert.h"
60 #include "fvm_morton.h"
61
62 #include "cs_all_to_all.h"
63 #include "cs_order.h"
64 #include "cs_parall.h"
65 #include "cs_sort_partition.h"
66
67 /*----------------------------------------------------------------------------
68 * Header for the current file
69 *----------------------------------------------------------------------------*/
70
71 #include "fvm_io_num.h"
72
73 /*----------------------------------------------------------------------------*/
74
75 BEGIN_C_DECLS
76
77 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
78
79 /*============================================================================
80 * Local structure definitions
81 *============================================================================*/
82
83 /*----------------------------------------------------------------------------
84 * Structure defining an I/O numbering scheme
85 *----------------------------------------------------------------------------*/
86
87 /*
88 * Notes:
89 *
90 * This structure currently only contains a global numbering array containing
91 * each entity's global number (a 1 to n index). In the future, it may
92 * also contain information relative to ghost zones, for I/O to file formats
93 * enabling domain splitting with (with multiple groups of processes writing
94 * different subsets, using ghost zones for entities on boundaries between
95 * mesh parts assigned to different process groups). In such a case, the
96 * main numbering would only be global as regards a process group, at least
97 * for use with formats such as that of EnSight Gold.
98 *
99 * In some cases, a global number may appear on multiple processes (for a face
100 * or vertex on a processor boundary). This means that multiple processes may
101 * update the corresponding value in a gather-type operation preceding or
102 * associated with I/O, in an undefined order. If the associated data is up to
103 * date on each process, it should be identical, so this should not be a
104 * problem. MPI-IO file writes or MPI-2 one-sided communication PUT operations
105 * also authorize this, though for the latter, the MPI-2 standard indicates
106 * that this is authorized if processes overlapping PUT operations should use
107 * the same predefined datatype, which seems to exclude similar indexed
108 * datatypes with different indexes. To avoid problems if we wish to use
109 * MPI-2 one-sided communication, one relatively simple solution would be to
110 * consider that for processes other than that of lowest rank in which an
111 * entity appears, it appears after all entities not occuring in processes
112 * of lower rank. In that case, we would have two array sizes:
113 * global_num_size_tot defining the array's full size, and global_num_size
114 * defining the size of the portion to use in gather type operations.
115 */
116
117 struct _fvm_io_num_t {
118
119 cs_gnum_t global_count; /* Global number of entities */
120 cs_lnum_t global_num_size; /* Local size of global numbering array */
121 const cs_gnum_t *global_num; /* Global (possibly shared) entity
122 numbers (1 to n) */
123 cs_gnum_t *_global_num; /* Global entity numbers if owner,
124 NULL otherwise */
125
126 };
127
128 /*=============================================================================
129 * Static global variables
130 *============================================================================*/
131
132 /* Names of space-filling curve types */
133
134 const char *fvm_io_num_sfc_type_name[] = {N_("Morton (in bounding box)"),
135 N_("Morton (in bounding cube)"),
136 N_("Hilbert (in bounding box)"),
137 N_("Hilbert (in bounding cube)")};
138
139 static const int _sampling_factors[4] = {1, /* OD */
140 2, /* 1D */
141 2, /* 2D */
142 4, /* 3D */};
143
144 /*=============================================================================
145 * Private function definitions
146 *============================================================================*/
147
148 /*----------------------------------------------------------------------------
149 * Function pointer for conversion of a double precision value in
150 * range [0, 1] to the same value.
151 *
152 * This is a trivial function, passed through a function pointer
153 * to cs_sort_partition_dest_rank_id.
154 *
155 * parameters:
156 * s <-- coordinate between 0 and 1
157 * elt --> pointer to element
158 * input <-- pointer to optional (untyped) value or structure.
159 */
160 /*----------------------------------------------------------------------------*/
161
162 static void
_s_to_real(double s,void * elt,const void * input)163 _s_to_real(double s,
164 void *elt,
165 const void *input)
166 {
167 CS_UNUSED(input);
168
169 cs_real_t *v = elt;
170 *v = s;
171 }
172
173 /*----------------------------------------------------------------------------*/
174 /*!
175 * \brief function pointer for comparison of 2
176 *
177 * This function is the same type as that used by qsort_r.
178 *
179 * \param[in] elt1 coordinate between 0 and 1
180 * \param[in] elt2 pointer to optional (untyped) value or structure.
181 * \param[in] input pointer to optional (untyped) value or structure.
182 *
183 * \return < 0 if elt1 < elt2, 0 if elt1 == elt2, > 0 if elt1 > elt2
184 */
185 /*----------------------------------------------------------------------------*/
186
187 static int
_s_compare(const void * elt1,const void * elt2,const void * input)188 _s_compare(const void *elt1,
189 const void *elt2,
190 const void *input)
191 {
192 CS_UNUSED(input);
193
194 int retval = 0;
195 if ( *(const cs_real_t *)elt1
196 < *(const cs_real_t *)elt2)
197 retval = -1;
198 else if ( *(const cs_real_t *)elt1
199 > *(const cs_real_t *)elt2)
200 retval = 1;
201
202 return retval;
203 }
204
205 /*----------------------------------------------------------------------------
206 * Use bubble sort on an expectedly short sequence of coordinates
207 * to ensure lexicographical ordering.
208 *
209 * parameters:
210 * dim <-- spatial dimension
211 * start_id <-- start id in array
212 * end_id <-- past-the-end id in array
213 * coords <-- pointer to entity coordinates (interlaced)
214 * order <-> ordering array base on Morton encoding, or
215 * lexicographical coordinate ordering for ties
216 *----------------------------------------------------------------------------*/
217
218 inline static void
_reorder_coords_lexicographic(int dim,size_t start_id,size_t end_id,const cs_coord_t coords[],cs_lnum_t order[])219 _reorder_coords_lexicographic(int dim,
220 size_t start_id,
221 size_t end_id,
222 const cs_coord_t coords[],
223 cs_lnum_t order[])
224 {
225 size_t i;
226 bool g_swap;
227
228 do {
229
230 g_swap = false;
231
232 for (i = start_id + 1; i < end_id; i++) {
233
234 size_t j_prev = order[i-1], j = order[i];
235 bool l_swap = false;
236
237 if (dim == 3) {
238 if (coords[j_prev*3] < coords[j*3])
239 continue;
240 else if (coords[j_prev*3] > coords[j*3])
241 l_swap = true;
242 else if (coords[j_prev*3 + 1] < coords[j*3 + 1])
243 continue;
244 else if ( coords[j_prev*3 + 1] > coords[j*3 + 1]
245 || coords[j_prev*3 + 2] > coords[j*3 + 2])
246 l_swap = true;
247 }
248 else if (dim == 2) {
249 if (coords[j_prev*2] < coords[j*2 + 1])
250 continue;
251 else if ( coords[j_prev*2] > coords[j*2]
252 || coords[j_prev*2 + 1] > coords[j*2 + 1])
253 l_swap = true;
254 }
255 else { /* if (dim == 1) */
256 if (coords[j_prev] > coords[j])
257 l_swap = true;
258 }
259
260 if (l_swap) {
261 cs_lnum_t o_save = order[i-1];
262 order[i-1] = order[i];
263 order[i] = o_save;
264 g_swap = true;
265 }
266 }
267
268 } while (g_swap);
269 }
270
271 /*----------------------------------------------------------------------------
272 * Creation of an I/O numbering structure based on coordinates.
273 *
274 * The ordering is based on a Morton code, and it is expected that
275 * entities are unique (i.e. not duplicated on 2 or more ranks).
276 * In the case that 2 entities have a same Morton code, their global
277 * number will be different, but their order is undetermined.
278 *
279 * parameters:
280 * dim <-- spatial dimension
281 * n_entities <-- number of entities considered
282 * coords <-- pointer to entity coordinates (interlaced)
283 * m_code <-- Morton code associated with each entity
284 * order <-> ordering array base on Morton encoding, or
285 * lexicographical coordinate ordering for ties
286 *
287 * returns:
288 * pointer to I/O numbering structure
289 *----------------------------------------------------------------------------*/
290
291 static void
_check_morton_ordering(int dim,size_t n_entities,const cs_coord_t coords[],const fvm_morton_code_t m_code[],cs_lnum_t order[])292 _check_morton_ordering(int dim,
293 size_t n_entities,
294 const cs_coord_t coords[],
295 const fvm_morton_code_t m_code[],
296 cs_lnum_t order[])
297 {
298 size_t i_prev = 0, i = 1;
299
300 if (n_entities == 0)
301 return;
302
303 /* Check ordering; if two entities have the same Morton codes,
304 use lexicographical coordinates ordering to ensure the
305 final order is deterministic. */
306
307 for (i = 1; i < n_entities; i++) {
308
309 size_t j_prev = order[i_prev], j = order[i];
310
311 if ( m_code[j_prev].X[0] != m_code[j].X[0]
312 || m_code[j_prev].X[1] != m_code[j].X[1]
313 || m_code[j_prev].X[2] != m_code[j].X[2]) {
314
315 /* If successive values have the same Morton code,
316 order them lexicographically */
317 if (i_prev < i - 1)
318 _reorder_coords_lexicographic(dim, i_prev, i-1, coords, order);
319
320 i_prev = i;
321 }
322 }
323
324 if (i_prev < n_entities - 1)
325 _reorder_coords_lexicographic(dim, i_prev, n_entities - 1, coords, order);
326 }
327
328 /*----------------------------------------------------------------------------
329 * Maximum local global number associated with an I/O numbering structure.
330 *
331 * This function is to be used for ordered global numberings.
332 *
333 * parameters:
334 * this_io_num <-- pointer to partially initialized I/O numbering structure.
335 *
336 * returns:
337 * maximum global number associated with the I/O numbering
338 *----------------------------------------------------------------------------*/
339
340 static cs_gnum_t
_fvm_io_num_local_max(const fvm_io_num_t * this_io_num)341 _fvm_io_num_local_max(const fvm_io_num_t *this_io_num)
342 {
343 cs_gnum_t local_max;
344
345 /* Get maximum global number value */
346
347 size_t n_ent = this_io_num->global_num_size;
348 if (n_ent > 0)
349 local_max = this_io_num->global_num[n_ent - 1];
350 else
351 local_max = 0;
352
353 return local_max;
354 }
355
356 /*----------------------------------------------------------------------------
357 * When sub-entities have been added to an I/O numbering, switch from
358 * a numbering on the initial entities (shifted by number of sub-entities) to
359 * a numbering on the final sub-entities.
360 *
361 * Also update whether the numbering may be shared, and discard copy if not
362 * needed.
363 *
364 * parameters:
365 * this_io_num <-> pointer to structure that should be ordered
366 * n_sub_entities <-- optional number of sub-entities per initial entity,
367 * or NULL if unused
368 * may_be_shared <-- indicate if structure may be shared at this stage
369 *----------------------------------------------------------------------------*/
370
371 static void
_fvm_io_num_order_finalize(fvm_io_num_t * this_io_num,const cs_lnum_t n_sub_entities[],bool may_be_shared)372 _fvm_io_num_order_finalize(fvm_io_num_t *this_io_num,
373 const cs_lnum_t n_sub_entities[],
374 bool may_be_shared)
375 {
376 if (n_sub_entities != NULL) {
377
378 cs_lnum_t i, j, k;
379 cs_gnum_t *_global_num;
380
381 for (i = 0, j = 0; i < this_io_num->global_num_size; i++)
382 j += n_sub_entities[i];
383
384 BFT_MALLOC(_global_num, j, cs_gnum_t);
385
386 for (i = 0, j = 0; i < this_io_num->global_num_size; i++) {
387 for (k = 0; k < n_sub_entities[i]; j++, k++)
388 _global_num[j] = this_io_num->_global_num[i] - n_sub_entities[i] + k + 1;
389 }
390
391 BFT_FREE(this_io_num->_global_num);
392 this_io_num->_global_num = _global_num;
393
394 if (this_io_num->global_num_size != (cs_lnum_t)j) {
395 this_io_num->global_num_size = j;
396 may_be_shared = false;
397 }
398
399 if (may_be_shared == false)
400 this_io_num->global_num = this_io_num->_global_num;
401 }
402
403 /* If numbering was initially shared, check if it was changed or if it
404 may remain shared (in which case the copy may be discarded) */
405
406 if (may_be_shared == true) {
407 cs_lnum_t i;
408 for (i = 0; i < this_io_num->global_num_size; i++)
409 if (this_io_num->_global_num[i] != this_io_num->global_num[i])
410 break;
411 if (i < this_io_num->global_num_size)
412 this_io_num->global_num = this_io_num->_global_num;
413 else
414 BFT_FREE(this_io_num->_global_num);
415 }
416 }
417
418 /*----------------------------------------------------------------------------
419 * Local ordering associated with an I/O numbering structure.
420 *
421 * The structure should contain an initial ordering, which should
422 * be sorted, but need not be contiguous. On output, the numbering
423 * will be contiguous.
424 *
425 * As an option, a number of sub-entities per initial entity may be
426 * given, in which case sub-entities of a same entity will have contiguous
427 * numbers in the final ordering.
428 *
429 * parameters:
430 * this_io_num <-> pointer to structure that should be ordered
431 * n_sub_entities <-- optional number of sub-entities per initial entity,
432 * or NULL if unused
433 *----------------------------------------------------------------------------*/
434
435 static void
_fvm_io_num_local_order(fvm_io_num_t * this_io_num,const cs_lnum_t n_sub_entities[])436 _fvm_io_num_local_order(fvm_io_num_t *this_io_num,
437 const cs_lnum_t n_sub_entities[])
438 {
439 cs_gnum_t num_prev, num_cur;
440
441 bool may_be_shared = false;
442
443 cs_gnum_t current_gnum = 0;
444
445 /* If numbering is shared, we will check later it was changed or if
446 it can remain shared (in which case the copy may be discarded) */
447
448 if (this_io_num->global_num != this_io_num->_global_num)
449 may_be_shared = true;
450 else
451 may_be_shared = false;
452
453 size_t n_ent = this_io_num->global_num_size;
454
455 if (n_ent > 0) {
456
457 cs_lnum_t *b_order;
458 cs_gnum_t *b_gnum = this_io_num->_global_num;
459 const cs_lnum_t *b_nsub = n_sub_entities;
460
461 BFT_MALLOC(b_order, n_ent, cs_lnum_t);
462
463 cs_order_gnum_allocated(NULL,
464 b_gnum,
465 b_order,
466 n_ent);
467
468 /* Determine global order; requires ordering to loop through buffer by
469 increasing number (blocks associated with each process are
470 already sorted, but the whole "gathered" block is not).
471 We build an initial global order based on the initial global numbering,
472 such that for each block, the global number of an entity is equal to
473 the cumulative number of sub-entities */
474
475 if (b_nsub != NULL) {
476
477 current_gnum = b_nsub[b_order[0]];
478 num_prev = b_gnum[b_order[0]];
479 b_gnum[b_order[0]] = current_gnum;
480
481 for (size_t i = 1; i < n_ent; i++) {
482 num_cur = b_gnum[b_order[i]];
483 if (num_cur > num_prev)
484 current_gnum += b_nsub[b_order[i]];
485 b_gnum[b_order[i]] = current_gnum;
486 num_prev = num_cur;
487 }
488
489 }
490 else { /* if (b_n_sub == NULL) */
491
492 current_gnum = 1;
493 num_prev = b_gnum[b_order[0]];
494 b_gnum[b_order[0]] = current_gnum;
495
496 for (size_t i = 1; i < n_ent; i++) {
497 num_cur = b_gnum[b_order[i]];
498 if (num_cur > num_prev)
499 current_gnum += 1;
500 b_gnum[b_order[i]] = current_gnum;
501 num_prev = num_cur;
502 }
503
504 }
505
506 BFT_FREE(b_order);
507
508 }
509
510 /* When sub-entities have been added, now switch from a numbering on
511 the initial entities (shifted by number of sub-entities) to
512 a numbering on the final sub-entities */
513
514 _fvm_io_num_order_finalize(this_io_num,
515 n_sub_entities,
516 may_be_shared);
517
518 /* Get final maximum global number value */
519
520 this_io_num->global_count = _fvm_io_num_local_max(this_io_num);
521 }
522
523 /*----------------------------------------------------------------------------
524 * Copy selected shared global ordering information to private ordering
525 * information for an I/O numbering structure.
526 *
527 * parameters:
528 * this_io_num <-- pointer to numbering structure
529 *----------------------------------------------------------------------------*/
530
531 static void
_fvm_io_num_copy_on_write(fvm_io_num_t * const this_io_num)532 _fvm_io_num_copy_on_write(fvm_io_num_t *const this_io_num)
533 {
534 if (this_io_num->_global_num == NULL) {
535 cs_lnum_t i;
536 BFT_MALLOC(this_io_num->_global_num,
537 this_io_num->global_num_size,
538 cs_gnum_t);
539 for (i = 0; i < this_io_num->global_num_size; i++)
540 this_io_num->_global_num[i] = this_io_num->global_num[i];
541 this_io_num->global_num = this_io_num->_global_num;
542 }
543 assert(this_io_num->global_num == this_io_num->_global_num);
544 }
545
546 /*----------------------------------------------------------------------------
547 * Copy selected shared global ordering information to private ordering
548 * information for an I/O numbering structure.
549 *
550 * parameters:
551 * this_io_num <-> pointer to numbering structure
552 * parent_global_number <-- pointer to shared list of global parent
553 * entity numbers
554 *----------------------------------------------------------------------------*/
555
556 static void
_fvm_io_num_try_to_set_shared(fvm_io_num_t * const this_io_num,const cs_gnum_t parent_global_number[])557 _fvm_io_num_try_to_set_shared(fvm_io_num_t *const this_io_num,
558 const cs_gnum_t parent_global_number[])
559 {
560 if (this_io_num->_global_num != NULL && parent_global_number != NULL) {
561 cs_lnum_t i;
562 for (i = 0; i < this_io_num->global_num_size; i++)
563 if (this_io_num->_global_num[i] != parent_global_number[i])
564 break;
565 if (i < this_io_num->global_num_size)
566 this_io_num->global_num = this_io_num->_global_num;
567 else {
568 this_io_num->global_num = parent_global_number;
569 BFT_FREE(this_io_num->_global_num);
570 }
571 }
572 }
573
574 #if defined(HAVE_MPI)
575
576 /*----------------------------------------------------------------------------
577 * Maximum global number associated with an I/O numbering structure.
578 *
579 * This function is to be used for ordered global numberings.
580 *
581 * parameters:
582 * this_io_num <-- pointer to partially initialized I/O numbering structure.
583 * comm <-- associated MPI communicator
584 *
585 * returns:
586 * maximum global number associated with the I/O numbering
587 *----------------------------------------------------------------------------*/
588
589 static cs_gnum_t
_fvm_io_num_global_max(const fvm_io_num_t * const this_io_num,const MPI_Comm comm)590 _fvm_io_num_global_max(const fvm_io_num_t *const this_io_num,
591 const MPI_Comm comm)
592 {
593 cs_gnum_t local_max = _fvm_io_num_local_max(this_io_num);
594 cs_gnum_t global_max = 0;
595
596 MPI_Allreduce(&local_max, &global_max, 1, CS_MPI_GNUM, MPI_MAX, comm);
597
598 return global_max;
599 }
600
601 /*----------------------------------------------------------------------------
602 * Maximum global number associated with an I/O numbering structure.
603 *
604 * This function may be used for ordered global numberings.
605 *
606 * parameters:
607 * this_io_num <-- pointer to partially initialized I/O numbering structure.
608 * comm <-- associated MPI communicator
609 *
610 * returns:
611 * maximum global number associated with the I/O numbering
612 *----------------------------------------------------------------------------*/
613
614 static cs_gnum_t
_fvm_io_num_global_max_unordered(const fvm_io_num_t * const this_io_num,const MPI_Comm comm)615 _fvm_io_num_global_max_unordered(const fvm_io_num_t *const this_io_num,
616 const MPI_Comm comm)
617 {
618 size_t i, n_ent;
619 cs_gnum_t local_max = 0, global_max = 0;;
620
621 /* Get maximum global number value */
622
623 n_ent = this_io_num->global_num_size;
624 for (i = 0; i < n_ent; i++) {
625 cs_gnum_t local_val = this_io_num->global_num[i];
626 if (local_val > local_max)
627 local_max = local_val;
628 }
629
630 MPI_Allreduce(&local_max, &global_max, 1, CS_MPI_GNUM, MPI_MAX, comm);
631
632 return global_max;
633 }
634
635 /*----------------------------------------------------------------------------
636 * Global ordering associated with an I/O numbering structure.
637 *
638 * The structure should contain an initial ordering, which should
639 * be sorted, but need not be contiguous. On output, the numbering
640 * will be contiguous.
641 *
642 * As an option, a number of sub-entities per initial entity may be
643 * given, in which case sub-entities of a same entity will have contiguous
644 * numbers in the final ordering.
645 *
646 * parameters:
647 * this_io_num <-> pointer to structure that should be ordered
648 * n_sub_entities <-- optional number of sub-entities per initial entity,
649 * or NULL if unused
650 * comm <-- associated MPI communicator
651 *----------------------------------------------------------------------------*/
652
653 static void
_fvm_io_num_global_order(fvm_io_num_t * this_io_num,const cs_lnum_t n_sub_entities[],MPI_Comm comm)654 _fvm_io_num_global_order(fvm_io_num_t *this_io_num,
655 const cs_lnum_t n_sub_entities[],
656 MPI_Comm comm)
657 {
658 cs_gnum_t num_prev, num_cur;
659
660 bool may_be_shared = false;
661
662 cs_lnum_t *b_nsub = NULL;
663 int have_sub_loc = 0, have_sub_glob = 0;
664
665 int local_rank, size;
666
667 cs_gnum_t current_gnum = 0, gnum_shift = 0;
668
669 /* Initialization */
670
671 MPI_Comm_rank(comm, &local_rank);
672 MPI_Comm_size(comm, &size);
673
674 num_prev = 0; /* true initialization later for block 0, */
675
676 /* If numbering is shared, we will check later it was changed or if
677 it can remain shared (in which case the copy may be discarded) */
678
679 if (this_io_num->global_num != this_io_num->_global_num)
680 may_be_shared = true;
681 else
682 may_be_shared = false;
683
684 /* Get temporary maximum global number value */
685
686 this_io_num->global_count = _fvm_io_num_global_max(this_io_num, comm);
687
688 cs_block_dist_info_t
689 bi = cs_block_dist_compute_sizes(local_rank,
690 size,
691 1,
692 0,
693 this_io_num->global_count);
694
695 cs_all_to_all_t
696 *d = cs_all_to_all_create_from_block(this_io_num->global_num_size,
697 0, /* flags */
698 this_io_num->global_num,
699 bi,
700 comm);
701
702 cs_gnum_t *b_gnum = cs_all_to_all_copy_array(d,
703 CS_GNUM_TYPE,
704 1,
705 false, /* reverse */
706 this_io_num->global_num,
707 NULL);
708
709 cs_lnum_t b_size = cs_all_to_all_n_elts_dest(d);
710
711 /* Do we have sub-entities ? */
712
713 if (n_sub_entities != NULL)
714 have_sub_loc = 1;
715
716 MPI_Allreduce(&have_sub_loc, &have_sub_glob, 1, MPI_INT, MPI_MAX, comm);
717
718 if (have_sub_glob > 0)
719 b_nsub = cs_all_to_all_copy_array(d,
720 CS_LNUM_TYPE,
721 1,
722 false, /* reverse */
723 n_sub_entities,
724 NULL);
725
726 if (b_size > 0) {
727
728 cs_lnum_t *b_order;
729
730 BFT_MALLOC(b_order, b_size, cs_lnum_t);
731
732 cs_order_gnum_allocated(NULL,
733 b_gnum,
734 b_order,
735 b_size);
736
737 /* Determine global order; requires ordering to loop through buffer by
738 increasing number (blocks associated with each process are
739 already sorted, but the whole "gathered" block is not).
740 We build an initial global order based on the initial global numbering,
741 such that for each block, the global number of an entity is equal to
742 the cumulative number of sub-entities */
743
744 if (have_sub_glob > 0) {
745
746 current_gnum = b_nsub[b_order[0]];
747 num_prev = b_gnum[b_order[0]];
748 b_gnum[b_order[0]] = current_gnum;
749
750 for (cs_lnum_t i = 1; i < b_size; i++) {
751 num_cur = b_gnum[b_order[i]];
752 if (num_cur > num_prev)
753 current_gnum += b_nsub[b_order[i]];
754 b_gnum[b_order[i]] = current_gnum;
755 num_prev = num_cur;
756 }
757
758 }
759 else { /* if (have_sub_glob == 0) */
760
761 current_gnum = 1;
762 num_prev = b_gnum[b_order[0]];
763 b_gnum[b_order[0]] = current_gnum;
764
765 for (cs_lnum_t i = 1; i < b_size; i++) {
766 num_cur = b_gnum[b_order[i]];
767 if (num_cur > num_prev)
768 current_gnum += 1;
769 b_gnum[b_order[i]] = current_gnum;
770 num_prev = num_cur;
771 }
772
773 }
774
775 BFT_FREE(b_order);
776
777 }
778
779 /* Partial clean-up */
780
781 BFT_FREE(b_nsub);
782
783 /* At this stage, b_gnum[] is valid for this process, and
784 current_gnum indicates the total number of entities handled
785 by this process; we must now shift global numberings on different
786 processes by the cumulative total number of entities handled by
787 each process */
788
789 MPI_Scan(¤t_gnum, &gnum_shift, 1, CS_MPI_GNUM,
790 MPI_SUM, comm);
791 gnum_shift -= current_gnum;
792
793 for (cs_lnum_t i = 0; i < b_size; i++)
794 b_gnum[i] += gnum_shift;
795
796 /* Return global order to all ranks */
797
798 cs_all_to_all_copy_array(d,
799 CS_GNUM_TYPE,
800 1,
801 true, /* reverse */
802 b_gnum,
803 this_io_num->_global_num);
804
805 /* Free memory */
806
807 BFT_FREE(b_gnum);
808
809 cs_all_to_all_destroy(&d);
810
811 /* When sub-entities have been added, now switch from a numbering on
812 the initial entities (shifted by number of sub-entities) to
813 a numbering on the final sub-entities */
814
815 _fvm_io_num_order_finalize(this_io_num,
816 n_sub_entities,
817 may_be_shared);
818
819 /* Get final maximum global number value */
820
821 this_io_num->global_count = _fvm_io_num_global_max(this_io_num, comm);
822 }
823
824 /*----------------------------------------------------------------------------
825 * Global ordering associated with an I/O numbering structure.
826 *
827 * The structure does not need to contain an initial ordering,
828 * though the array should be allocated.
829 *
830 * parameters:
831 * this_io_num <-> pointer to structure that should be ordered
832 * stride <-- values per entity
833 * global_num <-- global numbering array
834 * comm <-- associated MPI communicator
835 *----------------------------------------------------------------------------*/
836
837 static void
_fvm_io_num_global_order_s(fvm_io_num_t * this_io_num,size_t stride,cs_gnum_t global_num[],MPI_Comm comm)838 _fvm_io_num_global_order_s(fvm_io_num_t *this_io_num,
839 size_t stride,
840 cs_gnum_t global_num[],
841 MPI_Comm comm)
842 {
843 int local_rank, size;
844 cs_gnum_t current_gnum = 0, gnum_shift = 0;
845
846 cs_gnum_t *r_gnum = NULL;
847
848 /* Initialization */
849
850 MPI_Comm_rank(comm, &local_rank);
851 MPI_Comm_size(comm, &size);
852
853 /* Get maximum global number value for first value of each series
854 (does not need to be exact, simply used to define blocks) */
855
856 {
857 cs_gnum_t local_max = 0, global_max = 0;
858 size_t n_ent = this_io_num->global_num_size;
859
860 if (n_ent > 0)
861 local_max = global_num[(n_ent-1)*stride];
862 MPI_Allreduce(&local_max, &global_max, 1, CS_MPI_GNUM, MPI_MAX, comm);
863 this_io_num->global_count = global_max;
864 }
865
866 /* block_info */
867
868 cs_block_dist_info_t
869 bi = cs_block_dist_compute_sizes(local_rank,
870 size,
871 1,
872 0,
873 this_io_num->global_count);
874
875 const cs_gnum_t block_size = bi.block_size;
876
877 int *dest_rank;
878 BFT_MALLOC(dest_rank, this_io_num->global_num_size, int);
879
880 for (cs_lnum_t i = 0; i < this_io_num->global_num_size; i++) {
881 cs_gnum_t g_elt_id = global_num[stride*i] - 1;
882 dest_rank[i] = (g_elt_id / block_size)*bi.rank_step;
883 }
884
885 cs_all_to_all_t
886 *d = cs_all_to_all_create(this_io_num->global_num_size,
887 0, /* flags */
888 NULL, /* dest_id */
889 dest_rank,
890 comm);
891
892 cs_all_to_all_transfer_dest_rank(d, &dest_rank);
893
894 cs_gnum_t *b_gnum = cs_all_to_all_copy_array(d,
895 CS_GNUM_TYPE,
896 stride,
897 false, /* reverse */
898 global_num,
899 NULL);
900
901 cs_lnum_t b_size = cs_all_to_all_n_elts_dest(d);
902
903 /* Order received data based on global number sets */
904
905 if (b_size > 0) {
906
907 cs_lnum_t *b_order = NULL;
908
909 BFT_MALLOC(r_gnum, b_size, cs_gnum_t);
910 BFT_MALLOC(b_order, b_size, cs_lnum_t);
911
912 cs_order_gnum_allocated_s(NULL,
913 b_gnum,
914 stride,
915 b_order,
916 b_size);
917
918 /* We build an initial global order based on the initial global numbering,
919 such that for each block, the global number of an entity is equal to
920 the cumulative number of sub-entities */
921
922 current_gnum = 1;
923 cs_lnum_t prev_id = b_order[0];
924 r_gnum[b_order[0]] = current_gnum;
925
926 const cs_lnum_t _stride = stride;
927
928 for (cs_lnum_t i = 1; i < b_size; i++) {
929 bool greater_than_prev = false;
930 cs_lnum_t cur_id = b_order[i];
931 for (cs_lnum_t j = 0; j < _stride; j++) {
932 if ( b_gnum[cur_id*_stride + j]
933 > b_gnum[prev_id*_stride + j])
934 greater_than_prev = true;
935 }
936 if (greater_than_prev)
937 current_gnum += 1;
938 r_gnum[b_order[i]] = current_gnum;
939 prev_id = cur_id;
940 }
941
942 BFT_FREE(b_order);
943
944 }
945
946 BFT_FREE(b_gnum);
947
948 /* At this stage, r_gnum[] is valid for this process, and
949 current_gnum indicates the total number of entities handled
950 by this process; we must now shift global numberings on different
951 processes by the cumulative total number of entities handled by
952 each process */
953
954 MPI_Scan(¤t_gnum, &gnum_shift, 1, CS_MPI_GNUM,
955 MPI_SUM, comm);
956 gnum_shift -= current_gnum;
957
958 for (cs_lnum_t i = 0; i < b_size; i++)
959 r_gnum[i] += gnum_shift;
960
961 /* Return global order to all ranks */
962
963 cs_all_to_all_copy_array(d,
964 CS_GNUM_TYPE,
965 1,
966 true, /* reverse */
967 r_gnum,
968 this_io_num->_global_num);
969
970 /* Partial clean-up */
971
972 BFT_FREE(r_gnum);
973
974 cs_all_to_all_destroy(&d);
975
976 /* Get final maximum global number value */
977
978 this_io_num->global_count = _fvm_io_num_global_max(this_io_num, comm);
979 }
980
981 /*----------------------------------------------------------------------------
982 * Compare two elements in an indexed list and returns true if element in
983 * position i1 is strictly greater than element in position i2.
984 *
985 * parameters:
986 * i1 <-- position in index for the first element
987 * i2 <-- position in index for the second element
988 * index <-- number of values to compare for each entity
989 * number <-- pointer to numbers of entities that should be ordered.
990 * (if NULL, a default 1 to n numbering is considered)
991 *
992 * returns:
993 * true or false
994 *----------------------------------------------------------------------------*/
995
996 inline static bool
_indexed_is_greater(size_t i1,size_t i2,const cs_lnum_t index[],const cs_gnum_t number[])997 _indexed_is_greater(size_t i1,
998 size_t i2,
999 const cs_lnum_t index[],
1000 const cs_gnum_t number[])
1001 {
1002 int i;
1003
1004 cs_lnum_t i1_s = index[i1], i1_e = index[i1+1], s1 = i1_e - i1_s;
1005 cs_lnum_t i2_s = index[i2], i2_e = index[i2+1], s2 = i2_e - i2_s;
1006
1007 if (s1 > s2) {
1008
1009 for (i = 0; i < s2; i++) {
1010 if (number[i1_s + i] > number[i2_s + i])
1011 return true;
1012 else if (number[i1_s + i] < number[i2_s + i])
1013 return false;
1014 }
1015
1016 return true;
1017 }
1018 else { /* s1 <= s2 */
1019
1020 for (i = 0; i < s1; i++) {
1021 if (number[i1_s + i] > number[i2_s + i])
1022 return true;
1023 else if (number[i1_s + i] < number[i2_s + i])
1024 return false;
1025 }
1026
1027 return false;
1028 }
1029
1030 }
1031
1032 /*----------------------------------------------------------------------------
1033 * Global indexed ordering associated with an I/O numbering structure.
1034 *
1035 * The structure should contain an initial ordering, which should
1036 * be sorted, but need not be contiguous. On output, the numbering
1037 * will be contiguous.
1038 *
1039 * parameters:
1040 * this_io_num <-> pointer to structure that should be ordered
1041 * index <-- index on entities for global_num[]
1042 * global_num <-- global number associated with each entity
1043 * comm <-- associated MPI communicator
1044 *----------------------------------------------------------------------------*/
1045
1046 static void
_fvm_io_num_global_order_index(fvm_io_num_t * this_io_num,cs_lnum_t index[],cs_gnum_t global_num[],MPI_Comm comm)1047 _fvm_io_num_global_order_index(fvm_io_num_t *this_io_num,
1048 cs_lnum_t index[],
1049 cs_gnum_t global_num[],
1050 MPI_Comm comm)
1051 {
1052 int local_rank, size;
1053
1054 cs_gnum_t current_global_num = 0, global_num_shift = 0;
1055 cs_gnum_t *block_global_num = NULL, *recv_global_num = NULL;
1056
1057 size_t n_ent = this_io_num->global_num_size;
1058
1059 /* Initialization */
1060
1061 MPI_Comm_rank(comm, &local_rank);
1062 MPI_Comm_size(comm, &size);
1063
1064 /* Get maximum global number value for first value of each series
1065 (does not need to be exact, simply used to define blocks) */
1066
1067 {
1068 cs_gnum_t local_max = 0, global_max = 0;
1069
1070 if (n_ent > 0)
1071 local_max = global_num[index[n_ent-1]];
1072 MPI_Allreduce(&local_max, &global_max, 1, CS_MPI_GNUM, MPI_MAX, comm);
1073 this_io_num->global_count = global_max;
1074 }
1075
1076 /* block_size = ceil(this_io_num->global_count/size) */
1077
1078 cs_block_dist_info_t
1079 bi = cs_block_dist_compute_sizes(local_rank,
1080 size,
1081 1,
1082 0,
1083 this_io_num->global_count);
1084
1085 cs_gnum_t _block_size = bi.block_size;
1086
1087 /* Build for each block, a new ordered indexed list from the received
1088 elements */
1089
1090 int *dest_rank;
1091 BFT_MALLOC(dest_rank, this_io_num->global_num_size, int);
1092 for (size_t i = 0; i < n_ent; i++)
1093 dest_rank[i] = ((global_num[index[i]] - 1) / _block_size) * bi.rank_step;
1094
1095 int flags = CS_ALL_TO_ALL_ORDER_BY_SRC_RANK;
1096
1097 cs_all_to_all_t *d = cs_all_to_all_create(n_ent,
1098 flags,
1099 NULL,
1100 dest_rank,
1101 comm);
1102 cs_all_to_all_transfer_dest_rank(d, &dest_rank);
1103
1104 /* Get recv_index */
1105
1106 cs_lnum_t *recv_index
1107 = cs_all_to_all_copy_index(d,
1108 false, /* reverse */
1109 index,
1110 NULL);
1111
1112 cs_lnum_t n_ent_recv = cs_all_to_all_n_elts_dest(d);
1113
1114 /* Send indexed list to the destination ranks;
1115 TODO future optimization: when data is sorted by increasing base global
1116 numbering, such as here, indicating it to the cs_all_to_all distributor
1117 for possible optimisation (using a flag or a varianf function) could
1118 avoid an extra copy to an internal buffer. */
1119
1120 recv_global_num = cs_all_to_all_copy_indexed(d,
1121 CS_GNUM_TYPE,
1122 false, /* reverse */
1123 index,
1124 global_num,
1125 recv_index,
1126 NULL);
1127
1128 if (n_ent_recv > 0) { /* Order received elements of the indexed list */
1129
1130 size_t prev_id, cur_id;
1131
1132 cs_lnum_t *recv_order = NULL;
1133 BFT_MALLOC(recv_order, n_ent_recv, cs_lnum_t);
1134
1135 cs_order_gnum_allocated_i(NULL,
1136 recv_global_num,
1137 recv_index,
1138 recv_order,
1139 n_ent_recv);
1140
1141 /* Determine global order; requires ordering to loop through buffer by
1142 increasing number (blocks associated with each process are
1143 already sorted, but the whole "gathered" block is not).
1144 We build an initial global order based on the initial global numbering,
1145 such that for each block, the global number of an entity is equal to
1146 the cumulative number of elements */
1147
1148 BFT_MALLOC(block_global_num, n_ent_recv, cs_gnum_t);
1149
1150 current_global_num = 1;
1151 prev_id = recv_order[0];
1152 block_global_num[recv_order[0]] = current_global_num;
1153
1154 for (cs_lnum_t i = 1; i < n_ent_recv; i++) {
1155
1156 cur_id = recv_order[i];
1157
1158 if (_indexed_is_greater(cur_id, prev_id, recv_index, recv_global_num))
1159 current_global_num += 1;
1160
1161 block_global_num[recv_order[i]] = current_global_num;
1162 prev_id = cur_id;
1163
1164 }
1165
1166 BFT_FREE(recv_order);
1167
1168 } /* End if n_ent_recv > 0 */
1169
1170 /* Partial clean-up */
1171
1172 BFT_FREE(recv_index);
1173 BFT_FREE(recv_global_num);
1174
1175 /* At this stage, block_global_num[] is valid for this rank, and
1176 current_global_num indicates the total number of entities handled
1177 by this rank; we must now shift global numberings on different
1178 ranks by the cumulative total number of entities handled by
1179 each rank */
1180
1181 MPI_Scan(¤t_global_num, &global_num_shift, 1, CS_MPI_GNUM,
1182 MPI_SUM, comm);
1183 global_num_shift -= current_global_num;
1184
1185 for (cs_lnum_t i = 0; i < n_ent_recv; i++)
1186 block_global_num[i] += global_num_shift;
1187
1188 /* Return global order to all ranks */
1189
1190 cs_all_to_all_copy_array(d,
1191 CS_GNUM_TYPE,
1192 1,
1193 true, /* reverse */
1194 block_global_num,
1195 this_io_num->_global_num);
1196
1197 /* Free memory */
1198
1199 BFT_FREE(block_global_num);
1200
1201 cs_all_to_all_destroy(&d);
1202
1203 /* Get final maximum global number value */
1204
1205 this_io_num->global_count = _fvm_io_num_global_max(this_io_num, comm);
1206 }
1207
1208 /*----------------------------------------------------------------------------
1209 * Return the global number of sub-entities associated with an initial
1210 * entity whose global numbering is known, given the number of
1211 * sub-entities per initial entity.
1212 *
1213 * parameters:
1214 * this_io_num <-- pointer to base io numbering
1215 * n_sub_entities <-- number of sub-entities per initial entity
1216 * comm <-- associated MPI communicator
1217 *
1218 * returns:
1219 * global number of sub-entities
1220 *----------------------------------------------------------------------------*/
1221
1222 static cs_gnum_t
_fvm_io_num_global_sub_size(const fvm_io_num_t * this_io_num,const cs_lnum_t n_sub_entities[],MPI_Comm comm)1223 _fvm_io_num_global_sub_size(const fvm_io_num_t *this_io_num,
1224 const cs_lnum_t n_sub_entities[],
1225 MPI_Comm comm)
1226 {
1227 cs_gnum_t global_count, num_prev, num_cur;
1228
1229 cs_gnum_t *send_global_num = NULL;
1230 cs_lnum_t *recv_n_sub = NULL, *recv_order = NULL;
1231 int have_sub_loc = 0, have_sub_glob = 0;
1232
1233 int size, local_rank;
1234
1235 cs_gnum_t current_global_num = 0;
1236 cs_gnum_t retval = 0;
1237
1238 /* Initialization */
1239
1240 MPI_Comm_size(comm, &size);
1241 MPI_Comm_rank(comm, &local_rank);
1242
1243 num_prev = 0; /* true initialization later for block 0, */
1244
1245 /* Get temporary maximum global number value */
1246
1247 global_count = _fvm_io_num_global_max(this_io_num, comm);
1248
1249 /* block_size = ceil(this_io_num->global_count/size) */
1250
1251 cs_block_dist_info_t
1252 bi = cs_block_dist_compute_sizes(local_rank,
1253 size,
1254 1,
1255 0,
1256 global_count);
1257
1258 cs_all_to_all_t
1259 *d = cs_all_to_all_create_from_block(this_io_num->global_num_size,
1260 0, /* flags */
1261 this_io_num->global_num,
1262 bi,
1263 comm);
1264
1265 /* As data is sorted by increasing base global numbering, we do not
1266 need to build an extra array, but only to send the correct parts
1267 of the n_sub_entities[] array to the correct processors */
1268
1269 if (this_io_num->_global_num != NULL)
1270 send_global_num = this_io_num->_global_num;
1271 else {
1272 BFT_MALLOC(send_global_num,
1273 this_io_num->global_num_size,
1274 cs_gnum_t);
1275 memcpy(send_global_num,
1276 this_io_num->global_num,
1277 this_io_num->global_num_size * sizeof(cs_gnum_t));
1278 }
1279
1280 cs_gnum_t *recv_global_num = cs_all_to_all_copy_array(d,
1281 CS_GNUM_TYPE,
1282 1,
1283 false, /* reverse */
1284 send_global_num,
1285 NULL);
1286
1287 cs_lnum_t n_ent_recv = cs_all_to_all_n_elts_dest(d);
1288
1289 BFT_MALLOC(recv_order, n_ent_recv, cs_lnum_t);
1290
1291 if (send_global_num != this_io_num->_global_num)
1292 BFT_FREE(send_global_num);
1293
1294 /* Do we have sub-entities ? */
1295
1296 if (n_sub_entities != NULL)
1297 have_sub_loc = 1;
1298
1299 MPI_Allreduce(&have_sub_loc, &have_sub_glob, 1, MPI_INT, MPI_MAX, comm);
1300
1301 if (have_sub_glob > 0) {
1302
1303 cs_lnum_t *send_n_sub;
1304 BFT_MALLOC(send_n_sub, this_io_num->global_num_size, cs_lnum_t);
1305
1306 if (n_sub_entities != NULL) {
1307 for (cs_lnum_t i = 0; i < this_io_num->global_num_size; i++)
1308 send_n_sub[i] = n_sub_entities[i];
1309 }
1310 else {
1311 for (cs_lnum_t i = 0; i < this_io_num->global_num_size; i++)
1312 send_n_sub[i] = 1;
1313 }
1314
1315 recv_n_sub = cs_all_to_all_copy_array(d,
1316 CS_LNUM_TYPE,
1317 1,
1318 false, /* reverse */
1319 send_n_sub,
1320 NULL);
1321
1322 BFT_FREE(send_n_sub);
1323 }
1324
1325 if (n_ent_recv > 0) {
1326
1327 cs_order_gnum_allocated(NULL,
1328 recv_global_num,
1329 recv_order,
1330 n_ent_recv);
1331
1332 /* Determine global order; requires ordering to loop through buffer by
1333 increasing number (blocks associated with each process are
1334 already sorted, but the whole "gathered" block is not).
1335 We build an initial global order based on the initial global numbering,
1336 such that for each block, the global number of an entity is equal to
1337 the cumulative number of sub-entities */
1338
1339 current_global_num = recv_n_sub[recv_order[0]];
1340 num_prev = recv_global_num[recv_order[0]];
1341 recv_global_num[recv_order[0]] = current_global_num;
1342
1343 for (cs_lnum_t i = 1; i < n_ent_recv; i++) {
1344 num_cur = recv_global_num[recv_order[i]];
1345 if (num_cur > num_prev)
1346 current_global_num += recv_n_sub[recv_order[i]];
1347 num_prev = num_cur;
1348 }
1349
1350 }
1351
1352 /* Partial clean-up */
1353
1354 BFT_FREE(recv_n_sub);
1355 BFT_FREE(recv_order);
1356 BFT_FREE(recv_global_num);
1357
1358 cs_all_to_all_destroy(&d);
1359
1360 /* At this stage, current_global_num indicates the total number of
1361 entities handled by this process; we must now shift global
1362 numberings on different processes by the cumulative total
1363 number of entities handled by each process */
1364
1365 MPI_Allreduce(¤t_global_num, &retval, 1, CS_MPI_GNUM, MPI_SUM, comm);
1366
1367 return retval;
1368 }
1369
1370 #endif /* defined(HAVE_MPI) */
1371
1372 /*----------------------------------------------------------------------------
1373 * Stretch extents in some directions.
1374 *
1375 * If the multiplication factor for a given axis direction is less than 1,
1376 * the extents in that direction will be defined so that the extent box's
1377 * size in that direction is its maximum size.
1378 *
1379 * parameters:
1380 * extents <-> box extents
1381 * box_to_cube <-- if 1, transform bounding box to boundign cube
1382 *----------------------------------------------------------------------------*/
1383
1384 static void
_adjust_extents(cs_coord_t extents[6],int box_to_cube)1385 _adjust_extents(cs_coord_t extents[6],
1386 int box_to_cube)
1387 {
1388 size_t i;
1389 cs_coord_t max_width = 0.;
1390 const double epsilon = 1e-12;
1391
1392 for (i = 0; i < 3; i++) {
1393 double w = fabs(extents[i+3] - extents[i]);
1394 max_width = CS_MAX(max_width, w);
1395 }
1396
1397 for (i = 0; i < 3; i++) {
1398 double mult = 1.0;
1399 double m = (extents[i] + extents[i+3])*0.5;
1400 double w = fabs(extents[i+3] - extents[i]);
1401 if (box_to_cube > 0) {
1402 if (w > epsilon)
1403 mult = max_width / w;
1404 }
1405 if (mult < 1.0 + epsilon)
1406 mult= 1.0 + epsilon;
1407 extents[i] = m - (w*0.5*mult);
1408 extents[i+3] = m + (w*0.5*mult);
1409 }
1410 }
1411
1412 /*----------------------------------------------------------------------------
1413 * Creation of an I/O numbering structure based on coordinates.
1414 *
1415 * The ordering is based on a Morton code, and it is expected that
1416 * entities are unique (i.e. not duplicated on 2 or more ranks).
1417 * In the case that 2 entities have a same Morton code, their global
1418 * number will be determined by lexicographical ordering of coordinates.
1419 *
1420 * parameters:
1421 * coords <-- pointer to entity coordinates (interlaced)
1422 * dim <-- spatial dimension
1423 * n_entities <-- number of entities considered
1424 * box_to_cube <-- if 1, use bounding cube instead of bounding box
1425 *
1426 * returns:
1427 * pointer to I/O numbering structure
1428 *----------------------------------------------------------------------------*/
1429
1430 static fvm_io_num_t *
_create_from_coords_morton(const cs_coord_t coords[],int dim,size_t n_entities,int box_to_cube)1431 _create_from_coords_morton(const cs_coord_t coords[],
1432 int dim,
1433 size_t n_entities,
1434 int box_to_cube)
1435 {
1436 size_t i;
1437 cs_coord_t extents[6];
1438
1439 #if defined(HAVE_MPI)
1440 MPI_Comm comm = cs_glob_mpi_comm;
1441 #endif
1442
1443 const int level = sizeof(fvm_morton_int_t)*8 - 1;
1444 const int n_ranks = cs_glob_n_ranks;
1445
1446 fvm_io_num_t *this_io_num = NULL;
1447
1448 /* Create structure */
1449
1450 BFT_MALLOC(this_io_num, 1, fvm_io_num_t);
1451
1452 this_io_num->global_num_size = n_entities;
1453
1454 BFT_MALLOC(this_io_num->_global_num, n_entities, cs_gnum_t);
1455 this_io_num->global_num = this_io_num->_global_num;
1456
1457 /* Build Morton encoding and order it */
1458
1459 #if defined(HAVE_MPI)
1460 fvm_morton_get_coord_extents(dim, n_entities, coords, extents, comm);
1461 #else
1462 fvm_morton_get_coord_extents(dim, n_entities, coords, extents);
1463 #endif
1464
1465 _adjust_extents(extents, box_to_cube);
1466
1467 #if defined(HAVE_MPI)
1468
1469 if (n_ranks > 1) {
1470
1471 int *dest_rank = NULL;
1472 cs_lnum_t *order = NULL;
1473 fvm_morton_code_t *m_code = NULL;
1474 int input[1] = {dim};
1475
1476 BFT_MALLOC(m_code, n_entities, fvm_morton_code_t);
1477 BFT_MALLOC(order, n_entities, cs_lnum_t);
1478 BFT_MALLOC(dest_rank, n_entities, int);
1479
1480 fvm_morton_encode_coords(dim, level, extents, n_entities, coords, m_code);
1481 fvm_morton_local_order(n_entities, m_code, order);
1482
1483 cs_sort_partition_dest_rank_id(_sampling_factors[dim],
1484 sizeof(fvm_morton_code_t),
1485 n_entities,
1486 m_code,
1487 NULL, /* weight */
1488 order,
1489 dest_rank,
1490 fvm_morton_s_to_code,
1491 fvm_morton_compare_o,
1492 input,
1493 comm);
1494
1495 BFT_FREE(order);
1496 BFT_FREE(m_code);
1497
1498 cs_all_to_all_t
1499 *d = cs_all_to_all_create(this_io_num->global_num_size,
1500 0, /* flags */
1501 NULL, /* dest_id */
1502 dest_rank,
1503 comm);
1504
1505 cs_all_to_all_transfer_dest_rank(d, &dest_rank);
1506
1507 cs_real_t *b_coords
1508 = cs_all_to_all_copy_array(d,
1509 CS_REAL_TYPE,
1510 3,
1511 false, /* reverse */
1512 coords,
1513 NULL);
1514
1515 size_t b_size = cs_all_to_all_n_elts_dest(d);
1516
1517 /* Now re-build Morton codes on block distribution
1518 (exchanging and ordering Morton codes directly would be more
1519 efficient, but using coordinates allows breaking ties for
1520 possibly identical codes through lexicographical ordering). */
1521
1522 BFT_MALLOC(order, b_size, cs_lnum_t);
1523 BFT_MALLOC(m_code, b_size, fvm_morton_code_t);
1524
1525 fvm_morton_encode_coords(dim,
1526 level,
1527 extents,
1528 b_size,
1529 b_coords,
1530 m_code);
1531
1532 fvm_morton_local_order(b_size, m_code, order);
1533
1534 /* Check ordering; if two entities have the same Morton codes,
1535 use lexicographical coordinates ordering to ensure the
1536 final order is deterministic. */
1537
1538 _check_morton_ordering(dim, b_size, b_coords, m_code, order);
1539
1540 BFT_FREE(m_code);
1541 BFT_FREE(b_coords);
1542
1543 /* Determine global order; requires ordering to loop through buffer by
1544 increasing number (blocks associated with each process are
1545 already sorted, but the whole "gathered" block is not).
1546 We build an initial global order based on the initial global numbering,
1547 such that for each block, the global number of an entity is equal to
1548 the cumulative number of sub-entities */
1549
1550 cs_gnum_t *b_gnum;
1551 BFT_MALLOC(b_gnum, b_size, cs_gnum_t);
1552
1553 for (i = 0; i < b_size; i++)
1554 b_gnum[order[i]] = i+1;
1555
1556 BFT_FREE(order);
1557
1558 cs_gnum_t gnum_shift = 0, current_gnum = b_size;
1559
1560 /* At this stage, b_gnum[] is valid for this process, and
1561 current_gnum indicates the total number of entities handled
1562 by this process; we must now shift global numberings on different
1563 processes by the cumulative total number of entities handled by
1564 each process */
1565
1566 MPI_Scan(¤t_gnum, &gnum_shift, 1, CS_MPI_GNUM,
1567 MPI_SUM, comm);
1568 gnum_shift -= current_gnum;
1569
1570 for (i = 0; i < b_size; i++)
1571 b_gnum[i] += gnum_shift;
1572
1573 /* Return global order to all ranks */
1574
1575 cs_all_to_all_copy_array(d,
1576 CS_GNUM_TYPE,
1577 1,
1578 true, /* reverse */
1579 b_gnum,
1580 this_io_num->_global_num);
1581
1582 /* Free memory */
1583
1584 BFT_FREE(b_gnum);
1585
1586 cs_all_to_all_destroy(&d);
1587
1588 /* Get final maximum global number value */
1589
1590 this_io_num->global_count
1591 = _fvm_io_num_global_max_unordered(this_io_num, comm);
1592
1593 }
1594
1595 #endif /* HAVE_MPI */
1596
1597 if (n_ranks == 1) {
1598
1599 cs_lnum_t *order = NULL;
1600 fvm_morton_code_t *m_code = NULL;
1601
1602 BFT_MALLOC(m_code, n_entities, fvm_morton_code_t);
1603 BFT_MALLOC(order, n_entities, cs_lnum_t);
1604
1605 fvm_morton_encode_coords(dim, level, extents, n_entities, coords, m_code);
1606 fvm_morton_local_order(n_entities, m_code, order);
1607
1608 _check_morton_ordering(dim, n_entities, coords, m_code, order);
1609
1610 BFT_FREE(m_code);
1611
1612 for (i = 0; i < n_entities; i++)
1613 this_io_num->_global_num[order[i]] = i+1;
1614
1615 BFT_FREE(order);
1616
1617 this_io_num->global_count = n_entities;
1618
1619 }
1620
1621 return this_io_num;
1622 }
1623
1624 /*----------------------------------------------------------------------------
1625 * Creation of an I/O numbering structure based on coordinates.
1626 *
1627 * The ordering is based on a Hilbert curve, and it is expected that
1628 * entities are unique (i.e. not duplicated on 2 or more ranks).
1629 * In the case that 2 entities have a same Hilbert code, their global
1630 * number will be determined by lexicographical ordering of coordinates.
1631 *
1632 * parameters:
1633 * coords <-- pointer to entity coordinates (interlaced)
1634 * dim <-- spatial dimension
1635 * n_entities <-- number of entities considered
1636 * box_to_cube <-- if 1, use bounding cube instead of bounding box
1637 *
1638 * returns:
1639 * pointer to I/O numbering structure
1640 *----------------------------------------------------------------------------*/
1641
1642 static fvm_io_num_t *
_create_from_coords_hilbert(const cs_coord_t coords[],int dim,size_t n_entities,int box_to_cube)1643 _create_from_coords_hilbert(const cs_coord_t coords[],
1644 int dim,
1645 size_t n_entities,
1646 int box_to_cube)
1647 {
1648 size_t i;
1649 cs_coord_t extents[6];
1650
1651 #if defined(HAVE_MPI)
1652 MPI_Comm comm = cs_glob_mpi_comm;
1653 #endif
1654
1655 const int n_ranks = cs_glob_n_ranks;
1656
1657 fvm_io_num_t *this_io_num = NULL;
1658
1659 /* Create structure */
1660
1661 BFT_MALLOC(this_io_num, 1, fvm_io_num_t);
1662
1663 this_io_num->global_num_size = n_entities;
1664
1665 BFT_MALLOC(this_io_num->_global_num, n_entities, cs_gnum_t);
1666 this_io_num->global_num = this_io_num->_global_num;
1667
1668 /* Build Hilbert encoding and order it */
1669
1670 #if defined(HAVE_MPI)
1671 fvm_hilbert_get_coord_extents(dim, n_entities, coords, extents, comm);
1672 #else
1673 fvm_hilbert_get_coord_extents(dim, n_entities, coords, extents);
1674 #endif
1675
1676 _adjust_extents(extents, box_to_cube);
1677
1678 #if defined(HAVE_MPI)
1679
1680 if (n_ranks > 1) {
1681
1682 int *dest_rank = NULL;
1683 cs_lnum_t *order = NULL;
1684 fvm_hilbert_code_t *h_code = NULL;
1685
1686 BFT_MALLOC(h_code, n_entities, fvm_hilbert_code_t);
1687 BFT_MALLOC(order, n_entities, cs_lnum_t);
1688 BFT_MALLOC(dest_rank, n_entities, int);
1689
1690 fvm_hilbert_encode_coords(dim, extents, n_entities, coords, h_code);
1691 fvm_hilbert_local_order(n_entities, h_code, order);
1692
1693 cs_sort_partition_dest_rank_id(_sampling_factors[dim],
1694 sizeof(fvm_hilbert_code_t),
1695 n_entities,
1696 h_code,
1697 NULL, /* weight */
1698 order,
1699 dest_rank,
1700 fvm_hilbert_s_to_code,
1701 fvm_hilbert_compare,
1702 NULL,
1703 comm);
1704
1705 BFT_FREE(order);
1706 BFT_FREE(h_code);
1707
1708 cs_all_to_all_t
1709 *d = cs_all_to_all_create(this_io_num->global_num_size,
1710 0, /* flags */
1711 NULL, /* dest_id */
1712 dest_rank,
1713 comm);
1714
1715 cs_all_to_all_transfer_dest_rank(d, &dest_rank);
1716
1717 cs_real_t *b_coords
1718 = cs_all_to_all_copy_array(d,
1719 CS_REAL_TYPE,
1720 3,
1721 false, /* reverse */
1722 coords,
1723 NULL);
1724
1725 size_t b_size = cs_all_to_all_n_elts_dest(d);
1726
1727 /* Now re-build Hilbert coords on block distribution
1728 (exchanging and ordering Hilbert codes directly would be more
1729 efficient, but using coordinates allows breaking ties for
1730 possibly identical codes through lexicographical ordering). */
1731
1732 BFT_MALLOC(order, b_size, cs_lnum_t);
1733
1734 fvm_hilbert_local_order_coords(dim,
1735 extents,
1736 b_size,
1737 b_coords,
1738 order);
1739
1740 BFT_FREE(b_coords);
1741
1742 /* Determine global order; requires ordering to loop through buffer by
1743 increasing number (blocks associated with each process are
1744 already sorted, but the whole "gathered" block is not).
1745 We build an initial global order based on the initial global numbering,
1746 such that for each block, the global number of an entity is equal to
1747 the cumulative number of sub-entities */
1748
1749 cs_gnum_t *b_gnum;
1750 BFT_MALLOC(b_gnum, b_size, cs_gnum_t);
1751
1752 for (i = 0; i < b_size; i++)
1753 b_gnum[order[i]] = i+1;
1754
1755 BFT_FREE(order);
1756
1757 cs_gnum_t gnum_shift = 0, current_gnum = b_size;
1758
1759 /* At this stage, b_gnum[] is valid for this process, and
1760 current_gnum indicates the total number of entities handled
1761 by this process; we must now shift global numberings on different
1762 processes by the cumulative total number of entities handled by
1763 each process */
1764
1765 MPI_Scan(¤t_gnum, &gnum_shift, 1, CS_MPI_GNUM,
1766 MPI_SUM, comm);
1767 gnum_shift -= current_gnum;
1768
1769 for (i = 0; i < b_size; i++)
1770 b_gnum[i] += gnum_shift;
1771
1772 /* Return global order to all ranks */
1773
1774 cs_all_to_all_copy_array(d,
1775 CS_GNUM_TYPE,
1776 1,
1777 true, /* reverse */
1778 b_gnum,
1779 this_io_num->_global_num);
1780
1781 /* Free memory */
1782
1783 BFT_FREE(b_gnum);
1784
1785 cs_all_to_all_destroy(&d);
1786
1787 /* Get final maximum global number value */
1788
1789 this_io_num->global_count
1790 = _fvm_io_num_global_max_unordered(this_io_num, comm);
1791
1792 }
1793
1794 #endif /* HAVE_MPI */
1795
1796 if (n_ranks == 1) {
1797
1798 cs_lnum_t *order = NULL;
1799 BFT_MALLOC(order, n_entities, cs_lnum_t);
1800
1801 fvm_hilbert_local_order_coords(dim,
1802 extents,
1803 n_entities,
1804 coords,
1805 order);
1806
1807 for (i = 0; i < n_entities; i++)
1808 this_io_num->_global_num[order[i]] = i+1;
1809
1810 BFT_FREE(order);
1811
1812 this_io_num->global_count = n_entities;
1813
1814 }
1815
1816 return this_io_num;
1817 }
1818
1819 /*----------------------------------------------------------------------------
1820 * Test if an array of global numbers is lexicographically ordered.
1821 *
1822 * parameters:
1823 * parent_entity_id <-- optional list (0 to n-1 numbering) of selected
1824 * entities (or NULL if all n_entities are selected).
1825 * This list may contain element ids in any order
1826 * number <-- array of all entity numbers (number of entity i given
1827 * by number[i] or number[list[i] - 1]) if list exists
1828 * (if NULL, a default 1 to n numbering is considered)
1829 * stride <-- stride of number array (number of values to compare)
1830 * n_entities <-- number of entities considered
1831 *
1832 * returns:
1833 * 1 if ordered, 0 otherwise.
1834 *----------------------------------------------------------------------------*/
1835
1836 static int
_is_gnum_ordered_s(const cs_lnum_t parent_entity_id[],const cs_gnum_t number[],size_t stride,size_t n_entities)1837 _is_gnum_ordered_s(const cs_lnum_t parent_entity_id[],
1838 const cs_gnum_t number[],
1839 size_t stride,
1840 size_t n_entities)
1841 {
1842 size_t j;
1843 size_t i = 0;
1844
1845 /* If numbering is explicit */
1846
1847 if (number != NULL) {
1848
1849 if (parent_entity_id != NULL) {
1850 for (i = 1 ; i < n_entities ; i++) {
1851 size_t j_prev, k;
1852 bool unordered = false;
1853 j_prev = parent_entity_id[i-1];
1854 j = parent_entity_id[i];
1855 for (k = 0; k < stride; k++) {
1856 if (number[j_prev*stride + k] < number[j*stride + k])
1857 break;
1858 else if (number[j_prev*stride + k] > number[j*stride + k])
1859 unordered = true;
1860 }
1861 if (unordered == true)
1862 break;
1863 }
1864 }
1865 else {
1866 for (i = 1 ; i < n_entities ; i++) {
1867 size_t i_prev, k;
1868 bool unordered = false;
1869 i_prev = i-1;
1870 for (k = 0; k < stride; k++) {
1871 if (number[i_prev*stride + k] < number[i*stride + k])
1872 break;
1873 else if (number[i_prev*stride + k] > number[i*stride + k])
1874 unordered = true;
1875 }
1876 if (unordered == true)
1877 break;
1878 }
1879 }
1880
1881 /* If numbering is implicit */
1882
1883 }
1884 else {
1885
1886 if (parent_entity_id != NULL) {
1887 for (i = 1 ; i < n_entities ; i++) {
1888 if (parent_entity_id[i] < parent_entity_id[i-1])
1889 break;
1890 }
1891 }
1892 else
1893 i = n_entities;
1894 }
1895
1896 if (i == n_entities || n_entities == 0)
1897 return 1;
1898 else
1899 return 0;
1900 }
1901
1902 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1903
1904 /*=============================================================================
1905 * Public function definitions
1906 *============================================================================*/
1907
1908 /*----------------------------------------------------------------------------
1909 * Creation of an I/O numbering structure.
1910 *
1911 * This function is similar to fvm_io_num_create_from_select, albeit
1912 * using parent entity numbers (1 to n) instead of ids (0 to n-1).
1913 *
1914 * parameters:
1915 * parent_entity_number <-- pointer to list of selected entitie's parent's
1916 * numbers, or NULL if all first n_ent entities
1917 * are used
1918 * parent_global_number <-- pointer to list of global (i.e. domain splitting
1919 * independent) parent entity numbers
1920 * n_entities <-- number of entities considered
1921 * share_parent_global <-- if non zero, try to share parent_global_number
1922 * instead of using a local copy
1923 *
1924 * returns:
1925 * pointer to I/O numbering structure
1926 *----------------------------------------------------------------------------*/
1927
1928 fvm_io_num_t *
fvm_io_num_create(const cs_lnum_t parent_entity_number[],const cs_gnum_t parent_global_number[],size_t n_entities,int share_parent_global)1929 fvm_io_num_create(const cs_lnum_t parent_entity_number[],
1930 const cs_gnum_t parent_global_number[],
1931 size_t n_entities,
1932 int share_parent_global)
1933 {
1934 cs_lnum_t *parent_entity_id = NULL;
1935
1936 if (parent_entity_number != NULL) {
1937 BFT_MALLOC(parent_entity_id, n_entities, cs_lnum_t);
1938 for (size_t i = 0; i < n_entities; i++)
1939 parent_entity_id[i] = parent_entity_number[i] - 1;
1940 }
1941
1942 fvm_io_num_t *this_io_num
1943 = fvm_io_num_create_from_select(parent_entity_id,
1944 parent_global_number,
1945 n_entities,
1946 share_parent_global);
1947
1948 BFT_FREE(parent_entity_id);
1949
1950 return this_io_num;
1951 }
1952
1953 /*----------------------------------------------------------------------------
1954 * Creation of an I/O numbering structure based on a selection of entities.
1955 *
1956 * parameters:
1957 * parent_entity_id <-- pointer to list of selected entitie's parent's
1958 * ids, or NULL if all first n_ent entities
1959 * are used
1960 * parent_global_number <-- pointer to list of global (i.e. domain splitting
1961 * independent) parent entity numbers
1962 * n_entities <-- number of entities considered
1963 * share_parent_global <-- if non zero, try to share parent_global_number
1964 * instead of using a local copy
1965 *
1966 * returns:
1967 * pointer to I/O numbering structure
1968 *----------------------------------------------------------------------------*/
1969
1970 fvm_io_num_t *
fvm_io_num_create_from_select(const cs_lnum_t parent_entity_id[],const cs_gnum_t parent_global_number[],size_t n_entities,int share_parent_global)1971 fvm_io_num_create_from_select(const cs_lnum_t parent_entity_id[],
1972 const cs_gnum_t parent_global_number[],
1973 size_t n_entities,
1974 int share_parent_global)
1975 {
1976 size_t i;
1977 cs_lnum_t *order = NULL;
1978
1979 fvm_io_num_t *this_io_num = NULL;
1980
1981 /* Initial checks */
1982
1983 if (cs_glob_n_ranks < 2 && parent_global_number == NULL)
1984 return NULL;
1985
1986 /* Create structure */
1987
1988 BFT_MALLOC(this_io_num, 1, fvm_io_num_t);
1989
1990 this_io_num->global_num_size = n_entities;
1991
1992 BFT_MALLOC(this_io_num->_global_num, n_entities, cs_gnum_t);
1993 this_io_num->global_num = this_io_num->_global_num;
1994
1995 if (n_entities > 0) {
1996
1997 /* Assign initial global numbers */
1998
1999 if (parent_entity_id != NULL) {
2000 for (i = 0 ; i < n_entities ; i++)
2001 this_io_num->_global_num[i] = parent_global_number[parent_entity_id[i]];
2002 }
2003 else {
2004 for (i = 0 ; i < n_entities ; i++)
2005 this_io_num->_global_num[i] = parent_global_number[i];
2006 }
2007
2008 if (cs_order_gnum_test(NULL,
2009 this_io_num->_global_num,
2010 n_entities) == false) {
2011 cs_gnum_t *tmp_num;
2012 order = cs_order_gnum(NULL,
2013 this_io_num->_global_num,
2014 n_entities);
2015 BFT_MALLOC(tmp_num, n_entities, cs_gnum_t);
2016 for (i = 0; i < n_entities; i++)
2017 tmp_num[i] = this_io_num->_global_num[order[i]];
2018 memcpy(this_io_num->_global_num, tmp_num, n_entities*sizeof(cs_gnum_t));
2019 BFT_FREE(tmp_num);
2020 }
2021 }
2022
2023 this_io_num->global_count = n_entities;
2024
2025 _fvm_io_num_copy_on_write(this_io_num);
2026
2027 /* Order globally */
2028
2029 #if defined(HAVE_MPI)
2030
2031 if (cs_glob_n_ranks > 1)
2032 _fvm_io_num_global_order(this_io_num,
2033 NULL,
2034 cs_glob_mpi_comm);
2035
2036 #endif
2037
2038 if (cs_glob_n_ranks == 1)
2039 _fvm_io_num_local_order(this_io_num,
2040 NULL);
2041
2042 if (order != NULL) {
2043 cs_gnum_t *tmp_num;
2044 BFT_MALLOC(tmp_num, n_entities, cs_gnum_t);
2045 for (i = 0; i < n_entities; i++)
2046 tmp_num[order[i]] = this_io_num->_global_num[i];
2047 memcpy(this_io_num->_global_num, tmp_num, n_entities*sizeof(cs_gnum_t));
2048 BFT_FREE(tmp_num);
2049 BFT_FREE(order);
2050 }
2051
2052 if (share_parent_global != 0)
2053 _fvm_io_num_try_to_set_shared(this_io_num,
2054 parent_global_number);
2055
2056 return this_io_num;
2057 }
2058
2059 /*----------------------------------------------------------------------------
2060 * Creation of an I/O numbering structure,
2061 * sharing a given global numbering array.
2062 *
2063 * parameters:
2064 * global_number <-- pointer to list of global (i.e. domain splitting
2065 * independent) entity numbers
2066 * global_count <-- global number of entities
2067 * n_entities <-- number of local entities considered
2068 *
2069 * returns:
2070 * pointer to I/O numbering structure
2071 *----------------------------------------------------------------------------*/
2072
2073 fvm_io_num_t *
fvm_io_num_create_shared(const cs_gnum_t global_number[],cs_gnum_t global_count,size_t n_entities)2074 fvm_io_num_create_shared(const cs_gnum_t global_number[],
2075 cs_gnum_t global_count,
2076 size_t n_entities)
2077 {
2078 fvm_io_num_t *this_io_num;
2079
2080 /* Create structure */
2081
2082 BFT_MALLOC(this_io_num, 1, fvm_io_num_t);
2083
2084 this_io_num->global_count = global_count;
2085 this_io_num->global_num_size = n_entities;
2086
2087 this_io_num->global_num = global_number;
2088 this_io_num->_global_num = NULL;
2089
2090 return this_io_num;
2091 }
2092
2093 /*----------------------------------------------------------------------------
2094 * Creation of an I/O numbering structure based on an an initial
2095 * I/O numbering and a number of new entities per base entity.
2096 *
2097 * This is useful for example to create an I/O numbering for
2098 * triangles based on split polygons, whose I/O numbering is defined.
2099 *
2100 * parameters:
2101 * base_io_num <-- pointer to base I/O numbering structure
2102 * n_sub_entities <-- number of new entities per base entity
2103 *
2104 * returns:
2105 * pointer to I/O numbering structure
2106 *----------------------------------------------------------------------------*/
2107
2108 fvm_io_num_t *
fvm_io_num_create_from_sub(const fvm_io_num_t * base_io_num,const cs_lnum_t n_sub_entities[])2109 fvm_io_num_create_from_sub(const fvm_io_num_t *base_io_num,
2110 const cs_lnum_t n_sub_entities[])
2111 {
2112 fvm_io_num_t *this_io_num = NULL;
2113
2114 /* Initial checks */
2115
2116 if (base_io_num == NULL)
2117 return NULL;
2118
2119 /* Create structure */
2120
2121 BFT_MALLOC(this_io_num, 1, fvm_io_num_t);
2122
2123 cs_lnum_t n_ent = base_io_num->global_num_size;
2124
2125 BFT_MALLOC(this_io_num->_global_num, n_ent, cs_gnum_t);
2126 this_io_num->global_num = this_io_num->_global_num;
2127
2128 this_io_num->global_num_size = n_ent;
2129
2130 /* Assign initial global numbers */
2131
2132 for (cs_lnum_t i = 0 ; i < n_ent ; i++)
2133 this_io_num->_global_num[i] = base_io_num->global_num[i];
2134
2135 this_io_num->global_count = n_ent;
2136
2137 _fvm_io_num_copy_on_write(this_io_num);
2138
2139 /* Order globally */
2140
2141 #if defined(HAVE_MPI)
2142
2143 if (cs_glob_n_ranks > 1)
2144 _fvm_io_num_global_order(this_io_num,
2145 n_sub_entities,
2146 cs_glob_mpi_comm);
2147
2148 #endif
2149
2150 if (cs_glob_n_ranks == 1)
2151 _fvm_io_num_local_order(this_io_num,
2152 n_sub_entities);
2153
2154 return this_io_num;
2155 }
2156
2157 /*----------------------------------------------------------------------------
2158 * Creation of an I/O numbering structure based on a strided adjacency.
2159 *
2160 * The corresponding entities must be locally ordered.
2161 *
2162 * parameters:
2163 * parent_entity_id <-- pointer to list of selected entitie's parent's ids,
2164 * or NULL if all first n_ent entities are used
2165 * index <-- index on entities for adjacency
2166 * adjacency <-- entity adjacency (1 to n global numbering)
2167 * n_entities <-- number of entities considered
2168 *
2169 * returns:
2170 * pointer to I/O numbering structure
2171 *----------------------------------------------------------------------------*/
2172
2173 fvm_io_num_t *
fvm_io_num_create_from_adj_s(const cs_lnum_t parent_entity_id[],const cs_gnum_t adjacency[],size_t n_entities,size_t stride)2174 fvm_io_num_create_from_adj_s(const cs_lnum_t parent_entity_id[],
2175 const cs_gnum_t adjacency[],
2176 size_t n_entities,
2177 size_t stride)
2178 {
2179 fvm_io_num_t *this_io_num = NULL;
2180
2181 /* Initial checks */
2182
2183 if (cs_glob_n_ranks < 2)
2184 return NULL;
2185
2186 assert(_is_gnum_ordered_s(parent_entity_id,
2187 adjacency,
2188 stride,
2189 n_entities) == true);
2190
2191 #if defined(HAVE_MPI)
2192 {
2193 cs_gnum_t *_adjacency = NULL;
2194
2195 /* Create structure */
2196
2197 BFT_MALLOC(this_io_num, 1, fvm_io_num_t);
2198
2199 this_io_num->global_num_size = n_entities;
2200
2201 BFT_MALLOC(this_io_num->_global_num, n_entities, cs_gnum_t);
2202 this_io_num->global_num = this_io_num->_global_num;
2203
2204 if (n_entities > 0) {
2205
2206 size_t i, j;
2207
2208 /* Assign initial global numbers */
2209
2210 BFT_MALLOC(_adjacency, n_entities*stride, cs_gnum_t);
2211
2212 if (parent_entity_id != NULL) {
2213 for (i = 0 ; i < n_entities ; i++) {
2214 for (j = 0; j < stride; j++)
2215 _adjacency[i*stride + j]
2216 = adjacency[parent_entity_id[i]*stride + j];
2217 }
2218 }
2219 else
2220 memcpy(_adjacency, adjacency, n_entities*stride*sizeof(cs_gnum_t));
2221
2222 }
2223
2224 /* Order globally */
2225
2226 this_io_num->global_count = n_entities;
2227
2228 _fvm_io_num_global_order_s(this_io_num,
2229 stride,
2230 _adjacency,
2231 cs_glob_mpi_comm);
2232
2233 BFT_FREE(_adjacency);
2234 }
2235 #endif
2236
2237 return this_io_num;
2238 }
2239
2240 /*----------------------------------------------------------------------------
2241 * Creation of an I/O numbering structure based on an indexed adjacency.
2242 *
2243 * The corresponding entities do not need to be locally ordered.
2244 *
2245 * parameters:
2246 * parent_entity_id <-- pointer to list of selected entitie's parent's ids,
2247 * or NULL if all first n_ent entities are used
2248 * index <-- index on entities for adjacency
2249 * adjacency <-- entity adjacency (1 to n global numbering)
2250 * n_entities <-- number of entities considered
2251 *
2252 * returns:
2253 * pointer to I/O numbering structure
2254 *----------------------------------------------------------------------------*/
2255
2256 fvm_io_num_t *
fvm_io_num_create_from_adj_i(const cs_lnum_t parent_entity_id[],const cs_lnum_t index[],const cs_gnum_t adjacency[],cs_lnum_t n_entities)2257 fvm_io_num_create_from_adj_i(const cs_lnum_t parent_entity_id[],
2258 const cs_lnum_t index[],
2259 const cs_gnum_t adjacency[],
2260 cs_lnum_t n_entities)
2261 {
2262 fvm_io_num_t *this_io_num = NULL;
2263
2264 /* Initial checks */
2265
2266 if (cs_glob_n_ranks < 2)
2267 return NULL;
2268
2269 #if defined(HAVE_MPI)
2270 {
2271 cs_lnum_t *_index = NULL;
2272 cs_gnum_t *_adjacency = NULL;
2273
2274 #if defined(DEBUG) && !defined(NDEBUG)
2275 const char no_adjacent_elt_msg[]
2276 = " Error during the creation of a fvm_io_num_t structure\n"
2277 " for an indexed adjacency.\n\n"
2278 " At least one entity has no adjacent element, and\n"
2279 " this case is not currently handled.\n\n"
2280 " Check if this definition is correct.";
2281 #endif
2282
2283 /* Create structure */
2284
2285 BFT_MALLOC(this_io_num, 1, fvm_io_num_t);
2286
2287 this_io_num->global_num_size = n_entities;
2288
2289 BFT_MALLOC(this_io_num->_global_num, n_entities, cs_gnum_t);
2290 this_io_num->global_num = this_io_num->_global_num;
2291
2292 BFT_MALLOC(_index, n_entities + 1, cs_lnum_t);
2293 _index[0] = 0;
2294
2295 if (n_entities > 0) {
2296
2297 cs_lnum_t i, j, k, ent_id, _shift;
2298
2299 /* Assign initial global numbers */
2300
2301 if (parent_entity_id != NULL) {
2302
2303 #if defined(DEBUG) && !defined(NDEBUG)
2304 for (i = 0 ; i < n_entities ; i++) {
2305 ent_id = parent_entity_id[i];
2306 if ((index[ent_id+1] - index[ent_id]) == 0)
2307 bft_error(__FILE__, __LINE__, 0, no_adjacent_elt_msg);
2308 }
2309 #endif
2310
2311 /* Count reduced size */
2312
2313 for (i = 0 ; i < n_entities ; i++) {
2314 ent_id = parent_entity_id[i];
2315 _index[i+1] = index[ent_id+1] - index[ent_id];
2316 }
2317
2318 for (i = 0 ; i < n_entities ; i++)
2319 _index[i+1] += _index[i];
2320
2321 BFT_MALLOC(_adjacency, _index[n_entities], cs_gnum_t);
2322
2323 /* Define reduced index and adjacency */
2324
2325 for (i = 0 ; i < n_entities ; i++) {
2326
2327 ent_id = parent_entity_id[i];
2328 _shift = _index[i];
2329
2330 for (j = index[ent_id], k = 0; j < index[ent_id+1]; j++, k++)
2331 _adjacency[_shift + k] = adjacency[j];
2332
2333 }
2334
2335 }
2336 else {
2337
2338 #if defined(DEBUG) && !defined(NDEBUG)
2339 for (i = 0 ; i < n_entities ; i++) {
2340 if ((index[i+1] - index[i]) == 0)
2341 bft_error(__FILE__, __LINE__, 0, no_adjacent_elt_msg);
2342 }
2343 #endif
2344
2345 BFT_MALLOC(_adjacency, index[n_entities], cs_gnum_t);
2346
2347 memcpy(_index, index, (n_entities+1)*sizeof(cs_lnum_t));
2348 memcpy(_adjacency, adjacency, index[n_entities]*sizeof(cs_gnum_t));
2349
2350 }
2351
2352 }
2353
2354 /* Order globally */
2355
2356 this_io_num->global_count = n_entities;
2357
2358 _fvm_io_num_global_order_index(this_io_num,
2359 _index,
2360 _adjacency,
2361 cs_glob_mpi_comm);
2362
2363 if (_adjacency != NULL)
2364 BFT_FREE(_adjacency);
2365 if (_index != NULL)
2366 BFT_FREE(_index);
2367 }
2368 #endif
2369
2370 return this_io_num;
2371 }
2372
2373 /*----------------------------------------------------------------------------
2374 * Creation of an I/O numbering structure based on a space-filling curve.
2375 *
2376 * It is expected that entities are unique (i.e. not duplicated on 2 or
2377 * more ranks). If 2 entities have a same Morton or Hilbert code, their global
2378 * number will be determined by lexicographical ordering of coordinates.
2379 *
2380 * parameters:
2381 * coords <-- pointer to entity coordinates (interlaced)
2382 * dim <-- spatial dimension
2383 * n_entities <-- number of entities considered
2384 * sfc_type <-- type of space-filling curve (Morton or Hilbert)
2385 *
2386 * returns:
2387 * pointer to I/O numbering structure
2388 *----------------------------------------------------------------------------*/
2389
2390 fvm_io_num_t *
fvm_io_num_create_from_sfc(const cs_coord_t coords[],int dim,size_t n_entities,fvm_io_num_sfc_t sfc_type)2391 fvm_io_num_create_from_sfc(const cs_coord_t coords[],
2392 int dim,
2393 size_t n_entities,
2394 fvm_io_num_sfc_t sfc_type)
2395 {
2396 fvm_io_num_t *this_io_num = NULL;
2397
2398 switch(sfc_type) {
2399 case FVM_IO_NUM_SFC_MORTON_BOX:
2400 this_io_num = _create_from_coords_morton(coords, dim, n_entities, 0);
2401 break;
2402 case FVM_IO_NUM_SFC_MORTON_CUBE:
2403 this_io_num = _create_from_coords_morton(coords, dim, n_entities, 1);
2404 break;
2405 case FVM_IO_NUM_SFC_HILBERT_BOX:
2406 this_io_num = _create_from_coords_hilbert(coords, dim, n_entities, 0);
2407 break;
2408 case FVM_IO_NUM_SFC_HILBERT_CUBE:
2409 this_io_num = _create_from_coords_hilbert(coords, dim, n_entities, 1);
2410 break;
2411 default:
2412 assert(0);
2413 }
2414
2415 return this_io_num;
2416 }
2417
2418 /*----------------------------------------------------------------------------
2419 * Creation of an I/O numbering structure based on real values, assuming
2420 * ordering by increasing values.
2421 *
2422 * It is expected that entities are unique (i.e. not duplicated on 2 or
2423 * more ranks). If 2 entities have a same value, their global
2424 * number will be determined by their initial order.
2425 *
2426 * parameters:
2427 * val <-- pointer to real values
2428 * n_entities <-- number of entities considered
2429 *
2430 * returns:
2431 * pointer to I/O numbering structure
2432 *----------------------------------------------------------------------------*/
2433
2434 fvm_io_num_t *
fvm_io_num_create_from_real(const cs_real_t val[],size_t n_entities)2435 fvm_io_num_create_from_real(const cs_real_t val[],
2436 size_t n_entities)
2437 {
2438 size_t i;
2439
2440 #if defined(HAVE_MPI)
2441 MPI_Comm comm = cs_glob_mpi_comm;
2442 #endif
2443
2444 const int n_ranks = cs_glob_n_ranks;
2445
2446 fvm_io_num_t *this_io_num = NULL;
2447
2448 /* Create structure */
2449
2450 BFT_MALLOC(this_io_num, 1, fvm_io_num_t);
2451
2452 this_io_num->global_num_size = n_entities;
2453
2454 BFT_MALLOC(this_io_num->_global_num, n_entities, cs_gnum_t);
2455 this_io_num->global_num = this_io_num->_global_num;
2456
2457 /* Scale values */
2458
2459 double v_min = DBL_MAX;
2460 double v_max = -DBL_MAX;
2461
2462 for (i = 0; i < n_entities; i++) {
2463 if (val[i] < v_min) v_min = val[i];
2464 if (val[i] > v_max) v_max = val[i];
2465 }
2466
2467 cs_parall_min(1, CS_DOUBLE, &v_min);
2468 cs_parall_max(1, CS_DOUBLE, &v_max);
2469
2470 if (v_max <= v_min) {
2471 cs_gnum_t n_g_entities = n_entities;
2472 cs_parall_counter(&n_g_entities, 1);
2473 if (n_g_entities > 1)
2474 bft_error(__FILE__, __LINE__, 0,
2475 _("%s: point set contains identical values."),
2476 __func__);
2477 }
2478
2479 #if defined(HAVE_MPI)
2480
2481 if (n_ranks > 1) {
2482
2483 double scale;
2484 if (v_max > v_min)
2485 scale = (1.0 - 1.e-12) / (v_max - v_min);
2486 else
2487 scale = 0;
2488
2489 cs_real_t *s_val;
2490 BFT_MALLOC(s_val, n_entities, cs_real_t);
2491
2492 for (i = 0; i < n_entities; i++)
2493 s_val[i] = (val[i] - v_min)*scale;
2494
2495 int *dest_rank = NULL;
2496 cs_lnum_t *order = NULL;
2497
2498 BFT_MALLOC(order, n_entities, cs_lnum_t);
2499 BFT_MALLOC(dest_rank, n_entities, int);
2500
2501 cs_order_real_allocated(NULL, s_val, order, n_entities);
2502
2503 cs_sort_partition_dest_rank_id(_sampling_factors[1],
2504 sizeof(cs_real_t),
2505 n_entities,
2506 s_val,
2507 NULL, /* weight */
2508 order,
2509 dest_rank,
2510 _s_to_real,
2511 _s_compare,
2512 NULL,
2513 comm);
2514
2515 BFT_FREE(order);
2516
2517 cs_all_to_all_t
2518 *d = cs_all_to_all_create(this_io_num->global_num_size,
2519 0, /* flags */
2520 NULL, /* dest_id */
2521 dest_rank,
2522 comm);
2523
2524 cs_all_to_all_transfer_dest_rank(d, &dest_rank);
2525
2526 cs_real_t *b_val
2527 = cs_all_to_all_copy_array(d,
2528 CS_REAL_TYPE,
2529 1,
2530 false, /* reverse */
2531 s_val,
2532 NULL);
2533
2534 BFT_FREE(s_val);
2535
2536 size_t b_size = cs_all_to_all_n_elts_dest(d);
2537
2538 /* Now re-build order on block distribution. */
2539
2540 BFT_MALLOC(order, b_size, cs_lnum_t);
2541
2542 cs_order_real_allocated(NULL, b_val, order, b_size);
2543
2544 BFT_FREE(b_val);
2545
2546 /* Determine global order; requires ordering to loop through buffer by
2547 increasing number (blocks associated with each process are
2548 already sorted, but the whole "gathered" block is not).
2549 We build an initial global order based on the initial global numbering,
2550 such that for each block, the global number of an entity is equal to
2551 the cumulative number of sub-entities */
2552
2553 cs_gnum_t *b_gnum;
2554 BFT_MALLOC(b_gnum, b_size, cs_gnum_t);
2555
2556 for (i = 0; i < b_size; i++)
2557 b_gnum[order[i]] = i+1;
2558
2559 BFT_FREE(order);
2560
2561 cs_gnum_t gnum_shift = 0, current_gnum = b_size;
2562
2563 /* At this stage, b_gnum[] is valid for this process, and
2564 current_gnum indicates the total number of entities handled
2565 by this process; we must now shift global numberings on different
2566 processes by the cumulative total number of entities handled by
2567 each process */
2568
2569 MPI_Scan(¤t_gnum, &gnum_shift, 1, CS_MPI_GNUM,
2570 MPI_SUM, comm);
2571 gnum_shift -= current_gnum;
2572
2573 for (i = 0; i < b_size; i++)
2574 b_gnum[i] += gnum_shift;
2575
2576 /* Return global order to all ranks */
2577
2578 cs_all_to_all_copy_array(d,
2579 CS_GNUM_TYPE,
2580 1,
2581 true, /* reverse */
2582 b_gnum,
2583 this_io_num->_global_num);
2584
2585 /* Free memory */
2586
2587 BFT_FREE(b_gnum);
2588
2589 cs_all_to_all_destroy(&d);
2590
2591 /* Get final maximum global number value */
2592
2593 this_io_num->global_count
2594 = _fvm_io_num_global_max_unordered(this_io_num, comm);
2595
2596 }
2597
2598 #endif /* HAVE_MPI */
2599
2600 if (n_ranks == 1) {
2601
2602 cs_lnum_t *order = NULL;
2603 BFT_MALLOC(order, n_entities, cs_lnum_t);
2604
2605 cs_order_real_allocated(NULL, val, order, n_entities);
2606
2607 for (i = 0; i < n_entities; i++)
2608 this_io_num->_global_num[order[i]] = i+1;
2609
2610 BFT_FREE(order);
2611
2612 this_io_num->global_count = n_entities;
2613
2614 }
2615
2616 return this_io_num;
2617 }
2618
2619 /*----------------------------------------------------------------------------
2620 * Creation of an I/O numbering structure based on a simple accumulation
2621 * (i.e. scan) of counts on successive ranks.
2622 *
2623 * parameters:
2624 * n_entities <-- number of entities considered
2625 *
2626 * returns:
2627 * pointer to I/O numbering structure
2628 *----------------------------------------------------------------------------*/
2629
2630 fvm_io_num_t *
fvm_io_num_create_from_scan(size_t n_entities)2631 fvm_io_num_create_from_scan(size_t n_entities)
2632 {
2633 fvm_io_num_t *this_io_num = NULL;
2634
2635 /* Create structure */
2636
2637 BFT_MALLOC(this_io_num, 1, fvm_io_num_t);
2638
2639 BFT_MALLOC(this_io_num->_global_num, n_entities, cs_gnum_t);
2640 this_io_num->global_num = this_io_num->_global_num;
2641
2642 this_io_num->global_num_size = n_entities;
2643
2644 #if defined(HAVE_MPI)
2645 if (cs_glob_n_ranks > 1) {
2646 cs_gnum_t gnum_base = n_entities;
2647 cs_gnum_t gnum_sum = n_entities;
2648 cs_gnum_t gnum_shift = 0;
2649
2650 MPI_Comm comm = cs_glob_mpi_comm;
2651
2652 MPI_Scan(&gnum_base, &gnum_shift, 1, CS_MPI_GNUM, MPI_SUM, comm);
2653
2654 gnum_base = gnum_shift - gnum_base + 1;
2655
2656 for (size_t i = 0; i < n_entities; i++)
2657 this_io_num->_global_num[i] = gnum_base + i;
2658
2659 gnum_base = n_entities;
2660
2661 MPI_Allreduce(&gnum_base, &gnum_sum, 1, CS_MPI_GNUM, MPI_SUM, comm);
2662
2663 this_io_num->global_count = gnum_sum;
2664 }
2665 #endif
2666
2667 if (cs_glob_n_ranks < 2) {
2668 for (size_t i = 0; i < n_entities; i++)
2669 this_io_num->_global_num[i] = i+1;
2670
2671 this_io_num->global_count = n_entities;
2672 }
2673
2674 return this_io_num;
2675 }
2676
2677 /*----------------------------------------------------------------------------
2678 * Destruction of a I/O numbering structure.
2679 *
2680 * parameters:
2681 * this_io_num <-- pointer to structure that should be destroyed
2682 *
2683 * returns:
2684 * NULL pointer
2685 *----------------------------------------------------------------------------*/
2686
2687 fvm_io_num_t *
fvm_io_num_destroy(fvm_io_num_t * this_io_num)2688 fvm_io_num_destroy(fvm_io_num_t * this_io_num)
2689 {
2690 if (this_io_num != NULL) {
2691 BFT_FREE(this_io_num->_global_num);
2692 BFT_FREE(this_io_num);
2693 }
2694
2695 return this_io_num;
2696 }
2697
2698 /*----------------------------------------------------------------------------
2699 * Transfer ownership of global numbering array from IO numbering structure.
2700 *
2701 * parameters:
2702 * this_io_num <-> pointer to structure transferring array ownership.
2703 *
2704 * returns:
2705 * pointer to transferred array
2706 *----------------------------------------------------------------------------*/
2707
2708 cs_gnum_t *
fvm_io_num_transfer_global_num(fvm_io_num_t * this_io_num)2709 fvm_io_num_transfer_global_num(fvm_io_num_t * this_io_num)
2710 {
2711 cs_gnum_t *retval = NULL;
2712
2713 if (this_io_num != NULL) {
2714 retval = this_io_num->_global_num;
2715 this_io_num->_global_num = NULL;
2716 }
2717
2718 return retval;
2719 }
2720
2721 /*----------------------------------------------------------------------------
2722 * Return local number of entities associated with an I/O numbering
2723 * structure.
2724 *
2725 * parameters:
2726 * this_io_num <-- pointer to I/O/ numbering structure
2727 *
2728 * returns:
2729 * local number of associated entities
2730 *----------------------------------------------------------------------------*/
2731
2732 cs_lnum_t
fvm_io_num_get_local_count(const fvm_io_num_t * const this_io_num)2733 fvm_io_num_get_local_count(const fvm_io_num_t *const this_io_num)
2734 {
2735 assert(this_io_num != NULL);
2736
2737 return this_io_num->global_num_size;
2738 }
2739
2740 /*----------------------------------------------------------------------------
2741 * Return global number of entities associated with an I/O numbering
2742 * structure.
2743 *
2744 * parameters:
2745 * this_io_num <-- pointer to I/O/ numbering structure
2746 *
2747 * returns:
2748 * global number of associated entities
2749 *----------------------------------------------------------------------------*/
2750
2751 cs_gnum_t
fvm_io_num_get_global_count(const fvm_io_num_t * const this_io_num)2752 fvm_io_num_get_global_count(const fvm_io_num_t *const this_io_num)
2753 {
2754 assert(this_io_num != NULL);
2755
2756 return this_io_num->global_count;
2757 }
2758
2759 /*----------------------------------------------------------------------------
2760 * Return global numbering associated with an I/O numbering structure.
2761 *
2762 * parameters:
2763 * this_io_num <-- pointer to I/O/ numbering structure
2764 *
2765 * returns:
2766 * pointer to array of global numbers associated with local entities
2767 * (1 to n numbering)
2768 *----------------------------------------------------------------------------*/
2769
2770 const cs_gnum_t *
fvm_io_num_get_global_num(const fvm_io_num_t * const this_io_num)2771 fvm_io_num_get_global_num(const fvm_io_num_t *const this_io_num)
2772 {
2773 assert(this_io_num != NULL);
2774
2775 return this_io_num->global_num;
2776 }
2777
2778 /*----------------------------------------------------------------------------
2779 * Return the global number of sub-entities associated with an initial
2780 * entity whose global numbering is known, given the number of
2781 * sub-entities per initial entity.
2782 *
2783 * parameters:
2784 * this_io_num <-> pointer to base io numbering
2785 * n_sub_entities <-- number of sub-entities per initial entity
2786 * comm <-- associated MPI communicator
2787 *
2788 * returns:
2789 * global number of sub-entities
2790 *----------------------------------------------------------------------------*/
2791
2792 cs_gnum_t
fvm_io_num_global_sub_size(const fvm_io_num_t * this_io_num,const cs_lnum_t n_sub_entities[])2793 fvm_io_num_global_sub_size(const fvm_io_num_t *this_io_num,
2794 const cs_lnum_t n_sub_entities[])
2795 {
2796 cs_gnum_t retval = 0;
2797
2798 /* Initial checks */
2799
2800 if (this_io_num == NULL)
2801 return retval;
2802
2803 #if defined(HAVE_MPI)
2804 if (cs_glob_n_ranks > 1) {
2805 int have_sub_loc = 0, have_sub_glob = 0;
2806
2807 /* Caution: we may have sub-entities on some ranks and not on others */
2808
2809 if (n_sub_entities != NULL)
2810 have_sub_loc = 1;
2811
2812 MPI_Allreduce(&have_sub_loc, &have_sub_glob, 1, MPI_INT, MPI_MAX,
2813 cs_glob_mpi_comm);
2814
2815 if (have_sub_glob > 0)
2816 retval = _fvm_io_num_global_sub_size(this_io_num,
2817 n_sub_entities,
2818 cs_glob_mpi_comm);
2819 }
2820 #endif
2821
2822 if (cs_glob_n_ranks == 1 && n_sub_entities != NULL) {
2823 for (size_t i = 0; i < (size_t)(this_io_num->global_num_size); i++)
2824 retval += n_sub_entities[i];
2825 }
2826
2827 return retval;
2828 }
2829
2830 /*----------------------------------------------------------------------------
2831 * Dump printout of a I/O numbering structure.
2832 *
2833 * parameters:
2834 * this_io_num <-- pointer to structure that should be dumped
2835 *----------------------------------------------------------------------------*/
2836
2837 void
fvm_io_num_dump(const fvm_io_num_t * const this_io_num)2838 fvm_io_num_dump(const fvm_io_num_t *const this_io_num)
2839 {
2840 cs_lnum_t i;
2841
2842 if (this_io_num == NULL) {
2843 bft_printf(" global numbering: nil\n");
2844 return;
2845 }
2846
2847 bft_printf(" global numbering size: %u\n",
2848 (unsigned)this_io_num->global_num_size);
2849
2850 bft_printf("\n"
2851 " pointer to shareable array:\n"
2852 " global_num: %p\n",
2853 (const void *)this_io_num->global_num);
2854
2855 bft_printf("\n"
2856 " pointer to local array:\n"
2857 " _global_num: %p\n",
2858 (const void *)this_io_num->global_num);
2859
2860 if (this_io_num->global_num_size > 0) {
2861
2862 bft_printf("\n global number:\n\n");
2863 for (i = 0 ; i < this_io_num->global_num_size ; i++)
2864 bft_printf(" %10u : %10llu\n",
2865 (unsigned)i + 1,
2866 (unsigned long long)this_io_num->global_num[i]);
2867 }
2868 }
2869
2870 /*----------------------------------------------------------------------------*/
2871
2872 END_C_DECLS
2873