1
/*============================================================================
2
* Set of subroutines for:
3
* - merging equivalent vertices,
4
* - managing tolerance reduction
5
*===========================================================================*/
8
This file is part of Code_Saturne, a general-purpose CFD tool.
10
Copyright (C) 1998-2011 EDF S.A.
12
This program is free software; you can redistribute it and/or modify it under
13
the terms of the GNU General Public License as published by the Free Software
14
Foundation; either version 2 of the License, or (at your option) any later
17
This program is distributed in the hope that it will be useful, but WITHOUT
18
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
22
You should have received a copy of the GNU General Public License along with
23
this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
24
Street, Fifth Floor, Boston, MA 02110-1301, USA.
27
/*----------------------------------------------------------------------------*/
29
#if defined(HAVE_CONFIG_H)
30
#include "cs_config.h"
33
/*----------------------------------------------------------------------------
34
* Standard C library headers
35
*---------------------------------------------------------------------------*/
44
/*----------------------------------------------------------------------------
46
*---------------------------------------------------------------------------*/
49
#include <bft_timer.h>
50
#include <bft_printf.h>
52
/*----------------------------------------------------------------------------
54
*---------------------------------------------------------------------------*/
56
#include <fvm_io_num.h>
57
#include <fvm_order.h>
58
#include <fvm_parall.h>
60
/*----------------------------------------------------------------------------
62
*---------------------------------------------------------------------------*/
64
#include "cs_search.h"
65
#include "cs_join_post.h"
67
/*----------------------------------------------------------------------------
68
* Header for the current file
69
*---------------------------------------------------------------------------*/
71
#include "cs_join_merge.h"
73
/*---------------------------------------------------------------------------*/
77
/*============================================================================
78
* Structure and type definitions
79
*===========================================================================*/
81
/*============================================================================
83
*===========================================================================*/
85
/* Turn on (1) or off (0) the tolerance reduc. */
86
#define CS_JOIN_MERGE_TOL_REDUC 1
87
#define CS_JOIN_MERGE_INV_TOL 1
89
/*============================================================================
90
* Global variable definitions
91
*===========================================================================*/
93
/* Parameters to control the vertex merge */
97
CS_JOIN_MERGE_MAX_GLOB_ITERS = 50, /* Max. number of glob. iter. for finding
98
equivalent vertices */
99
CS_JOIN_MERGE_MAX_LOC_ITERS = 100 /* Max. number of loc. iter. for finding
100
equivalent vertices */
103
/* Coefficient to deal with rounding approximations */
105
static const double cs_join_tol_eps_coef2 = 1.0001*1.001;
107
/* Counter on the number of loops useful to converge for the merge operation */
109
static int _glob_merge_counter = 0, _loc_merge_counter = 0;
111
/*============================================================================
112
* Private function definitions
113
*===========================================================================*/
115
/*----------------------------------------------------------------------------
116
* Initialize counter for the merge operation
117
*---------------------------------------------------------------------------*/
120
_initialize_merge_counter(void)
122
_glob_merge_counter = 0;
123
_loc_merge_counter = 0;
126
#if 0 && defined(DEBUG) && !defined(NDEBUG)
128
/*----------------------------------------------------------------------------
129
* Dump an cs_join_eset_t structure on vertices.
132
* e_set <-- cs_join_eset_t structure to dump
133
* mesh <-- cs_join_mesh_t structure associated
134
* logfile <-- handle to log file
135
*---------------------------------------------------------------------------*/
138
_dump_vtx_eset(const cs_join_eset_t *e_set,
139
const cs_join_mesh_t *mesh,
144
fprintf(logfile, "\n Dump an cs_join_eset_t structure (%p)\n",
145
(const void *)e_set);
146
fprintf(logfile, " n_max_equiv: %10d\n", e_set->n_max_equiv);
147
fprintf(logfile, " n_equiv : %10d\n\n", e_set->n_equiv);
149
for (i = 0; i < e_set->n_equiv; i++) {
151
cs_int_t v1_num = e_set->equiv_couple[2*i];
152
cs_int_t v2_num = e_set->equiv_couple[2*i+1];
155
" %10d - local: (%9d, %9d) - global: (%10llu, %10llu)\n",
157
(unsigned long long)(mesh->vertices[v1_num-1]).gnum,
158
(unsigned long long)(mesh->vertices[v2_num-1]).gnum);
164
#endif /* Only in debug mode */
166
/*----------------------------------------------------------------------------
167
* Compute the length of a segment between two vertices.
170
* v1 <-- cs_join_vertex_t structure for the first vertex of the segment
171
* v2 <-- cs_join_vertex_t structure for the second vertex of the segment
174
* length of the segment
175
*---------------------------------------------------------------------------*/
177
inline static cs_real_t
178
_compute_length(cs_join_vertex_t v1,
182
cs_real_t len = 0.0, d2 = 0.0;
184
for (k = 0; k < 3; k++) {
185
cs_real_t d = v1.coord[k] - v2.coord[k];
193
/*----------------------------------------------------------------------------
194
* Compute a new cs_join_vertex_t structure.
197
* curv_abs <-- curvilinear abscissa of the intersection
198
* gnum <-- global number associated to the new
199
* cs_join_vertex_t structure
200
* vtx_couple <-- couple of vertex numbers defining the current edge
201
* work <-- local cs_join_mesh_t structure
204
* a new cs_join_vertex_t structure
205
*---------------------------------------------------------------------------*/
207
static cs_join_vertex_t
208
_get_new_vertex(float curv_abs,
210
const cs_int_t vtx_couple[],
211
const cs_join_mesh_t *work)
214
cs_join_vertex_t new_vtx_data;
216
cs_join_vertex_t v1 = work->vertices[vtx_couple[0]-1];
217
cs_join_vertex_t v2 = work->vertices[vtx_couple[1]-1];
219
assert(curv_abs >= 0.0);
220
assert(curv_abs <= 1.0);
222
/* New vertex features */
224
new_vtx_data.state = CS_JOIN_STATE_NEW;
225
new_vtx_data.gnum = gnum;
226
new_vtx_data.tolerance = (1-curv_abs)*v1.tolerance + curv_abs*v2.tolerance;
228
for (k = 0; k < 3; k++)
229
new_vtx_data.coord[k] = (1-curv_abs)*v1.coord[k] + curv_abs*v2.coord[k];
234
/*----------------------------------------------------------------------------
235
* Define a tag (3 values) to globally order intersections.
238
* tag <-> tag to fill
239
* e1_gnum <-- global number for the first edge
240
* e2_gnum <-- global number for the second edge
241
* link_vtx_gnum <-- global number of the vertex associated to the current
243
*---------------------------------------------------------------------------*/
246
_define_inter_tag(fvm_gnum_t tag[],
249
fvm_gnum_t link_vtx_gnum)
251
if (e1_gnum < e2_gnum) {
260
tag[2] = link_vtx_gnum;
263
/*----------------------------------------------------------------------------
264
* Creation of new vertices.
266
* Update list of equivalent vertices.
269
* work <-- pointer to a cs_join_mesh_t structure
270
* edges <-- list of edges
271
* inter_set <-- structure including data on edge intersections
272
* init_max_vtx_gnum <-- initial max. global numbering for vertices
273
* n_iwm_vertices <-- initial local number of vertices (work struct)
274
* n_new_vertices <-- local number of new vertices to define
275
* p_n_g_new_vertices <-> pointer to the global number of new vertices
276
* p_new_vtx_gnum <-> pointer to the global numbering array for the
278
*---------------------------------------------------------------------------*/
281
_compute_new_vertex_gnum(const cs_join_mesh_t *work,
282
const cs_join_edges_t *edges,
283
const cs_join_inter_set_t *inter_set,
284
fvm_gnum_t init_max_vtx_gnum,
285
cs_int_t n_iwm_vertices,
286
cs_int_t n_new_vertices,
287
fvm_gnum_t *p_n_g_new_vertices,
288
fvm_gnum_t *p_new_vtx_gnum[])
292
fvm_gnum_t n_g_new_vertices = 0;
293
cs_int_t n_new_vertices_save = n_new_vertices;
294
fvm_lnum_t *order = NULL;
295
fvm_gnum_t *inter_tag = NULL, *adjacency = NULL, *new_vtx_gnum = NULL;
296
fvm_io_num_t *new_vtx_io_num = NULL;
298
/* Define a fvm_io_num_t structure to get the global numbering
299
for the new vertices.
300
First, build a tag associated to each intersection */
302
BFT_MALLOC(new_vtx_gnum, n_new_vertices, fvm_gnum_t);
303
BFT_MALLOC(inter_tag, 3*n_new_vertices, fvm_gnum_t);
307
for (i = 0; i < inter_set->n_inter; i++) {
309
cs_join_inter_t inter1 = inter_set->inter_lst[2*i];
310
cs_join_inter_t inter2 = inter_set->inter_lst[2*i+1];
311
fvm_gnum_t e1_gnum = edges->gnum[inter1.edge_id];
312
fvm_gnum_t e2_gnum = edges->gnum[inter2.edge_id];
314
if (inter1.vtx_id + 1 > n_iwm_vertices) {
316
if (inter2.vtx_id + 1 > n_iwm_vertices)
317
_define_inter_tag(&(inter_tag[3*n_new_vertices]),
321
_define_inter_tag(&(inter_tag[3*n_new_vertices]),
323
(work->vertices[inter2.vtx_id]).gnum);
327
} /* New vertices for this intersection */
329
if (inter2.vtx_id + 1 > n_iwm_vertices) {
331
if (inter1.vtx_id + 1 > n_iwm_vertices)
332
_define_inter_tag(&(inter_tag[3*n_new_vertices]),
334
init_max_vtx_gnum + 1);
336
_define_inter_tag(&(inter_tag[3*n_new_vertices]),
338
(work->vertices[inter1.vtx_id]).gnum);
342
} /* New vertices for this intersection */
344
} /* End of loop on intersections */
346
if (n_new_vertices != n_new_vertices_save)
347
bft_error(__FILE__, __LINE__, 0,
348
_(" The number of new vertices to create is not consistent.\n"
349
" Previous number: %10d\n"
350
" Current number: %10d\n\n"),
351
n_new_vertices_save, n_new_vertices);
353
/* Create a new fvm_io_num_t structure */
355
BFT_MALLOC(order, n_new_vertices, fvm_lnum_t);
357
fvm_order_local_allocated_s(NULL, inter_tag, 3, order, n_new_vertices);
359
BFT_MALLOC(adjacency, 3*n_new_vertices, fvm_gnum_t);
361
for (i = 0; i < n_new_vertices; i++) {
363
cs_int_t o_id = order[i];
365
adjacency[3*i] = inter_tag[3*o_id];
366
adjacency[3*i+1] = inter_tag[3*o_id+1];
367
adjacency[3*i+2] = inter_tag[3*o_id+2];
373
if (cs_glob_n_ranks > 1) {
375
const fvm_gnum_t *global_num = NULL;
378
fvm_io_num_create_from_adj_s(NULL, adjacency, n_new_vertices, 3);
380
n_g_new_vertices = fvm_io_num_get_global_count(new_vtx_io_num);
381
global_num = fvm_io_num_get_global_num(new_vtx_io_num);
383
for (i = 0; i < n_new_vertices; i++)
384
new_vtx_gnum[order[i]] = global_num[i] + init_max_vtx_gnum;
386
fvm_io_num_destroy(new_vtx_io_num);
388
} /* End of parallel treatment */
392
if (n_new_vertices > 0) {
394
fvm_gnum_t new_gnum = init_max_vtx_gnum + 1;
396
new_vtx_gnum[order[0]] = new_gnum;
398
for (i = 1; i < n_new_vertices; i++) {
400
if (adjacency[3*i] != adjacency[3*(i-1)])
403
if (adjacency[3*i+1] != adjacency[3*(i-1)+1])
406
if (adjacency[3*i+2] != adjacency[3*(i-1)+2])
410
new_vtx_gnum[order[i]] = new_gnum;
414
} /* End if n_new_vertices > 0 */
416
n_g_new_vertices = n_new_vertices;
418
} /* End of serial treatment */
427
*p_n_g_new_vertices = n_g_new_vertices;
428
*p_new_vtx_gnum = new_vtx_gnum;
432
/*----------------------------------------------------------------------------
433
* Get vertex id associated to the current intersection.
435
* Create a new vertex id if needed. Update n_new_vertices in this case.
438
* inter <-- a inter_t structure
439
* vtx_couple <-- couple of vertex numbers defining the current edge
440
* n_init_vertices <-- initial number of vertices
441
* n_new_vertices <-- number of new vertices created
444
* vertex id associated to the current intersection.
445
*---------------------------------------------------------------------------*/
448
_get_vtx_id(cs_join_inter_t inter,
449
const cs_int_t vtx_couple[],
450
cs_int_t n_init_vertices,
451
cs_int_t *p_n_new_vertices)
453
cs_int_t vtx_id = -1;
454
cs_int_t n_new_vertices = *p_n_new_vertices;
456
assert(inter.curv_abs >= 0.0);
457
assert(inter.curv_abs <= 1.0);
459
if (inter.curv_abs <= 0.0)
460
vtx_id = vtx_couple[0] - 1;
462
else if (inter.curv_abs >= 1.0)
463
vtx_id = vtx_couple[1] - 1;
467
assert(inter.curv_abs > 0 && inter.curv_abs < 1.0);
468
vtx_id = n_init_vertices + n_new_vertices;
473
assert(vtx_id != -1);
475
*p_n_new_vertices = n_new_vertices;
481
/*----------------------------------------------------------------------------
482
* Test if we have to continue to spread the tag associate to each vertex
485
* n_vertices <-- local number of vertices
486
* prev_vtx_tag <-- previous tag for each vertex
487
* vtx_tag <-- tag for each vertex
491
*---------------------------------------------------------------------------*/
494
_is_spread_not_converged(cs_int_t n_vertices,
495
const fvm_gnum_t prev_vtx_tag[],
496
const fvm_gnum_t vtx_tag[])
500
cs_bool_t have_to_continue = true;
502
for (i = 0; i < n_vertices; i++)
503
if (vtx_tag[i] != prev_vtx_tag[i])
507
have_to_continue = false;
509
return have_to_continue;
512
/*----------------------------------------------------------------------------
513
* Spread the tag associated to each vertex according the rule:
514
* Between two equivalent vertices, the tag associated to each considered
515
* vertex is equal to the minimal global number.
518
* n_vertices <-- local number of vertices
519
* vtx_eset <-- structure dealing with vertices equivalences
520
* vtx_tag <-> tag for each vertex
521
*---------------------------------------------------------------------------*/
524
_spread_tag(cs_int_t n_vertices,
525
const cs_join_eset_t *vtx_eset,
526
fvm_gnum_t vtx_tag[])
528
cs_int_t i, v1_id, v2_id;
529
fvm_gnum_t v1_gnum, v2_gnum;
530
cs_int_t *equiv_lst = vtx_eset->equiv_couple;
532
for (i = 0; i < vtx_eset->n_equiv; i++) {
534
v1_id = equiv_lst[2*i] - 1, v2_id = equiv_lst[2*i+1] - 1;
535
assert(v1_id < n_vertices);
536
assert(v1_id < n_vertices);
537
v1_gnum = vtx_tag[v1_id], v2_gnum = vtx_tag[v2_id];
539
if (v1_gnum != v2_gnum) {
541
fvm_gnum_t min_gnum = CS_MIN(v1_gnum, v2_gnum);
543
vtx_tag[v1_id] = min_gnum;
544
vtx_tag[v2_id] = min_gnum;
547
} /* End of loop on vertex equivalences */
550
/*----------------------------------------------------------------------------
551
* Define an array wich keeps the new vertex id of each vertex.
553
* If two vertices have the same vertex id, they should merge.
556
* vtx_eset <-- structure dealing with vertex equivalences
557
* n_vertices <-- local number of vertices
558
* prev_vtx_tag <-> previous tag for each vertex
559
* vtx_tag <-> tag for each vertex
560
*---------------------------------------------------------------------------*/
563
_local_spread(const cs_join_eset_t *vtx_eset,
565
fvm_gnum_t prev_vtx_tag[],
566
fvm_gnum_t vtx_tag[])
570
_loc_merge_counter++;
572
_spread_tag(n_vertices, vtx_eset, vtx_tag);
574
while (_is_spread_not_converged(n_vertices, prev_vtx_tag, vtx_tag)) {
576
_loc_merge_counter++;
578
if (_loc_merge_counter > CS_JOIN_MERGE_MAX_LOC_ITERS)
579
bft_error(__FILE__, __LINE__, 0,
580
_("\n The authorized maximum number of iterations "
581
" for the merge of vertices has been reached.\n"
582
" Local counter on iteration : %d (MAX =%d)\n"
583
" Check the fraction parameter.\n"),
584
_loc_merge_counter, CS_JOIN_MERGE_MAX_LOC_ITERS);
586
for (i = 0; i < n_vertices; i++)
587
prev_vtx_tag[i] = vtx_tag[i];
589
_spread_tag(n_vertices, vtx_eset, vtx_tag);
593
#if defined(HAVE_MPI)
595
/*----------------------------------------------------------------------------
596
* Exchange local vtx_tag buffer over the ranks and update global vtx_tag
597
* buffers. Apply modifications observed on the global vtx_tag to the local
601
* block_size <-- size of block for the current rank
602
* work <-- local cs_join_mesh_t structure which has initial
604
* vtx_tag <-> local vtx_tag for the local vertices
605
* glob_vtx_tag <-> global vtx_tag affected to the local rank
607
* prev_glob_vtx_tag <-> same but for the previous iteration
608
* recv2glob <-> buffer used to place correctly receive elements
609
* send_count <-> buffer used to count the number of elts to send
610
* send_shift <-> index on ranks of the elements to send
611
* send_glob_buffer <-> buffer used to save elements to send
612
* recv_count <-> buffer used to count the number of elts to receive
613
* recv_shift <-> index on ranks of the elements to receive
614
* recv_glob_buffer <-> buffer used to save elements to receive
617
* true if we have to continue the spread, false otherwise.
618
*---------------------------------------------------------------------------*/
621
_global_spread(cs_int_t block_size,
622
const cs_join_mesh_t *work,
623
fvm_gnum_t vtx_tag[],
624
fvm_gnum_t glob_vtx_tag[],
625
fvm_gnum_t prev_glob_vtx_tag[],
626
fvm_gnum_t recv2glob[],
627
cs_int_t send_count[],
628
cs_int_t send_shift[],
629
fvm_gnum_t send_glob_buffer[],
630
cs_int_t recv_count[],
631
cs_int_t recv_shift[],
632
fvm_gnum_t recv_glob_buffer[])
635
cs_int_t i, local_value, global_value;
637
cs_int_t n_vertices = work->n_vertices;
638
int n_ranks = cs_glob_n_ranks;
639
MPI_Comm mpi_comm = cs_glob_mpi_comm;
641
_glob_merge_counter++;
643
/* Push modifications in local vtx_tag to the global vtx_tag */
645
for (i = 0; i < n_ranks; i++)
648
for (i = 0; i < n_vertices; i++) {
650
int rank = (work->vertices[i].gnum - 1) % n_ranks;
651
cs_int_t shift = send_shift[rank] + send_count[rank];
653
send_glob_buffer[shift] = vtx_tag[i];
654
send_count[rank] += 1;
658
MPI_Alltoallv(send_glob_buffer, send_count, send_shift, FVM_MPI_GNUM,
659
recv_glob_buffer, recv_count, recv_shift, FVM_MPI_GNUM,
662
/* Apply update to glob_vtx_tag */
664
for (i = 0; i < recv_shift[n_ranks]; i++) {
665
cs_int_t cur_id = recv2glob[i];
666
glob_vtx_tag[cur_id] = CS_MIN(glob_vtx_tag[cur_id], recv_glob_buffer[i]);
669
ret_value = _is_spread_not_converged(block_size,
673
if (ret_value == false)
678
MPI_Allreduce(&local_value, &global_value, 1, MPI_INT, MPI_SUM, mpi_comm);
680
if (global_value > 0) { /* Store the current state as the previous one
681
Update local vtx_tag */
683
if (_glob_merge_counter > CS_JOIN_MERGE_MAX_GLOB_ITERS)
684
bft_error(__FILE__, __LINE__, 0,
685
_("\n The authorized maximum number of iterations "
686
" for the merge of vertices has been reached.\n"
687
" Global counter on iteration : %d (MAX =%d)\n"
688
" Check the fraction parameter.\n"),
689
_glob_merge_counter, CS_JOIN_MERGE_MAX_GLOB_ITERS);
691
for (i = 0; i < block_size; i++)
692
prev_glob_vtx_tag[i] = glob_vtx_tag[i];
694
for (i = 0; i < recv_shift[n_ranks]; i++)
695
recv_glob_buffer[i] = glob_vtx_tag[recv2glob[i]];
697
MPI_Alltoallv(recv_glob_buffer, recv_count, recv_shift, FVM_MPI_GNUM,
698
send_glob_buffer, send_count, send_shift, FVM_MPI_GNUM,
703
for (i = 0; i < n_ranks; i++)
706
assert(send_shift[n_ranks] == n_vertices);
708
for (i = 0; i < n_vertices; i++) {
710
int rank = (work->vertices[i].gnum - 1) % n_ranks;
711
cs_int_t shift = send_shift[rank] + send_count[rank];
713
vtx_tag[i] = CS_MIN(send_glob_buffer[shift], vtx_tag[i]);
714
send_count[rank] += 1;
720
} /* End if prev_glob_vtx_tag != glob_vtx_tag */
723
return false; /* No need to continue */
726
/*----------------------------------------------------------------------------
727
* Initialize and allocate buffers for the tag operation in parallel mode.
730
* n_g_vertices_to_treat <-- global number of vertices to consider for the
731
* merge operation (existing + created vertices)
732
* work <-- local cs_join_mesh_t structure which has
733
* initial vertex data
734
* p_block_size <-> size of block for the current rank
735
* p_send_count <-> buf. for counting the number of elts to send
736
* p_send_shift <-> index on ranks of the elements to send
737
* p_send_glob_buf <-> buf. for saving elements to send
738
* p_recv_count <-> buf. for counting the number of elts to receive
739
* p_recv_shift <-> index on ranks of the elements to receive
740
* p_recv_glob_buf <-> buf. for storing elements to receive
741
* p_recv2glob <-> buf. for putting correctly received elements
742
* p_glob_vtx_tag <-> vtx_tag locally treated (size = block_size)
743
* p_prev_glob_vtx_tag <-> idem but for the previous iteration
744
*---------------------------------------------------------------------------*/
747
_parall_tag_init(fvm_gnum_t n_g_vertices_to_treat,
748
const cs_join_mesh_t *work,
749
cs_int_t *p_block_size,
750
cs_int_t *p_send_count[],
751
cs_int_t *p_send_shift[],
752
fvm_gnum_t *p_send_glob_buf[],
753
cs_int_t *p_recv_count[],
754
cs_int_t *p_recv_shift[],
755
fvm_gnum_t *p_recv_glob_buf[],
756
fvm_gnum_t *p_recv2glob[],
757
fvm_gnum_t *p_glob_vtx_tag[],
758
fvm_gnum_t *p_prev_glob_vtx_tag[])
762
cs_int_t n_vertices = work->n_vertices;
763
cs_int_t block_size = 0, left_over = 0;
764
cs_int_t *send_count = NULL, *recv_count = NULL;
765
cs_int_t *send_shift = NULL, *recv_shift = NULL;
766
fvm_gnum_t *recv2glob = NULL;
767
fvm_gnum_t *recv_glob_buffer = NULL, *send_glob_buffer = NULL;
768
fvm_gnum_t *glob_vtx_tag = NULL, *prev_glob_vtx_tag = NULL;
769
MPI_Comm mpi_comm = cs_glob_mpi_comm;
771
const int n_ranks = cs_glob_n_ranks;
772
const int local_rank = CS_MAX(cs_glob_rank_id, 0);
774
/* Allocate and intialize vtx_tag associated to the local rank */
776
block_size = n_g_vertices_to_treat / n_ranks;
777
left_over = n_g_vertices_to_treat % n_ranks;
778
if (local_rank < left_over)
781
BFT_MALLOC(glob_vtx_tag, block_size, fvm_gnum_t);
782
BFT_MALLOC(prev_glob_vtx_tag, block_size, fvm_gnum_t);
784
for (i = 0; i < block_size; i++) {
785
prev_glob_vtx_tag[i] = i*n_ranks + local_rank + 1;
786
glob_vtx_tag[i] = i*n_ranks + local_rank + 1;
789
/* Allocate and define send/recv_count/shift */
791
BFT_MALLOC(send_count, n_ranks, cs_int_t);
792
BFT_MALLOC(recv_count, n_ranks, cs_int_t);
793
BFT_MALLOC(send_shift, n_ranks + 1, cs_int_t);
794
BFT_MALLOC(recv_shift, n_ranks + 1, cs_int_t);
799
for (i = 0; i < n_ranks; i++)
802
for (i = 0; i < n_vertices; i++) {
803
int rank = (work->vertices[i].gnum - 1) % n_ranks;
804
send_count[rank] += 1;
807
MPI_Alltoall(&(send_count[0]), 1, FVM_MPI_LNUM,
808
&(recv_count[0]), 1, FVM_MPI_LNUM, mpi_comm);
812
for (i = 0; i < n_ranks; i++) {
813
send_shift[i+1] = send_shift[i] + send_count[i];
814
recv_shift[i+1] = recv_shift[i] + recv_count[i];
817
assert(send_shift[n_ranks] == n_vertices);
819
/* Allocate and define recv2glob */
821
BFT_MALLOC(send_glob_buffer, send_shift[n_ranks], fvm_gnum_t);
822
BFT_MALLOC(recv2glob, recv_shift[n_ranks], fvm_gnum_t);
823
BFT_MALLOC(recv_glob_buffer, recv_shift[n_ranks], fvm_gnum_t);
825
for (i = 0; i < n_ranks; i++)
828
for (i = 0; i < n_vertices; i++) {
830
int rank = (work->vertices[i].gnum - 1) % n_ranks;
831
cs_int_t shift = send_shift[rank] + send_count[rank];
833
send_glob_buffer[shift] = (work->vertices[i].gnum - 1) / n_ranks;
834
send_count[rank] += 1;
838
MPI_Alltoallv(send_glob_buffer, send_count, send_shift, FVM_MPI_GNUM,
839
recv2glob, recv_count, recv_shift, FVM_MPI_GNUM,
842
/* Return pointers */
844
*p_block_size = block_size;
845
*p_send_count = send_count;
846
*p_recv_count = recv_count;
847
*p_send_shift = send_shift;
848
*p_recv_shift = recv_shift;
849
*p_send_glob_buf = send_glob_buffer;
850
*p_recv_glob_buf = recv_glob_buffer;
851
*p_recv2glob = recv2glob;
852
*p_glob_vtx_tag = glob_vtx_tag;
853
*p_prev_glob_vtx_tag = prev_glob_vtx_tag;
857
#endif /* HAVE_MPI */
859
/*----------------------------------------------------------------------------
860
* Tag with the same number all the vertices which might be merged together
863
* n_g_vertices_tot <-- global number of vertices to consider for the
864
* merge operation (existing + created vertices)
865
* vtx_eset <-- structure dealing with vertex equivalences
866
* work <-- local cs_join_mesh_t structure which has initial
868
* verbosity <-- level of detail in information display
869
* p_vtx_tag --> pointer to the vtx_tag for the local vertices
870
*---------------------------------------------------------------------------*/
873
_tag_equiv_vertices(fvm_gnum_t n_g_vertices_tot,
874
const cs_join_eset_t *vtx_eset,
875
const cs_join_mesh_t *work,
877
fvm_gnum_t *p_vtx_tag[])
881
fvm_gnum_t *vtx_tag = NULL;
882
fvm_gnum_t *prev_vtx_tag = NULL;
883
FILE *logfile = cs_glob_join_log;
885
const cs_int_t n_vertices = work->n_vertices;
886
const int n_ranks = cs_glob_n_ranks;
888
/* Local initialization : we tag each vertex by its global number */
890
BFT_MALLOC(prev_vtx_tag, n_vertices, fvm_gnum_t);
891
BFT_MALLOC(vtx_tag, n_vertices, fvm_gnum_t);
893
for (i = 0; i < work->n_vertices; i++) {
895
fvm_gnum_t v_gnum = work->vertices[i].gnum;
898
prev_vtx_tag[i] = v_gnum;
902
#if 0 && defined(DEBUG) && !defined(NDEBUG)
903
for (i = 0; i < n_vertices; i++)
904
fprintf(logfile, " Initial vtx_tag[%6d] = %9llu\n",
905
i, (unsigned long long)vtx_tag[i]);
909
/* Compute vtx_tag */
911
_local_spread(vtx_eset, n_vertices, prev_vtx_tag, vtx_tag);
913
if (n_ranks > 1) { /* Parallel treatment */
917
cs_int_t block_size = 0;
918
cs_int_t *send_count = NULL, *recv_count = NULL;
919
cs_int_t *send_shift = NULL, *recv_shift = NULL;
920
fvm_gnum_t *recv2glob = NULL;
921
fvm_gnum_t *recv_glob_buffer = NULL, *send_glob_buffer = NULL;
922
fvm_gnum_t *glob_vtx_tag = NULL, *prev_glob_vtx_tag = NULL;
924
#if defined(HAVE_MPI)
925
_parall_tag_init(n_g_vertices_tot,
928
&send_count, &send_shift, &send_glob_buffer,
929
&recv_count, &recv_shift, &recv_glob_buffer,
934
go_on = _global_spread(block_size,
940
send_count, send_shift, send_glob_buffer,
941
recv_count, recv_shift, recv_glob_buffer);
943
while (go_on == true) {
945
/* Local convergence of vtx_tag */
947
_local_spread(vtx_eset, n_vertices, prev_vtx_tag, vtx_tag);
949
/* Global update and test to continue */
951
go_on = _global_spread(block_size,
957
send_count, send_shift, send_glob_buffer,
958
recv_count, recv_shift, recv_glob_buffer);
964
BFT_FREE(glob_vtx_tag);
965
BFT_FREE(prev_glob_vtx_tag);
966
BFT_FREE(send_count);
967
BFT_FREE(send_shift);
968
BFT_FREE(send_glob_buffer);
969
BFT_FREE(recv_count);
970
BFT_FREE(recv_shift);
972
BFT_FREE(recv_glob_buffer);
975
} /* End of parallel treatment */
977
BFT_FREE(prev_vtx_tag);
981
"\n Number of local iterations to converge on vertex"
982
" equivalences: %3d\n", _loc_merge_counter);
985
" Number of global iterations to converge on vertex"
986
" equivalences: %3d\n\n", _glob_merge_counter);
990
#if 0 && defined(DEBUG) && !defined(NDEBUG)
991
if (logfile != NULL) {
992
for (i = 0; i < n_vertices; i++)
993
fprintf(logfile, " Final vtx_tag[%6d] = %9llu\n",
994
i, (unsigned long long)vtx_tag[i]);
999
/* Returns pointer */
1001
*p_vtx_tag = vtx_tag;
1004
#if defined(HAVE_MPI)
1006
/*----------------------------------------------------------------------------
1007
* Build in parallel a cs_join_gset_t structure to store all the potential
1008
* merges between vertices and its associated cs_join_vertex_t structure.
1011
* work <-- local cs_join_mesh_t structure which
1012
* has initial vertex data
1013
* vtx_tag <-- local vtx_tag for the local vertices
1014
* send_count <-> buffer used to count the number of elts to send
1015
* send_shift <-> index on ranks of the elements to send
1016
* recv_count <-> buffer used to count the number of elts to receive
1017
* recv_shift <-> index on ranks of the elements to receive
1018
* p_vtx_merge_data <-> a pointer to a cs_join_vertex_t structure which
1019
* stores data about merged vertices
1020
* p_merge_set <-> pointer to a cs_join_gset_t struct. storing the
1021
* evolution of each global vtx number
1022
*---------------------------------------------------------------------------*/
1025
_build_parall_merge_structures(const cs_join_mesh_t *work,
1026
const fvm_gnum_t vtx_tag[],
1027
cs_int_t send_count[],
1028
cs_int_t send_shift[],
1029
cs_int_t recv_count[],
1030
cs_int_t recv_shift[],
1031
cs_join_vertex_t *p_vtx_merge_data[],
1032
cs_join_gset_t **p_merge_set)
1036
cs_int_t n_vertices = work->n_vertices;
1037
fvm_gnum_t *recv_gbuf = NULL, *send_gbuf = NULL;
1038
cs_join_vertex_t *send_vtx_data = NULL, *recv_vtx_data = NULL;
1039
cs_join_gset_t *merge_set = NULL;
1041
MPI_Datatype CS_MPI_JOIN_VERTEX = cs_join_mesh_create_vtx_datatype();
1042
MPI_Comm mpi_comm = cs_glob_mpi_comm;
1044
const int n_ranks = cs_glob_n_ranks;
1046
for (i = 0; i < n_ranks; i++)
1049
for (i = 0; i < n_vertices; i++) {
1050
int rank = (vtx_tag[i] - 1) % n_ranks;
1051
send_count[rank] += 1;
1054
MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, mpi_comm);
1061
for (i = 0; i < n_ranks; i++) {
1062
send_shift[i+1] = send_shift[i] + send_count[i];
1063
recv_shift[i+1] = recv_shift[i] + recv_count[i];
1066
assert(send_shift[n_ranks] == n_vertices);
1068
/* Allocate and define recv_gbuf and send_gbuf */
1070
BFT_MALLOC(send_gbuf, send_shift[n_ranks], fvm_gnum_t);
1071
BFT_MALLOC(recv_gbuf, recv_shift[n_ranks], fvm_gnum_t);
1073
for (i = 0; i < n_ranks; i++)
1076
for (i = 0; i < n_vertices; i++) {
1078
int rank = (vtx_tag[i] - 1) % n_ranks;
1079
cs_int_t shift = send_shift[rank] + send_count[rank];
1081
send_gbuf[shift] = vtx_tag[i];
1082
send_count[rank] += 1;
1086
MPI_Alltoallv(send_gbuf, send_count, send_shift, FVM_MPI_GNUM,
1087
recv_gbuf, recv_count, recv_shift, FVM_MPI_GNUM,
1090
/* Allocate and build send_vtx_data, receive recv_vtx_data. */
1092
BFT_MALLOC(recv_vtx_data, recv_shift[n_ranks], cs_join_vertex_t);
1093
BFT_MALLOC(send_vtx_data, send_shift[n_ranks], cs_join_vertex_t);
1095
for (i = 0; i < n_ranks; i++)
1098
for (i = 0; i < n_vertices; i++) {
1100
int rank = (vtx_tag[i] - 1) % n_ranks;
1101
cs_int_t shift = send_shift[rank] + send_count[rank];
1103
send_vtx_data[shift] = work->vertices[i];
1104
send_count[rank] += 1;
1108
MPI_Alltoallv(send_vtx_data, send_count, send_shift, CS_MPI_JOIN_VERTEX,
1109
recv_vtx_data, recv_count, recv_shift, CS_MPI_JOIN_VERTEX,
1112
/* Partial free memory */
1114
BFT_FREE(send_vtx_data);
1115
BFT_FREE(send_gbuf);
1116
MPI_Type_free(&CS_MPI_JOIN_VERTEX);
1118
/* Build merge set */
1120
merge_set = cs_join_gset_create_from_tag(recv_shift[n_ranks], recv_gbuf);
1122
cs_join_gset_sort_sublist(merge_set);
1126
BFT_FREE(recv_gbuf);
1128
#if 0 && defined(DEBUG) && !defined(NDEBUG)
1129
if (cs_glob_join_log != NULL) {
1130
FILE *logfile = cs_glob_join_log;
1132
"\n Number of vertices to treat for the merge step: %d\n",
1133
recv_shift[n_ranks]);
1135
" List of vertices to treat:\n");
1136
for (i = 0; i < recv_shift[n_ranks]; i++) {
1137
fprintf(logfile, " %9d - ", i);
1138
cs_join_mesh_dump_vertex(logfile, recv_vtx_data[i]);
1144
/* Set return pointers */
1146
*p_merge_set = merge_set;
1147
*p_vtx_merge_data = recv_vtx_data;
1150
/*----------------------------------------------------------------------------
1151
* Exchange the updated cs_join_vertex_t array over the ranks.
1154
* work <-- local cs_join_mesh_t structure which
1155
* has initial vertex data
1156
* vtx_tag <-- local vtx_tag for the local vertices
1157
* send_count <-- buffer used to count the number of elts to send
1158
* send_shift <-- index on ranks of the elements to send
1159
* recv_count <-- buffer used to count the number of elts to receive
1160
* recv_shift <-- index on ranks of the elements to receive
1161
* vtx_merge_data <-- pointer to a cs_join_vertex_t structure
1162
*---------------------------------------------------------------------------*/
1165
_exchange_merged_vertices(const cs_join_mesh_t *work,
1166
const fvm_gnum_t vtx_tag[],
1167
cs_int_t send_count[],
1168
cs_int_t send_shift[],
1169
cs_int_t recv_count[],
1170
cs_int_t recv_shift[],
1171
cs_join_vertex_t vtx_merge_data[])
1175
cs_join_vertex_t *updated_vtx_data = NULL;
1176
int n_ranks = cs_glob_n_ranks;
1177
MPI_Datatype cs_mpi_join_vertex = cs_join_mesh_create_vtx_datatype();
1178
MPI_Comm mpi_comm = cs_glob_mpi_comm;
1180
/* Allocate send_vtx_data and exchange vtx_merge_data */
1182
BFT_MALLOC(updated_vtx_data, send_shift[n_ranks], cs_join_vertex_t);
1184
MPI_Alltoallv(vtx_merge_data, recv_count, recv_shift, cs_mpi_join_vertex,
1185
updated_vtx_data, send_count, send_shift, cs_mpi_join_vertex,
1188
/* Replace work->vertices by the updated structure after merge
1191
assert(send_shift[n_ranks] == work->n_vertices);
1193
for (i = 0; i < n_ranks; i++)
1196
for (i = 0; i < work->n_vertices; i++) {
1198
int rank = (vtx_tag[i] - 1) % n_ranks;
1199
cs_int_t shift = send_shift[rank] + send_count[rank];
1201
work->vertices[i] = updated_vtx_data[shift];
1202
send_count[rank] += 1;
1208
MPI_Type_free(&cs_mpi_join_vertex);
1209
BFT_FREE(updated_vtx_data);
1212
#endif /* HAVE_MPI */
1214
/*----------------------------------------------------------------------------
1215
* Get the resulting cs_join_vertex_t structure after the merge of a set
1219
* n_elts <-- size of the set
1220
* set <-- set of vertices
1223
* a cs_join_vertex_t structure for the resulting vertex
1224
*---------------------------------------------------------------------------*/
1226
static cs_join_vertex_t
1227
_compute_merged_vertex(cs_int_t n_elts,
1228
const cs_join_vertex_t set[])
1232
cs_join_vertex_t mvtx;
1234
cs_real_t denum = 0.0;
1238
/* Initialize cs_join_vertex_t structure */
1240
mvtx.state = CS_JOIN_STATE_UNDEF;
1241
mvtx.gnum = set[0].gnum;
1242
mvtx.tolerance = set[0].tolerance;
1244
for (k = 0; k < 3; k++)
1245
mvtx.coord[k] = 0.0;
1247
/* Compute the resulting vertex */
1249
for (i = 0; i < n_elts; i++) {
1251
mvtx.tolerance = CS_MIN(set[i].tolerance, mvtx.tolerance);
1252
mvtx.gnum = CS_MIN(set[i].gnum, mvtx.gnum);
1253
mvtx.state = CS_MAX(set[i].state, mvtx.state);
1255
/* Compute the resulting coordinates of the merged vertices */
1257
#if CS_JOIN_MERGE_INV_TOL
1258
w = 1.0/set[i].tolerance;
1264
for (k = 0; k < 3; k++)
1265
mvtx.coord[k] += w * set[i].coord[k];
1269
for (k = 0; k < 3; k++)
1270
mvtx.coord[k] /= denum;
1272
if (mvtx.state == CS_JOIN_STATE_ORIGIN)
1273
mvtx.state = CS_JOIN_STATE_MERGE;
1274
else if (mvtx.state == CS_JOIN_STATE_PERIO)
1275
mvtx.state = CS_JOIN_STATE_PERIO_MERGE;
1280
/*----------------------------------------------------------------------------
1281
* Merge between identical vertices.
1283
* Only the vertex numbering and the related tolerance may be different.
1284
* Store new data associated to the merged vertices in vertices array.
1287
* param <-- set of user-defined parameters
1288
* merge_set <-> a pointer to a cs_join_vertex_t structure which
1289
* stores data about merged vertices
1290
* n_vertices <-- number of vertices in vertices array
1291
* vertices <-> array of cs_join_vertex_t structures
1292
* equiv_gnum --> equivalence between id in vertices (same global number
1293
* initially or identical vertices: same coordinates)
1294
*---------------------------------------------------------------------------*/
1297
_pre_merge(cs_join_param_t param,
1298
cs_join_gset_t *merge_set,
1299
cs_join_vertex_t vertices[],
1300
cs_join_gset_t **p_equiv_gnum)
1302
cs_int_t i, j, j1, j2, k, k1, k2, n_sub_elts;
1303
cs_real_t deltad, deltat, limit, min_tol;
1304
cs_join_vertex_t mvtx, coupled_vertices[2];
1306
cs_int_t max_n_sub_elts = 0, n_local_pre_merge = 0;
1307
cs_int_t *merge_index = merge_set->index;
1308
fvm_gnum_t *merge_list = merge_set->g_list;
1309
fvm_gnum_t *sub_list = NULL, *init_list = NULL;
1310
cs_join_gset_t *equiv_gnum = NULL;
1312
const cs_real_t pmf = param.pre_merge_factor;
1314
cs_join_gset_sort_sublist(merge_set);
1316
#if 0 && defined(DEBUG) && !defined(NDEBUG)
1319
FILE *dbg_file = NULL;
1320
char *filename = NULL;
1322
len = strlen("JoinDBG_InitMergeSet.dat")+1+2+4;
1323
BFT_MALLOC(filename, len, char);
1324
sprintf(filename, "Join%02dDBG_InitMergeSet%04d.dat",
1325
param.num, CS_MAX(cs_glob_rank_id, 0));
1326
dbg_file = fopen(filename, "w");
1328
cs_join_gset_dump(dbg_file, merge_set);
1336
/* Compute the max. size of a sub list */
1338
for (i = 0; i < merge_set->n_elts; i++)
1339
max_n_sub_elts = CS_MAX(max_n_sub_elts,
1340
merge_index[i+1] - merge_index[i]);
1342
BFT_MALLOC(sub_list, max_n_sub_elts, fvm_gnum_t);
1344
/* Store initial merge list */
1346
BFT_MALLOC(init_list, merge_index[merge_set->n_elts], fvm_gnum_t);
1348
for (i = 0; i < merge_index[merge_set->n_elts]; i++)
1349
init_list[i] = merge_list[i];
1353
for (i = 0; i < merge_set->n_elts; i++) {
1355
cs_int_t f_s = merge_index[i];
1356
cs_int_t f_e = merge_index[i+1];
1358
n_sub_elts = f_e - f_s;
1360
for (j = f_s, k = 0; j < f_e; j++, k++)
1361
sub_list[k] = merge_list[j];
1363
for (j1 = 0; j1 < n_sub_elts - 1; j1++) {
1365
cs_int_t v1_id = sub_list[j1];
1366
cs_join_vertex_t v1 = vertices[v1_id];
1368
for (j2 = j1 + 1; j2 < n_sub_elts; j2++) {
1370
cs_int_t v2_id = sub_list[j2];
1371
cs_join_vertex_t v2 = vertices[v2_id];
1373
if (v1.gnum == v2.gnum) { /* Possible if n_ranks > 1 */
1375
if (sub_list[j1] < sub_list[j2])
1380
for (k = 0; k < n_sub_elts; k++)
1381
if (sub_list[k] == sub_list[k2])
1382
sub_list[k] = sub_list[k1];
1387
min_tol = CS_MIN(v1.tolerance, v2.tolerance);
1388
limit = min_tol * pmf;
1389
deltat = CS_ABS(v1.tolerance - v2.tolerance);
1391
if (deltat < limit) {
1393
deltad = _compute_length(v1, v2);
1395
if (deltad < limit) { /* Do a pre-merge */
1397
n_local_pre_merge++;
1399
if (v1.gnum < v2.gnum)
1404
for (k = 0; k < n_sub_elts; k++)
1405
if (sub_list[k] == sub_list[k2])
1406
sub_list[k] = sub_list[k1];
1408
coupled_vertices[0] = v1, coupled_vertices[1] = v2;
1409
mvtx = _compute_merged_vertex(2, coupled_vertices);
1410
vertices[v1_id] = mvtx;
1411
vertices[v2_id] = mvtx;
1413
} /* deltad < limit */
1415
} /* deltat < limit */
1417
} /* v1.gnum != v2.gnum */
1419
} /* End of loop on j2 */
1420
} /* End of loop on j1 */
1422
/* Update vertices */
1424
for (j = f_s, k = 0; j < f_e; j++, k++)
1425
vertices[merge_list[j]] = vertices[sub_list[k]];
1427
/* Update merge list */
1429
for (j = f_s, k = 0; j < f_e; j++, k++)
1430
merge_list[j] = sub_list[k];
1432
} /* End of loop on merge_set elements */
1434
/* Keep equivalences between identical vertices in equiv_gnum */
1436
equiv_gnum = cs_join_gset_create_by_equiv(merge_set, init_list);
1438
/* Clean merge set */
1440
cs_join_gset_clean(merge_set);
1442
/* Display information about the joining */
1444
if (param.verbosity > 0) {
1446
fvm_gnum_t n_g_counter = n_local_pre_merge;
1447
fvm_parall_counter(&n_g_counter, 1);
1449
bft_printf(_("\n Pre-merge for %llu global element couples.\n"),
1450
(unsigned long long)n_g_counter);
1452
if (param.verbosity > 2) {
1453
fprintf(cs_glob_join_log, "\n Local number of pre-merges: %d\n",
1461
BFT_FREE(init_list);
1463
/* Return pointer */
1465
*p_equiv_gnum = equiv_gnum;
1468
/*----------------------------------------------------------------------------
1469
* Check if all vertices in the set include the ref_vertex in their tolerance.
1472
* set_size <-- size of set of vertices
1473
* vertices <-- set of vertices to check
1474
* ref_vertex <-- ref. vertex
1477
* true if all vertices have ref_vertex in their tolerance, false otherwise
1478
*---------------------------------------------------------------------------*/
1481
_is_in_tolerance(cs_int_t set_size,
1482
const cs_join_vertex_t set[],
1483
cs_join_vertex_t ref_vertex)
1487
for (i = 0; i < set_size; i++) {
1489
cs_real_t d2ref = _compute_length(set[i], ref_vertex);
1490
cs_real_t tolerance = set[i].tolerance * cs_join_tol_eps_coef2;
1492
if (d2ref > tolerance)
1500
/*----------------------------------------------------------------------------
1501
* Test if we have to continue to the subset building
1504
* set_size <-- size of set
1505
* prev_num <-> array used to store previous subset_num
1506
* new_num <-> number associated to each vertices of the set
1510
*---------------------------------------------------------------------------*/
1513
_continue_subset_building(int set_size,
1514
const cs_int_t prev_num[],
1515
const cs_int_t new_num[])
1519
for (i = 0; i < set_size; i++)
1520
if (new_num[i] != prev_num[i])
1526
/*----------------------------------------------------------------------------
1527
* Define subsets of vertices.
1530
* set_size <-- size of set
1531
* state <-- array keeping the state of the link
1532
* subset_num <-> number associated to each vertices of the set
1533
*---------------------------------------------------------------------------*/
1536
_iter_subset_building(cs_int_t set_size,
1537
const cs_int_t state[],
1538
cs_int_t subset_num[])
1542
for (k = 0, i1 = 0; i1 < set_size-1; i1++) {
1543
for (i2 = i1 + 1; i2 < set_size; i2++, k++) {
1545
if (state[k] == 1) { /* v1 - v2 are in tolerance each other */
1547
int _min = CS_MIN(subset_num[i1], subset_num[i2]);
1549
subset_num[i1] = _min;
1550
subset_num[i2] = _min;
1559
/*----------------------------------------------------------------------------
1560
* Define subsets of vertices.
1563
* set_size <-- size of set
1564
* state <-- array keeping the state of the link
1565
* prev_num <-> array used to store previous subset_num
1566
* subset_num <-> number associated to each vertices of the set
1567
*---------------------------------------------------------------------------*/
1570
_build_subsets(cs_int_t set_size,
1571
const cs_int_t state[],
1572
cs_int_t prev_num[],
1573
cs_int_t subset_num[])
1576
cs_int_t n_loops = 0;
1580
for (i = 0; i < set_size; i++) {
1581
subset_num[i] = i+1;
1582
prev_num[i] = subset_num[i];
1585
_iter_subset_building(set_size, state, subset_num);
1587
while ( _continue_subset_building(set_size, prev_num, subset_num)
1588
&& n_loops < CS_JOIN_MERGE_MAX_LOC_ITERS ) {
1592
for (i = 0; i < set_size; i++)
1593
prev_num[i] = subset_num[i];
1595
_iter_subset_building(set_size, state, subset_num);
1599
#if 0 && defined(DEBUG) && !defined(NDEBUG)
1600
if (cs_glob_join_log != NULL && n_loops >= CS_JOIN_MERGE_MAX_LOC_ITERS)
1601
fprintf(cs_glob_join_log,
1602
"WARNING max sub_loops to build subset reached\n");
1607
/*----------------------------------------------------------------------------
1608
* Check if each subset is consistent with tolerance of vertices
1609
* If a transitivity is found, modify the state of the link
1610
* state = 1 (each other in their tolerance)
1611
* = 0 (not each other in their tolerance)
1614
* set_size <-- size of set
1615
* set <-- pointer to a set of vertices
1616
* state <-> array keeping the state of the link
1617
* subset_num <-> number associated to each vertices of the set
1618
* issues <-> numbering of inconsistent subsets
1619
* verbosity <-- level of verbosity
1622
* number of subsets not consistent
1623
*---------------------------------------------------------------------------*/
1626
_check_tol_consistency(cs_int_t set_size,
1627
const cs_join_vertex_t set[],
1629
cs_int_t subset_num[],
1633
cs_int_t i1, i2, j, k;
1635
cs_int_t n_issues = 0;
1636
FILE *logfile = cs_glob_join_log;
1638
for (k = 0, i1 = 0; i1 < set_size-1; i1++) {
1639
for (i2 = i1 + 1; i2 < set_size; i2++, k++) {
1641
if (state[k] == 0) {
1643
if (subset_num[i1] == subset_num[i2]) {
1647
" Transitivity detected between (%llu, %llu)\n",
1648
(unsigned long long)set[i1].gnum,
1649
(unsigned long long)set[i2].gnum);
1651
for (j = 0; j < n_issues; j++)
1652
if (issues[j] == subset_num[i1])
1655
issues[n_issues++] = subset_num[i1];
1660
} /* End of loop on i2 */
1661
} /* ENd of loop on i1 */
1663
return n_issues; /* Not a subset number */
1666
/*----------------------------------------------------------------------------
1667
* Check if the merged vertex related to a subset is consistent with tolerance
1668
* of each vertex of the subset.
1671
* set_size <-- size of set
1672
* subset_num <-- number associated to each vertices of the set
1673
* set <-- pointer to a set of vertices
1674
* merge_set <-> merged vertex related to each subset
1675
* work_set <-> work array of vertices
1678
* true if all subsets are consistent otherwise false
1679
*---------------------------------------------------------------------------*/
1682
_check_subset_consistency(cs_int_t set_size,
1683
const cs_int_t subset_num[],
1684
const cs_join_vertex_t set[],
1685
cs_join_vertex_t merge_set[],
1686
cs_join_vertex_t work_set[])
1688
cs_int_t i, set_id, subset_size;
1690
cs_bool_t is_consistent = true;
1692
/* Apply merged to each subset */
1694
for (set_id = 0; set_id < set_size; set_id++) {
1697
for (i = 0; i < set_size; i++)
1698
if (subset_num[i] == set_id+1)
1699
work_set[subset_size++] = set[i];
1701
if (subset_size > 0) {
1703
merge_set[set_id] = _compute_merged_vertex(subset_size, work_set);
1705
if (!_is_in_tolerance(subset_size, work_set, merge_set[set_id]))
1706
is_consistent = false;
1710
} /* End of loop on subsets */
1712
return is_consistent;
1715
/*----------------------------------------------------------------------------
1716
* Get position of the link between vertices i1 and i2.
1719
* i1 <-- id in set for vertex 1
1720
* i2 <-- id in set for vertex 2
1721
* idx <-- array of positions
1724
* position in an array like distances or state
1725
*---------------------------------------------------------------------------*/
1727
inline static cs_int_t
1728
_get_pos(cs_int_t i1,
1730
const cs_int_t idx[])
1735
pos = idx[i1] + i2-i1-1;
1738
pos = idx[i2] + i1-i2-1;
1744
/*----------------------------------------------------------------------------
1745
* Break equivalences for vertices implied in transitivity issue
1748
* param <-- parameters used to manage the joining algorithm
1749
* set_size <-- size of set
1750
* set <-- pointer to a set of vertices
1751
* state <-> array keeping the state of the link
1752
* n_issues <-- number of detected transitivity issues
1753
* issues <-- subset numbers of subset with a transitivity issue
1754
* idx <-- position of vertices couple in array like distances
1755
* subset_num <-- array of subset numbers
1756
* distances <-- array keeping the distances between vertices
1757
*---------------------------------------------------------------------------*/
1760
_break_equivalence(cs_join_param_t param,
1762
const cs_join_vertex_t set[],
1767
const int subset_num[],
1768
const double distances[])
1770
cs_int_t i, i1, i2, k;
1772
for (i = 0; i < n_issues; i++) {
1774
/* Find the weakest equivalence and break it.
1775
Purpose: Have the minimal number of equivalences to break
1776
for each subset where an inconsistency has been detected */
1779
double rtf = -1.0, dist_save = 0.0;
1781
for (k = 0, i1 = 0; i1 < set_size-1; i1++) {
1782
for (i2 = i1 + 1; i2 < set_size; i2++, k++) {
1784
if (state[k] == 1) { /* v1, v2 are equivalent */
1786
if ( subset_num[i1] == issues[i]
1787
&& subset_num[i2] == issues[i]) {
1789
/* Vertices belongs to a subset where an inconsistency
1792
double rtf12 = distances[k]/set[i1].tolerance;
1793
double rtf21 = distances[k]/set[i2].tolerance;
1795
assert(rtf12 < 1.0); /* Because they are equivalent */
1796
assert(rtf21 < 1.0);
1798
if (rtf12 >= rtf21) {
1802
dist_save = distances[k];
1809
dist_save = distances[k];
1816
} /* End of loop on i1 */
1817
} /* End of loop on i2 */
1821
/* Break equivalence between i_save and all vertices linked to
1822
i_save with a distance to i_save >= dist_save */
1824
for (i2 = 0; i2 < set_size; i2++) {
1828
k = _get_pos(i_save, i2, idx);
1829
if (distances[k] >= dist_save && state[k] == 1) {
1831
state[k] = 0; /* Break equivalence */
1833
if (param.verbosity > 3)
1834
fprintf(cs_glob_join_log,
1835
" %2d - Break equivalence between [%llu, %llu]"
1836
" (dist_ref: %6.4e)\n",
1838
(unsigned long long)set[i_save].gnum,
1839
(unsigned long long)set[i2].gnum, dist_save);
1844
} /* End of loop on vertices */
1848
} /* End of loop on issues */
1852
/*----------------------------------------------------------------------------
1853
* Break equivalences between vertices until each vertex of the list has
1854
* the resulting vertex of the merge under its tolerance.
1857
* param <-- set of user-defined parameters
1858
* set_size <-- size of the set of vertices
1859
* set <-> set of vertices
1860
* vbuf <-> tmp buffer
1861
* rbuf <-> tmp buffer
1862
* ibuf <-> tmp buffer
1865
* number of loops necessary to build consistent subsets
1866
*---------------------------------------------------------------------------*/
1869
_solve_transitivity(cs_join_param_t param,
1871
cs_join_vertex_t set[],
1872
cs_join_vertex_t vbuf[],
1876
cs_int_t i1, i2, k, n_issues;
1878
cs_int_t n_loops = 0;
1879
cs_bool_t is_end = false;
1880
cs_int_t *subset_num = NULL, *state = NULL, *prev_num = NULL;
1881
cs_int_t *subset_issues = NULL, *idx = NULL;
1882
cs_real_t *distances = NULL;
1883
cs_join_vertex_t *merge_set = NULL, *work_set = NULL;
1887
assert(set_size > 0);
1889
/* Define temporary buffers */
1891
subset_num = &(ibuf[0]);
1892
prev_num = &(ibuf[set_size]);
1893
subset_issues = &(ibuf[2*set_size]);
1894
idx = &(ibuf[3*set_size]);
1895
state = &(ibuf[4*set_size]);
1896
distances = &(rbuf[0]);
1897
merge_set = &(vbuf[0]);
1898
work_set = &(vbuf[set_size]);
1900
/* Compute distances between each couple of vertices among the set */
1902
for (k = 0, i1 = 0; i1 < set_size-1; i1++)
1903
for (i2 = i1 + 1; i2 < set_size; i2++, k++)
1904
distances[k] = _compute_length(set[i1], set[i2]);
1906
/* Compute initial state of equivalences between vertices */
1908
for (k = 0, i1 = 0; i1 < set_size-1; i1++) {
1909
for (i2 = i1 + 1; i2 < set_size; i2++, k++) {
1910
if ( set[i1].tolerance < distances[k]
1911
|| set[i2].tolerance < distances[k])
1919
for (k = 1; k < set_size - 1; k++)
1920
idx[k] = set_size - k + idx[k-1];
1922
#if 0 && defined(DEBUG) && !defined(NDEBUG)
1923
if (cs_glob_join_log != NULL) {
1924
cs_join_dump_array(cs_glob_join_log, "double", "\nDist",
1925
set_size*(set_size-1)/2, distances);
1926
cs_join_dump_array(cs_glob_join_log, "int", "\nInit. State",
1927
set_size*(set_size-1)/2, state);
1931
_build_subsets(set_size, state, prev_num, subset_num);
1933
while (is_end == false && n_loops < param.n_max_equiv_breaks) {
1937
n_issues = _check_tol_consistency(set_size,
1945
_break_equivalence(param,
1955
_build_subsets(set_size, state, prev_num, subset_num);
1957
is_end = _check_subset_consistency(set_size,
1963
} /* End of while */
1965
if (param.verbosity > 3) {
1967
fprintf(cs_glob_join_log,
1968
" Number of tolerance reductions: %4d\n", n_loops);
1970
#if 0 && defined(DEBUG) && !defined(NDEBUG)
1971
cs_join_dump_array(cs_glob_join_log, "int", "\nFinal Subset",
1972
set_size, subset_num);
1977
/* Apply merged to each subset */
1979
for (k = 0; k < set_size; k++)
1980
set[k] = merge_set[subset_num[k]-1];
1985
/*----------------------------------------------------------------------------
1986
* Merge between vertices. Store new data associated to the merged vertices
1990
* param <-- set of user-defined parameters
1991
* merge_set <-> a pointer to a cs_join_vertex_t structure which
1992
* stores data about merged vertices
1993
* n_vertices <-- number of vertices in vertices array
1994
* vertices <-> array of cs_join_vertex_t structures
1995
*---------------------------------------------------------------------------*/
1998
_merge_vertices(cs_join_param_t param,
1999
cs_join_gset_t *merge_set,
2000
cs_int_t n_vertices,
2001
cs_join_vertex_t vertices[])
2003
cs_int_t i, j, k, list_size;
2004
cs_join_vertex_t merged_vertex;
2007
cs_int_t max_list_size = 0, vv_max_list_size = 0;
2008
cs_int_t n_loops = 0, n_max_loops = 0, n_transitivity = 0;
2010
cs_join_gset_t *equiv_gnum = NULL;
2011
cs_real_t *rbuf = NULL;
2012
cs_int_t *merge_index = NULL, *ibuf = NULL;
2013
fvm_gnum_t *merge_list = NULL, *merge_ref_elts = NULL;
2014
fvm_gnum_t *list = NULL;
2015
cs_join_vertex_t *set = NULL, *vbuf = NULL;
2016
FILE *logfile = cs_glob_join_log;
2018
const int verbosity = param.verbosity;
2022
assert(param.merge_tol_coef >= 0.0);
2024
/* Pre-merge of identical vertices */
2026
_pre_merge(param, merge_set, vertices, &equiv_gnum);
2028
#if 0 && defined(DEBUG) && !defined(NDEBUG)
2031
FILE *dbg_file = NULL;
2032
char *filename = NULL;
2034
len = strlen("JoinDBG_MergeSet.dat")+1+2+4;
2035
BFT_MALLOC(filename, len, char);
2036
sprintf(filename, "Join%02dDBG_MergeSet%04d.dat",
2037
param.num, CS_MAX(cs_glob_rank_id, 0));
2038
dbg_file = fopen(filename, "w");
2040
cs_join_gset_dump(dbg_file, merge_set);
2046
#endif /* defined(DEBUG) && !defined(NDEBUG) */
2048
/* Modify the tolerance for the merge operation if needed */
2050
if (fabs(param.merge_tol_coef - 1.0) > 1e-30) {
2051
for (i = 0; i < n_vertices; i++)
2052
vertices[i].tolerance *= param.merge_tol_coef;
2055
/* Compute the max. size of a sub list */
2057
merge_index = merge_set->index;
2058
merge_list = merge_set->g_list;
2059
merge_ref_elts = merge_set->g_elts;
2061
for (i = 0; i < merge_set->n_elts; i++) {
2062
list_size = merge_index[i+1] - merge_index[i];
2063
max_list_size = CS_MAX(max_list_size, list_size);
2065
vv_max_list_size = ((max_list_size-1)*max_list_size)/2;
2067
if (verbosity > 0) { /* Display information */
2069
fvm_lnum_t g_max_list_size = max_list_size;
2070
fvm_parall_counter_max(&g_max_list_size, 1);
2072
if (g_max_list_size < 2) {
2073
cs_join_gset_destroy(&equiv_gnum);
2074
bft_printf(_("\n No need to merge vertices.\n"));
2078
bft_printf(_("\n Max size of a merge set of vertices: %llu\n"),
2079
(unsigned long long)g_max_list_size);
2082
/* Temporary buffers allocation */
2084
BFT_MALLOC(ibuf, 4*max_list_size + vv_max_list_size, cs_int_t);
2085
BFT_MALLOC(rbuf, vv_max_list_size, cs_real_t);
2086
BFT_MALLOC(vbuf, 2*max_list_size, cs_join_vertex_t);
2087
BFT_MALLOC(list, max_list_size, fvm_gnum_t);
2088
BFT_MALLOC(set, max_list_size, cs_join_vertex_t);
2090
/* Merge set of vertices */
2092
for (i = 0; i < merge_set->n_elts; i++) {
2094
list_size = merge_index[i+1] - merge_index[i];
2096
if (list_size > 1) {
2098
for (j = 0, k = merge_index[i]; k < merge_index[i+1]; k++, j++) {
2099
list[j] = merge_list[k];
2100
set[j] = vertices[list[j]];
2103
/* Define the resulting cs_join_vertex_t structure of the merge */
2105
merged_vertex = _compute_merged_vertex(list_size, set);
2107
/* Check if the vertex resulting of the merge is in the tolerance
2108
for each vertex of the list */
2110
ok = _is_in_tolerance(list_size, set, merged_vertex);
2112
#if CS_JOIN_MERGE_TOL_REDUC
2113
if (ok == false) { /*
2114
The merged vertex is not in the tolerance of
2115
each vertex. This is a transitivity problem.
2116
We have to split the initial set into several
2122
/* Display information on vertices to merge */
2123
if (verbosity > 3) {
2125
"\n Begin merge for ref. elt: %llu - list_size: %d\n",
2126
(unsigned long long)merge_ref_elts[i],
2127
merge_index[i+1] - merge_index[i]);
2128
for (j = 0; j < list_size; j++) {
2129
fprintf(logfile, "%9llu -", (unsigned long long)list[j]);
2130
cs_join_mesh_dump_vertex(logfile, set[j]);
2132
fprintf(logfile, "\nMerged vertex rejected:\n");
2133
cs_join_mesh_dump_vertex(logfile, merged_vertex);
2136
n_loops = _solve_transitivity(param,
2143
for (j = 0; j < list_size; j++)
2144
vertices[list[j]] = set[j];
2146
n_max_loops = CS_MAX(n_max_loops, n_loops);
2148
if (verbosity > 3) { /* Display information */
2149
fprintf(logfile, "\n %3d loop(s) to get consistent subsets\n",
2151
fprintf(logfile, "\n End merge for ref. elt: %llu - list_size: %d\n",
2152
(unsigned long long)merge_ref_elts[i],
2153
merge_index[i+1] - merge_index[i]);
2154
for (j = 0; j < list_size; j++) {
2155
fprintf(logfile, "%7llu -", (unsigned long long)list[j]);
2156
cs_join_mesh_dump_vertex(logfile, vertices[list[j]]);
2158
fprintf(logfile, "\n");
2162
else /* New vertex data for the sub-elements */
2164
#endif /* CS_JOIN_MERGE_TOL_REDUC */
2166
for (j = 0; j < list_size; j++)
2167
vertices[list[j]] = merged_vertex;
2169
} /* list_size > 1 */
2171
} /* End of loop on potential merges */
2173
/* Apply merge to vertex initially identical */
2175
if (equiv_gnum != NULL) {
2177
#if 0 && defined(DEBUG) && !defined(NDEBUG)
2180
FILE *dbg_file = NULL;
2181
char *filename = NULL;
2183
len = strlen("JoinDBG_EquivMerge.dat")+1+2+4;
2184
BFT_MALLOC(filename, len, char);
2185
sprintf(filename, "Join%02dDBG_EquivMerge%04d.dat",
2186
param.num, CS_MAX(cs_glob_rank_id, 0));
2187
dbg_file = fopen(filename, "w");
2189
cs_join_gset_dump(dbg_file, equiv_gnum);
2195
#endif /* defined(DEBUG) && !defined(NDEBUG) */
2197
for (i = 0; i < equiv_gnum->n_elts; i++) {
2199
cs_int_t start = equiv_gnum->index[i];
2200
cs_int_t end = equiv_gnum->index[i+1];
2201
cs_int_t ref_id = equiv_gnum->g_elts[i];
2203
for (j = start; j < end; j++)
2204
vertices[equiv_gnum->g_list[j]] = vertices[ref_id];
2209
if (verbosity > 0) {
2211
fvm_gnum_t n_g_counter = n_transitivity;
2212
fvm_parall_counter(&n_g_counter, 1);
2214
bft_printf(_("\n Excessive transitivity for %llu set(s) of vertices.\n"),
2215
(unsigned long long)n_g_counter);
2217
if (verbosity > 1) {
2218
fvm_lnum_t g_n_max_loops = n_max_loops;
2219
fvm_parall_counter_max(&g_n_max_loops, 1);
2220
bft_printf(_("\n Max. number of iterations to solve transitivity"
2221
" excess: %llu\n"), (unsigned long long)g_n_max_loops);
2233
cs_join_gset_destroy(&equiv_gnum);
2236
/*----------------------------------------------------------------------------
2237
* Keep an history of the evolution of each vertex id before/after the merge
2241
* n_iwm_vertices <-- number of vertices before intersection for the
2242
* work cs_join_mesh_t structure
2243
* iwm_vtx_gnum <-- initial global vertex num. (work mesh struct.)
2244
* init_max_vtx_gnum <-- initial max. global numbering for vertices
2245
* n_vertices <-- number of vertices before merge/after intersection
2246
* vertices <-- array of cs_join_vertex_t structures
2247
* p_o2n_vtx_gnum --> distributed array by block on the new global vertex
2248
* numbering for the initial vertices (before inter.)
2249
*---------------------------------------------------------------------------*/
2252
_keep_global_vtx_evolution(cs_int_t n_iwm_vertices,
2253
const fvm_gnum_t iwm_vtx_gnum[],
2254
fvm_gnum_t init_max_vtx_gnum,
2255
cs_int_t n_vertices,
2256
const cs_join_vertex_t vertices[],
2257
fvm_gnum_t *p_o2n_vtx_gnum[])
2261
fvm_gnum_t *o2n_vtx_gnum = NULL;
2263
const int n_ranks = cs_glob_n_ranks;
2264
const int local_rank = CS_MAX(cs_glob_rank_id, 0);
2266
assert(n_iwm_vertices <= n_vertices); /* after inter. >= init */
2270
BFT_MALLOC(o2n_vtx_gnum, n_iwm_vertices, fvm_gnum_t);
2272
for (i = 0; i < n_iwm_vertices; i++)
2273
o2n_vtx_gnum[i] = vertices[i].gnum;
2277
#if defined(HAVE_MPI) /* Parallel treatment */
2281
cs_int_t shift, rank, n_recv_elts;
2283
cs_int_t *send_shift = NULL, *recv_shift = NULL;
2284
cs_int_t *send_count = NULL, *recv_count = NULL;
2285
fvm_gnum_t *send_glist = NULL, *recv_glist = NULL;
2287
cs_join_block_info_t block_info = cs_join_get_block_info(init_max_vtx_gnum,
2291
MPI_Comm mpi_comm = cs_glob_mpi_comm;
2293
/* Initialize o2n_vtx_gnum */
2295
BFT_MALLOC(o2n_vtx_gnum, block_info.local_size, fvm_gnum_t);
2297
for (ii = 0; ii < block_info.local_size; ii++)
2298
o2n_vtx_gnum[ii] = block_info.first_gnum + ii;
2300
/* Send new vtx global number to the related rank = the good block */
2302
BFT_MALLOC(send_count, n_ranks, cs_int_t);
2303
BFT_MALLOC(recv_count, n_ranks, cs_int_t);
2305
for (i = 0; i < n_ranks; i++)
2308
for (i = 0; i < n_iwm_vertices; i++) {
2309
rank = (iwm_vtx_gnum[i] - 1)/block_info.size;
2310
send_count[rank] += 2;
2313
MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, mpi_comm);
2315
BFT_MALLOC(send_shift, n_ranks + 1, cs_int_t);
2316
BFT_MALLOC(recv_shift, n_ranks + 1, cs_int_t);
2321
for (rank = 0; rank < n_ranks; rank++) {
2322
send_shift[rank + 1] = send_shift[rank] + send_count[rank];
2323
recv_shift[rank + 1] = recv_shift[rank] + recv_count[rank];
2326
assert(send_shift[n_ranks] == 2*n_iwm_vertices);
2328
/* Build send_list */
2330
BFT_MALLOC(send_glist, send_shift[n_ranks], fvm_gnum_t);
2331
BFT_MALLOC(recv_glist, recv_shift[n_ranks], fvm_gnum_t);
2333
for (i = 0; i < n_ranks; i++)
2336
for (i = 0; i < n_iwm_vertices; i++) {
2338
rank = (iwm_vtx_gnum[i] - 1)/block_info.size;
2339
shift = send_shift[rank] + send_count[rank];
2341
send_glist[shift] = iwm_vtx_gnum[i]; /* Old global number */
2342
send_glist[shift+1] = vertices[i].gnum; /* New global number */
2343
send_count[rank] += 2;
2347
MPI_Alltoallv(send_glist, send_count, send_shift, FVM_MPI_GNUM,
2348
recv_glist, recv_count, recv_shift, FVM_MPI_GNUM,
2351
n_recv_elts = recv_shift[n_ranks]/2;
2353
BFT_FREE(send_count);
2354
BFT_FREE(send_shift);
2355
BFT_FREE(send_glist);
2356
BFT_FREE(recv_count);
2358
/* Update o2n_vtx_gnum */
2360
for (rank = 0; rank < n_ranks; rank++) {
2362
for (i = recv_shift[rank]; i < recv_shift[rank+1]; i+=2) {
2364
fvm_gnum_t o_gnum = recv_glist[i];
2365
fvm_gnum_t n_gnum = recv_glist[i+1];
2366
cs_int_t id = o_gnum - block_info.first_gnum;
2368
#if 0 && defined(DEBUG) && !defined(NDEBUG)
2369
if (o2n_vtx_gnum[id] != block_info.first_gnum + id)
2370
assert(o2n_vtx_gnum[id] == n_gnum);
2373
o2n_vtx_gnum[id] = n_gnum;
2377
} /* End of loop on ranks */
2379
BFT_FREE(recv_shift);
2380
BFT_FREE(recv_glist);
2383
#endif /* HAVE_MPI */
2385
/* Set return pointer */
2387
*p_o2n_vtx_gnum = o2n_vtx_gnum;
2390
/*----------------------------------------------------------------------------
2391
* Keep a history of the evolution of each vertex id before/after the merge
2392
* operation for the current mesh (local point of view).
2395
* n_vertices <-- number of vertices before merge/after intersection
2396
* vertices <-- array of cs_join_vertex_t structures
2397
* p_n_am_vertices --> number of vertices after the merge step
2398
* p_o2n_vtx_id --> array keeping the evolution of the vertex ids
2399
*---------------------------------------------------------------------------*/
2402
_keep_local_vtx_evolution(cs_int_t n_vertices,
2403
const cs_join_vertex_t vertices[],
2404
cs_int_t *p_n_am_vertices,
2405
cs_int_t *p_o2n_vtx_id[])
2410
cs_int_t n_am_vertices = 0;
2411
cs_int_t *o2n_vtx_id = NULL;
2412
fvm_lnum_t *order = NULL;
2413
fvm_gnum_t *vtx_gnum = NULL;
2415
if (n_vertices == 0)
2418
BFT_MALLOC(vtx_gnum, n_vertices, fvm_gnum_t);
2420
for (i = 0; i < n_vertices; i++)
2421
vtx_gnum[i] = vertices[i].gnum;
2423
/* Order vertices according to their global numbering */
2425
BFT_MALLOC(order, n_vertices, fvm_lnum_t);
2427
fvm_order_local_allocated(NULL, vtx_gnum, order, n_vertices);
2429
/* Delete vertices sharing the same global number. Keep only one */
2431
BFT_MALLOC(o2n_vtx_id, n_vertices, cs_int_t);
2433
prev = vtx_gnum[order[0]];
2434
o2n_vtx_id[order[0]] = n_am_vertices;
2436
for (i = 1; i < n_vertices; i++) {
2438
cs_int_t o_id = order[i];
2439
fvm_gnum_t cur = vtx_gnum[o_id];
2444
o2n_vtx_id[o_id] = n_am_vertices;
2447
o2n_vtx_id[o_id] = n_am_vertices;
2449
} /* End of loop on vertices */
2451
/* n_am_vertices is an id */
2454
assert(n_am_vertices <= n_vertices); /* after merge <= after inter. */
2461
/* Set return pointers */
2463
*p_n_am_vertices = n_am_vertices;
2464
*p_o2n_vtx_id = o2n_vtx_id;
2467
/*----------------------------------------------------------------------------
2468
* Search for new elements to add to the definition of the current edge
2469
* These new sub-elements come from initial edges which are now (after the
2470
* merge step) sub-edge of the current edge
2474
* edge_id <-- id of the edge to deal with
2475
* inter_edges <-- structure keeping data on edge intersections
2476
* edges <-- edges definition
2477
* n_iwm_vertices <-- initial number of vertices in work_mesh
2478
* n_new_sub_elts --> number of new elements to add in the edge definition
2481
* number of new sub-elements related to this edge
2482
*---------------------------------------------------------------------------*/
2485
_count_new_sub_edge_elts(cs_int_t edge_id,
2486
const cs_join_inter_edges_t *inter_edges,
2487
const cs_join_edges_t *edges,
2488
cs_int_t n_iwm_vertices)
2490
cs_int_t j, k, j1, j2, sub_edge_id;
2491
cs_int_t start, end, _start, _end, v1_num, v2_num;
2494
cs_int_t n_new_sub_elts = 0;
2496
start = inter_edges->index[edge_id];
2497
end = inter_edges->index[edge_id+1];
2499
for (j1 = start; j1 < end-1; j1++) {
2501
v1_num = inter_edges->vtx_lst[j1];
2503
if (v1_num <= n_iwm_vertices) { /* v1 is an initial vertex */
2504
for (j2 = j1+1; j2 < end; j2++) {
2506
v2_num = inter_edges->vtx_lst[j2];
2508
if (v2_num <= n_iwm_vertices) { /* (v1,v2) is an initial edge */
2510
sub_edge_id = CS_ABS(cs_join_mesh_get_edge(v1_num,
2513
assert(sub_edge_id != -1);
2514
_start = inter_edges->index[sub_edge_id];
2515
_end = inter_edges->index[sub_edge_id+1];
2517
for (j = _start; j < _end; j++) {
2520
for (k = j1+1; k < j2; k++)
2521
if (inter_edges->vtx_glst[k] == inter_edges->vtx_glst[j])
2525
n_new_sub_elts += 1;
2527
} /* End of loop on sub_edge_id definition */
2531
} /* End of loop on j2 */
2534
} /* End of loop on j1 */
2536
return n_new_sub_elts;
2539
/*----------------------------------------------------------------------------
2540
* Update a cs_join_inter_edges_t structure after the merge operation.
2541
* cs_join_inter_edges_t structure should be not NULL.
2544
* param <-- user-defined parameters for the joining algorithm
2545
* n_iwm_vertices <-- initial number of vertices in work_mesh
2546
* o2n_vtx_id <-- array keeping the evolution of the vertex ids
2547
* edges <-- edges definition
2548
* p_inter_edges <-> pointer to the structure keeping data on
2549
* edge intersections
2550
*---------------------------------------------------------------------------*/
2553
_update_inter_edges_after_merge(cs_join_param_t param,
2554
cs_int_t n_iwm_vertices,
2555
const cs_int_t o2n_vtx_id[],
2556
const cs_join_edges_t *edges,
2557
const cs_join_mesh_t *mesh,
2558
cs_join_inter_edges_t **p_inter_edges)
2560
cs_int_t i, j,k, j1, j2, start_shift, idx_shift;
2561
cs_int_t save, _start, _end, start, end;
2562
cs_int_t v1_num, v2_num, v1_id, v2_id, sub_edge_id;
2563
fvm_gnum_t v1_gnum, v2_gnum, new_gnum, prev_gnum;
2566
cs_int_t n_adds = 0;
2568
cs_join_inter_edges_t *inter_edges = *p_inter_edges;
2569
cs_join_inter_edges_t *new_inter_edges = NULL;
2570
cs_int_t n_edges = inter_edges->n_edges;
2571
cs_int_t init_list_size = inter_edges->index[n_edges];
2572
FILE *logfile = cs_glob_join_log;
2574
assert(n_edges == edges->n_edges);
2576
/* Define vtx_glst to compare global vertex numbering */
2578
if (inter_edges->vtx_glst == NULL)
2579
BFT_MALLOC(inter_edges->vtx_glst, inter_edges->index[n_edges], fvm_gnum_t);
2581
for (i = 0; i < inter_edges->index[n_edges]; i++) {
2582
v1_id = inter_edges->vtx_lst[i] - 1;
2583
inter_edges->vtx_glst[i] = mesh->vertices[v1_id].gnum;
2586
/* Delete redundancies of vertices sharing the same global numbering
2587
after the merge step and define a new index */
2590
save = inter_edges->index[0];
2592
for (i = 0; i < n_edges; i++) {
2595
end = inter_edges->index[i+1];
2597
if (end - start > 0) {
2599
start_shift = start;
2600
v1_id = edges->def[2*i] - 1;
2601
v2_id = edges->def[2*i+1] - 1;
2602
v1_gnum = mesh->vertices[v1_id].gnum;
2603
v2_gnum = mesh->vertices[v2_id].gnum;
2604
prev_gnum = inter_edges->vtx_glst[start_shift];
2606
/* Don't take into account vertices with the same number as the
2607
first edge element */
2609
while (prev_gnum == v1_gnum && start_shift + 1 < end)
2610
prev_gnum = inter_edges->vtx_glst[++start_shift];
2612
if (prev_gnum != v1_gnum && start_shift < end) {
2614
inter_edges->vtx_lst[idx_shift] = inter_edges->vtx_lst[start_shift];
2615
inter_edges->abs_lst[idx_shift] = inter_edges->abs_lst[start_shift];
2616
inter_edges->vtx_glst[idx_shift] = inter_edges->vtx_glst[start_shift];
2619
for (j = start_shift + 1; j < end; j++) {
2621
new_gnum = inter_edges->vtx_glst[j];
2623
/* Don't take into account redundancies and vertices with the same
2624
number as the second edge element */
2626
if (prev_gnum != new_gnum && new_gnum != v2_gnum) {
2627
prev_gnum = new_gnum;
2628
inter_edges->vtx_lst[idx_shift] = inter_edges->vtx_lst[j];
2629
inter_edges->abs_lst[idx_shift] = inter_edges->abs_lst[j];
2630
inter_edges->vtx_glst[idx_shift] = inter_edges->vtx_glst[j];
2636
} /* If start_shift < end */
2638
} /* end - start > 0 */
2640
save = inter_edges->index[i+1];
2641
inter_edges->index[i+1] = idx_shift;
2643
} /* End of loop on edge intersections */
2645
inter_edges->max_sub_size = 0;
2647
for (i = 0; i < n_edges; i++)
2648
inter_edges->max_sub_size =
2649
CS_MAX(inter_edges->max_sub_size,
2650
inter_edges->index[i+1] - inter_edges->index[i]);
2652
assert(inter_edges->index[n_edges] <= init_list_size);
2654
BFT_REALLOC(inter_edges->vtx_lst, inter_edges->index[n_edges], cs_int_t);
2655
BFT_REALLOC(inter_edges->abs_lst, inter_edges->index[n_edges], float);
2657
#if 0 && defined(DEBUG) && !defined(NDEBUG) /* Dump local structures */
2658
fprintf(logfile, " AFTER REDUNDANCIES CLEAN\n");
2659
cs_join_inter_edges_dump(logfile, inter_edges, edges, mesh);
2662
/* Add new vertices from initial edges which are now sub-edges */
2664
for (i = 0; i < n_edges; i++)
2665
n_adds += _count_new_sub_edge_elts(i, inter_edges, edges, n_iwm_vertices);
2667
if (param.verbosity > 2)
2669
" Number of sub-elements to add to edge definition: %8d\n",
2672
if (n_adds > 0) { /* Define a new inter_edges structure */
2674
new_inter_edges = cs_join_inter_edges_create(n_edges);
2676
BFT_MALLOC(new_inter_edges->vtx_lst,
2677
inter_edges->index[n_edges] + n_adds, cs_int_t);
2678
BFT_MALLOC(new_inter_edges->abs_lst,
2679
inter_edges->index[n_edges] + n_adds, float);
2683
new_inter_edges->index[0] = 0;
2685
for (i = 0; i < n_edges; i++) {
2687
new_inter_edges->edge_gnum[i] = inter_edges->edge_gnum[i];
2688
start = inter_edges->index[i];
2689
end = inter_edges->index[i+1];
2691
if (start - end > 0) {
2693
for (j1 = start; j1 < end-1; j1++) {
2695
v1_num = inter_edges->vtx_lst[j1];
2696
new_inter_edges->vtx_lst[idx_shift] = v1_num;
2697
new_inter_edges->abs_lst[idx_shift] = inter_edges->abs_lst[j1];
2700
if (v1_num <= n_iwm_vertices) { /* v1 is an initial vertex */
2701
for (j2 = j1+1; j2 < end; j2++) {
2703
v2_num = inter_edges->vtx_lst[j2];
2705
if (v2_num <= n_iwm_vertices) { /* (v1,v2) is an initial edge */
2708
CS_ABS(cs_join_mesh_get_edge(v1_num, v2_num, edges))-1;
2709
assert(sub_edge_id != -1);
2711
_start = inter_edges->index[sub_edge_id];
2712
_end = inter_edges->index[sub_edge_id+1];
2714
for (j = _start; j < _end; j++) {
2717
for (k = j1+1; k < j2; k++)
2718
if (inter_edges->vtx_glst[k] == inter_edges->vtx_glst[j])
2721
if (found == false) {
2723
new_inter_edges->vtx_lst[idx_shift] =
2724
inter_edges->vtx_lst[j];
2725
new_inter_edges->abs_lst[idx_shift] =
2726
inter_edges->abs_lst[j];
2731
} /* End of loop on sub_edge_id definition */
2735
} /* End of loop on j2 */
2738
} /* End of loop on j1 */
2740
/* Add last vertex in the previous edge definition */
2742
new_inter_edges->vtx_lst[idx_shift] = inter_edges->vtx_lst[end-1];
2743
new_inter_edges->abs_lst[idx_shift] = inter_edges->abs_lst[end-1];
2746
} /* If end - start > 0 */
2748
new_inter_edges->index[i+1] = idx_shift;
2750
} /* End of loop on edges */
2752
cs_join_inter_edges_destroy(&inter_edges);
2753
inter_edges = new_inter_edges;
2755
inter_edges->max_sub_size = 0;
2757
for (i = 0; i < n_edges; i++)
2758
inter_edges->max_sub_size = CS_MAX(inter_edges->max_sub_size,
2759
inter_edges->index[i+1]);
2761
} /* End if n_adds > 0 */
2763
#if 0 && defined(DEBUG) && !defined(NDEBUG) /* Dump local structures */
2764
if (logfile != NULL) {
2765
fprintf(logfile, " AFTER SUB ELTS ADD\n");
2766
cs_join_inter_edges_dump(logfile, inter_edges, edges, mesh);
2770
/* Update cs_join_inter_edges_t structure */
2772
for (i = 0; i < n_edges; i++) {
2774
start = inter_edges->index[i];
2775
end = inter_edges->index[i+1];
2777
for (j = start; j < end; j++) {
2779
cs_int_t old_id = inter_edges->vtx_lst[j] - 1;
2781
inter_edges->vtx_lst[j] = o2n_vtx_id[old_id] + 1;
2787
/* Return pointer */
2789
*p_inter_edges = inter_edges;
2792
#if defined(HAVE_MPI)
2794
/*----------------------------------------------------------------------------
2795
* Define send_rank_index and send_faces to prepare the exchange of new faces
2796
* between mesh structures.
2799
* n_faces <-- number of faces to send
2800
* n_g_faces <-- global number of faces to be joined
2801
* face_gnum <-- global face number
2802
* gnum_rank_index <-- index on ranks for the init. global face numbering
2803
* p_send_rank_index --> index on ranks for sending face
2804
* p_send_faces --> list of face ids to send
2805
*---------------------------------------------------------------------------*/
2808
_get_faces_to_send(cs_int_t n_faces,
2809
fvm_gnum_t n_g_faces,
2810
const fvm_gnum_t face_gnum[],
2811
const fvm_gnum_t gnum_rank_index[],
2812
cs_int_t *p_send_rank_index[],
2813
cs_int_t *p_send_faces[])
2815
cs_int_t i, rank, shift, reduce_rank;
2816
fvm_gnum_t start_gnum, end_gnum;
2817
cs_join_block_info_t block_info;
2819
cs_int_t reduce_size = 0;
2820
cs_int_t *send_rank_index = NULL, *send_faces = NULL;
2821
cs_int_t *reduce_ids = NULL, *count = NULL;
2822
fvm_gnum_t *reduce_index = NULL;
2824
const int local_rank = CS_MAX(cs_glob_rank_id, 0);
2825
const int n_ranks = cs_glob_n_ranks;
2829
assert(gnum_rank_index != NULL);
2830
assert(n_ranks > 1);
2832
/* Compute block_size */
2834
block_info = cs_join_get_block_info(n_g_faces, n_ranks, local_rank);
2835
start_gnum = block_info.first_gnum;
2836
end_gnum = block_info.first_gnum + block_info.local_size;
2838
/* Compact init. global face distribution. Remove ranks without face
2841
for (i = 0; i < n_ranks; i++)
2842
if (gnum_rank_index[i] < gnum_rank_index[i+1])
2845
BFT_MALLOC(reduce_index, reduce_size+1, fvm_gnum_t);
2846
BFT_MALLOC(reduce_ids, reduce_size, cs_int_t);
2849
reduce_index[0] = gnum_rank_index[0] + 1;
2851
for (i = 0; i < n_ranks; i++) {
2853
/* Add +1 to gnum_rank_index because it's an id and we work on numbers */
2855
if (gnum_rank_index[i] < gnum_rank_index[i+1]) {
2856
reduce_index[reduce_size+1] = gnum_rank_index[i+1] + 1;
2857
reduce_ids[reduce_size++] = i;
2862
BFT_MALLOC(send_rank_index, n_ranks + 1, cs_int_t);
2864
for (i = 0; i < n_ranks + 1; i++)
2865
send_rank_index[i] = 0;
2867
/* Count number of ranks associated to each new face */
2869
for (i = 0; i < n_faces; i++) {
2871
if (face_gnum[i] >= start_gnum && face_gnum[i] < end_gnum) {
2873
/* The current face is a "main" face for the local rank */
2875
reduce_rank = cs_search_gindex_binary(reduce_size,
2879
assert(reduce_rank > -1);
2880
assert(reduce_rank < reduce_size);
2882
rank = reduce_ids[reduce_rank];
2883
send_rank_index[rank+1] += 1;
2889
for (i = 0; i < n_ranks; i++)
2890
send_rank_index[i+1] += send_rank_index[i];
2892
BFT_MALLOC(send_faces, send_rank_index[n_ranks], cs_int_t);
2893
BFT_MALLOC(count, n_ranks, cs_int_t);
2895
for (i = 0; i < n_ranks; i++)
2898
/* Fill the list of ranks */
2900
for (i = 0; i < n_faces; i++) {
2902
if (face_gnum[i] >= start_gnum && face_gnum[i] < end_gnum) {
2904
/* The current face is a "main" face for the local rank */
2906
reduce_rank = cs_search_gindex_binary(reduce_size,
2910
assert(reduce_rank > -1);
2911
assert(reduce_rank < reduce_size);
2913
rank = reduce_ids[reduce_rank];
2914
shift = send_rank_index[rank] + count[rank];
2915
send_faces[shift] = i;
2918
} /* End of loop on initial faces */
2924
#if 0 && defined(DEBUG) && !defined(NDEBUG)
2925
if (cs_glob_join_log != NULL) {
2926
FILE *logfile = cs_glob_join_log;
2927
for (rank = 0; rank < n_ranks; rank++) {
2928
fprintf(logfile, " rank %d: ", rank);
2929
for (i = send_rank_index[rank]; i < send_rank_index[rank+1]; i++)
2930
fprintf(logfile, " %d (%llu)",
2931
send_faces[i], (unsigned long long)face_gnum[send_faces[i]]);
2932
fprintf(logfile, "\n");
2939
BFT_FREE(reduce_ids);
2940
BFT_FREE(reduce_index);
2942
/* Set return pointers */
2944
*p_send_rank_index = send_rank_index;
2945
*p_send_faces = send_faces;
2948
#endif /* defined(HAVE_MPI) */
2950
/*----------------------------------------------------------------------------
2951
* Update local_mesh by redistributing mesh.
2952
* Send back to the original rank the new face and vertex description.
2955
* gnum_rank_index <-- index on ranks for the old global face numbering
2956
* send_mesh <-- distributed mesh on faces to join
2957
* p_recv_mesh <-> mesh on local selected faces to be joined
2958
*---------------------------------------------------------------------------*/
2961
_redistribute_mesh(const fvm_gnum_t gnum_rank_index[],
2962
const cs_join_mesh_t *send_mesh,
2963
cs_join_mesh_t **p_recv_mesh)
2965
cs_join_mesh_t *recv_mesh = *p_recv_mesh;
2967
const int n_ranks = cs_glob_n_ranks;
2971
assert(send_mesh != NULL);
2972
assert(recv_mesh != NULL);
2975
cs_join_mesh_copy(&recv_mesh, send_mesh);
2977
#if defined(HAVE_MPI)
2978
if (n_ranks > 1) { /* Parallel mode */
2980
cs_int_t *send_rank_index = NULL, *send_faces = NULL;
2982
MPI_Comm mpi_comm = cs_glob_mpi_comm;
2984
/* Free some structures of the mesh */
2986
cs_join_mesh_reset(recv_mesh);
2988
_get_faces_to_send(send_mesh->n_faces,
2989
send_mesh->n_g_faces,
2990
send_mesh->face_gnum,
2995
assert(send_rank_index[n_ranks] <= send_mesh->n_faces);
2997
/* Get the new face connectivity from the distributed send_mesh */
2999
cs_join_mesh_exchange(n_ranks,
3006
BFT_FREE(send_faces);
3007
BFT_FREE(send_rank_index);
3012
/* Return pointers */
3014
*p_recv_mesh = recv_mesh;
3018
/*============================================================================
3019
* Public function definitions
3020
*===========================================================================*/
3022
/*----------------------------------------------------------------------------
3023
* Creation of new vertices.
3025
* Update list of equivalent vertices, and assign a vertex (existing or
3026
* newly created) to each intersection.
3029
* verbosity <-- verbosity level
3030
* edges <-- list of edges
3031
* work <-> joining mesh maintaining initial vertex data
3032
* inter_set <-> cs_join_inter_set_t structure including
3033
* data on edge-edge intersections
3034
* init_max_vtx_gnum <-- initial max. global numbering for vertices
3035
* p_n_g_new_vertices <-> pointer to the global number of new vertices
3036
* p_vtx_eset <-> pointer to a structure dealing with vertex
3038
*---------------------------------------------------------------------------*/
3041
cs_join_create_new_vertices(int verbosity,
3042
const cs_join_edges_t *edges,
3043
cs_join_mesh_t *work,
3044
cs_join_inter_set_t *inter_set,
3045
fvm_gnum_t init_max_vtx_gnum,
3046
fvm_gnum_t *p_n_g_new_vertices,
3047
cs_join_eset_t **p_vtx_eset)
3051
cs_join_vertex_t v1, v2;
3053
cs_int_t n_new_vertices = 0;
3054
fvm_gnum_t n_g_new_vertices = 0;
3055
fvm_gnum_t *new_vtx_gnum = NULL;
3056
cs_int_t n_iwm_vertices = work->n_vertices;
3057
cs_join_eset_t *vtx_equiv = *p_vtx_eset;
3059
/* Count the number of new vertices. Update cs_join_inter_set_t struct. */
3061
for (i = 0; i < inter_set->n_inter; i++) {
3063
cs_join_inter_t inter1 = inter_set->inter_lst[2*i];
3064
cs_join_inter_t inter2 = inter_set->inter_lst[2*i+1];
3066
inter1.vtx_id = _get_vtx_id(inter1,
3067
&(edges->def[2*inter1.edge_id]),
3071
inter2.vtx_id = _get_vtx_id(inter2,
3072
&(edges->def[2*inter2.edge_id]),
3076
inter_set->inter_lst[2*i] = inter1;
3077
inter_set->inter_lst[2*i+1] = inter2;
3079
} /* End of loop on intersections */
3081
/* Compute the global numbering for the new vertices (Take into account
3082
potential redundancies) */
3084
_compute_new_vertex_gnum(work,
3094
bft_printf(_("\n Global number of new vertices to create: %10llu\n"),
3095
(unsigned long long)n_g_new_vertices);
3097
/* Define new vertices */
3099
work->n_vertices += n_new_vertices;
3100
work->n_g_vertices += n_g_new_vertices;
3102
BFT_REALLOC(work->vertices, work->n_vertices, cs_join_vertex_t);
3104
#if defined(DEBUG) && !defined(NDEBUG) /* Prepare sanity checks */
3106
cs_join_vertex_t incoherency;
3108
/* Initialize to incoherent values new vertices structures */
3110
incoherency.gnum = 0;
3111
incoherency.coord[0] = -9999.9999;
3112
incoherency.coord[1] = -9999.9999;
3113
incoherency.coord[2] = -9999.9999;
3114
incoherency.tolerance = -1.0;
3115
incoherency.state = CS_JOIN_STATE_UNDEF;
3117
for (i = 0; i < n_new_vertices; i++)
3118
work->vertices[n_iwm_vertices + i] = incoherency;
3123
/* Fill vertices structure with new vertex definitions */
3125
for (i = 0; i < inter_set->n_inter; i++) {
3127
cs_join_inter_t inter1 = inter_set->inter_lst[2*i];
3128
cs_join_inter_t inter2 = inter_set->inter_lst[2*i+1];
3129
cs_int_t v1_num = inter1.vtx_id + 1;
3130
cs_int_t v2_num = inter2.vtx_id + 1;
3131
cs_int_t equiv_id = vtx_equiv->n_equiv;
3133
assert(inter1.vtx_id < work->n_vertices);
3134
assert(inter2.vtx_id < work->n_vertices);
3136
/* Create new vertices if needed */
3138
if (v1_num > n_iwm_vertices) {
3140
shift = inter1.vtx_id - n_iwm_vertices;
3141
v1 = _get_new_vertex(inter1.curv_abs,
3142
new_vtx_gnum[shift],
3143
&(edges->def[2*inter1.edge_id]),
3145
tol_min = v1.tolerance;
3149
tol_min = work->vertices[v1_num-1].tolerance;
3151
if (v2_num > n_iwm_vertices) {
3153
shift = inter2.vtx_id - n_iwm_vertices;
3154
v2 = _get_new_vertex(inter2.curv_abs,
3155
new_vtx_gnum[shift],
3156
&(edges->def[2*inter2.edge_id]),
3158
tol_min = CS_MIN(tol_min, v2.tolerance);
3162
tol_min = CS_MIN(tol_min, work->vertices[v2_num-1].tolerance);
3164
/* A new vertex has a tolerance equal to the minimal tolerance
3165
between the two vertices implied in the intersection */
3167
if (v1_num > n_iwm_vertices) {
3168
v1.tolerance = tol_min;
3169
work->vertices[inter1.vtx_id] = v1;
3171
if (v2_num > n_iwm_vertices) {
3172
v2.tolerance = tol_min;
3173
work->vertices[inter2.vtx_id] = v2;
3176
/* Add equivalence between the two current vertices */
3178
cs_join_eset_check_size(equiv_id, &vtx_equiv);
3180
if (v1_num < v2_num) {
3181
vtx_equiv->equiv_couple[2*equiv_id] = v1_num;
3182
vtx_equiv->equiv_couple[2*equiv_id+1] = v2_num;
3185
vtx_equiv->equiv_couple[2*equiv_id] = v2_num;
3186
vtx_equiv->equiv_couple[2*equiv_id+1] = v1_num;
3189
vtx_equiv->n_equiv += 1;
3191
} /* End of loop on intersections */
3195
BFT_FREE(new_vtx_gnum);
3197
#if defined(DEBUG) && !defined(NDEBUG) /* Sanity checks */
3198
for (i = 0; i < work->n_vertices; i++) {
3200
cs_join_vertex_t vtx = work->vertices[i];
3202
if (vtx.gnum == 0 || vtx.tolerance < -0.99)
3203
bft_error(__FILE__, __LINE__, 0,
3204
_(" Inconsistent value found in cs_join_vertex_t struct.:\n"
3205
" Vertex %d is defined by:\n"
3206
" %u - [%7.4le, %7.4le, %7.4le] - %lg\n"),
3207
i, vtx.gnum, vtx.coord[0], vtx.coord[1], vtx.coord[2],
3210
} /* End of loop on vertices */
3213
_dump_vtx_eset(vtx_equiv, work, cs_glob_join_log);
3218
/* Set return pointers */
3220
*p_n_g_new_vertices = n_g_new_vertices;
3221
*p_vtx_eset = vtx_equiv;
3224
/*----------------------------------------------------------------------------
3225
* Merge of equivalent vertices (and tolerance reduction if necessary)
3227
* Define a new cs_join_vertex_t structure (stored in "work" structure).
3228
* Returns an updated cs_join_mesh_t and cs_join_edges_t structures.
3231
* param <-- set of user-defined parameters for the joining
3232
* n_g_vertices_tot <-- global number of vertices (initial parent mesh)
3233
* work <-> pointer to a cs_join_mesh_t structure
3234
* vtx_eset <-- structure storing equivalences between vertices
3235
* (two vertices are equivalent if they are within
3236
* each other's tolerance)
3237
*---------------------------------------------------------------------------*/
3240
cs_join_merge_vertices(cs_join_param_t param,
3241
fvm_gnum_t n_g_vertices_tot,
3242
cs_join_mesh_t *work,
3243
const cs_join_eset_t *vtx_eset)
3245
double clock_start, clock_end, cpu_start, cpu_end;
3247
fvm_gnum_t *vtx_tags = NULL;
3248
cs_join_gset_t *merge_set = NULL;
3250
const int n_ranks = cs_glob_n_ranks;
3252
/* Initialize counters for the merge operation */
3254
_initialize_merge_counter();
3256
#if 0 && defined(DEBUG) && !defined(NDEBUG) /* Dump local structures */
3257
_dump_vtx_eset(vtx_eset, work, cs_glob_join_log);
3260
if (param.verbosity > 2) {
3261
fvm_gnum_t g_n_equiv = vtx_eset->n_equiv;
3262
fvm_parall_counter(&g_n_equiv, 1);
3263
fprintf(cs_glob_join_log,
3265
" Final number of equiv. between vertices; local: %9d\n"
3267
vtx_eset->n_equiv, (unsigned long long)g_n_equiv);
3270
/* Operate merge between equivalent vertices.
3271
Manage reduction of tolerance if necessary */
3273
clock_start = bft_timer_wtime();
3274
cpu_start = bft_timer_cpu_time();
3276
/* Tag with the same number all the vertices which might be merged together */
3278
_tag_equiv_vertices(n_g_vertices_tot,
3284
if (n_ranks == 1) { /* Serial mode */
3286
/* Build a merge list */
3288
merge_set = cs_join_gset_create_from_tag(work->n_vertices, vtx_tags);
3290
/* Merge of equivalent vertices */
3292
_merge_vertices(param,
3299
#if defined(HAVE_MPI)
3300
if (n_ranks > 1) { /* Parallel mode: we work by block */
3302
cs_int_t *send_count = NULL, *recv_count = NULL;
3303
cs_int_t *send_shift = NULL, *recv_shift = NULL;
3304
cs_join_vertex_t *vtx_merge_data = NULL;
3306
BFT_MALLOC(send_count, n_ranks, cs_int_t);
3307
BFT_MALLOC(recv_count, n_ranks, cs_int_t);
3308
BFT_MALLOC(send_shift, n_ranks+1, cs_int_t);
3309
BFT_MALLOC(recv_shift, n_ranks+1, cs_int_t);
3311
/* Build a merge list in parallel */
3313
_build_parall_merge_structures(work,
3315
send_count, send_shift,
3316
recv_count, recv_shift,
3320
/* Merge of equivalent vertices for the current block */
3322
_merge_vertices(param,
3324
recv_shift[n_ranks],
3327
_exchange_merged_vertices(work,
3329
send_count, send_shift,
3330
recv_count, recv_shift,
3333
BFT_FREE(send_count);
3334
BFT_FREE(send_shift);
3335
BFT_FREE(recv_count);
3336
BFT_FREE(recv_shift);
3337
BFT_FREE(vtx_merge_data);
3340
#endif /* HAVE_MPI */
3346
clock_end = bft_timer_wtime();
3347
cpu_end = bft_timer_cpu_time();
3349
cs_join_gset_destroy(&merge_set);
3351
if (param.verbosity > 1)
3353
" Vertex merge (only)\n"
3354
" wall clock time: %10.3g\n"
3355
" cpu time: %10.3g\n"),
3356
clock_end - clock_start, cpu_end - cpu_start);
3359
/*----------------------------------------------------------------------------
3360
* Merge of equivalent vertices (and reduction of tolerance if necessary)
3362
* Define a new cs_join_vertex_t structure (stored in "work" structure)
3363
* Returns an updated cs_join_mesh_t and cs_join_edges_t structures.
3366
* param <-- set of user-defined parameters for the joining
3367
* n_iwm_vertices <-- initial number of vertices (work mesh struct.)
3368
* iwm_vtx_gnum <-- initial global vertex num. (work mesh struct)
3369
* init_max_vtx_gnum <-- initial max. global numbering for vertices
3370
* rank_face_gnum_index <-- index on face global numbering to determine
3372
* p_mesh <-> pointer to cs_join_mesh_t structure
3373
* p_edges <-> pointer to cs_join_edges_t structure
3374
* p_inter_edges <-> pointer to a cs_join_inter_edges_t struct.
3375
* p_local_mesh <-> pointer to a cs_join_mesh_t structure
3376
* p_o2n_vtx_gnum --> array on blocks on the new global vertex
3377
* numbering for the init. vertices (before inter.)
3378
*---------------------------------------------------------------------------*/
3381
cs_join_merge_update_struct(cs_join_param_t param,
3382
cs_int_t n_iwm_vertices,
3383
const fvm_gnum_t iwm_vtx_gnum[],
3384
fvm_gnum_t init_max_vtx_gnum,
3385
const fvm_gnum_t rank_face_gnum_index[],
3386
cs_join_mesh_t **p_mesh,
3387
cs_join_edges_t **p_edges,
3388
cs_join_inter_edges_t **p_inter_edges,
3389
cs_join_mesh_t **p_local_mesh,
3390
fvm_gnum_t *p_o2n_vtx_gnum[])
3392
cs_int_t n_am_vertices = 0; /* new number of vertices after merge */
3393
cs_int_t *o2n_vtx_id = NULL;
3394
fvm_gnum_t *o2n_vtx_gnum = NULL;
3395
cs_join_mesh_t *mesh = *p_mesh;
3396
cs_join_mesh_t *local_mesh = *p_local_mesh;
3397
cs_join_edges_t *edges = *p_edges;
3398
cs_join_inter_edges_t *inter_edges = *p_inter_edges;
3400
/* Keep an history of the evolution of each vertex */
3402
_keep_global_vtx_evolution(n_iwm_vertices, /* n_vertices before inter */
3405
mesh->n_vertices, /* n_vertices after inter */
3407
&o2n_vtx_gnum); /* defined by block in // */
3409
_keep_local_vtx_evolution(mesh->n_vertices, /* n_vertices after inter */
3411
&n_am_vertices, /* n_vertices after merge */
3414
/* Update all structures which keeps data about vertices */
3416
if (inter_edges != NULL) { /* The join type is not conform */
3418
/* Update inter_edges structure */
3420
_update_inter_edges_after_merge(param,
3422
o2n_vtx_id, /* size of mesh->n_vertices */
3427
assert(edges->n_edges == inter_edges->n_edges); /* Else: problem for
3430
/* Update cs_join_mesh_t structure after the merge of vertices
3431
numbering of the old vertices + add new vertices */
3433
cs_join_mesh_update(mesh,
3436
inter_edges->vtx_lst,
3440
} /* End if inter_edges != NULL */
3443
/* Update cs_join_mesh_t structure after the merge of vertices
3444
numbering of the old vertices + add new vertices */
3446
cs_join_mesh_update(mesh,
3453
BFT_FREE(o2n_vtx_id);
3455
/* Update local_mesh by redistributing mesh */
3457
_redistribute_mesh(rank_face_gnum_index,
3461
/* Set return pointers */
3465
*p_inter_edges = inter_edges;
3466
*p_o2n_vtx_gnum = o2n_vtx_gnum;
3467
*p_local_mesh = local_mesh;
3470
/*---------------------------------------------------------------------------*/
added
removed
src/mesh/cs_join_merge.c
