1 /*============================================================================
2 * Set of subroutines for:
3 * - merging equivalent vertices,
4 * - managing tolerance reduction
5 *===========================================================================*/
6
7 /*
8 This file is part of Code_Saturne, a general-purpose CFD tool.
9
10 Copyright (C) 1998-2021 EDF S.A.
11
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
15 version.
16
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
20 details.
21
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.
25 */
26
27 /*----------------------------------------------------------------------------*/
28
29 #include "cs_defs.h"
30
31 /*----------------------------------------------------------------------------
32 * Standard C library headers
33 *---------------------------------------------------------------------------*/
34
35 #include <stdlib.h>
36 #include <stdio.h>
37 #include <string.h>
38 #include <math.h>
39 #include <float.h>
40 #include <assert.h>
41
42 /*----------------------------------------------------------------------------
43 * Local headers
44 *---------------------------------------------------------------------------*/
45
46 #include "bft_mem.h"
47 #include "bft_printf.h"
48
49 #include "fvm_io_num.h"
50
51 #include "cs_all_to_all.h"
52 #include "cs_block_dist.h"
53 #include "cs_log.h"
54 #include "cs_order.h"
55 #include "cs_search.h"
56 #include "cs_join_post.h"
57 #include "cs_parall.h"
58
59 /*----------------------------------------------------------------------------
60 * Header for the current file
61 *---------------------------------------------------------------------------*/
62
63 #include "cs_join_merge.h"
64
65 /*---------------------------------------------------------------------------*/
66
67 BEGIN_C_DECLS
68
69 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
70
71 /*============================================================================
72 * Local macro definitions
73 *===========================================================================*/
74
75 /* Turn on (1) or off (0) the tolerance reduc. */
76 #define CS_JOIN_MERGE_TOL_REDUC 1
77 #define CS_JOIN_MERGE_INV_TOL 1
78
79 /*============================================================================
80 * Local structure and type definitions
81 *===========================================================================*/
82
83 /*============================================================================
84 * Global variable definitions
85 *===========================================================================*/
86
87 /* Parameters to control the vertex merge */
88
89 enum {
90
91 CS_JOIN_MERGE_MAX_GLOB_ITERS = 50, /* Max. number of glob. iter. for finding
92 equivalent vertices */
93 CS_JOIN_MERGE_MAX_LOC_ITERS = 100 /* Max. number of loc. iter. for finding
94 equivalent vertices */
95 };
96
97 /* Coefficient to deal with rounding approximations */
98
99 static const double cs_join_tol_eps_coef2 = 1.0001*1.001;
100
101 /* Counter on the number of loops useful to converge for the merge operation */
102
103 static int _glob_merge_counter = 0, _loc_merge_counter = 0;
104
105 /*============================================================================
106 * Private function definitions
107 *===========================================================================*/
108
109 /*----------------------------------------------------------------------------
110 * Initialize counter for the merge operation
111 *---------------------------------------------------------------------------*/
112
113 static void
_initialize_merge_counter(void)114 _initialize_merge_counter(void)
115 {
116 _glob_merge_counter = 0;
117 _loc_merge_counter = 0;
118 }
119
120 #if 0 && defined(DEBUG) && !defined(NDEBUG)
121
122 /*----------------------------------------------------------------------------
123 * Dump an cs_join_eset_t structure on vertices.
124 *
125 * parameters:
126 * e_set <-- cs_join_eset_t structure to dump
127 * mesh <-- cs_join_mesh_t structure associated
128 * logfile <-- handle to log file
129 *---------------------------------------------------------------------------*/
130
131 static void
132 _dump_vtx_eset(const cs_join_eset_t *e_set,
133 const cs_join_mesh_t *mesh,
134 FILE *logfile)
135 {
136 int i;
137
138 fprintf(logfile, "\n Dump an cs_join_eset_t structure (%p)\n",
139 (const void *)e_set);
140 fprintf(logfile, " n_max_equiv: %10d\n", e_set->n_max_equiv);
141 fprintf(logfile, " n_equiv : %10d\n\n", e_set->n_equiv);
142
143 for (i = 0; i < e_set->n_equiv; i++) {
144
145 cs_lnum_t v1_num = e_set->equiv_couple[2*i];
146 cs_lnum_t v2_num = e_set->equiv_couple[2*i+1];
147
148 fprintf(logfile,
149 " %10d - local: (%9d, %9d) - global: (%10llu, %10llu)\n",
150 i, v1_num, v2_num,
151 (unsigned long long)(mesh->vertices[v1_num-1]).gnum,
152 (unsigned long long)(mesh->vertices[v2_num-1]).gnum);
153
154 }
155 fflush(logfile);
156 }
157
158 #endif /* Only in debug mode */
159
160 /*----------------------------------------------------------------------------
161 * Compute the length of a segment between two vertices.
162 *
163 * parameters:
164 * v1 <-- cs_join_vertex_t structure for the first vertex of the segment
165 * v2 <-- cs_join_vertex_t structure for the second vertex of the segment
166 *
167 * returns:
168 * length of the segment
169 *---------------------------------------------------------------------------*/
170
171 inline static cs_real_t
_compute_length(cs_join_vertex_t v1,cs_join_vertex_t v2)172 _compute_length(cs_join_vertex_t v1,
173 cs_join_vertex_t v2)
174 {
175 cs_lnum_t k;
176 cs_real_t len = 0.0, d2 = 0.0;
177
178 for (k = 0; k < 3; k++) {
179 cs_real_t d = v1.coord[k] - v2.coord[k];
180 d2 += d * d;
181 }
182 len = sqrt(d2);
183
184 return len;
185 }
186
187 /*----------------------------------------------------------------------------
188 * Compute a new cs_join_vertex_t structure.
189 *
190 * parameters:
191 * curv_abs <-- curvilinear abscissa of the intersection
192 * gnum <-- global number associated to the new
193 * cs_join_vertex_t structure
194 * vtx_couple <-- couple of vertex numbers defining the current edge
195 * work <-- local cs_join_mesh_t structure
196 *
197 * returns:
198 * a new cs_join_vertex_t structure
199 *---------------------------------------------------------------------------*/
200
201 static cs_join_vertex_t
_get_new_vertex(cs_coord_t curv_abs,cs_gnum_t gnum,const cs_lnum_t vtx_couple[],const cs_join_mesh_t * work)202 _get_new_vertex(cs_coord_t curv_abs,
203 cs_gnum_t gnum,
204 const cs_lnum_t vtx_couple[],
205 const cs_join_mesh_t *work)
206 {
207 cs_lnum_t k;
208 cs_join_vertex_t new_vtx_data;
209
210 #if defined(DEBUG) && !defined(NDEBUG)
211 /* Avoid Valgrind warnings in byte copies due to padding */
212 memset(&new_vtx_data, 0, sizeof(cs_join_vertex_t));
213 #endif
214
215 cs_join_vertex_t v1 = work->vertices[vtx_couple[0]-1];
216 cs_join_vertex_t v2 = work->vertices[vtx_couple[1]-1];
217
218 assert(curv_abs >= 0.0);
219 assert(curv_abs <= 1.0);
220
221 /* New vertex features */
222
223 new_vtx_data.state = CS_JOIN_STATE_NEW;
224 new_vtx_data.gnum = gnum;
225 new_vtx_data.tolerance = (1-curv_abs)*v1.tolerance + curv_abs*v2.tolerance;
226
227 for (k = 0; k < 3; k++)
228 new_vtx_data.coord[k] = (1-curv_abs)*v1.coord[k] + curv_abs*v2.coord[k];
229
230 return new_vtx_data;
231 }
232
233 /*----------------------------------------------------------------------------
234 * Define a tag (3 values) to globally order intersections.
235 *
236 * parameters:
237 * tag <-> tag to fill
238 * e1_gnum <-- global number for the first edge
239 * e2_gnum <-- global number for the second edge
240 * link_vtx_gnum <-- global number of the vertex associated to the current
241 * intersection
242 *---------------------------------------------------------------------------*/
243
244 static void
_define_inter_tag(cs_gnum_t tag[],cs_gnum_t e1_gnum,cs_gnum_t e2_gnum,cs_gnum_t link_vtx_gnum)245 _define_inter_tag(cs_gnum_t tag[],
246 cs_gnum_t e1_gnum,
247 cs_gnum_t e2_gnum,
248 cs_gnum_t link_vtx_gnum)
249 {
250 if (e1_gnum < e2_gnum) {
251 tag[0] = e1_gnum;
252 tag[1] = e2_gnum;
253 }
254 else {
255 tag[0] = e2_gnum;
256 tag[1] = e1_gnum;
257 }
258
259 tag[2] = link_vtx_gnum;
260 }
261
262 /*----------------------------------------------------------------------------
263 * Creation of new vertices.
264 *
265 * Update list of equivalent vertices.
266 *
267 * parameters:
268 * work <-- pointer to a cs_join_mesh_t structure
269 * edges <-- list of edges
270 * inter_set <-- structure including data on edge intersections
271 * init_max_vtx_gnum <-- initial max. global numbering for vertices
272 * n_iwm_vertices <-- initial local number of vertices (work struct)
273 * n_new_vertices <-- local number of new vertices to define
274 * p_n_g_new_vertices <-> pointer to the global number of new vertices
275 * p_new_vtx_gnum <-> pointer to the global numbering array for the
276 * new vertices
277 *---------------------------------------------------------------------------*/
278
279 static void
_compute_new_vertex_gnum(const cs_join_mesh_t * work,const cs_join_edges_t * edges,const cs_join_inter_set_t * inter_set,cs_gnum_t init_max_vtx_gnum,cs_lnum_t n_iwm_vertices,cs_lnum_t n_new_vertices,cs_gnum_t * p_n_g_new_vertices,cs_gnum_t * p_new_vtx_gnum[])280 _compute_new_vertex_gnum(const cs_join_mesh_t *work,
281 const cs_join_edges_t *edges,
282 const cs_join_inter_set_t *inter_set,
283 cs_gnum_t init_max_vtx_gnum,
284 cs_lnum_t n_iwm_vertices,
285 cs_lnum_t n_new_vertices,
286 cs_gnum_t *p_n_g_new_vertices,
287 cs_gnum_t *p_new_vtx_gnum[])
288 {
289 cs_lnum_t i;
290
291 cs_gnum_t n_g_new_vertices = 0;
292 cs_lnum_t n_new_vertices_save = n_new_vertices;
293 cs_lnum_t *order = NULL;
294 cs_gnum_t *inter_tag = NULL, *adjacency = NULL, *new_vtx_gnum = NULL;
295 fvm_io_num_t *new_vtx_io_num = NULL;
296
297 /* Define a fvm_io_num_t structure to get the global numbering
298 for the new vertices.
299 First, build a tag associated to each intersection */
300
301 BFT_MALLOC(new_vtx_gnum, n_new_vertices, cs_gnum_t);
302 BFT_MALLOC(inter_tag, 3*n_new_vertices, cs_gnum_t);
303
304 n_new_vertices = 0;
305
306 for (i = 0; i < inter_set->n_inter; i++) {
307
308 cs_join_inter_t inter1 = inter_set->inter_lst[2*i];
309 cs_join_inter_t inter2 = inter_set->inter_lst[2*i+1];
310 cs_gnum_t e1_gnum = edges->gnum[inter1.edge_id];
311 cs_gnum_t e2_gnum = edges->gnum[inter2.edge_id];
312
313 if (inter1.vtx_id + 1 > n_iwm_vertices) {
314
315 if (inter2.vtx_id + 1 > n_iwm_vertices)
316 _define_inter_tag(&(inter_tag[3*n_new_vertices]),
317 e1_gnum, e2_gnum,
318 0);
319 else
320 _define_inter_tag(&(inter_tag[3*n_new_vertices]),
321 e1_gnum, e2_gnum,
322 (work->vertices[inter2.vtx_id]).gnum);
323
324 n_new_vertices++;
325
326 } /* New vertices for this intersection */
327
328 if (inter2.vtx_id + 1 > n_iwm_vertices) {
329
330 if (inter1.vtx_id + 1 > n_iwm_vertices)
331 _define_inter_tag(&(inter_tag[3*n_new_vertices]),
332 e1_gnum, e2_gnum,
333 init_max_vtx_gnum + 1);
334 else
335 _define_inter_tag(&(inter_tag[3*n_new_vertices]),
336 e1_gnum, e2_gnum,
337 (work->vertices[inter1.vtx_id]).gnum);
338
339 n_new_vertices++;
340
341 } /* New vertices for this intersection */
342
343 } /* End of loop on intersections */
344
345 if (n_new_vertices != n_new_vertices_save)
346 bft_error(__FILE__, __LINE__, 0,
347 _(" The number of new vertices to create is not consistent.\n"
348 " Previous number: %10ld\n"
349 " Current number: %10ld\n\n"),
350 (long)n_new_vertices_save, (long)n_new_vertices);
351
352 /* Create a new fvm_io_num_t structure */
353
354 BFT_MALLOC(order, n_new_vertices, cs_lnum_t);
355
356 cs_order_gnum_allocated_s(NULL, inter_tag, 3, order, n_new_vertices);
357
358 BFT_MALLOC(adjacency, 3*n_new_vertices, cs_gnum_t);
359
360 for (i = 0; i < n_new_vertices; i++) {
361
362 cs_lnum_t o_id = order[i];
363
364 adjacency[3*i] = inter_tag[3*o_id];
365 adjacency[3*i+1] = inter_tag[3*o_id+1];
366 adjacency[3*i+2] = inter_tag[3*o_id+2];
367
368 }
369
370 BFT_FREE(inter_tag);
371
372 if (cs_glob_n_ranks > 1) {
373
374 const cs_gnum_t *global_num = NULL;
375
376 new_vtx_io_num =
377 fvm_io_num_create_from_adj_s(NULL, adjacency, n_new_vertices, 3);
378
379 n_g_new_vertices = fvm_io_num_get_global_count(new_vtx_io_num);
380 global_num = fvm_io_num_get_global_num(new_vtx_io_num);
381
382 for (i = 0; i < n_new_vertices; i++)
383 new_vtx_gnum[order[i]] = global_num[i] + init_max_vtx_gnum;
384
385 fvm_io_num_destroy(new_vtx_io_num);
386
387 } /* End of parallel treatment */
388
389 else {
390
391 if (n_new_vertices > 0) {
392
393 cs_gnum_t new_gnum = init_max_vtx_gnum + 1;
394
395 new_vtx_gnum[order[0]] = new_gnum;
396
397 for (i = 1; i < n_new_vertices; i++) {
398
399 if (adjacency[3*i] != adjacency[3*(i-1)])
400 new_gnum += 1;
401 else {
402 if (adjacency[3*i+1] != adjacency[3*(i-1)+1])
403 new_gnum += 1;
404 else
405 if (adjacency[3*i+2] != adjacency[3*(i-1)+2])
406 new_gnum += 1;
407 }
408
409 new_vtx_gnum[order[i]] = new_gnum;
410
411 }
412
413 } /* End if n_new_vertices > 0 */
414
415 n_g_new_vertices = n_new_vertices;
416
417 } /* End of serial treatment */
418
419 /* Free memory */
420
421 BFT_FREE(order);
422 BFT_FREE(adjacency);
423
424 /* Return pointer */
425
426 *p_n_g_new_vertices = n_g_new_vertices;
427 *p_new_vtx_gnum = new_vtx_gnum;
428
429 }
430
431 /*----------------------------------------------------------------------------
432 * Get vertex id associated to the current intersection.
433 *
434 * Create a new vertex id if needed. Update n_new_vertices in this case.
435 *
436 * parameters:
437 * inter <-- a inter_t structure
438 * vtx_couple <-- couple of vertex numbers defining the current edge
439 * n_init_vertices <-- initial number of vertices
440 * n_new_vertices <-- number of new vertices created
441 *
442 * returns:
443 * vertex id associated to the current intersection.
444 *---------------------------------------------------------------------------*/
445
446 static cs_lnum_t
_get_vtx_id(cs_join_inter_t inter,const cs_lnum_t vtx_couple[],cs_lnum_t n_init_vertices,cs_lnum_t * p_n_new_vertices)447 _get_vtx_id(cs_join_inter_t inter,
448 const cs_lnum_t vtx_couple[],
449 cs_lnum_t n_init_vertices,
450 cs_lnum_t *p_n_new_vertices)
451 {
452 cs_lnum_t vtx_id = -1;
453 cs_lnum_t n_new_vertices = *p_n_new_vertices;
454
455 assert(inter.curv_abs >= 0.0);
456 assert(inter.curv_abs <= 1.0);
457
458 if (inter.curv_abs <= 0.0)
459 vtx_id = vtx_couple[0] - 1;
460
461 else if (inter.curv_abs >= 1.0)
462 vtx_id = vtx_couple[1] - 1;
463
464 else {
465
466 assert(inter.curv_abs > 0 && inter.curv_abs < 1.0);
467 vtx_id = n_init_vertices + n_new_vertices;
468 n_new_vertices++;
469
470 }
471
472 assert(vtx_id != -1);
473
474 *p_n_new_vertices = n_new_vertices;
475
476 return vtx_id;
477 }
478
479
480 /*----------------------------------------------------------------------------
481 * Test if we have to continue to spread the tag associate to each vertex
482 *
483 * parameters:
484 * n_vertices <-- local number of vertices
485 * prev_vtx_tag <-- previous tag for each vertex
486 * vtx_tag <-- tag for each vertex
487 *
488 * returns:
489 * 1 for true, 0 for false
490 *---------------------------------------------------------------------------*/
491
492 static int
_is_spread_not_converged(cs_lnum_t n_vertices,const cs_gnum_t prev_vtx_tag[],const cs_gnum_t vtx_tag[])493 _is_spread_not_converged(cs_lnum_t n_vertices,
494 const cs_gnum_t prev_vtx_tag[],
495 const cs_gnum_t vtx_tag[])
496 {
497 int have_to_continue = 0;
498
499 for (cs_lnum_t i = 0; i < n_vertices; i++) {
500 if (vtx_tag[i] != prev_vtx_tag[i]) {
501 have_to_continue = 1;
502 break;
503 }
504 }
505
506 return have_to_continue;
507 }
508
509 /*----------------------------------------------------------------------------
510 * Spread the tag associated to each vertex according the rule:
511 * Between two equivalent vertices, the tag associated to each considered
512 * vertex is equal to the minimal global number.
513 *
514 * parameters:
515 * n_vertices <-- local number of vertices
516 * vtx_eset <-- structure dealing with vertices equivalences
517 * vtx_tag <-> tag for each vertex
518 *---------------------------------------------------------------------------*/
519
520 static void
_spread_tag(cs_lnum_t n_vertices,const cs_join_eset_t * vtx_eset,cs_gnum_t vtx_tag[])521 _spread_tag(cs_lnum_t n_vertices,
522 const cs_join_eset_t *vtx_eset,
523 cs_gnum_t vtx_tag[])
524 {
525 cs_lnum_t i, v1_id, v2_id;
526 cs_gnum_t v1_gnum, v2_gnum;
527 cs_lnum_t *equiv_lst = vtx_eset->equiv_couple;
528
529 for (i = 0; i < vtx_eset->n_equiv; i++) {
530
531 v1_id = equiv_lst[2*i] - 1, v2_id = equiv_lst[2*i+1] - 1;
532 assert(v1_id < n_vertices);
533 assert(v1_id < n_vertices);
534 v1_gnum = vtx_tag[v1_id], v2_gnum = vtx_tag[v2_id];
535
536 if (v1_gnum != v2_gnum) {
537
538 cs_gnum_t min_gnum = CS_MIN(v1_gnum, v2_gnum);
539
540 vtx_tag[v1_id] = min_gnum;
541 vtx_tag[v2_id] = min_gnum;
542 }
543
544 } /* End of loop on vertex equivalences */
545 }
546
547 /*----------------------------------------------------------------------------
548 * Define an array wich keeps the new vertex id of each vertex.
549 *
550 * If two vertices have the same vertex id, they should merge.
551 *
552 * parameters:
553 * vtx_eset <-- structure dealing with vertex equivalences
554 * n_vertices <-- local number of vertices
555 * prev_vtx_tag <-> previous tag for each vertex
556 * vtx_tag <-> tag for each vertex
557 *---------------------------------------------------------------------------*/
558
559 static void
_local_spread(const cs_join_eset_t * vtx_eset,cs_lnum_t n_vertices,cs_gnum_t prev_vtx_tag[],cs_gnum_t vtx_tag[])560 _local_spread(const cs_join_eset_t *vtx_eset,
561 cs_lnum_t n_vertices,
562 cs_gnum_t prev_vtx_tag[],
563 cs_gnum_t vtx_tag[])
564 {
565 cs_lnum_t i;
566
567 _loc_merge_counter++;
568
569 _spread_tag(n_vertices, vtx_eset, vtx_tag);
570
571 while (_is_spread_not_converged(n_vertices, prev_vtx_tag, vtx_tag)) {
572
573 _loc_merge_counter++;
574
575 if (_loc_merge_counter > CS_JOIN_MERGE_MAX_LOC_ITERS)
576 bft_error(__FILE__, __LINE__, 0,
577 _("\n The authorized maximum number of iterations "
578 " for the merge of vertices has been reached.\n"
579 " Local counter on iteration : %d (MAX =%d)\n"
580 " Check the fraction parameter.\n"),
581 _loc_merge_counter, CS_JOIN_MERGE_MAX_LOC_ITERS);
582
583 for (i = 0; i < n_vertices; i++)
584 prev_vtx_tag[i] = vtx_tag[i];
585
586 _spread_tag(n_vertices, vtx_eset, vtx_tag);
587 }
588 }
589
590 #if defined(HAVE_MPI)
591
592 /*----------------------------------------------------------------------------
593 * Exchange local vtx_tag buffer over the ranks and update global vtx_tag
594 * buffers. Apply modifications observed on the global vtx_tag to the local
595 * vtx_tag.
596 *
597 * parameters:
598 * block_size <-- size of block for the current rank
599 * d <-- all to all distributor
600 * work <-- local cs_join_mesh_t structure which has initial
601 * vertex data
602 * vtx_tag <-> local vtx_tag for the local vertices
603 * glob_vtx_tag <-> global vtx_tag affected to the local rank
604 * (size: block_size)
605 * prev_glob_vtx_tag <-> same but for the previous iteration
606 * recv2glob <-> buffer used to place correctly receive elements
607 * send_count <-> buffer used to count the number of elts to send
608 * send_shift <-> index on ranks of the elements to send
609 * send_glob_buffer <-> buffer used to save elements to send
610 * recv_count <-> buffer used to count the number of elts to receive
611 * recv_shift <-> index on ranks of the elements to receive
612 * recv_glob_buffer <-> buffer used to save elements to receive
613 *
614 * returns:
615 * true if we have to continue the spread, false otherwise.
616 *---------------------------------------------------------------------------*/
617
618 static bool
_global_spread(cs_lnum_t block_size,cs_all_to_all_t * d,const cs_join_mesh_t * work,cs_gnum_t vtx_tag[],cs_gnum_t glob_vtx_tag[],cs_gnum_t prev_glob_vtx_tag[],cs_gnum_t recv2glob[],cs_gnum_t send_glob_buffer[],cs_gnum_t recv_glob_buffer[])619 _global_spread(cs_lnum_t block_size,
620 cs_all_to_all_t *d,
621 const cs_join_mesh_t *work,
622 cs_gnum_t vtx_tag[],
623 cs_gnum_t glob_vtx_tag[],
624 cs_gnum_t prev_glob_vtx_tag[],
625 cs_gnum_t recv2glob[],
626 cs_gnum_t send_glob_buffer[],
627 cs_gnum_t recv_glob_buffer[])
628 {
629 int global_value;
630
631 cs_lnum_t n_vertices = work->n_vertices;
632 MPI_Comm mpi_comm = cs_glob_mpi_comm;
633
634 _glob_merge_counter++;
635
636 /* Push modifications in local vtx_tag to the global vtx_tag */
637
638 cs_all_to_all_copy_array(d,
639 CS_GNUM_TYPE,
640 1,
641 false, /* reverse */
642 vtx_tag,
643 recv_glob_buffer);
644
645 /* Apply update to glob_vtx_tag */
646
647 cs_lnum_t n_recv = cs_all_to_all_n_elts_dest(d);
648
649 for (cs_lnum_t i = 0; i < n_recv; i++) {
650 cs_lnum_t cur_id = recv2glob[i];
651 glob_vtx_tag[cur_id] = CS_MIN(glob_vtx_tag[cur_id], recv_glob_buffer[i]);
652 }
653
654 int local_value = _is_spread_not_converged(block_size,
655 prev_glob_vtx_tag,
656 glob_vtx_tag);
657
658 MPI_Allreduce(&local_value, &global_value, 1, MPI_INT, MPI_SUM, mpi_comm);
659
660 if (global_value > 0) { /* Store the current state as the previous one
661 Update local vtx_tag */
662
663 if (_glob_merge_counter > CS_JOIN_MERGE_MAX_GLOB_ITERS)
664 bft_error(__FILE__, __LINE__, 0,
665 _("\n The authorized maximum number of iterations "
666 " for the merge of vertices has been reached.\n"
667 " Global counter on iteration : %d (MAX =%d)\n"
668 " Check the fraction parameter.\n"),
669 _glob_merge_counter, CS_JOIN_MERGE_MAX_GLOB_ITERS);
670
671 for (cs_lnum_t i = 0; i < block_size; i++)
672 prev_glob_vtx_tag[i] = glob_vtx_tag[i];
673
674 for (cs_lnum_t i = 0; i < n_recv; i++)
675 recv_glob_buffer[i] = glob_vtx_tag[recv2glob[i]];
676
677 cs_all_to_all_copy_array(d,
678 CS_GNUM_TYPE,
679 1,
680 true, /* reverse */
681 recv_glob_buffer,
682 send_glob_buffer);
683
684 /* Update vtx_tag */
685
686 for (cs_lnum_t i = 0; i < n_vertices; i++)
687 vtx_tag[i] = CS_MIN(send_glob_buffer[i], vtx_tag[i]);
688
689 return true;
690
691 } /* End if prev_glob_vtx_tag != glob_vtx_tag */
692
693 else
694 return false; /* No need to continue */
695 }
696
697 /*----------------------------------------------------------------------------
698 * Initialize and allocate buffers for the tag operation in parallel mode.
699 *
700 * parameters:
701 * bi <-- block distribution info
702 * work <-- local cs_join_mesh_t structure which has
703 * initial vertex data
704 * p_all_to_all_d <-> pointer to all to all distributor
705 * p_recv2glob <-> buf. for putting correctly received elements
706 * p_glob_vtx_tag <-> vtx_tag locally treated (size = block_size)
707 * p_prev_glob_vtx_tag <-> idem but for the previous iteration
708 *---------------------------------------------------------------------------*/
709
710 static void
_parall_tag_init(cs_block_dist_info_t bi,const cs_join_mesh_t * work,cs_all_to_all_t ** p_all_to_all_d,cs_gnum_t * p_recv2glob[],cs_gnum_t * p_glob_vtx_tag[],cs_gnum_t * p_prev_glob_vtx_tag[])711 _parall_tag_init(cs_block_dist_info_t bi,
712 const cs_join_mesh_t *work,
713 cs_all_to_all_t **p_all_to_all_d,
714 cs_gnum_t *p_recv2glob[],
715 cs_gnum_t *p_glob_vtx_tag[],
716 cs_gnum_t *p_prev_glob_vtx_tag[])
717 {
718 cs_lnum_t n_vertices = work->n_vertices;
719 cs_gnum_t *glob_vtx_tag = NULL, *prev_glob_vtx_tag = NULL;
720 MPI_Comm mpi_comm = cs_glob_mpi_comm;
721
722 const int n_ranks = cs_glob_n_ranks;
723 const int local_rank = CS_MAX(cs_glob_rank_id, 0);
724 const cs_gnum_t _n_ranks = n_ranks, _local_rank = local_rank;
725
726 /* Allocate and initialize vtx_tag associated to the local rank */
727
728 BFT_MALLOC(glob_vtx_tag, bi.block_size, cs_gnum_t);
729 BFT_MALLOC(prev_glob_vtx_tag, bi.block_size, cs_gnum_t);
730
731 for (cs_lnum_t i = 0; i < bi.block_size; i++) {
732 cs_gnum_t gi = i;
733 prev_glob_vtx_tag[i] = gi*_n_ranks + _local_rank + 1;
734 glob_vtx_tag[i] = gi*_n_ranks + _local_rank + 1;
735 }
736
737 /* Create all to all distributor */
738
739 int *dest_rank;
740 BFT_MALLOC(dest_rank, n_vertices, int);
741
742 cs_gnum_t *wv_gnum;
743 BFT_MALLOC(wv_gnum, n_vertices, cs_gnum_t);
744
745 for (cs_lnum_t i = 0; i < n_vertices; i++) {
746 dest_rank[i] = (work->vertices[i].gnum - 1) % _n_ranks;
747 wv_gnum[i] = (work->vertices[i].gnum - 1) / _n_ranks;
748 }
749
750 cs_all_to_all_t *d
751 = cs_all_to_all_create(n_vertices,
752 0, /* flags */
753 NULL,
754 dest_rank,
755 mpi_comm);
756
757 cs_all_to_all_transfer_dest_rank(d, &dest_rank);
758
759 /* Allocate and define recv2glob */
760
761 cs_gnum_t *recv2glob = cs_all_to_all_copy_array(d,
762 CS_GNUM_TYPE,
763 1,
764 false, /* reverse */
765 wv_gnum,
766 NULL);
767
768 BFT_FREE(wv_gnum);
769
770 /* Return pointers */
771
772 *p_all_to_all_d = d;
773 *p_recv2glob = recv2glob;
774 *p_glob_vtx_tag = glob_vtx_tag;
775 *p_prev_glob_vtx_tag = prev_glob_vtx_tag;
776
777 }
778
779 #endif /* HAVE_MPI */
780
781 /*----------------------------------------------------------------------------
782 * Tag with the same number all the vertices which might be merged together
783 *
784 * parameters:
785 * n_g_vertices_tot <-- global number of vertices to consider for the
786 * merge operation (existing + created vertices)
787 * vtx_eset <-- structure dealing with vertex equivalences
788 * work <-- local cs_join_mesh_t structure which has initial
789 * vertex data
790 * verbosity <-- level of detail in information display
791 * p_vtx_tag --> pointer to the vtx_tag for the local vertices
792 *---------------------------------------------------------------------------*/
793
794 static void
_tag_equiv_vertices(cs_gnum_t n_g_vertices_tot,const cs_join_eset_t * vtx_eset,const cs_join_mesh_t * work,int verbosity,cs_gnum_t * p_vtx_tag[])795 _tag_equiv_vertices(cs_gnum_t n_g_vertices_tot,
796 const cs_join_eset_t *vtx_eset,
797 const cs_join_mesh_t *work,
798 int verbosity,
799 cs_gnum_t *p_vtx_tag[])
800 {
801 cs_lnum_t i;
802
803 cs_gnum_t *vtx_tag = NULL;
804 cs_gnum_t *prev_vtx_tag = NULL;
805 FILE *logfile = cs_glob_join_log;
806
807 const cs_lnum_t n_vertices = work->n_vertices;
808 const int n_ranks = cs_glob_n_ranks;
809
810 /* Local initialization : we tag each vertex by its global number */
811
812 BFT_MALLOC(prev_vtx_tag, n_vertices, cs_gnum_t);
813 BFT_MALLOC(vtx_tag, n_vertices, cs_gnum_t);
814
815 for (i = 0; i < work->n_vertices; i++) {
816
817 cs_gnum_t v_gnum = work->vertices[i].gnum;
818
819 vtx_tag[i] = v_gnum;
820 prev_vtx_tag[i] = v_gnum;
821
822 }
823
824 #if 0 && defined(DEBUG) && !defined(NDEBUG)
825 for (i = 0; i < n_vertices; i++)
826 fprintf(logfile, " Initial vtx_tag[%6d] = %9llu\n",
827 i, (unsigned long long)vtx_tag[i]);
828 fflush(logfile);
829 #endif
830
831 /* Compute vtx_tag */
832
833 _local_spread(vtx_eset, n_vertices, prev_vtx_tag, vtx_tag);
834
835 #if defined(HAVE_MPI)
836
837 if (n_ranks > 1) { /* Parallel treatment */
838
839 bool go_on;
840
841 cs_gnum_t *glob_vtx_tag = NULL, *prev_glob_vtx_tag = NULL;
842 cs_gnum_t *recv2glob;
843
844 const int local_rank = CS_MAX(cs_glob_rank_id, 0);
845
846 cs_block_dist_info_t bi = cs_block_dist_compute_sizes(local_rank,
847 n_ranks,
848 1,
849 0,
850 n_g_vertices_tot);
851 cs_all_to_all_t *d = NULL;
852
853 _parall_tag_init(bi,
854 work,
855 &d,
856 &recv2glob,
857 &glob_vtx_tag,
858 &prev_glob_vtx_tag);
859
860 cs_lnum_t n_recv = cs_all_to_all_n_elts_dest(d);
861 cs_gnum_t *send_glob_buffer, *recv_glob_buffer;
862 BFT_MALLOC(send_glob_buffer, n_vertices, cs_gnum_t);
863 BFT_MALLOC(recv_glob_buffer, n_recv, cs_gnum_t);
864
865 go_on = _global_spread(bi.block_size,
866 d,
867 work,
868 vtx_tag,
869 glob_vtx_tag,
870 prev_glob_vtx_tag,
871 recv2glob,
872 send_glob_buffer,
873 recv_glob_buffer);
874
875 while (go_on == true) {
876
877 /* Local convergence of vtx_tag */
878
879 _local_spread(vtx_eset, n_vertices, prev_vtx_tag, vtx_tag);
880
881 /* Global update and test to continue */
882
883 go_on = _global_spread(bi.block_size,
884 d,
885 work,
886 vtx_tag,
887 glob_vtx_tag,
888 prev_glob_vtx_tag,
889 recv2glob,
890 send_glob_buffer,
891 recv_glob_buffer);
892
893 }
894
895 /* Partial free */
896
897 BFT_FREE(glob_vtx_tag);
898 BFT_FREE(prev_glob_vtx_tag);
899 BFT_FREE(send_glob_buffer);
900 BFT_FREE(recv2glob);
901 BFT_FREE(recv_glob_buffer);
902
903 cs_all_to_all_destroy(&d);
904
905 } /* End of parallel treatment */
906
907 #endif
908
909 BFT_FREE(prev_vtx_tag);
910
911 if (verbosity > 3) {
912 fprintf(logfile,
913 "\n Number of local iterations to converge on vertex"
914 " equivalences: %3d\n", _loc_merge_counter);
915 if (n_ranks > 1)
916 fprintf(logfile,
917 " Number of global iterations to converge on vertex"
918 " equivalences: %3d\n\n", _glob_merge_counter);
919 fflush(logfile);
920 }
921
922 #if 0 && defined(DEBUG) && !defined(NDEBUG)
923 if (logfile != NULL) {
924 for (i = 0; i < n_vertices; i++)
925 fprintf(logfile, " Final vtx_tag[%6d] = %9llu\n",
926 i, (unsigned long long)vtx_tag[i]);
927 fflush(logfile);
928 }
929 #endif
930
931 /* Returns pointer */
932
933 *p_vtx_tag = vtx_tag;
934 }
935
936 #if defined(HAVE_MPI)
937
938 /*----------------------------------------------------------------------------
939 * Build in parallel a cs_join_gset_t structure to store all the potential
940 * merges between vertices and its associated cs_join_vertex_t structure.
941 *
942 * parameters:
943 * work <-- local cs_join_mesh_t structure which
944 * has initial vertex data
945 * vtx_tag <-- local vtx_tag for the local vertices
946 * d <-> all to all distributor
947 * p_vtx_merge_data <-> a pointer to a cs_join_vertex_t structure which
948 * stores data about merged vertices
949 * p_merge_set <-> pointer to a cs_join_gset_t struct. storing the
950 * evolution of each global vtx number
951 *---------------------------------------------------------------------------*/
952
953 static void
_build_parall_merge_structures(const cs_join_mesh_t * work,const cs_gnum_t vtx_tag[],cs_all_to_all_t * d,cs_join_vertex_t * p_vtx_merge_data[],cs_join_gset_t ** p_merge_set)954 _build_parall_merge_structures(const cs_join_mesh_t *work,
955 const cs_gnum_t vtx_tag[],
956 cs_all_to_all_t *d,
957 cs_join_vertex_t *p_vtx_merge_data[],
958 cs_join_gset_t **p_merge_set)
959 {
960 /* Distribute vertex tags */
961
962 cs_gnum_t *recv_gbuf = cs_all_to_all_copy_array(d,
963 CS_GNUM_TYPE,
964 1,
965 false, /* reverse */
966 vtx_tag,
967 NULL);
968
969 /* Allocate and build send_vtx_data, receive recv_vtx_data. */
970
971 cs_join_vertex_t *recv_vtx_data
972 = cs_all_to_all_copy_array(d,
973 CS_CHAR,
974 sizeof(cs_join_vertex_t),
975 false, /* reverse */
976 work->vertices,
977 NULL);
978
979 /* Build merge set */
980
981 const cs_lnum_t n_recv = cs_all_to_all_n_elts_dest(d);
982
983 cs_join_gset_t *merge_set
984 = cs_join_gset_create_from_tag(n_recv, recv_gbuf);
985
986 cs_join_gset_sort_sublist(merge_set);
987
988 /* Free memory */
989
990 BFT_FREE(recv_gbuf);
991
992 #if 0 && defined(DEBUG) && !defined(NDEBUG)
993 if (cs_glob_join_log != NULL) {
994 FILE *logfile = cs_glob_join_log;
995 fprintf(logfile,
996 "\n Number of vertices to treat for the merge step: %d\n",
997 recv_shift[n_ranks]);
998 fprintf(logfile,
999 " List of vertices to treat:\n");
1000 for (i = 0; i < recv_shift[n_ranks]; i++) {
1001 fprintf(logfile, " %9d - ", i);
1002 cs_join_mesh_dump_vertex(logfile, recv_vtx_data[i]);
1003 }
1004 fflush(logfile);
1005 }
1006 #endif
1007
1008 /* Set return pointers */
1009
1010 *p_merge_set = merge_set;
1011 *p_vtx_merge_data = recv_vtx_data;
1012 }
1013
1014 #endif /* HAVE_MPI */
1015
1016 /*----------------------------------------------------------------------------
1017 * Get the resulting cs_join_vertex_t structure after the merge of a set
1018 * of vertices.
1019 *
1020 * parameters:
1021 * n_elts <-- size of the set
1022 * set <-- set of vertices
1023 *
1024 * returns:
1025 * a cs_join_vertex_t structure for the resulting vertex
1026 *---------------------------------------------------------------------------*/
1027
1028 static cs_join_vertex_t
_compute_merged_vertex(cs_lnum_t n_elts,const cs_join_vertex_t set[])1029 _compute_merged_vertex(cs_lnum_t n_elts,
1030 const cs_join_vertex_t set[])
1031 {
1032 cs_lnum_t i, k;
1033 cs_real_t w;
1034 cs_join_vertex_t mvtx;
1035
1036 cs_real_t denum = 0.0;
1037
1038 #if defined(DEBUG) && !defined(NDEBUG)
1039 /* Avoid Valgrind warnings in byte copies due to padding */
1040 memset(&mvtx, 0, sizeof(cs_join_vertex_t));
1041 #endif
1042
1043 assert(n_elts > 0);
1044
1045 /* Initialize cs_join_vertex_t structure */
1046
1047 mvtx.state = CS_JOIN_STATE_UNDEF;
1048 mvtx.gnum = set[0].gnum;
1049 mvtx.tolerance = set[0].tolerance;
1050
1051 for (k = 0; k < 3; k++)
1052 mvtx.coord[k] = 0.0;
1053
1054 /* Compute the resulting vertex */
1055
1056 for (i = 0; i < n_elts; i++) {
1057
1058 mvtx.tolerance = CS_MIN(set[i].tolerance, mvtx.tolerance);
1059 mvtx.gnum = CS_MIN(set[i].gnum, mvtx.gnum);
1060 mvtx.state = CS_MAX(set[i].state, mvtx.state);
1061
1062 /* Compute the resulting coordinates of the merged vertices */
1063
1064 #if CS_JOIN_MERGE_INV_TOL
1065 w = 1.0/set[i].tolerance;
1066 #else
1067 w = 1.0;
1068 #endif
1069 denum += w;
1070
1071 for (k = 0; k < 3; k++)
1072 mvtx.coord[k] += w * set[i].coord[k];
1073
1074 }
1075
1076 for (k = 0; k < 3; k++)
1077 mvtx.coord[k] /= denum;
1078
1079 if (mvtx.state == CS_JOIN_STATE_ORIGIN)
1080 mvtx.state = CS_JOIN_STATE_MERGE;
1081 else if (mvtx.state == CS_JOIN_STATE_PERIO)
1082 mvtx.state = CS_JOIN_STATE_PERIO_MERGE;
1083
1084 return mvtx;
1085 }
1086
1087 /*----------------------------------------------------------------------------
1088 * Merge between identical vertices.
1089 *
1090 * Only the vertex numbering and the related tolerance may be different.
1091 * Store new data associated to the merged vertices in vertices array.
1092 *
1093 * parameters:
1094 * param <-- set of user-defined parameters
1095 * merge_set <-> a pointer to a cs_join_vertex_t structure which
1096 * stores data about merged vertices
1097 * n_vertices <-- number of vertices in vertices array
1098 * vertices <-> array of cs_join_vertex_t structures
1099 * equiv_gnum --> equivalence between id in vertices (same global number
1100 * initially or identical vertices: same coordinates)
1101 *---------------------------------------------------------------------------*/
1102
1103 static void
_pre_merge(cs_join_param_t param,cs_join_gset_t * merge_set,cs_join_vertex_t vertices[],cs_join_gset_t ** p_equiv_gnum)1104 _pre_merge(cs_join_param_t param,
1105 cs_join_gset_t *merge_set,
1106 cs_join_vertex_t vertices[],
1107 cs_join_gset_t **p_equiv_gnum)
1108 {
1109 cs_lnum_t i, j, j1, j2, k, k1, k2, n_sub_elts;
1110 cs_real_t deltad, deltat, limit, min_tol;
1111 cs_join_vertex_t mvtx, coupled_vertices[2];
1112
1113 cs_lnum_t max_n_sub_elts = 0, n_local_pre_merge = 0;
1114 cs_lnum_t *merge_index = merge_set->index;
1115 cs_gnum_t *merge_list = merge_set->g_list;
1116 cs_gnum_t *sub_list = NULL, *init_list = NULL;
1117 cs_join_gset_t *equiv_gnum = NULL;
1118
1119 const cs_real_t pmf = param.pre_merge_factor;
1120
1121 cs_join_gset_sort_sublist(merge_set);
1122
1123 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1124 {
1125 int len;
1126 FILE *dbg_file = NULL;
1127 char *filename = NULL;
1128
1129 len = strlen("JoinDBG_InitMergeSet.dat")+1+2+4;
1130 BFT_MALLOC(filename, len, char);
1131 sprintf(filename, "Join%02dDBG_InitMergeSet%04d.dat",
1132 param.num, CS_MAX(cs_glob_rank_id, 0));
1133 dbg_file = fopen(filename, "w");
1134
1135 cs_join_gset_dump(dbg_file, merge_set);
1136
1137 fflush(dbg_file);
1138 BFT_FREE(filename);
1139 fclose(dbg_file);
1140 }
1141 #endif
1142
1143 /* Compute the max. size of a sub list */
1144
1145 for (i = 0; i < merge_set->n_elts; i++)
1146 max_n_sub_elts = CS_MAX(max_n_sub_elts,
1147 merge_index[i+1] - merge_index[i]);
1148
1149 BFT_MALLOC(sub_list, max_n_sub_elts, cs_gnum_t);
1150
1151 /* Store initial merge list */
1152
1153 BFT_MALLOC(init_list, merge_index[merge_set->n_elts], cs_gnum_t);
1154
1155 for (i = 0; i < merge_index[merge_set->n_elts]; i++)
1156 init_list[i] = merge_list[i];
1157
1158 /* Apply merge */
1159
1160 for (i = 0; i < merge_set->n_elts; i++) {
1161
1162 cs_lnum_t f_s = merge_index[i];
1163 cs_lnum_t f_e = merge_index[i+1];
1164
1165 n_sub_elts = f_e - f_s;
1166
1167 for (j = f_s, k = 0; j < f_e; j++, k++)
1168 sub_list[k] = merge_list[j];
1169
1170 for (j1 = 0; j1 < n_sub_elts - 1; j1++) {
1171
1172 cs_lnum_t v1_id = sub_list[j1];
1173 cs_join_vertex_t v1 = vertices[v1_id];
1174
1175 for (j2 = j1 + 1; j2 < n_sub_elts; j2++) {
1176
1177 cs_lnum_t v2_id = sub_list[j2];
1178 cs_join_vertex_t v2 = vertices[v2_id];
1179
1180 if (v1.gnum == v2.gnum) { /* Possible if n_ranks > 1 */
1181
1182 if (sub_list[j1] < sub_list[j2])
1183 k1 = j1, k2 = j2;
1184 else
1185 k1 = j2, k2 = j1;
1186
1187 for (k = 0; k < n_sub_elts; k++)
1188 if (sub_list[k] == sub_list[k2])
1189 sub_list[k] = sub_list[k1];
1190
1191 }
1192 else {
1193
1194 min_tol = CS_MIN(v1.tolerance, v2.tolerance);
1195 limit = min_tol * pmf;
1196 deltat = CS_ABS(v1.tolerance - v2.tolerance);
1197
1198 if (deltat < limit) {
1199
1200 deltad = _compute_length(v1, v2);
1201
1202 if (deltad < limit) { /* Do a pre-merge */
1203
1204 n_local_pre_merge++;
1205
1206 if (v1.gnum < v2.gnum)
1207 k1 = j1, k2 = j2;
1208 else
1209 k1 = j2, k2 = j1;
1210
1211 for (k = 0; k < n_sub_elts; k++)
1212 if (sub_list[k] == sub_list[k2])
1213 sub_list[k] = sub_list[k1];
1214
1215 coupled_vertices[0] = v1, coupled_vertices[1] = v2;
1216 mvtx = _compute_merged_vertex(2, coupled_vertices);
1217 vertices[v1_id] = mvtx;
1218 vertices[v2_id] = mvtx;
1219
1220 } /* deltad < limit */
1221
1222 } /* deltat < limit */
1223
1224 } /* v1.gnum != v2.gnum */
1225
1226 } /* End of loop on j2 */
1227 } /* End of loop on j1 */
1228
1229 /* Update vertices */
1230
1231 for (j = f_s, k = 0; j < f_e; j++, k++)
1232 vertices[merge_list[j]] = vertices[sub_list[k]];
1233
1234 /* Update merge list */
1235
1236 for (j = f_s, k = 0; j < f_e; j++, k++)
1237 merge_list[j] = sub_list[k];
1238
1239 } /* End of loop on merge_set elements */
1240
1241 /* Keep equivalences between identical vertices in equiv_gnum */
1242
1243 equiv_gnum = cs_join_gset_create_by_equiv(merge_set, init_list);
1244
1245 /* Clean merge set */
1246
1247 cs_join_gset_clean(merge_set);
1248
1249 /* Display information about the joining */
1250
1251 if (param.verbosity > 0) {
1252
1253 cs_gnum_t n_g_counter = n_local_pre_merge;
1254 cs_parall_counter(&n_g_counter, 1);
1255
1256 bft_printf(_("\n Pre-merge for %llu global element couples.\n"),
1257 (unsigned long long)n_g_counter);
1258
1259 if (param.verbosity > 2) {
1260 fprintf(cs_glob_join_log, "\n Local number of pre-merges: %ld\n",
1261 (long)n_local_pre_merge);
1262 }
1263 }
1264
1265 /* Free memory */
1266
1267 BFT_FREE(sub_list);
1268 BFT_FREE(init_list);
1269
1270 /* Return pointer */
1271
1272 *p_equiv_gnum = equiv_gnum;
1273 }
1274
1275 /*----------------------------------------------------------------------------
1276 * Check if all vertices in the set include the ref_vertex in their tolerance.
1277 *
1278 * parameters:
1279 * set_size <-- size of set of vertices
1280 * vertices <-- set of vertices to check
1281 * ref_vertex <-- ref. vertex
1282 *
1283 * returns:
1284 * true if all vertices have ref_vertex in their tolerance, false otherwise
1285 *---------------------------------------------------------------------------*/
1286
1287 static bool
_is_in_tolerance(cs_lnum_t set_size,const cs_join_vertex_t set[],cs_join_vertex_t ref_vertex)1288 _is_in_tolerance(cs_lnum_t set_size,
1289 const cs_join_vertex_t set[],
1290 cs_join_vertex_t ref_vertex)
1291 {
1292 cs_lnum_t i;
1293
1294 for (i = 0; i < set_size; i++) {
1295
1296 cs_real_t d2ref = _compute_length(set[i], ref_vertex);
1297 cs_real_t tolerance = set[i].tolerance * cs_join_tol_eps_coef2;
1298
1299 if (d2ref > tolerance)
1300 return false;
1301
1302 }
1303
1304 return true;
1305 }
1306
1307 /*----------------------------------------------------------------------------
1308 * Test if we have to continue to the subset building
1309 *
1310 * parameters:
1311 * set_size <-- size of set
1312 * prev_num <-> array used to store previous subset_num
1313 * new_num <-> number associated to each vertices of the set
1314 *
1315 * returns:
1316 * true or false
1317 *---------------------------------------------------------------------------*/
1318
1319 static bool
_continue_subset_building(int set_size,const cs_lnum_t prev_num[],const cs_lnum_t new_num[])1320 _continue_subset_building(int set_size,
1321 const cs_lnum_t prev_num[],
1322 const cs_lnum_t new_num[])
1323 {
1324 int i;
1325
1326 for (i = 0; i < set_size; i++)
1327 if (new_num[i] != prev_num[i])
1328 return true;
1329
1330 return false;
1331 }
1332
1333 /*----------------------------------------------------------------------------
1334 * Define subsets of vertices.
1335 *
1336 * parameters:
1337 * set_size <-- size of set
1338 * state <-- array keeping the state of the link
1339 * subset_num <-> number associated to each vertices of the set
1340 *---------------------------------------------------------------------------*/
1341
1342 static void
_iter_subset_building(cs_lnum_t set_size,const cs_lnum_t state[],cs_lnum_t subset_num[])1343 _iter_subset_building(cs_lnum_t set_size,
1344 const cs_lnum_t state[],
1345 cs_lnum_t subset_num[])
1346 {
1347 cs_lnum_t i1, i2, k;
1348
1349 for (k = 0, i1 = 0; i1 < set_size-1; i1++) {
1350 for (i2 = i1 + 1; i2 < set_size; i2++, k++) {
1351
1352 if (state[k] == 1) { /* v1 - v2 are in tolerance each other */
1353
1354 int _min = CS_MIN(subset_num[i1], subset_num[i2]);
1355
1356 subset_num[i1] = _min;
1357 subset_num[i2] = _min;
1358
1359 }
1360
1361 }
1362 }
1363
1364 }
1365
1366 /*----------------------------------------------------------------------------
1367 * Define subsets of vertices.
1368 *
1369 * parameters:
1370 * set_size <-- size of set
1371 * state <-- array keeping the state of the link
1372 * prev_num <-> array used to store previous subset_num
1373 * subset_num <-> number associated to each vertices of the set
1374 *---------------------------------------------------------------------------*/
1375
1376 static void
_build_subsets(cs_lnum_t set_size,const cs_lnum_t state[],cs_lnum_t prev_num[],cs_lnum_t subset_num[])1377 _build_subsets(cs_lnum_t set_size,
1378 const cs_lnum_t state[],
1379 cs_lnum_t prev_num[],
1380 cs_lnum_t subset_num[])
1381 {
1382 int i;
1383 cs_lnum_t n_loops = 0;
1384
1385 /* Initialize */
1386
1387 for (i = 0; i < set_size; i++) {
1388 subset_num[i] = i+1;
1389 prev_num[i] = subset_num[i];
1390 }
1391
1392 _iter_subset_building(set_size, state, subset_num);
1393
1394 while ( _continue_subset_building(set_size, prev_num, subset_num)
1395 && n_loops < CS_JOIN_MERGE_MAX_LOC_ITERS ) {
1396
1397 n_loops++;
1398
1399 for (i = 0; i < set_size; i++)
1400 prev_num[i] = subset_num[i];
1401
1402 _iter_subset_building(set_size, state, subset_num);
1403
1404 }
1405
1406 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1407 if (cs_glob_join_log != NULL && n_loops >= CS_JOIN_MERGE_MAX_LOC_ITERS)
1408 fprintf(cs_glob_join_log,
1409 "WARNING max sub_loops to build subset reached\n");
1410 #endif
1411
1412 }
1413
1414 /*----------------------------------------------------------------------------
1415 * Check if each subset is consistent with tolerance of vertices
1416 * If a transitivity is found, modify the state of the link
1417 * state = 1 (each other in their tolerance)
1418 * = 0 (not each other in their tolerance)
1419 *
1420 * parameters:
1421 * set_size <-- size of set
1422 * set <-- pointer to a set of vertices
1423 * state <-> array keeping the state of the link
1424 * subset_num <-> number associated to each vertices of the set
1425 * issues <-> numbering of inconsistent subsets
1426 * verbosity <-- level of verbosity
1427 *
1428 * returns:
1429 * number of subsets not consistent
1430 *---------------------------------------------------------------------------*/
1431
1432 static cs_lnum_t
_check_tol_consistency(cs_lnum_t set_size,const cs_join_vertex_t set[],cs_lnum_t state[],cs_lnum_t subset_num[],cs_lnum_t issues[],cs_lnum_t verbosity)1433 _check_tol_consistency(cs_lnum_t set_size,
1434 const cs_join_vertex_t set[],
1435 cs_lnum_t state[],
1436 cs_lnum_t subset_num[],
1437 cs_lnum_t issues[],
1438 cs_lnum_t verbosity)
1439 {
1440 cs_lnum_t i1, i2, j, k;
1441
1442 cs_lnum_t n_issues = 0;
1443 FILE *logfile = cs_glob_join_log;
1444
1445 for (k = 0, i1 = 0; i1 < set_size-1; i1++) {
1446 for (i2 = i1 + 1; i2 < set_size; i2++, k++) {
1447
1448 if (state[k] == 0) {
1449
1450 if (subset_num[i1] == subset_num[i2]) {
1451
1452 if (verbosity > 4)
1453 fprintf(logfile,
1454 " Transitivity detected between (%llu, %llu)\n",
1455 (unsigned long long)set[i1].gnum,
1456 (unsigned long long)set[i2].gnum);
1457
1458 for (j = 0; j < n_issues; j++)
1459 if (issues[j] == subset_num[i1])
1460 break;
1461 if (j == n_issues)
1462 issues[n_issues++] = subset_num[i1];
1463
1464 }
1465 }
1466
1467 } /* End of loop on i2 */
1468 } /* ENd of loop on i1 */
1469
1470 return n_issues; /* Not a subset number */
1471 }
1472
1473 /*----------------------------------------------------------------------------
1474 * Check if the merged vertex related to a subset is consistent with tolerance
1475 * of each vertex of the subset.
1476 *
1477 * parameters:
1478 * set_size <-- size of set
1479 * subset_num <-- number associated to each vertices of the set
1480 * set <-- pointer to a set of vertices
1481 * merge_set <-> merged vertex related to each subset
1482 * work_set <-> work array of vertices
1483 *
1484 * returns:
1485 * true if all subsets are consistent otherwise false
1486 *---------------------------------------------------------------------------*/
1487
1488 static bool
_check_subset_consistency(cs_lnum_t set_size,const cs_lnum_t subset_num[],const cs_join_vertex_t set[],cs_join_vertex_t merge_set[],cs_join_vertex_t work_set[])1489 _check_subset_consistency(cs_lnum_t set_size,
1490 const cs_lnum_t subset_num[],
1491 const cs_join_vertex_t set[],
1492 cs_join_vertex_t merge_set[],
1493 cs_join_vertex_t work_set[])
1494 {
1495 cs_lnum_t i, set_id, subset_size;
1496
1497 bool is_consistent = true;
1498
1499 /* Apply merged to each subset */
1500
1501 for (set_id = 0; set_id < set_size; set_id++) {
1502
1503 subset_size = 0;
1504 for (i = 0; i < set_size; i++)
1505 if (subset_num[i] == set_id+1)
1506 work_set[subset_size++] = set[i];
1507
1508 if (subset_size > 0) {
1509
1510 merge_set[set_id] = _compute_merged_vertex(subset_size, work_set);
1511
1512 if (!_is_in_tolerance(subset_size, work_set, merge_set[set_id]))
1513 is_consistent = false;
1514
1515 }
1516
1517 } /* End of loop on subsets */
1518
1519 return is_consistent;
1520 }
1521
1522 /*----------------------------------------------------------------------------
1523 * Get position of the link between vertices i1 and i2.
1524 *
1525 * parameters:
1526 * i1 <-- id in set for vertex 1
1527 * i2 <-- id in set for vertex 2
1528 * idx <-- array of positions
1529 *
1530 * returns:
1531 * position in an array like distances or state
1532 *---------------------------------------------------------------------------*/
1533
1534 inline static cs_lnum_t
_get_pos(cs_lnum_t i1,cs_lnum_t i2,const cs_lnum_t idx[])1535 _get_pos(cs_lnum_t i1,
1536 cs_lnum_t i2,
1537 const cs_lnum_t idx[])
1538 {
1539 cs_lnum_t pos = -1;
1540
1541 if (i1 < i2)
1542 pos = idx[i1] + i2-i1-1;
1543 else {
1544 assert(i1 != i2);
1545 pos = idx[i2] + i1-i2-1;
1546 }
1547
1548 return pos;
1549 }
1550
1551 /*----------------------------------------------------------------------------
1552 * Break equivalences for vertices implied in transitivity issue
1553 *
1554 * parameters:
1555 * param <-- parameters used to manage the joining algorithm
1556 * set_size <-- size of set
1557 * set <-- pointer to a set of vertices
1558 * state <-> array keeping the state of the link
1559 * n_issues <-- number of detected transitivity issues
1560 * issues <-- subset numbers of subset with a transitivity issue
1561 * idx <-- position of vertices couple in array like distances
1562 * subset_num <-- array of subset numbers
1563 * distances <-- array keeping the distances between vertices
1564 *---------------------------------------------------------------------------*/
1565
1566 static void
_break_equivalence(cs_join_param_t param,cs_lnum_t set_size,const cs_join_vertex_t set[],cs_lnum_t state[],cs_lnum_t n_issues,const cs_lnum_t issues[],const cs_lnum_t idx[],const cs_lnum_t subset_num[],const double distances[])1567 _break_equivalence(cs_join_param_t param,
1568 cs_lnum_t set_size,
1569 const cs_join_vertex_t set[],
1570 cs_lnum_t state[],
1571 cs_lnum_t n_issues,
1572 const cs_lnum_t issues[],
1573 const cs_lnum_t idx[],
1574 const cs_lnum_t subset_num[],
1575 const double distances[])
1576 {
1577 cs_lnum_t i, i1, i2, k;
1578
1579 for (i = 0; i < n_issues; i++) {
1580
1581 /* Find the weakest equivalence and break it.
1582 Purpose: Have the minimal number of equivalences to break
1583 for each subset where an inconsistency has been detected */
1584
1585 int i_save = 0;
1586 double rtf = -1.0, dist_save = 0.0;
1587
1588 for (k = 0, i1 = 0; i1 < set_size-1; i1++) {
1589 for (i2 = i1 + 1; i2 < set_size; i2++, k++) {
1590
1591 if (state[k] == 1) { /* v1, v2 are equivalent */
1592
1593 if ( subset_num[i1] == issues[i]
1594 && subset_num[i2] == issues[i]) {
1595
1596 /* Vertices belongs to a subset where an inconsistency
1597 has been found */
1598
1599 double rtf12 = distances[k]/set[i1].tolerance;
1600 double rtf21 = distances[k]/set[i2].tolerance;
1601
1602 assert(rtf12 < 1.0); /* Because they are equivalent */
1603 assert(rtf21 < 1.0);
1604
1605 if (rtf12 >= rtf21) {
1606 if (rtf12 > rtf) {
1607 rtf = rtf12;
1608 i_save = i1;
1609 dist_save = distances[k];
1610 }
1611 }
1612 else {
1613 if (rtf21 > rtf) {
1614 rtf = rtf21;
1615 i_save = i2;
1616 dist_save = distances[k];
1617 }
1618 }
1619
1620 }
1621 }
1622
1623 } /* End of loop on i1 */
1624 } /* End of loop on i2 */
1625
1626 if (rtf > 0.0) {
1627
1628 /* Break equivalence between i_save and all vertices linked to
1629 i_save with a distance to i_save >= dist_save */
1630
1631 for (i2 = 0; i2 < set_size; i2++) {
1632
1633 if (i2 != i_save) {
1634
1635 k = _get_pos(i_save, i2, idx);
1636 if (distances[k] >= dist_save && state[k] == 1) {
1637
1638 state[k] = 0; /* Break equivalence */
1639
1640 if (param.verbosity > 3)
1641 fprintf(cs_glob_join_log,
1642 " %2ld - Break equivalence between [%llu, %llu]"
1643 " (dist_ref: %6.4e)\n",
1644 (long)issues[i],
1645 (unsigned long long)set[i_save].gnum,
1646 (unsigned long long)set[i2].gnum, dist_save);
1647
1648 }
1649 }
1650
1651 } /* End of loop on vertices */
1652
1653 } /* rtf > 0.0 */
1654
1655 } /* End of loop on issues */
1656
1657 }
1658
1659 /*----------------------------------------------------------------------------
1660 * Break equivalences between vertices until each vertex of the list has
1661 * the resulting vertex of the merge under its tolerance.
1662 *
1663 * parameters:
1664 * param <-- set of user-defined parameters
1665 * set_size <-- size of the set of vertices
1666 * set <-> set of vertices
1667 * vbuf <-> tmp buffer
1668 * rbuf <-> tmp buffer
1669 * ibuf <-> tmp buffer
1670 *
1671 * returns:
1672 * number of loops necessary to build consistent subsets
1673 *---------------------------------------------------------------------------*/
1674
1675 static cs_lnum_t
_solve_transitivity(cs_join_param_t param,cs_lnum_t set_size,cs_join_vertex_t set[],cs_join_vertex_t vbuf[],cs_real_t rbuf[],cs_lnum_t ibuf[])1676 _solve_transitivity(cs_join_param_t param,
1677 cs_lnum_t set_size,
1678 cs_join_vertex_t set[],
1679 cs_join_vertex_t vbuf[],
1680 cs_real_t rbuf[],
1681 cs_lnum_t ibuf[])
1682 {
1683 cs_lnum_t i1, i2, k, n_issues;
1684
1685 int n_loops = 0;
1686 bool is_end = false;
1687 cs_lnum_t *subset_num = NULL, *state = NULL, *prev_num = NULL;
1688 cs_lnum_t *subset_issues = NULL, *idx = NULL;
1689 cs_real_t *distances = NULL;
1690 cs_join_vertex_t *merge_set = NULL, *work_set = NULL;
1691
1692 /* Sanity checks */
1693
1694 assert(set_size > 0);
1695
1696 /* Define temporary buffers */
1697
1698 subset_num = &(ibuf[0]);
1699 prev_num = &(ibuf[set_size]);
1700 subset_issues = &(ibuf[2*set_size]);
1701 idx = &(ibuf[3*set_size]);
1702 state = &(ibuf[4*set_size]);
1703 distances = &(rbuf[0]);
1704 merge_set = &(vbuf[0]);
1705 work_set = &(vbuf[set_size]);
1706
1707 /* Compute distances between each couple of vertices among the set */
1708
1709 for (k = 0, i1 = 0; i1 < set_size-1; i1++)
1710 for (i2 = i1 + 1; i2 < set_size; i2++, k++)
1711 distances[k] = _compute_length(set[i1], set[i2]);
1712
1713 /* Compute initial state of equivalences between vertices */
1714
1715 for (k = 0, i1 = 0; i1 < set_size-1; i1++) {
1716 for (i2 = i1 + 1; i2 < set_size; i2++, k++) {
1717 if ( set[i1].tolerance < distances[k]
1718 || set[i2].tolerance < distances[k])
1719 state[k] = 0;
1720 else
1721 state[k] = 1;
1722 }
1723 }
1724
1725 idx[0] = 0;
1726 for (k = 1; k < set_size - 1; k++)
1727 idx[k] = set_size - k + idx[k-1];
1728
1729 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1730 if (cs_glob_join_log != NULL) {
1731 cs_join_dump_array(cs_glob_join_log, "double", "\nDist",
1732 set_size*(set_size-1)/2, distances);
1733 cs_join_dump_array(cs_glob_join_log, "int", "\nInit. State",
1734 set_size*(set_size-1)/2, state);
1735 }
1736 #endif
1737
1738 _build_subsets(set_size, state, prev_num, subset_num);
1739
1740 while (is_end == false && n_loops < param.n_max_equiv_breaks) {
1741
1742 n_loops++;
1743
1744 n_issues = _check_tol_consistency(set_size,
1745 set,
1746 state,
1747 subset_num,
1748 subset_issues,
1749 param.verbosity);
1750
1751 if (n_issues > 0)
1752 _break_equivalence(param,
1753 set_size,
1754 set,
1755 state,
1756 n_issues,
1757 subset_issues,
1758 idx,
1759 subset_num,
1760 distances);
1761
1762 _build_subsets(set_size, state, prev_num, subset_num);
1763
1764 is_end = _check_subset_consistency(set_size,
1765 subset_num,
1766 set,
1767 merge_set,
1768 work_set);
1769
1770 } /* End of while */
1771
1772 if (param.verbosity > 3) {
1773
1774 fprintf(cs_glob_join_log,
1775 " Number of tolerance reductions: %4d\n", n_loops);
1776
1777 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1778 cs_join_dump_array(cs_glob_join_log, "int", "\nFinal Subset",
1779 set_size, subset_num);
1780 #endif
1781
1782 }
1783
1784 /* Apply merged to each subset */
1785
1786 for (k = 0; k < set_size; k++)
1787 set[k] = merge_set[subset_num[k]-1];
1788
1789 return n_loops;
1790 }
1791
1792 /*----------------------------------------------------------------------------
1793 * Merge between vertices. Store new data associated to the merged vertices
1794 * in vertices.
1795 *
1796 * parameters:
1797 * param <-- set of user-defined parameters
1798 * merge_set <-> a pointer to a cs_join_vertex_t structure which
1799 * stores data about merged vertices
1800 * n_vertices <-- number of vertices in vertices array
1801 * vertices <-> array of cs_join_vertex_t structures
1802 *---------------------------------------------------------------------------*/
1803
1804 static void
_merge_vertices(cs_join_param_t param,cs_join_gset_t * merge_set,cs_lnum_t n_vertices,cs_join_vertex_t vertices[])1805 _merge_vertices(cs_join_param_t param,
1806 cs_join_gset_t *merge_set,
1807 cs_lnum_t n_vertices,
1808 cs_join_vertex_t vertices[])
1809 {
1810 cs_lnum_t i, j, k, list_size;
1811 cs_join_vertex_t merged_vertex;
1812 bool ok;
1813
1814 cs_lnum_t max_list_size = 0, vv_max_list_size = 0;
1815 cs_lnum_t n_transitivity = 0;
1816 int n_loops = 0, n_max_loops = 0;
1817
1818 cs_join_gset_t *equiv_gnum = NULL;
1819 cs_real_t *rbuf = NULL;
1820 cs_lnum_t *merge_index = NULL, *ibuf = NULL;
1821 cs_gnum_t *merge_list = NULL, *merge_ref_elts = NULL;
1822 cs_gnum_t *list = NULL;
1823 cs_join_vertex_t *set = NULL, *vbuf = NULL;
1824 FILE *logfile = cs_glob_join_log;
1825
1826 const int verbosity = param.verbosity;
1827
1828 /* Sanity check */
1829
1830 assert(param.merge_tol_coef >= 0.0);
1831
1832 /* Pre-merge of identical vertices */
1833
1834 _pre_merge(param, merge_set, vertices, &equiv_gnum);
1835
1836 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1837 {
1838 int len;
1839 FILE *dbg_file = NULL;
1840 char *filename = NULL;
1841
1842 len = strlen("JoinDBG_MergeSet.dat")+1+2+4;
1843 BFT_MALLOC(filename, len, char);
1844 sprintf(filename, "Join%02dDBG_MergeSet%04d.dat",
1845 param.num, CS_MAX(cs_glob_rank_id, 0));
1846 dbg_file = fopen(filename, "w");
1847
1848 cs_join_gset_dump(dbg_file, merge_set);
1849
1850 fflush(dbg_file);
1851 BFT_FREE(filename);
1852 fclose(dbg_file);
1853 }
1854 #endif /* defined(DEBUG) && !defined(NDEBUG) */
1855
1856 /* Modify the tolerance for the merge operation if needed */
1857
1858 if (fabs(param.merge_tol_coef - 1.0) > 1e-30) {
1859 for (i = 0; i < n_vertices; i++)
1860 vertices[i].tolerance *= param.merge_tol_coef;
1861 }
1862
1863 /* Compute the max. size of a sub list */
1864
1865 merge_index = merge_set->index;
1866 merge_list = merge_set->g_list;
1867 merge_ref_elts = merge_set->g_elts;
1868
1869 for (i = 0; i < merge_set->n_elts; i++) {
1870 list_size = merge_index[i+1] - merge_index[i];
1871 max_list_size = CS_MAX(max_list_size, list_size);
1872 }
1873 vv_max_list_size = ((max_list_size-1)*max_list_size)/2;
1874
1875 if (verbosity > 0) { /* Display information */
1876
1877 cs_lnum_t g_max_list_size = max_list_size;
1878 cs_parall_counter_max(&g_max_list_size, 1);
1879
1880 if (g_max_list_size < 2) {
1881 cs_join_gset_destroy(&equiv_gnum);
1882 bft_printf(_("\n No need to merge vertices.\n"));
1883 return;
1884 }
1885 else
1886 bft_printf(_("\n Max size of a merge set of vertices: %llu\n"),
1887 (unsigned long long)g_max_list_size);
1888 }
1889
1890 /* Temporary buffers allocation */
1891
1892 BFT_MALLOC(ibuf, 4*max_list_size + vv_max_list_size, cs_lnum_t);
1893 BFT_MALLOC(rbuf, vv_max_list_size, cs_real_t);
1894 BFT_MALLOC(vbuf, 2*max_list_size, cs_join_vertex_t);
1895 BFT_MALLOC(list, max_list_size, cs_gnum_t);
1896 BFT_MALLOC(set, max_list_size, cs_join_vertex_t);
1897
1898 /* Merge set of vertices */
1899
1900 for (i = 0; i < merge_set->n_elts; i++) {
1901
1902 list_size = merge_index[i+1] - merge_index[i];
1903
1904 if (list_size > 1) {
1905
1906 for (j = 0, k = merge_index[i]; k < merge_index[i+1]; k++, j++) {
1907 list[j] = merge_list[k];
1908 set[j] = vertices[list[j]];
1909 }
1910
1911 /* Define the resulting cs_join_vertex_t structure of the merge */
1912
1913 merged_vertex = _compute_merged_vertex(list_size, set);
1914
1915 /* Check if the vertex resulting of the merge is in the tolerance
1916 for each vertex of the list */
1917
1918 ok = _is_in_tolerance(list_size, set, merged_vertex);
1919
1920 #if CS_JOIN_MERGE_TOL_REDUC
1921 if (ok == false) { /*
1922 The merged vertex is not in the tolerance of
1923 each vertex. This is a transitivity problem.
1924 We have to split the initial set into several
1925 subsets.
1926 */
1927
1928 n_transitivity++;
1929
1930 /* Display information on vertices to merge */
1931 if (verbosity > 3) {
1932 fprintf(logfile,
1933 "\n Begin merge for ref. elt: %llu - list_size: %ld\n",
1934 (unsigned long long)merge_ref_elts[i],
1935 (long)(merge_index[i+1] - merge_index[i]));
1936 for (j = 0; j < list_size; j++) {
1937 fprintf(logfile, "%9llu -", (unsigned long long)list[j]);
1938 cs_join_mesh_dump_vertex(logfile, set[j]);
1939 }
1940 fprintf(logfile, "\nMerged vertex rejected:\n");
1941 cs_join_mesh_dump_vertex(logfile, merged_vertex);
1942 }
1943
1944 n_loops = _solve_transitivity(param,
1945 list_size,
1946 set,
1947 vbuf,
1948 rbuf,
1949 ibuf);
1950
1951 for (j = 0; j < list_size; j++)
1952 vertices[list[j]] = set[j];
1953
1954 n_max_loops = CS_MAX(n_max_loops, n_loops);
1955
1956 if (verbosity > 3) { /* Display information */
1957 fprintf(logfile, "\n %3d loop(s) to get consistent subsets\n",
1958 n_loops);
1959 fprintf(logfile, "\n End merge for ref. elt: %llu - list_size: %ld\n",
1960 (unsigned long long)merge_ref_elts[i],
1961 (long)(merge_index[i+1] - merge_index[i]));
1962 for (j = 0; j < list_size; j++) {
1963 fprintf(logfile, "%7llu -", (unsigned long long)list[j]);
1964 cs_join_mesh_dump_vertex(logfile, vertices[list[j]]);
1965 }
1966 fprintf(logfile, "\n");
1967 }
1968
1969 }
1970 else /* New vertex data for the sub-elements */
1971
1972 #endif /* CS_JOIN_MERGE_TOL_REDUC */
1973
1974 for (j = 0; j < list_size; j++)
1975 vertices[list[j]] = merged_vertex;
1976
1977 } /* list_size > 1 */
1978
1979 } /* End of loop on potential merges */
1980
1981 /* Apply merge to vertex initially identical */
1982
1983 if (equiv_gnum != NULL) {
1984
1985 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1986 {
1987 int len;
1988 FILE *dbg_file = NULL;
1989 char *filename = NULL;
1990
1991 len = strlen("JoinDBG_EquivMerge.dat")+1+2+4;
1992 BFT_MALLOC(filename, len, char);
1993 sprintf(filename, "Join%02dDBG_EquivMerge%04d.dat",
1994 param.num, CS_MAX(cs_glob_rank_id, 0));
1995 dbg_file = fopen(filename, "w");
1996
1997 cs_join_gset_dump(dbg_file, equiv_gnum);
1998
1999 fflush(dbg_file);
2000 BFT_FREE(filename);
2001 fclose(dbg_file);
2002 }
2003 #endif /* defined(DEBUG) && !defined(NDEBUG) */
2004
2005 for (i = 0; i < equiv_gnum->n_elts; i++) {
2006
2007 cs_lnum_t start = equiv_gnum->index[i];
2008 cs_lnum_t end = equiv_gnum->index[i+1];
2009 cs_lnum_t ref_id = equiv_gnum->g_elts[i];
2010
2011 for (j = start; j < end; j++)
2012 vertices[equiv_gnum->g_list[j]] = vertices[ref_id];
2013
2014 }
2015 }
2016
2017 if (verbosity > 0) {
2018
2019 cs_gnum_t n_g_counter = n_transitivity;
2020 cs_parall_counter(&n_g_counter, 1);
2021
2022 bft_printf(_("\n Excessive transitivity for %llu set(s) of vertices.\n"),
2023 (unsigned long long)n_g_counter);
2024
2025 if (verbosity > 1) {
2026 cs_lnum_t g_n_max_loops = n_max_loops;
2027 cs_parall_counter_max(&g_n_max_loops, 1);
2028 bft_printf(_("\n Max. number of iterations to solve transitivity"
2029 " excess: %llu\n"), (unsigned long long)g_n_max_loops);
2030 }
2031 }
2032
2033 /* Free memory */
2034
2035 BFT_FREE(ibuf);
2036 BFT_FREE(vbuf);
2037 BFT_FREE(rbuf);
2038 BFT_FREE(set);
2039 BFT_FREE(list);
2040
2041 cs_join_gset_destroy(&equiv_gnum);
2042 }
2043
2044 /*----------------------------------------------------------------------------
2045 * Keep an history of the evolution of each vertex id before/after the merge
2046 * operation.
2047 *
2048 * parameters:
2049 * n_iwm_vertices <-- number of vertices before intersection for the
2050 * work cs_join_mesh_t structure
2051 * iwm_vtx_gnum <-- initial global vertex num. (work mesh struct.)
2052 * init_max_vtx_gnum <-- initial max. global numbering for vertices
2053 * n_vertices <-- number of vertices before merge/after intersection
2054 * vertices <-- array of cs_join_vertex_t structures
2055 * p_o2n_vtx_gnum --> distributed array by block on the new global vertex
2056 * numbering for the initial vertices (before inter.)
2057 *---------------------------------------------------------------------------*/
2058
2059 static void
_keep_global_vtx_evolution(cs_lnum_t n_iwm_vertices,const cs_gnum_t iwm_vtx_gnum[],cs_gnum_t init_max_vtx_gnum,cs_lnum_t n_vertices,const cs_join_vertex_t vertices[],cs_gnum_t * p_o2n_vtx_gnum[])2060 _keep_global_vtx_evolution(cs_lnum_t n_iwm_vertices,
2061 const cs_gnum_t iwm_vtx_gnum[],
2062 cs_gnum_t init_max_vtx_gnum,
2063 cs_lnum_t n_vertices,
2064 const cs_join_vertex_t vertices[],
2065 cs_gnum_t *p_o2n_vtx_gnum[])
2066 {
2067 cs_gnum_t *o2n_vtx_gnum = NULL;
2068
2069 const int n_ranks = cs_glob_n_ranks;
2070
2071 assert(n_iwm_vertices <= n_vertices); /* after inter. >= init */
2072
2073 if (n_ranks == 1) {
2074
2075 BFT_MALLOC(o2n_vtx_gnum, n_iwm_vertices, cs_gnum_t);
2076
2077 for (cs_lnum_t i = 0; i < n_iwm_vertices; i++)
2078 o2n_vtx_gnum[i] = vertices[i].gnum;
2079
2080 }
2081
2082 #if defined(HAVE_MPI) /* Parallel treatment */
2083
2084 if (n_ranks > 1) {
2085
2086 cs_lnum_t block_size = 0;
2087
2088 const int local_rank = CS_MAX(cs_glob_rank_id, 0);
2089
2090 cs_block_dist_info_t bi = cs_block_dist_compute_sizes(local_rank,
2091 n_ranks,
2092 1,
2093 0,
2094 init_max_vtx_gnum);
2095
2096 MPI_Comm mpi_comm = cs_glob_mpi_comm;
2097
2098 if (bi.gnum_range[1] > bi.gnum_range[0])
2099 block_size = bi.gnum_range[1] - bi.gnum_range[0];
2100
2101 /* Initialize o2n_vtx_gnum */
2102
2103 BFT_MALLOC(o2n_vtx_gnum, block_size, cs_gnum_t);
2104
2105 for (cs_lnum_t i = 0; i < block_size; i++) {
2106 cs_gnum_t g_id = i;
2107 o2n_vtx_gnum[i] = bi.gnum_range[0] + g_id;
2108 }
2109
2110 /* Send new vtx global number to the related rank = the good block */
2111
2112 cs_all_to_all_t *d
2113 = cs_all_to_all_create_from_block(n_iwm_vertices,
2114 0, /* flags */
2115 iwm_vtx_gnum,
2116 bi,
2117 mpi_comm);
2118
2119 /* Build send_list */
2120
2121 cs_gnum_t *send_glist = NULL;
2122 BFT_MALLOC(send_glist, n_iwm_vertices*2, cs_gnum_t);
2123
2124 for (cs_lnum_t i = 0; i < n_iwm_vertices; i++) {
2125 send_glist[i*2] = iwm_vtx_gnum[i]; /* Old global number */
2126 send_glist[i*2+1] = vertices[i].gnum; /* New global number */
2127 }
2128
2129 cs_gnum_t *recv_glist
2130 = cs_all_to_all_copy_array(d,
2131 CS_GNUM_TYPE,
2132 2,
2133 false, /* reverse, */
2134 send_glist,
2135 NULL);
2136
2137 BFT_FREE(send_glist);
2138
2139 /* Update o2n_vtx_gnum */
2140
2141 const cs_lnum_t n_recv = cs_all_to_all_n_elts_dest(d);
2142
2143 for (cs_lnum_t i = 0; i < n_recv; i++) {
2144
2145 cs_gnum_t o_gnum = recv_glist[i*2];
2146 cs_gnum_t n_gnum = recv_glist[i*2+1];
2147 cs_lnum_t id = o_gnum - bi.gnum_range[0];
2148
2149 #if 0 && defined(DEBUG) && !defined(NDEBUG)
2150 if (o2n_vtx_gnum[id] != bi.gnum_range[0] + id)
2151 assert(o2n_vtx_gnum[id] == n_gnum);
2152 #endif
2153
2154 o2n_vtx_gnum[id] = n_gnum;
2155
2156 }
2157
2158 BFT_FREE(recv_glist);
2159
2160 cs_all_to_all_destroy(&d);
2161
2162 }
2163 #endif /* HAVE_MPI */
2164
2165 /* Set return pointer */
2166
2167 *p_o2n_vtx_gnum = o2n_vtx_gnum;
2168 }
2169
2170 /*----------------------------------------------------------------------------
2171 * Keep a history of the evolution of each vertex id before/after the merge
2172 * operation for the current mesh (local point of view).
2173 *
2174 * parameters:
2175 * n_vertices <-- number of vertices before merge/after intersection
2176 * vertices <-- array of cs_join_vertex_t structures
2177 * p_n_am_vertices --> number of vertices after the merge step
2178 * p_o2n_vtx_id --> array keeping the evolution of the vertex ids
2179 *---------------------------------------------------------------------------*/
2180
2181 static void
_keep_local_vtx_evolution(cs_lnum_t n_vertices,const cs_join_vertex_t vertices[],cs_lnum_t * p_n_am_vertices,cs_lnum_t * p_o2n_vtx_id[])2182 _keep_local_vtx_evolution(cs_lnum_t n_vertices,
2183 const cs_join_vertex_t vertices[],
2184 cs_lnum_t *p_n_am_vertices,
2185 cs_lnum_t *p_o2n_vtx_id[])
2186 {
2187 cs_lnum_t i;
2188 cs_gnum_t prev;
2189
2190 cs_lnum_t n_am_vertices = 0;
2191 cs_lnum_t *o2n_vtx_id = NULL;
2192 cs_lnum_t *order = NULL;
2193 cs_gnum_t *vtx_gnum = NULL;
2194
2195 if (n_vertices == 0)
2196 return;
2197
2198 BFT_MALLOC(vtx_gnum, n_vertices, cs_gnum_t);
2199
2200 for (i = 0; i < n_vertices; i++)
2201 vtx_gnum[i] = vertices[i].gnum;
2202
2203 /* Order vertices according to their global numbering */
2204
2205 BFT_MALLOC(order, n_vertices, cs_lnum_t);
2206
2207 cs_order_gnum_allocated(NULL, vtx_gnum, order, n_vertices);
2208
2209 /* Delete vertices sharing the same global number. Keep only one */
2210
2211 BFT_MALLOC(o2n_vtx_id, n_vertices, cs_lnum_t);
2212
2213 prev = vtx_gnum[order[0]];
2214 o2n_vtx_id[order[0]] = n_am_vertices;
2215
2216 for (i = 1; i < n_vertices; i++) {
2217
2218 cs_lnum_t o_id = order[i];
2219 cs_gnum_t cur = vtx_gnum[o_id];
2220
2221 if (cur != prev) {
2222 prev = cur;
2223 n_am_vertices++;
2224 o2n_vtx_id[o_id] = n_am_vertices;
2225 }
2226 else
2227 o2n_vtx_id[o_id] = n_am_vertices;
2228
2229 } /* End of loop on vertices */
2230
2231 /* n_am_vertices is an id */
2232 n_am_vertices += 1;
2233
2234 assert(n_am_vertices <= n_vertices); /* after merge <= after inter. */
2235
2236 /* Free memory */
2237
2238 BFT_FREE(order);
2239 BFT_FREE(vtx_gnum);
2240
2241 /* Set return pointers */
2242
2243 *p_n_am_vertices = n_am_vertices;
2244 *p_o2n_vtx_id = o2n_vtx_id;
2245 }
2246
2247 /*----------------------------------------------------------------------------
2248 * Search for new elements to add to the definition of the current edge
2249 * These new sub-elements come from initial edges which are now (after the
2250 * merge step) sub-edge of the current edge
2251 * Count step
2252 *
2253 * parameters:
2254 * edge_id <-- id of the edge to deal with
2255 * inter_edges <-- structure keeping data on edge intersections
2256 * edges <-- edges definition
2257 * n_iwm_vertices <-- initial number of vertices in work_mesh
2258 * n_new_sub_elts --> number of new elements to add in the edge definition
2259 *
2260 * returns:
2261 * number of new sub-elements related to this edge
2262 *---------------------------------------------------------------------------*/
2263
2264 static cs_lnum_t
_count_new_sub_edge_elts(cs_lnum_t edge_id,const cs_join_inter_edges_t * inter_edges,const cs_join_edges_t * edges,cs_lnum_t n_iwm_vertices)2265 _count_new_sub_edge_elts(cs_lnum_t edge_id,
2266 const cs_join_inter_edges_t *inter_edges,
2267 const cs_join_edges_t *edges,
2268 cs_lnum_t n_iwm_vertices)
2269 {
2270 cs_lnum_t j, k, j1, j2, sub_edge_id;
2271 cs_lnum_t start, end, _start, _end, v1_num, v2_num;
2272 bool found;
2273
2274 cs_lnum_t n_new_sub_elts = 0;
2275
2276 start = inter_edges->index[edge_id];
2277 end = inter_edges->index[edge_id+1];
2278
2279 for (j1 = start; j1 < end-1; j1++) {
2280
2281 v1_num = inter_edges->vtx_lst[j1];
2282
2283 if (v1_num <= n_iwm_vertices) { /* v1 is an initial vertex */
2284 for (j2 = j1+1; j2 < end; j2++) {
2285
2286 v2_num = inter_edges->vtx_lst[j2];
2287
2288 if (v2_num <= n_iwm_vertices) { /* (v1,v2) is an initial edge */
2289
2290 sub_edge_id = CS_ABS(cs_join_mesh_get_edge(v1_num,
2291 v2_num,
2292 edges))-1;
2293 assert(sub_edge_id != -1);
2294 _start = inter_edges->index[sub_edge_id];
2295 _end = inter_edges->index[sub_edge_id+1];
2296
2297 for (j = _start; j < _end; j++) {
2298
2299 found = false;
2300 for (k = j1+1; k < j2; k++)
2301 if (inter_edges->vtx_glst[k] == inter_edges->vtx_glst[j])
2302 found = true;
2303
2304 if (found == false)
2305 n_new_sub_elts += 1;
2306
2307 } /* End of loop on sub_edge_id definition */
2308
2309 }
2310
2311 } /* End of loop on j2 */
2312 }
2313
2314 } /* End of loop on j1 */
2315
2316 return n_new_sub_elts;
2317 }
2318
2319 /*----------------------------------------------------------------------------
2320 * Update a cs_join_inter_edges_t structure after the merge operation.
2321 * cs_join_inter_edges_t structure should be not NULL.
2322 *
2323 * parameters:
2324 * param <-- user-defined parameters for the joining algorithm
2325 * n_iwm_vertices <-- initial number of vertices in work_mesh
2326 * o2n_vtx_id <-- array keeping the evolution of the vertex ids
2327 * edges <-- edges definition
2328 * p_inter_edges <-> pointer to the structure keeping data on
2329 * edge intersections
2330 *---------------------------------------------------------------------------*/
2331
2332 static void
_update_inter_edges_after_merge(cs_join_param_t param,cs_lnum_t n_iwm_vertices,const cs_lnum_t o2n_vtx_id[],const cs_join_edges_t * edges,const cs_join_mesh_t * mesh,cs_join_inter_edges_t ** p_inter_edges)2333 _update_inter_edges_after_merge(cs_join_param_t param,
2334 cs_lnum_t n_iwm_vertices,
2335 const cs_lnum_t o2n_vtx_id[],
2336 const cs_join_edges_t *edges,
2337 const cs_join_mesh_t *mesh,
2338 cs_join_inter_edges_t **p_inter_edges)
2339 {
2340 cs_lnum_t i, j,k, j1, j2, start_shift, idx_shift;
2341 cs_lnum_t save, _start, _end, start, end;
2342 cs_lnum_t v1_num, v2_num, v1_id, v2_id, sub_edge_id;
2343 cs_gnum_t v1_gnum, v2_gnum, new_gnum, prev_gnum;
2344 bool found;
2345
2346 cs_lnum_t n_adds = 0;
2347
2348 cs_join_inter_edges_t *inter_edges = *p_inter_edges;
2349 cs_join_inter_edges_t *new_inter_edges = NULL;
2350 cs_lnum_t n_edges = inter_edges->n_edges;
2351 cs_lnum_t init_list_size = inter_edges->index[n_edges];
2352 FILE *logfile = cs_glob_join_log;
2353
2354 assert(n_edges == edges->n_edges);
2355
2356 /* Define vtx_glst to compare global vertex numbering */
2357
2358 if (inter_edges->vtx_glst == NULL)
2359 BFT_MALLOC(inter_edges->vtx_glst, inter_edges->index[n_edges], cs_gnum_t);
2360
2361 for (i = 0; i < inter_edges->index[n_edges]; i++) {
2362 v1_id = inter_edges->vtx_lst[i] - 1;
2363 inter_edges->vtx_glst[i] = mesh->vertices[v1_id].gnum;
2364 }
2365
2366 /* Delete redundancies of vertices sharing the same global numbering
2367 after the merge step and define a new index */
2368
2369 idx_shift = 0;
2370 save = inter_edges->index[0];
2371
2372 for (i = 0; i < n_edges; i++) {
2373
2374 start = save;
2375 end = inter_edges->index[i+1];
2376
2377 if (end - start > 0) {
2378
2379 start_shift = start;
2380 v1_id = edges->def[2*i] - 1;
2381 v2_id = edges->def[2*i+1] - 1;
2382 v1_gnum = mesh->vertices[v1_id].gnum;
2383 v2_gnum = mesh->vertices[v2_id].gnum;
2384 prev_gnum = inter_edges->vtx_glst[start_shift];
2385
2386 /* Don't take into account vertices with the same number as the
2387 first edge element */
2388
2389 while (prev_gnum == v1_gnum && start_shift + 1 < end)
2390 prev_gnum = inter_edges->vtx_glst[++start_shift];
2391
2392 if (prev_gnum != v1_gnum && start_shift < end) {
2393
2394 inter_edges->vtx_lst[idx_shift] = inter_edges->vtx_lst[start_shift];
2395 inter_edges->abs_lst[idx_shift] = inter_edges->abs_lst[start_shift];
2396 inter_edges->vtx_glst[idx_shift] = inter_edges->vtx_glst[start_shift];
2397 idx_shift += 1;
2398
2399 for (j = start_shift + 1; j < end; j++) {
2400
2401 new_gnum = inter_edges->vtx_glst[j];
2402
2403 /* Don't take into account redundancies and vertices with the same
2404 number as the second edge element */
2405
2406 if (prev_gnum != new_gnum && new_gnum != v2_gnum) {
2407 prev_gnum = new_gnum;
2408 inter_edges->vtx_lst[idx_shift] = inter_edges->vtx_lst[j];
2409 inter_edges->abs_lst[idx_shift] = inter_edges->abs_lst[j];
2410 inter_edges->vtx_glst[idx_shift] = inter_edges->vtx_glst[j];
2411 idx_shift += 1;
2412 }
2413
2414 }
2415
2416 } /* If start_shift < end */
2417
2418 } /* end - start > 0 */
2419
2420 save = inter_edges->index[i+1];
2421 inter_edges->index[i+1] = idx_shift;
2422
2423 } /* End of loop on edge intersections */
2424
2425 inter_edges->max_sub_size = 0;
2426
2427 for (i = 0; i < n_edges; i++)
2428 inter_edges->max_sub_size =
2429 CS_MAX(inter_edges->max_sub_size,
2430 inter_edges->index[i+1] - inter_edges->index[i]);
2431
2432 assert(inter_edges->index[n_edges] <= init_list_size);
2433
2434 BFT_REALLOC(inter_edges->vtx_lst, inter_edges->index[n_edges], cs_lnum_t);
2435 BFT_REALLOC(inter_edges->abs_lst, inter_edges->index[n_edges], cs_coord_t);
2436
2437 #if 0 && defined(DEBUG) && !defined(NDEBUG) /* Dump local structures */
2438 fprintf(logfile, " AFTER REDUNDANCIES CLEAN\n");
2439 cs_join_inter_edges_dump(logfile, inter_edges, edges, mesh);
2440 #endif
2441
2442 /* Add new vertices from initial edges which are now sub-edges */
2443
2444 for (i = 0; i < n_edges; i++)
2445 n_adds += _count_new_sub_edge_elts(i, inter_edges, edges, n_iwm_vertices);
2446
2447 if (param.verbosity > 2)
2448 fprintf(logfile,
2449 " Number of sub-elements to add to edge definition: %8ld\n",
2450 (long)n_adds);
2451
2452 if (n_adds > 0) { /* Define a new inter_edges structure */
2453
2454 new_inter_edges = cs_join_inter_edges_create(n_edges);
2455
2456 BFT_MALLOC(new_inter_edges->vtx_lst,
2457 inter_edges->index[n_edges] + n_adds, cs_lnum_t);
2458 BFT_MALLOC(new_inter_edges->abs_lst,
2459 inter_edges->index[n_edges] + n_adds, cs_coord_t);
2460
2461 n_adds = 0;
2462 idx_shift = 0;
2463 new_inter_edges->index[0] = 0;
2464
2465 for (i = 0; i < n_edges; i++) {
2466
2467 new_inter_edges->edge_gnum[i] = inter_edges->edge_gnum[i];
2468 start = inter_edges->index[i];
2469 end = inter_edges->index[i+1];
2470
2471 if (start - end > 0) {
2472
2473 for (j1 = start; j1 < end-1; j1++) {
2474
2475 v1_num = inter_edges->vtx_lst[j1];
2476 new_inter_edges->vtx_lst[idx_shift] = v1_num;
2477 new_inter_edges->abs_lst[idx_shift] = inter_edges->abs_lst[j1];
2478 idx_shift++;
2479
2480 if (v1_num <= n_iwm_vertices) { /* v1 is an initial vertex */
2481 for (j2 = j1+1; j2 < end; j2++) {
2482
2483 v2_num = inter_edges->vtx_lst[j2];
2484
2485 if (v2_num <= n_iwm_vertices) { /* (v1,v2) is an initial edge */
2486
2487 sub_edge_id =
2488 CS_ABS(cs_join_mesh_get_edge(v1_num, v2_num, edges))-1;
2489 assert(sub_edge_id != -1);
2490
2491 _start = inter_edges->index[sub_edge_id];
2492 _end = inter_edges->index[sub_edge_id+1];
2493
2494 for (j = _start; j < _end; j++) {
2495
2496 found = false;
2497 for (k = j1+1; k < j2; k++)
2498 if (inter_edges->vtx_glst[k] == inter_edges->vtx_glst[j])
2499 found = true;
2500
2501 if (found == false) {
2502
2503 new_inter_edges->vtx_lst[idx_shift] =
2504 inter_edges->vtx_lst[j];
2505 new_inter_edges->abs_lst[idx_shift] =
2506 inter_edges->abs_lst[j];
2507 idx_shift++;
2508
2509 }
2510
2511 } /* End of loop on sub_edge_id definition */
2512
2513 }
2514
2515 } /* End of loop on j2 */
2516 }
2517
2518 } /* End of loop on j1 */
2519
2520 /* Add last vertex in the previous edge definition */
2521
2522 new_inter_edges->vtx_lst[idx_shift] = inter_edges->vtx_lst[end-1];
2523 new_inter_edges->abs_lst[idx_shift] = inter_edges->abs_lst[end-1];
2524 idx_shift++;
2525
2526 } /* If end - start > 0 */
2527
2528 new_inter_edges->index[i+1] = idx_shift;
2529
2530 } /* End of loop on edges */
2531
2532 cs_join_inter_edges_destroy(&inter_edges);
2533 inter_edges = new_inter_edges;
2534
2535 inter_edges->max_sub_size = 0;
2536
2537 for (i = 0; i < n_edges; i++)
2538 inter_edges->max_sub_size = CS_MAX(inter_edges->max_sub_size,
2539 inter_edges->index[i+1]);
2540
2541 } /* End if n_adds > 0 */
2542
2543 #if 0 && defined(DEBUG) && !defined(NDEBUG) /* Dump local structures */
2544 if (logfile != NULL) {
2545 fprintf(logfile, " AFTER SUB ELTS ADD\n");
2546 cs_join_inter_edges_dump(logfile, inter_edges, edges, mesh);
2547 }
2548 #endif
2549
2550 /* Update cs_join_inter_edges_t structure */
2551
2552 for (i = 0; i < n_edges; i++) {
2553
2554 start = inter_edges->index[i];
2555 end = inter_edges->index[i+1];
2556
2557 for (j = start; j < end; j++) {
2558
2559 cs_lnum_t old_id = inter_edges->vtx_lst[j] - 1;
2560
2561 inter_edges->vtx_lst[j] = o2n_vtx_id[old_id] + 1;
2562
2563 }
2564
2565 }
2566
2567 /* Return pointer */
2568
2569 *p_inter_edges = inter_edges;
2570 }
2571
2572 #if defined(HAVE_MPI)
2573
2574 /*----------------------------------------------------------------------------
2575 * Define send_rank_index and send_faces to prepare the exchange of new faces
2576 * between mesh structures.
2577 *
2578 * parameters:
2579 * n_send <-- number of faces to send
2580 * n_g_faces <-- global number of faces to be joined
2581 * face_gnum <-- global face number
2582 * gnum_rank_index <-- index on ranks for the init. global face numbering
2583 * p_n_send --> number of face/rank couples to send
2584 * p_send_rank --> rank ids to which to send
2585 * p_send_faces --> list of face ids to send
2586 *---------------------------------------------------------------------------*/
2587
2588 static void
_get_faces_to_send(cs_lnum_t n_faces,cs_gnum_t n_g_faces,const cs_gnum_t face_gnum[],const cs_gnum_t gnum_rank_index[],cs_lnum_t * p_n_send,int * p_send_rank[],cs_lnum_t * p_send_faces[])2589 _get_faces_to_send(cs_lnum_t n_faces,
2590 cs_gnum_t n_g_faces,
2591 const cs_gnum_t face_gnum[],
2592 const cs_gnum_t gnum_rank_index[],
2593 cs_lnum_t *p_n_send,
2594 int *p_send_rank[],
2595 cs_lnum_t *p_send_faces[])
2596 {
2597 cs_lnum_t i, rank, reduce_rank;
2598 cs_block_dist_info_t bi;
2599
2600 cs_lnum_t reduce_size = 0, n_send = 0;
2601 int *send_rank = NULL;
2602 cs_lnum_t *send_faces = NULL;
2603 cs_lnum_t *reduce_ids = NULL, *count = NULL;
2604 cs_gnum_t *reduce_index = NULL;
2605
2606 const int local_rank = CS_MAX(cs_glob_rank_id, 0);
2607 const int n_ranks = cs_glob_n_ranks;
2608
2609 /* Sanity checks */
2610
2611 assert(gnum_rank_index != NULL);
2612 assert(n_ranks > 1);
2613
2614 /* Compute block_size */
2615
2616 bi = cs_block_dist_compute_sizes(local_rank,
2617 n_ranks,
2618 1,
2619 0,
2620 n_g_faces);
2621
2622 /* Compact init. global face distribution. Remove ranks without face
2623 at the begining */
2624
2625 for (i = 0; i < n_ranks; i++)
2626 if (gnum_rank_index[i] < gnum_rank_index[i+1])
2627 reduce_size++;
2628
2629 BFT_MALLOC(reduce_index, reduce_size+1, cs_gnum_t);
2630 BFT_MALLOC(reduce_ids, reduce_size, cs_lnum_t);
2631
2632 reduce_size = 0;
2633 reduce_index[0] = gnum_rank_index[0] + 1;
2634
2635 for (i = 0; i < n_ranks; i++) {
2636
2637 /* Add +1 to gnum_rank_index because it's an id and we work on numbers */
2638
2639 if (gnum_rank_index[i] < gnum_rank_index[i+1]) {
2640 reduce_index[reduce_size+1] = gnum_rank_index[i+1] + 1;
2641 reduce_ids[reduce_size++] = i;
2642 }
2643
2644 }
2645
2646 BFT_MALLOC(send_rank, n_faces, int);
2647 BFT_MALLOC(send_faces, n_faces, cs_lnum_t);
2648
2649 /* Fill the list of ranks */
2650
2651 n_send = 0;
2652
2653 for (i = 0; i < n_faces; i++) {
2654
2655 if (face_gnum[i] >= bi.gnum_range[0] && face_gnum[i] < bi.gnum_range[1]) {
2656
2657 /* The current face is a "main" face for the local rank */
2658
2659 reduce_rank = cs_search_gindex_binary(reduce_size,
2660 face_gnum[i],
2661 reduce_index);
2662
2663 assert(reduce_rank > -1);
2664 assert(reduce_rank < reduce_size);
2665
2666 rank = reduce_ids[reduce_rank];
2667 send_rank[n_send] = rank;
2668 send_faces[n_send] = i;
2669 n_send += 1;
2670
2671 } /* End of loop on initial faces */
2672
2673 }
2674
2675 BFT_REALLOC(send_rank, n_send, int);
2676 BFT_REALLOC(send_faces, n_send, cs_lnum_t);
2677
2678 /* Free memory */
2679
2680 #if 0 && defined(DEBUG) && !defined(NDEBUG)
2681 if (cs_glob_join_log != NULL) {
2682 FILE *logfile = cs_glob_join_log;
2683 for (i = 0; i < n_send; i++)
2684 fprintf(logfile, " %d (%llu) to rank %d\n",
2685 send_faces[i], (unsigned long long)face_gnum[send_faces[i]],
2686 send_rank[i]);
2687 fflush(logfile);
2688 }
2689 #endif
2690
2691 BFT_FREE(count);
2692 BFT_FREE(reduce_ids);
2693 BFT_FREE(reduce_index);
2694
2695 /* Set return pointers */
2696
2697 *p_n_send = n_send;
2698 *p_send_rank = send_rank;
2699 *p_send_faces = send_faces;
2700 }
2701
2702 #endif /* defined(HAVE_MPI) */
2703
2704 /*----------------------------------------------------------------------------
2705 * Update local_mesh by redistributing mesh.
2706 * Send back to the original rank the new face and vertex description.
2707 *
2708 * parameters:
2709 * gnum_rank_index <-- index on ranks for the old global face numbering
2710 * send_mesh <-- distributed mesh on faces to join
2711 * p_recv_mesh <-> mesh on local selected faces to be joined
2712 *---------------------------------------------------------------------------*/
2713
2714 static void
_redistribute_mesh(const cs_gnum_t gnum_rank_index[],const cs_join_mesh_t * send_mesh,cs_join_mesh_t ** p_recv_mesh)2715 _redistribute_mesh(const cs_gnum_t gnum_rank_index[],
2716 const cs_join_mesh_t *send_mesh,
2717 cs_join_mesh_t **p_recv_mesh)
2718 {
2719 cs_join_mesh_t *recv_mesh = *p_recv_mesh;
2720
2721 const int n_ranks = cs_glob_n_ranks;
2722
2723 /* sanity checks */
2724
2725 assert(send_mesh != NULL);
2726 assert(recv_mesh != NULL);
2727
2728 if (n_ranks == 1)
2729 cs_join_mesh_copy(&recv_mesh, send_mesh);
2730
2731 #if defined(HAVE_MPI)
2732 if (n_ranks > 1) { /* Parallel mode */
2733
2734 cs_lnum_t n_send = 0;
2735 int *send_rank = NULL;
2736 cs_lnum_t *send_faces = NULL;
2737
2738 MPI_Comm mpi_comm = cs_glob_mpi_comm;
2739
2740 /* Free some structures of the mesh */
2741
2742 cs_join_mesh_reset(recv_mesh);
2743
2744 _get_faces_to_send(send_mesh->n_faces,
2745 send_mesh->n_g_faces,
2746 send_mesh->face_gnum,
2747 gnum_rank_index,
2748 &n_send,
2749 &send_rank,
2750 &send_faces);
2751
2752 assert(n_send <= send_mesh->n_faces);
2753
2754 /* Get the new face connectivity from the distributed send_mesh */
2755
2756 cs_join_mesh_exchange(n_send,
2757 send_rank,
2758 send_faces,
2759 send_mesh,
2760 recv_mesh,
2761 mpi_comm);
2762
2763 BFT_FREE(send_faces);
2764 BFT_FREE(send_rank);
2765
2766 }
2767 #endif
2768
2769 /* Return pointers */
2770
2771 *p_recv_mesh = recv_mesh;
2772
2773 }
2774
2775 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
2776
2777 /*============================================================================
2778 * Public function definitions
2779 *===========================================================================*/
2780
2781 /*----------------------------------------------------------------------------
2782 * Creation of new vertices.
2783 *
2784 * Update list of equivalent vertices, and assign a vertex (existing or
2785 * newly created) to each intersection.
2786 *
2787 * parameters:
2788 * verbosity <-- verbosity level
2789 * edges <-- list of edges
2790 * work <-> joining mesh maintaining initial vertex data
2791 * inter_set <-> cs_join_inter_set_t structure including
2792 * data on edge-edge intersections
2793 * init_max_vtx_gnum <-- initial max. global numbering for vertices
2794 * p_n_g_new_vertices <-> pointer to the global number of new vertices
2795 * p_vtx_eset <-> pointer to a structure dealing with vertex
2796 * equivalences
2797 *---------------------------------------------------------------------------*/
2798
2799 void
cs_join_create_new_vertices(int verbosity,const cs_join_edges_t * edges,cs_join_mesh_t * work,cs_join_inter_set_t * inter_set,cs_gnum_t init_max_vtx_gnum,cs_gnum_t * p_n_g_new_vertices,cs_join_eset_t ** p_vtx_eset)2800 cs_join_create_new_vertices(int verbosity,
2801 const cs_join_edges_t *edges,
2802 cs_join_mesh_t *work,
2803 cs_join_inter_set_t *inter_set,
2804 cs_gnum_t init_max_vtx_gnum,
2805 cs_gnum_t *p_n_g_new_vertices,
2806 cs_join_eset_t **p_vtx_eset)
2807 {
2808 cs_lnum_t i, shift;
2809 double tol_min;
2810 cs_join_vertex_t v1, v2;
2811
2812 cs_lnum_t n_new_vertices = 0;
2813 cs_gnum_t n_g_new_vertices = 0;
2814 cs_gnum_t *new_vtx_gnum = NULL;
2815 cs_lnum_t n_iwm_vertices = work->n_vertices;
2816 cs_join_eset_t *vtx_equiv = *p_vtx_eset;
2817
2818 /* Count the number of new vertices. Update cs_join_inter_set_t struct. */
2819
2820 for (i = 0; i < inter_set->n_inter; i++) {
2821
2822 cs_join_inter_t inter1 = inter_set->inter_lst[2*i];
2823 cs_join_inter_t inter2 = inter_set->inter_lst[2*i+1];
2824
2825 inter1.vtx_id = _get_vtx_id(inter1,
2826 &(edges->def[2*inter1.edge_id]),
2827 n_iwm_vertices,
2828 &n_new_vertices);
2829
2830 inter2.vtx_id = _get_vtx_id(inter2,
2831 &(edges->def[2*inter2.edge_id]),
2832 n_iwm_vertices,
2833 &n_new_vertices);
2834
2835 inter_set->inter_lst[2*i] = inter1;
2836 inter_set->inter_lst[2*i+1] = inter2;
2837
2838 } /* End of loop on intersections */
2839
2840 /* Compute the global numbering for the new vertices (Take into account
2841 potential redundancies) */
2842
2843 _compute_new_vertex_gnum(work,
2844 edges,
2845 inter_set,
2846 init_max_vtx_gnum,
2847 n_iwm_vertices,
2848 n_new_vertices,
2849 &n_g_new_vertices,
2850 &new_vtx_gnum);
2851
2852 if (verbosity > 0)
2853 bft_printf(_("\n Global number of new vertices to create: %10llu\n"),
2854 (unsigned long long)n_g_new_vertices);
2855
2856 /* Define new vertices */
2857
2858 work->n_vertices += n_new_vertices;
2859 work->n_g_vertices += n_g_new_vertices;
2860
2861 BFT_REALLOC(work->vertices, work->n_vertices, cs_join_vertex_t);
2862
2863 #if defined(DEBUG) && !defined(NDEBUG) /* Prepare sanity checks */
2864 {
2865 cs_join_vertex_t incoherency;
2866
2867 /* Initialize to incoherent values new vertices structures */
2868
2869 incoherency.gnum = 0;
2870 incoherency.coord[0] = -9999.9999;
2871 incoherency.coord[1] = -9999.9999;
2872 incoherency.coord[2] = -9999.9999;
2873 incoherency.tolerance = -1.0;
2874 incoherency.state = CS_JOIN_STATE_UNDEF;
2875
2876 for (i = 0; i < n_new_vertices; i++)
2877 work->vertices[n_iwm_vertices + i] = incoherency;
2878
2879 }
2880 #endif
2881
2882 /* Fill vertices structure with new vertex definitions */
2883
2884 for (i = 0; i < inter_set->n_inter; i++) {
2885
2886 cs_join_inter_t inter1 = inter_set->inter_lst[2*i];
2887 cs_join_inter_t inter2 = inter_set->inter_lst[2*i+1];
2888 cs_lnum_t v1_num = inter1.vtx_id + 1;
2889 cs_lnum_t v2_num = inter2.vtx_id + 1;
2890 cs_lnum_t equiv_id = vtx_equiv->n_equiv;
2891
2892 assert(inter1.vtx_id < work->n_vertices);
2893 assert(inter2.vtx_id < work->n_vertices);
2894
2895 /* Create new vertices if needed */
2896
2897 if (v1_num > n_iwm_vertices) {
2898
2899 shift = inter1.vtx_id - n_iwm_vertices;
2900 v1 = _get_new_vertex(inter1.curv_abs,
2901 new_vtx_gnum[shift],
2902 &(edges->def[2*inter1.edge_id]),
2903 work);
2904 tol_min = v1.tolerance;
2905
2906 }
2907 else
2908 tol_min = work->vertices[v1_num-1].tolerance;
2909
2910 if (v2_num > n_iwm_vertices) {
2911
2912 shift = inter2.vtx_id - n_iwm_vertices;
2913 v2 = _get_new_vertex(inter2.curv_abs,
2914 new_vtx_gnum[shift],
2915 &(edges->def[2*inter2.edge_id]),
2916 work);
2917 tol_min = CS_MIN(tol_min, v2.tolerance);
2918
2919 }
2920 else
2921 tol_min = CS_MIN(tol_min, work->vertices[v2_num-1].tolerance);
2922
2923 /* A new vertex has a tolerance equal to the minimal tolerance
2924 between the two vertices implied in the intersection */
2925
2926 if (v1_num > n_iwm_vertices) {
2927 v1.tolerance = tol_min;
2928 work->vertices[inter1.vtx_id] = v1;
2929 }
2930 if (v2_num > n_iwm_vertices) {
2931 v2.tolerance = tol_min;
2932 work->vertices[inter2.vtx_id] = v2;
2933 }
2934
2935 /* Add equivalence between the two current vertices */
2936
2937 cs_join_eset_check_size(equiv_id, &vtx_equiv);
2938
2939 if (v1_num < v2_num) {
2940 vtx_equiv->equiv_couple[2*equiv_id] = v1_num;
2941 vtx_equiv->equiv_couple[2*equiv_id+1] = v2_num;
2942 }
2943 else {
2944 vtx_equiv->equiv_couple[2*equiv_id] = v2_num;
2945 vtx_equiv->equiv_couple[2*equiv_id+1] = v1_num;
2946 }
2947
2948 vtx_equiv->n_equiv += 1;
2949
2950 } /* End of loop on intersections */
2951
2952 /* Free memory */
2953
2954 BFT_FREE(new_vtx_gnum);
2955
2956 #if defined(DEBUG) && !defined(NDEBUG) /* Sanity checks */
2957 for (i = 0; i < work->n_vertices; i++) {
2958
2959 cs_join_vertex_t vtx = work->vertices[i];
2960
2961 if (vtx.gnum == 0 || vtx.tolerance < -0.99)
2962 bft_error(__FILE__, __LINE__, 0,
2963 _(" Inconsistent value found in cs_join_vertex_t struct.:\n"
2964 " Vertex %ld is defined by:\n"
2965 " %llu - [%7.4le, %7.4le, %7.4le] - %lg\n"),
2966 (long)i, (unsigned long long)vtx.gnum,
2967 vtx.coord[0], vtx.coord[1], vtx.coord[2],
2968 vtx.tolerance);
2969
2970 } /* End of loop on vertices */
2971
2972 #if 0
2973 _dump_vtx_eset(vtx_equiv, work, cs_glob_join_log);
2974 #endif
2975
2976 #endif
2977
2978 /* Set return pointers */
2979
2980 *p_n_g_new_vertices = n_g_new_vertices;
2981 *p_vtx_eset = vtx_equiv;
2982 }
2983
2984 /*----------------------------------------------------------------------------
2985 * Merge of equivalent vertices (and tolerance reduction if necessary)
2986 *
2987 * Define a new cs_join_vertex_t structure (stored in "work" structure).
2988 * Returns an updated cs_join_mesh_t and cs_join_edges_t structures.
2989 *
2990 * parameters:
2991 * param <-- set of user-defined parameters for the joining
2992 * n_g_vertices_tot <-- global number of vertices (initial parent mesh)
2993 * work <-> pointer to a cs_join_mesh_t structure
2994 * vtx_eset <-- structure storing equivalences between vertices
2995 * (two vertices are equivalent if they are within
2996 * each other's tolerance)
2997 *---------------------------------------------------------------------------*/
2998
2999 void
cs_join_merge_vertices(cs_join_param_t param,cs_gnum_t n_g_vertices_tot,cs_join_mesh_t * work,const cs_join_eset_t * vtx_eset)3000 cs_join_merge_vertices(cs_join_param_t param,
3001 cs_gnum_t n_g_vertices_tot,
3002 cs_join_mesh_t *work,
3003 const cs_join_eset_t *vtx_eset)
3004 {
3005 cs_gnum_t *vtx_tags = NULL;
3006 cs_join_gset_t *merge_set = NULL;
3007
3008 const int n_ranks = cs_glob_n_ranks;
3009
3010 /* Initialize counters for the merge operation */
3011
3012 _initialize_merge_counter();
3013
3014 #if 0 && defined(DEBUG) && !defined(NDEBUG) /* Dump local structures */
3015 _dump_vtx_eset(vtx_eset, work, cs_glob_join_log);
3016 #endif
3017
3018 if (param.verbosity > 2) {
3019 cs_gnum_t g_n_equiv = vtx_eset->n_equiv;
3020 cs_parall_counter(&g_n_equiv, 1);
3021 fprintf(cs_glob_join_log,
3022 "\n"
3023 " Final number of equiv. between vertices; local: %9ld\n"
3024 " global: %9llu\n",
3025 (long)vtx_eset->n_equiv, (unsigned long long)g_n_equiv);
3026 }
3027
3028 /* Operate merge between equivalent vertices.
3029 Manage reduction of tolerance if necessary */
3030
3031 /* Tag with the same number all the vertices which might be merged together */
3032
3033 _tag_equiv_vertices(n_g_vertices_tot,
3034 vtx_eset,
3035 work,
3036 param.verbosity,
3037 &vtx_tags);
3038
3039 if (n_ranks == 1) { /* Serial mode */
3040
3041 /* Build a merge list */
3042
3043 merge_set = cs_join_gset_create_from_tag(work->n_vertices, vtx_tags);
3044
3045 /* Merge of equivalent vertices */
3046
3047 _merge_vertices(param,
3048 merge_set,
3049 work->n_vertices,
3050 work->vertices);
3051
3052 }
3053
3054 #if defined(HAVE_MPI)
3055 if (n_ranks > 1) { /* Parallel mode: we work by block */
3056
3057 MPI_Comm mpi_comm = cs_glob_mpi_comm;
3058
3059 const cs_lnum_t n_vertices = work->n_vertices;
3060 const cs_gnum_t _n_ranks = n_ranks;
3061
3062 int *dest_rank = NULL;
3063 BFT_MALLOC(dest_rank, n_vertices, int);
3064
3065 for (cs_lnum_t i = 0; i < n_vertices; i++)
3066 dest_rank[i] = (vtx_tags[i] - 1) % _n_ranks;
3067
3068 cs_all_to_all_t *d
3069 = cs_all_to_all_create(n_vertices,
3070 0, /* flags */
3071 NULL,
3072 dest_rank,
3073 mpi_comm);
3074
3075 cs_all_to_all_transfer_dest_rank(d, &dest_rank);
3076
3077 /* Build a merge list in parallel */
3078
3079 cs_join_vertex_t *vtx_merge_data = NULL;
3080
3081 _build_parall_merge_structures(work,
3082 vtx_tags,
3083 d,
3084 &vtx_merge_data,
3085 &merge_set);
3086
3087 /* Merge of equivalent vertices for the current block */
3088
3089 const cs_lnum_t n_recv = cs_all_to_all_n_elts_dest(d);
3090
3091 _merge_vertices(param,
3092 merge_set,
3093 n_recv,
3094 vtx_merge_data);
3095
3096 /* Allocate send_vtx_data and exchange vtx_merge_data */
3097
3098 cs_all_to_all_copy_array(d,
3099 CS_CHAR,
3100 sizeof(cs_join_vertex_t),
3101 true, /* reverse */
3102 vtx_merge_data,
3103 work->vertices);
3104
3105 BFT_FREE(vtx_merge_data);
3106
3107 cs_all_to_all_destroy(&d);
3108
3109 }
3110 #endif /* HAVE_MPI */
3111
3112 /* Free memory */
3113
3114 BFT_FREE(vtx_tags);
3115
3116 cs_join_gset_destroy(&merge_set);
3117
3118 if (param.verbosity > 1)
3119 bft_printf(_("\n"
3120 " Merging of equivalent vertices done.\n"));
3121 }
3122
3123 /*----------------------------------------------------------------------------
3124 * Merge of equivalent vertices (and reduction of tolerance if necessary)
3125 *
3126 * Define a new cs_join_vertex_t structure (stored in "work" structure)
3127 * Returns an updated cs_join_mesh_t and cs_join_edges_t structures.
3128 *
3129 * parameters:
3130 * param <-- set of user-defined parameters for the joining
3131 * n_iwm_vertices <-- initial number of vertices (work mesh struct.)
3132 * iwm_vtx_gnum <-- initial global vertex num. (work mesh struct)
3133 * init_max_vtx_gnum <-- initial max. global numbering for vertices
3134 * rank_face_gnum_index <-- index on face global numbering to determine
3135 * the related rank
3136 * p_mesh <-> pointer to cs_join_mesh_t structure
3137 * p_edges <-> pointer to cs_join_edges_t structure
3138 * p_inter_edges <-> pointer to a cs_join_inter_edges_t struct.
3139 * p_local_mesh <-> pointer to a cs_join_mesh_t structure
3140 * p_o2n_vtx_gnum --> array on blocks on the new global vertex
3141 * numbering for the init. vertices (before inter.)
3142 *---------------------------------------------------------------------------*/
3143
3144 void
cs_join_merge_update_struct(cs_join_param_t param,cs_lnum_t n_iwm_vertices,const cs_gnum_t iwm_vtx_gnum[],cs_gnum_t init_max_vtx_gnum,const cs_gnum_t rank_face_gnum_index[],cs_join_mesh_t ** p_mesh,cs_join_edges_t ** p_edges,cs_join_inter_edges_t ** p_inter_edges,cs_join_mesh_t ** p_local_mesh,cs_gnum_t * p_o2n_vtx_gnum[])3145 cs_join_merge_update_struct(cs_join_param_t param,
3146 cs_lnum_t n_iwm_vertices,
3147 const cs_gnum_t iwm_vtx_gnum[],
3148 cs_gnum_t init_max_vtx_gnum,
3149 const cs_gnum_t rank_face_gnum_index[],
3150 cs_join_mesh_t **p_mesh,
3151 cs_join_edges_t **p_edges,
3152 cs_join_inter_edges_t **p_inter_edges,
3153 cs_join_mesh_t **p_local_mesh,
3154 cs_gnum_t *p_o2n_vtx_gnum[])
3155 {
3156 cs_lnum_t n_am_vertices = 0; /* new number of vertices after merge */
3157 cs_lnum_t *o2n_vtx_id = NULL;
3158 cs_gnum_t *o2n_vtx_gnum = NULL;
3159 cs_join_mesh_t *mesh = *p_mesh;
3160 cs_join_mesh_t *local_mesh = *p_local_mesh;
3161 cs_join_edges_t *edges = *p_edges;
3162 cs_join_inter_edges_t *inter_edges = *p_inter_edges;
3163
3164 /* Keep an history of the evolution of each vertex */
3165
3166 _keep_global_vtx_evolution(n_iwm_vertices, /* n_vertices before inter */
3167 iwm_vtx_gnum,
3168 init_max_vtx_gnum,
3169 mesh->n_vertices, /* n_vertices after inter */
3170 mesh->vertices,
3171 &o2n_vtx_gnum); /* defined by block in // */
3172
3173 _keep_local_vtx_evolution(mesh->n_vertices, /* n_vertices after inter */
3174 mesh->vertices,
3175 &n_am_vertices, /* n_vertices after merge */
3176 &o2n_vtx_id);
3177
3178 /* Update all structures which keeps data about vertices */
3179
3180 if (inter_edges != NULL) { /* The join type is not conform */
3181
3182 /* Update inter_edges structure */
3183
3184 _update_inter_edges_after_merge(param,
3185 n_iwm_vertices,
3186 o2n_vtx_id, /* size of mesh->n_vertices */
3187 edges,
3188 mesh,
3189 &inter_edges);
3190
3191 assert(edges->n_edges == inter_edges->n_edges); /* Else: problem for
3192 future synchro. */
3193
3194 /* Update cs_join_mesh_t structure after the merge of vertices
3195 numbering of the old vertices + add new vertices */
3196
3197 cs_join_mesh_update(mesh,
3198 edges,
3199 inter_edges->index,
3200 inter_edges->vtx_lst,
3201 n_am_vertices,
3202 o2n_vtx_id);
3203
3204 } /* End if inter_edges != NULL */
3205
3206 else
3207 /* Update cs_join_mesh_t structure after the merge of vertices
3208 numbering of the old vertices + add new vertices */
3209
3210 cs_join_mesh_update(mesh,
3211 edges,
3212 NULL,
3213 NULL,
3214 n_am_vertices,
3215 o2n_vtx_id);
3216
3217 BFT_FREE(o2n_vtx_id);
3218
3219 /* Update local_mesh by redistributing mesh */
3220
3221 _redistribute_mesh(rank_face_gnum_index,
3222 mesh,
3223 &local_mesh);
3224
3225 /* Set return pointers */
3226
3227 *p_mesh = mesh;
3228 *p_edges = edges;
3229 *p_inter_edges = inter_edges;
3230 *p_o2n_vtx_gnum = o2n_vtx_gnum;
3231 *p_local_mesh = local_mesh;
3232 }
3233
3234 /*---------------------------------------------------------------------------*/
3235
3236 END_C_DECLS
3237