1 /*============================================================================
2 * Management of conforming and non-conforming joining
3 *===========================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *---------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <errno.h>
35 #include <assert.h>
36 #include <math.h>
37 #include <stdio.h>
38 #include <stdlib.h>
39 #include <string.h>
40
41 /*----------------------------------------------------------------------------
42 * Local headers
43 *---------------------------------------------------------------------------*/
44
45 #include "bft_mem.h"
46 #include "bft_error.h"
47 #include "bft_printf.h"
48
49 #include "fvm_periodicity.h"
50
51 #include "cs_all_to_all.h"
52 #include "cs_block_dist.h"
53 #include "cs_join_intersect.h"
54 #include "cs_join_merge.h"
55 #include "cs_join_mesh.h"
56 #include "cs_join_post.h"
57 #include "cs_join_perio.h"
58 #include "cs_join_set.h"
59 #include "cs_join_split.h"
60 #include "cs_join_update.h"
61 #include "cs_join_util.h"
62 #include "cs_log.h"
63 #include "cs_mesh_quantities.h"
64 #include "cs_parall.h"
65 #include "cs_post.h"
66 #include "cs_timer.h"
67
68 /*----------------------------------------------------------------------------
69 * Header for the current file
70 *---------------------------------------------------------------------------*/
71
72 #include "cs_join.h"
73
74 /*---------------------------------------------------------------------------*/
75
76 BEGIN_C_DECLS
77
78 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
79
80 /*============================================================================
81 * Private function definitions
82 *===========================================================================*/
83
84 #if 0 && defined(DEBUG) && !defined(NDEBUG)
85
86 /*----------------------------------------------------------------------------
87 * Dump a cs_mesh_t structure into a file for debugging purpose
88 *
89 * parameters:
90 * join_num <-- join number
91 * basename <-- string
92 * mesh <-- pointer to cs_mesh_t structure
93 *---------------------------------------------------------------------------*/
94
95 static void
96 _dump_mesh(const int join_num,
97 const char basename[],
98 const cs_mesh_t *mesh)
99 {
100 int base_len, len;
101 FILE *dbg_file = NULL;
102 char *filename = NULL;
103
104 base_len = strlen(basename);
105 len = strlen("log/JoinDBG__.dat")+1+4+2+base_len;
106 BFT_MALLOC(filename, len, char);
107 sprintf(filename, "log%cJoin%02dDBG_%s_%04d.dat", CS_DIR_SEPARATOR,
108 join_num, basename, cs_glob_rank_id);
109 dbg_file = fopen(filename, "w");
110
111 cs_mesh_dump_file(dbg_file, mesh);
112
113 fflush(dbg_file);
114 BFT_FREE(filename);
115 fclose(dbg_file);
116 }
117
118 /*----------------------------------------------------------------------------
119 * Dump a cs_join_gset_t structure into a file for debugging purpose
120 *
121 * parameters:
122 * join_num <-- join number
123 * basename <-- string
124 * set <-- pointer to cs_join_gset_t structure
125 *---------------------------------------------------------------------------*/
126
127 static void
128 _dump_gset(const int join_num,
129 const char basename[],
130 const cs_join_gset_t *set)
131 {
132 int base_len, len;
133 FILE *dbg_file = NULL;
134 char *filename = NULL;
135
136 base_len = strlen(basename);
137 len = strlen("log/JoinDBG__.dat")+1+4+2+base_len;
138 BFT_MALLOC(filename, len, char);
139 sprintf(filename, "log%cJoin%02dDBG_%s_%04d.dat", CS_DIR_SEPARATOR,
140 join_num, basename, cs_glob_rank_id);
141 dbg_file = fopen(filename, "w");
142
143 cs_join_gset_dump(dbg_file, set);
144
145 fflush(dbg_file);
146 BFT_FREE(filename);
147 fclose(dbg_file);
148 }
149
150 #endif
151
152 /*----------------------------------------------------------------------------
153 * Build a structure keeping data about entities selection and modify mesh
154 * in case of periodicity.
155 *
156 * parameters:
157 * this_join <-- pointer to a cs_join_t structure
158 * mesh <-> pointer to cs_mesh_t structure
159 *---------------------------------------------------------------------------*/
160
161 static void
_select_entities(cs_join_t * this_join,cs_mesh_t * mesh)162 _select_entities(cs_join_t *this_join,
163 cs_mesh_t *mesh)
164 {
165 cs_real_t *b_face_cog = NULL, *b_face_normal = NULL;
166 cs_join_param_t param = this_join->param;
167
168 const char *selection_criteria = this_join->criteria;
169
170 cs_mesh_init_group_classes(mesh);
171
172 cs_mesh_quantities_b_faces(mesh, &b_face_cog, &b_face_normal);
173
174 cs_glob_mesh->select_b_faces = fvm_selector_create(mesh->dim,
175 mesh->n_b_faces,
176 mesh->class_defs,
177 mesh->b_face_family,
178 1,
179 b_face_cog,
180 b_face_normal);
181
182 /* Get selected faces for this joining and define the related
183 cs_join_face_select_t structure.
184 - Compute the global number of selected faces
185 - Get the adjacent faces, ... */
186
187 this_join->selection = cs_join_select_create(selection_criteria,
188 param.perio_type,
189 param.verbosity);
190
191 /* Free arrays and structures needed for selection */
192
193 BFT_FREE(b_face_cog);
194 BFT_FREE(b_face_normal);
195
196 mesh->class_defs = fvm_group_class_set_destroy(mesh->class_defs);
197
198 if (mesh->select_b_faces != NULL)
199 mesh->select_b_faces = fvm_selector_destroy(mesh->select_b_faces);
200 if (mesh->class_defs != NULL)
201 mesh->class_defs = fvm_group_class_set_destroy(mesh->class_defs);
202
203 /* Return selection struct. */
204
205 if (mesh->verbosity > 0) {
206 bft_printf(_("\n Element selection successfully done.\n"));
207 bft_printf_flush();
208 }
209 }
210
211 /*----------------------------------------------------------------------------
212 * Define a cs_join_mesh_t structure on only faces which will be
213 * potentially modified by the joining operation.
214 *
215 * In serial mode, this is a subset of the local join mesh.
216 * In parallel mode, this is a distributed subset of the global join mesh.
217 *
218 * Subset is a restriction on faces which could intersect each other.
219 *
220 * Distribution is made so that ther is a well-balanced number of faces
221 * on each rank and so that faces in the mesh are spatially coherent
222 * to insure no problem for differents joining operations.
223 *
224 * Get the associated edges and the list of possible intersections
225 * between these edges.
226 *
227 * parameters:
228 * param <-- set of user-defined parameter
229 * rank_face_gnum_index <-- index on face global numering to determine
230 * the related rank
231 * local_mesh <-- pointer to a cs_join_mesh_t structure
232 * p_work_mesh --> pointer to the work cs_join_mesh_t structure
233 * p_work_edges --> pointer to the cs_join_edges_t structure
234 * p_work_face_normal --> pointer to the normal of faces defined in
235 * the work mesh struture.
236 * p_edge_edge_vis --> pointer to a cs_join_gset_t structure storing
237 * the visibility between edges
238 * stats <-> joining statistics
239 *---------------------------------------------------------------------------*/
240
241 static void
_get_work_struct(cs_join_param_t param,const cs_gnum_t rank_face_gnum_index[],const cs_join_mesh_t * local_mesh,cs_join_mesh_t ** p_work_mesh,cs_join_edges_t ** p_work_edges,cs_real_t * p_work_face_normal[],cs_join_gset_t ** p_edge_edge_vis,cs_join_stats_t * stats)242 _get_work_struct(cs_join_param_t param,
243 const cs_gnum_t rank_face_gnum_index[],
244 const cs_join_mesh_t *local_mesh,
245 cs_join_mesh_t **p_work_mesh,
246 cs_join_edges_t **p_work_edges,
247 cs_real_t *p_work_face_normal[],
248 cs_join_gset_t **p_edge_edge_vis,
249 cs_join_stats_t *stats)
250 {
251 cs_lnum_t n_inter_faces = 0;
252 char *mesh_name = NULL;
253 cs_real_t *face_normal = NULL;
254 cs_gnum_t *intersect_face_gnum = NULL;
255 cs_join_gset_t *face_face_vis = NULL, *edge_edge_vis = NULL;
256 cs_join_mesh_t *work_mesh = NULL;
257 cs_join_edges_t *work_edges = NULL;
258
259 const int n_ranks = cs_glob_n_ranks;
260 const int local_rank = CS_MAX(cs_glob_rank_id, 0);
261
262 /*
263 Build a bounding box for each selected face.
264 Find intersections between bounding boxes for the whole selected mesh
265 and retrieve a list (cs_join_gset_t structure) which has been
266 distributed over the ranks containing this information.
267 */
268
269 face_face_vis = cs_join_intersect_faces(param, local_mesh, stats);
270
271 /* Define an ordered list of all implied faces without redundancy */
272
273 /* TODO: check if this is necessary after cleanup done in face_face_vis */
274
275 cs_timer_t t0 = cs_timer_time();
276
277 cs_join_gset_single_order(face_face_vis,
278 &n_inter_faces,
279 &intersect_face_gnum);
280
281 cs_timer_t t1 = cs_timer_time();
282
283 if (param.verbosity > 1)
284 bft_printf(_("\n Sorted possible intersections between faces.\n"));
285
286 cs_timer_counter_add_diff(&(stats->t_inter_sort), &t0, &t1);
287
288 /* Define a distributed cs_join_mesh_t structure to store the connectivity
289 of the intersecting faces associated to their bounding boxes in
290 face_inter list */
291
292 if (n_ranks > 1) {
293
294 BFT_MALLOC(mesh_name, strlen("WorkMesh_j_n") + 2 + 5 + 1, char);
295 sprintf(mesh_name,"%s%02d%s%05d",
296 "WorkMesh_j", param.num, "_n", local_rank);
297
298 }
299 else {
300
301 BFT_MALLOC(mesh_name, strlen("WorkMesh_j") + 2 + 1, char);
302 sprintf(mesh_name,"%s%02d", "WorkMesh_j", param.num);
303
304 }
305
306 work_mesh = cs_join_mesh_create_from_glob_sel(mesh_name,
307 n_inter_faces,
308 intersect_face_gnum,
309 rank_face_gnum_index,
310 local_mesh);
311
312 /* Define a cs_join_edges_t structure associated to a cs_join_mesh_t
313 structure on which we work */
314
315 work_edges = cs_join_mesh_define_edges(work_mesh);
316
317 /* Transform face_inter into edge_inter */
318
319 edge_edge_vis = cs_join_intersect_face_to_edge(work_mesh,
320 work_edges,
321 face_face_vis);
322
323 /* Define the normal vector for each selected face before any modification */
324
325 face_normal = cs_join_mesh_get_face_normal(work_mesh);
326
327 /* Free memory */
328
329 BFT_FREE(mesh_name);
330 BFT_FREE(intersect_face_gnum);
331
332 cs_join_gset_destroy(&face_face_vis);
333
334 /* Return pointers */
335
336 *p_work_mesh = work_mesh;
337 *p_work_edges = work_edges;
338 *p_edge_edge_vis = edge_edge_vis;
339 *p_work_face_normal = face_normal;
340 }
341
342 /*----------------------------------------------------------------------------
343 * Build several structures useful to join faces.
344 *
345 * parameters:
346 * this_join <-- pointer to a cs_join_t structure
347 * mesh <-> pointer of pointer to cs_mesh_t structure
348 * p_loc_jmesh --> local cs_join_mesh_t structure based on local face
349 * selection
350 * p_work_jmesh --> distributed and balanced cs_join_mesh_t structure
351 * based on the global face selection
352 * p_work_join_edges --> edges definition related to work_jmesh
353 * p_work_face_normal --> unitary normal for the faces of work_jmesh
354 * p_edge_edge_vis --> list of all potential intersections between edges
355 *---------------------------------------------------------------------------*/
356
357 static void
_build_join_structures(cs_join_t * this_join,cs_mesh_t * mesh,cs_join_mesh_t ** p_loc_jmesh,cs_join_mesh_t ** p_work_jmesh,cs_join_edges_t ** p_work_join_edges,cs_real_t * p_work_face_normal[],cs_join_gset_t ** p_edge_edge_vis)358 _build_join_structures(cs_join_t *this_join,
359 cs_mesh_t *mesh,
360 cs_join_mesh_t **p_loc_jmesh,
361 cs_join_mesh_t **p_work_jmesh,
362 cs_join_edges_t **p_work_join_edges,
363 cs_real_t *p_work_face_normal[],
364 cs_join_gset_t **p_edge_edge_vis)
365 {
366 char *mesh_name = NULL;
367 cs_real_t *work_face_normal = NULL;
368 cs_join_gset_t *edge_edge_vis = NULL;
369 cs_join_mesh_t *loc_jmesh = NULL, *work_jmesh = NULL;
370 cs_join_edges_t *work_edges = NULL;
371 cs_join_param_t param = this_join->param;
372 cs_join_select_t *selection = this_join->selection;
373
374 cs_timer_t t0 = cs_timer_time();
375
376 /* Define a cs_join_mesh_structure from the selected connectivity */
377
378 if (cs_glob_n_ranks > 1) {
379 BFT_MALLOC(mesh_name, strlen("LocalMesh_n") + 5 + 1, char);
380 sprintf(mesh_name,"%s%05d", "LocalMesh_n", CS_MAX(cs_glob_rank_id, 0));
381 }
382 else {
383 BFT_MALLOC(mesh_name, strlen("LocalMesh") + 1, char);
384 sprintf(mesh_name,"%s", "LocalMesh");
385 }
386
387 loc_jmesh = cs_join_mesh_create_from_select(mesh_name,
388 param,
389 selection,
390 mesh->b_face_vtx_idx,
391 mesh->b_face_vtx_lst,
392 mesh->i_face_vtx_idx,
393 mesh->i_face_vtx_lst,
394 mesh->n_vertices,
395 mesh->vtx_coord,
396 mesh->global_vtx_num);
397
398 BFT_FREE(mesh_name);
399
400 if (param.perio_type != FVM_PERIODICITY_NULL)
401 cs_join_perio_apply(this_join, loc_jmesh, mesh);
402
403 if (param.verbosity > 0)
404 cs_join_mesh_minmax_tol(param, loc_jmesh);
405
406 cs_timer_t t1 = cs_timer_time();
407
408 if (param.visualization > 2)
409 cs_join_post_dump_mesh("LocalMesh", loc_jmesh, param);
410
411 /*
412 Define a cs_join_mesh_t structure on only faces which will be
413 potentially modified by the joining operation.
414
415 In serial mode, this is a subset of the local join mesh.
416 In parallel mode, this is a distributed subset of the global join mesh.
417
418 Subset is a restriction on faces which could be intersected each other.
419
420 Distribution is made so that there is a well-balanced number of faces
421 on each rank and so that faces in the mesh are spatially coherent
422 to insure no problem during the different joining operations.
423
424 Get the associated edges and the list of potential intersections
425 between these edges through an edge-edge visibility.
426 */
427
428 _get_work_struct(param,
429 selection->compact_rank_index,
430 loc_jmesh,
431 &work_jmesh,
432 &work_edges,
433 &work_face_normal,
434 &edge_edge_vis,
435 &(this_join->stats));
436
437 /* log performance info of previous step here only to simplify
438 "pretty printing". */
439
440 cs_timer_counter_add_diff(&(this_join->stats.t_l_join_mesh), &t0, &t1);
441
442 /* Return pointers */
443
444 *p_loc_jmesh = loc_jmesh;
445 *p_work_jmesh = work_jmesh;
446 *p_work_join_edges = work_edges;
447 *p_edge_edge_vis = edge_edge_vis;
448 *p_work_face_normal = work_face_normal;
449 }
450
451 /*----------------------------------------------------------------------------
452 * From real intersection between edges, define new vertices and/or
453 * update old vertices.
454 * Keep the relation between two intersecting edges through an equivalence
455 * between the vertex of each edge.
456 * Store also the new description of initial edges through a
457 * cs_join_inter_edges_t structure and synchronize it to get all the
458 * possible equivalences between vertices.
459 *
460 * parameters:
461 * this_join <-- pointer to a cs_join_t structure
462 * work_jmesh <-> pointer to a cs_join_mesh_t structure
463 * work_join_edges <-- pointer to a cs_join_edges_t structure
464 * p_edge_edge_vis <-> pointer to a cs_join_glist_t structure
465 * (freed here)
466 * init_max_vtx_gnum <-- initial max. global numbering for vertices
467 * p_n_g_new_vertices --> global number of vertices created during the
468 * intersection of edges
469 * p_vtx_eset --> structure storing equivalences between vertices
470 * Two vertices are equivalent if they are each
471 * other in their tolerance
472 * p_inter_edges --> structure storing the definition of new vertices
473 * on initial edges
474 *
475 * returns:
476 * type of joining detected
477 *---------------------------------------------------------------------------*/
478
479 static cs_join_type_t
_intersect_edges(cs_join_t * this_join,cs_join_mesh_t * work_jmesh,const cs_join_edges_t * work_join_edges,cs_join_gset_t ** p_edge_edge_vis,cs_gnum_t init_max_vtx_gnum,cs_gnum_t * p_n_g_new_vertices,cs_join_eset_t ** p_vtx_eset,cs_join_inter_edges_t ** p_inter_edges)480 _intersect_edges(cs_join_t *this_join,
481 cs_join_mesh_t *work_jmesh,
482 const cs_join_edges_t *work_join_edges,
483 cs_join_gset_t **p_edge_edge_vis,
484 cs_gnum_t init_max_vtx_gnum,
485 cs_gnum_t *p_n_g_new_vertices,
486 cs_join_eset_t **p_vtx_eset,
487 cs_join_inter_edges_t **p_inter_edges)
488 {
489 cs_join_type_t join_type = CS_JOIN_TYPE_NULL;
490
491 cs_gnum_t n_g_new_vertices = 0;
492 cs_join_inter_edges_t *inter_edges = NULL;
493 cs_join_eset_t *vtx_eset = NULL;
494 cs_join_inter_set_t *inter_set = NULL;
495 cs_join_param_t param = this_join->param;
496
497 const int n_ranks = cs_glob_n_ranks;
498
499 cs_timer_t t0 = cs_timer_time();
500
501 /*
502 Compute the intersections between edges.
503 Store the output in two data structures:
504 - a cs_join_eset_t struct. to store equiv. between vertices
505 issued from the same intersection
506 - a cs_join_inter_set_t struct. to store detected intersections
507 Return the type of the joining operation: conform or not.
508 */
509
510 join_type = cs_join_intersect_edges(param,
511 *p_edge_edge_vis,
512 work_join_edges,
513 work_jmesh,
514 &vtx_eset,
515 &inter_set);
516
517 cs_timer_t t1 = cs_timer_time();
518 cs_timer_counter_add_diff(&(this_join->stats.t_edge_inter), &t0, &t1);
519 t0 = t1;
520
521 cs_join_gset_destroy(p_edge_edge_vis);
522
523 #if 0 && defined(DEBUG) && !defined(NDEBUG) /* Dump structures after inter. */
524 cs_join_inter_set_dump(inter_set, work_join_edges, work_jmesh);
525 #endif
526
527 if (join_type == CS_JOIN_TYPE_NULL) {
528
529 if (cs_glob_mesh->verbosity > 0) {
530 bft_printf(_("\n Joining operation did modify any vertices.\n"));
531 bft_printf_flush();
532 }
533
534 cs_join_inter_set_destroy(&inter_set);
535
536 }
537 else if (join_type == CS_JOIN_TYPE_CONFORMING) {
538
539 if (cs_glob_mesh->verbosity > 0) {
540 bft_printf(_("\n Joining operation is conforming.\n"));
541 bft_printf_flush();
542 }
543
544 cs_join_inter_set_destroy(&inter_set);
545
546 }
547 else {
548
549 assert(join_type == CS_JOIN_TYPE_NON_CONFORMING);
550
551 if (cs_glob_mesh->verbosity > 0) {
552 bft_printf(_("\n Joining operation is non-conforming.\n"));
553 bft_printf_flush();
554 }
555
556 /* Creation of new vertices. Update list of equivalent vertices.
557 Associate to each intersection a vertex (old or created) */
558
559 cs_join_create_new_vertices(param.verbosity,
560 work_join_edges,
561 work_jmesh,
562 inter_set,
563 init_max_vtx_gnum,
564 &n_g_new_vertices,
565 &vtx_eset);
566
567 inter_edges = cs_join_inter_edges_define(work_join_edges, inter_set);
568 cs_join_inter_set_destroy(&inter_set);
569
570 #if 0 && defined(DEBUG) && !defined(NDEBUG) /* Dump structures after inter. */
571 cs_join_inter_edges_dump(inter_edges, work_join_edges, work_jmesh);
572 #endif
573
574 /* Synchronize inter_edges structure definition */
575
576 #if defined(HAVE_MPI)
577
578 if (n_ranks > 1 ) {
579
580 cs_join_inter_edges_t *sync_block = NULL;
581
582 sync_block = cs_join_inter_edges_part_to_block(work_jmesh,
583 work_join_edges,
584 inter_edges);
585
586 cs_join_inter_edges_block_to_part(work_join_edges->n_g_edges,
587 sync_block,
588 inter_edges);
589
590 cs_join_inter_edges_destroy(&sync_block);
591
592 /* Add new vertices to edge and mesh description if necessary */
593
594 cs_join_intersect_update_struct(param.verbosity,
595 work_join_edges,
596 work_jmesh,
597 &inter_edges);
598
599 cs_join_mesh_sync_vertices(work_jmesh);
600
601 }
602 #endif
603
604 /* Find if there are new equivalences between vertices on a same edge */
605
606 cs_join_add_equiv_from_edges(param,
607 work_jmesh,
608 work_join_edges,
609 inter_edges,
610 vtx_eset);
611
612 } /* non conforming joining operation */
613
614 /* Order and delete redundant equivalences */
615
616 cs_join_eset_clean(&vtx_eset);
617
618 /* Memory management: final state for vtx_eset (no more equiv. to get) */
619
620 vtx_eset->n_max_equiv = vtx_eset->n_equiv;
621 BFT_REALLOC(vtx_eset->equiv_couple, 2*vtx_eset->n_equiv, cs_lnum_t);
622
623 t1 = cs_timer_time();
624
625 cs_timer_counter_add_diff(&(this_join->stats.t_new_vtx), &t0, &t1);
626
627 if (param.verbosity > 0)
628 bft_printf(_("\n"
629 " Edge intersections and vertex creation done.\n"));
630
631 /* Returns pointers */
632
633 *p_vtx_eset = vtx_eset;
634 *p_inter_edges = inter_edges;
635 *p_n_g_new_vertices = n_g_new_vertices;
636
637 #if 0 && defined(DEBUG) && !defined(NDEBUG) /* Dump structures after inter. */
638 cs_join_inter_edges_dump(inter_edges, work_join_edges, work_jmesh);
639 #endif
640
641 return join_type;
642 }
643
644 #if defined(HAVE_MPI)
645
646 /*----------------------------------------------------------------------------
647 * Retrieve the local new global numbering for the initial vertices from
648 * the new vertex global numbering defined by block.
649 *
650 * parameters:
651 * param <-- set of user-defined parameter
652 * select <-- pointer to a cs_join_select_t structure
653 * mesh <-- pointer of pointer to cs_mesh_t structure
654 * init_max_vtx_gnum <-- initial max. global numbering for vertices
655 * p_o2n_vtx_gnum <-> in: array on blocks on the new global vertex
656 * out: local array on the new global vertex
657 *---------------------------------------------------------------------------*/
658
659 static void
_get_local_o2n_vtx_gnum(cs_join_param_t param,cs_join_select_t * select,cs_mesh_t * mesh,cs_gnum_t init_max_vtx_gnum,cs_gnum_t * p_o2n_vtx_gnum[])660 _get_local_o2n_vtx_gnum(cs_join_param_t param,
661 cs_join_select_t *select,
662 cs_mesh_t *mesh,
663 cs_gnum_t init_max_vtx_gnum,
664 cs_gnum_t *p_o2n_vtx_gnum[])
665 {
666 cs_lnum_t n_tot_vertices = mesh->n_vertices;
667
668 int *dest_rank = NULL;
669 cs_gnum_t *new_gnum_by_block = *p_o2n_vtx_gnum;
670 cs_gnum_t *new_local_gnum = NULL;
671 cs_all_to_all_t *d = NULL;
672
673 const int n_ranks = cs_glob_n_ranks;
674 const int local_rank = cs_glob_rank_id;
675
676 cs_block_dist_info_t bi = cs_block_dist_compute_sizes(local_rank,
677 n_ranks,
678 1,
679 0,
680 init_max_vtx_gnum);
681
682 MPI_Comm comm = cs_glob_mpi_comm;
683
684 if (param.perio_type != FVM_PERIODICITY_NULL)
685 n_tot_vertices = mesh->n_vertices + select->n_vertices;
686
687 /* Initialize new vertex gnum with old one */
688
689 BFT_MALLOC(new_local_gnum, n_tot_vertices, cs_gnum_t);
690 BFT_MALLOC(dest_rank, n_tot_vertices, int);
691
692 for (cs_lnum_t i = 0; i < mesh->n_vertices; i++) {
693 dest_rank[i] = (mesh->global_vtx_num[i] - 1)/(cs_gnum_t)(bi.block_size);
694 assert(dest_rank[i] >= 0 && dest_rank[i] < n_ranks);
695 new_local_gnum[i] = mesh->global_vtx_num[i];
696 }
697
698 if (param.perio_type != FVM_PERIODICITY_NULL) {
699 for (cs_lnum_t i = 0; i < select->n_vertices; i++) {
700 const cs_lnum_t j = mesh->n_vertices + i;
701 dest_rank[j] = (select->per_v_couples[2*i+1] - 1)
702 / (cs_gnum_t)(bi.block_size);
703 assert(dest_rank[j] >= 0 && dest_rank[j] < n_ranks);
704 new_local_gnum[j] = select->per_v_couples[2*i+1];
705 }
706 }
707
708 if (bi.rank_step > 1) {
709 for (cs_lnum_t i = 0; i < n_tot_vertices; i++)
710 dest_rank[i] *= bi.rank_step;
711 }
712
713 /* Send old vtx gnum to matching block */
714
715 d = cs_all_to_all_create(n_tot_vertices,
716 0, /* flags */
717 NULL, /* dest_id */
718 dest_rank,
719 comm);
720
721 cs_all_to_all_transfer_dest_rank(d, &dest_rank);
722
723 cs_gnum_t *b_data = cs_all_to_all_copy_array(d,
724 CS_GNUM_TYPE,
725 1,
726 false, /* reverse */
727 new_local_gnum,
728 NULL);
729
730 /* Request the new vtx gnum related to the initial vtx gnum */
731
732 cs_lnum_t n_b_vtx = cs_all_to_all_n_elts_dest(d);
733
734 for (cs_lnum_t i = 0; i < n_b_vtx; i++) {
735
736 /* Transform old to new vertex number */
737 cs_gnum_t o_shift = b_data[i] - bi.gnum_range[0];
738 b_data[i] = new_gnum_by_block[o_shift];
739
740 }
741
742 /* Send data back */
743
744 cs_all_to_all_copy_array(d,
745 CS_GNUM_TYPE,
746 1,
747 true, /* reverse */
748 b_data,
749 new_local_gnum);
750
751 BFT_FREE(b_data);
752
753 cs_all_to_all_destroy(&d);
754
755 /* Return pointer */
756
757 BFT_FREE(new_gnum_by_block);
758 *p_o2n_vtx_gnum = new_local_gnum;
759 }
760
761 #endif
762
763 /*----------------------------------------------------------------------------
764 * Define o2n_vtx_gnum for the current rank and in case of periodic apply
765 * changes from periodic faces to original faces. Update also per_v_couples
766 *
767 * parameters:
768 * this_join <-- pointer to a cs_join_t structure
769 * mesh <-- pointer to a cs_mesh_t structure
770 * local_jmesh <-> pointer to a local cs_join_mesh_t structure
771 * p_work_jmesh <-> distributed join mesh struct. on which operations
772 * take place
773 * p_work_edges <-> join edges struct. related to work_jmesh
774 * n_g_new_vertices <-- global number of vertices created with the
775 * intersection of edges
776 * init_max_vtx_gnum <-- initial max. global numbering for vertices
777 * o2n_vtx_gnum <-> in: array on blocks on the new global vertex
778 * out: local array on the new global vertex
779 *---------------------------------------------------------------------------*/
780 static void
_prepare_update_after_merge(cs_join_t * this_join,cs_mesh_t * mesh,cs_join_mesh_t * local_jmesh,cs_join_mesh_t ** p_work_jmesh,cs_join_edges_t ** p_work_edges,cs_gnum_t n_g_new_vertices,cs_gnum_t init_max_vtx_gnum,cs_gnum_t * p_o2n_vtx_gnum[])781 _prepare_update_after_merge(cs_join_t *this_join,
782 cs_mesh_t *mesh,
783 cs_join_mesh_t *local_jmesh,
784 cs_join_mesh_t **p_work_jmesh,
785 cs_join_edges_t **p_work_edges,
786 cs_gnum_t n_g_new_vertices,
787 cs_gnum_t init_max_vtx_gnum,
788 cs_gnum_t *p_o2n_vtx_gnum[])
789 {
790 int i, select_id, shift;
791
792 cs_gnum_t *o2n_vtx_gnum = *p_o2n_vtx_gnum;
793 cs_join_mesh_t *work_jmesh = *p_work_jmesh;
794 cs_join_edges_t *work_edges = *p_work_edges;
795 cs_join_select_t *selection = this_join->selection;
796 cs_join_param_t param = this_join->param;
797
798 const cs_lnum_t n_ranks = cs_glob_n_ranks;
799
800 /* Build an array keeping relation between old/new global vertex num. */
801
802 if (n_ranks == 1) {
803
804 cs_gnum_t *loc_vtx_gnum = NULL;
805
806 BFT_MALLOC(loc_vtx_gnum, mesh->n_vertices, cs_gnum_t);
807
808 /* Initialize array */
809
810 for (i = 0; i < mesh->n_vertices; i++)
811 loc_vtx_gnum[i] = i+1;
812
813 /* Update value for selected vertices */
814
815 for (i = 0, select_id = 0;
816 i < mesh->n_vertices && select_id < selection->n_vertices; i++) {
817
818 if (i + 1 == selection->vertices[select_id]) /* Is a selected vertex */
819 loc_vtx_gnum[i] = o2n_vtx_gnum[select_id++];
820
821 }
822
823 if (param.perio_type != FVM_PERIODICITY_NULL) { /* Update per_v_couples */
824
825 assert(selection->n_vertices == selection->n_couples);
826
827 for (i = 0; i < selection->n_vertices; i++) {
828 shift = selection->n_vertices + i;
829 selection->per_v_couples[2*i] = o2n_vtx_gnum[i];
830 selection->per_v_couples[2*i+1] = o2n_vtx_gnum[shift];
831 }
832
833 }
834
835 BFT_FREE(o2n_vtx_gnum);
836 o2n_vtx_gnum = loc_vtx_gnum; /* Without periodic vertices and for
837 all vertices in cs_mesh_t structure */
838
839 } /* End if serial mode */
840
841 #if defined(HAVE_MPI)
842 if (n_ranks > 1) {
843
844 _get_local_o2n_vtx_gnum(param,
845 selection,
846 mesh,
847 init_max_vtx_gnum,
848 &o2n_vtx_gnum);
849
850 if (param.perio_type != FVM_PERIODICITY_NULL) {
851
852 for (i = 0; i < selection->n_vertices; i++) {
853 select_id = selection->vertices[i] - 1;
854 shift = i + mesh->n_vertices;
855 selection->per_v_couples[2*i] = o2n_vtx_gnum[select_id];
856 selection->per_v_couples[2*i+1] = o2n_vtx_gnum[shift];
857 }
858
859 BFT_REALLOC(o2n_vtx_gnum, mesh->n_vertices, cs_gnum_t);
860
861 }
862
863 }
864 #endif
865
866 if (param.perio_type != FVM_PERIODICITY_NULL)
867 cs_join_perio_merge_back(this_join,
868 local_jmesh,
869 mesh,
870 &work_jmesh,
871 &work_edges,
872 init_max_vtx_gnum,
873 n_g_new_vertices);
874
875 /* Return pointer */
876
877 *p_o2n_vtx_gnum = o2n_vtx_gnum;
878 *p_work_jmesh = work_jmesh;
879 *p_work_edges = work_edges;
880
881 }
882
883 /*----------------------------------------------------------------------------
884 * Merge vertices from equivalences found between vertices.
885 * Update local and work structures after the merge step.
886 *
887 * parameters:
888 * this_join <-- pointer to a cs_join_t structure
889 * n_iwm_vertices <-- initial number of vertices in work struct.
890 * init_max_vtx_gnum <-- initial max. global numbering for vertices
891 * n_g_new_vertices <-- global number of vertices created with the
892 * intersection of edges
893 * vtx_eset <-> structure storing equivalences between vertices
894 * Two vertices are equivalent if they are each
895 * other in their tolerance; freed here after use
896 * inter_edges <-> structure storing the definition of new vertices
897 * on initial edges; freed here after use
898 * p_work_jmesh <-> pointer to a cs_join_mesh_t structure
899 * p_work_join_edges <-> pointer to a cs_join_edges_t structure
900 * p_local_jmesh <-> pointer to a cs_join_mesh_t structure
901 * mesh <-> pointer to a cs_mesh_t struct. to update
902 *---------------------------------------------------------------------------*/
903
904 static void
_merge_vertices(cs_join_t * this_join,cs_lnum_t n_iwm_vertices,cs_gnum_t init_max_vtx_gnum,cs_gnum_t n_g_new_vertices,cs_join_eset_t ** vtx_eset,cs_join_inter_edges_t ** inter_edges,cs_join_mesh_t ** p_work_jmesh,cs_join_edges_t ** p_work_join_edges,cs_join_mesh_t ** p_local_jmesh,cs_mesh_t * mesh)905 _merge_vertices(cs_join_t *this_join,
906 cs_lnum_t n_iwm_vertices,
907 cs_gnum_t init_max_vtx_gnum,
908 cs_gnum_t n_g_new_vertices,
909 cs_join_eset_t **vtx_eset,
910 cs_join_inter_edges_t **inter_edges,
911 cs_join_mesh_t **p_work_jmesh,
912 cs_join_edges_t **p_work_join_edges,
913 cs_join_mesh_t **p_local_jmesh,
914 cs_mesh_t *mesh)
915 {
916 int i;
917
918 cs_gnum_t new_max_vtx_gnum = init_max_vtx_gnum + n_g_new_vertices;
919 cs_gnum_t *iwm_vtx_gnum = NULL;
920 cs_gnum_t *o2n_vtx_gnum = NULL;
921 cs_join_mesh_t *local_jmesh = *p_local_jmesh;
922 cs_join_mesh_t *work_jmesh = *p_work_jmesh;
923 cs_join_edges_t *work_join_edges = *p_work_join_edges;
924 cs_join_param_t param = this_join->param;
925 cs_join_select_t *selection = this_join->selection;
926 cs_gnum_t *rank_face_gnum_index = selection->compact_rank_index;
927
928 assert(local_jmesh != NULL);
929 assert(work_jmesh != NULL);
930 assert(work_join_edges != NULL);
931
932 cs_timer_t t0 = cs_timer_time();
933
934 /*
935 Store the initial global vertex numbering
936 Initial vertices are between [0, n_init_vertices[
937 Added vertices from inter. are between [n_init_vertices, n_vertices]
938 */
939
940 BFT_MALLOC(iwm_vtx_gnum, n_iwm_vertices, cs_gnum_t);
941
942 for (i = 0; i < n_iwm_vertices; i++)
943 iwm_vtx_gnum[i] = (work_jmesh->vertices[i]).gnum;
944
945 /* Merge vertices */
946
947 cs_timer_t t1 = cs_timer_time();
948 cs_timer_counter_add_diff(&(this_join->stats.t_u_merge_vtx), &t0, &t1);
949 t0 = t1;
950
951 cs_join_merge_vertices(param,
952 new_max_vtx_gnum,
953 work_jmesh,
954 *vtx_eset);
955
956 t1 = cs_timer_time();
957 cs_timer_counter_add_diff(&(this_join->stats.t_merge_vtx), &t0, &t1);
958 t0 = t1;
959
960 cs_join_eset_destroy(vtx_eset);
961
962 /* Keep the evolution of vertex global numbering.
963 Update work and local structures after vertex merge */
964
965 cs_join_merge_update_struct(param,
966 n_iwm_vertices,
967 iwm_vtx_gnum,
968 init_max_vtx_gnum,
969 rank_face_gnum_index,
970 &work_jmesh,
971 &work_join_edges,
972 inter_edges,
973 &local_jmesh,
974 &o2n_vtx_gnum); /* Defined by slice in // */
975
976 /* Free memory */
977
978 BFT_FREE(iwm_vtx_gnum);
979
980 cs_join_inter_edges_destroy(inter_edges);
981
982 /* Post if required and level of verbosity is reached */
983
984 if (param.visualization > 3)
985 cs_join_post_dump_mesh("MergeBeforePerioWorkMesh", work_jmesh, param);
986
987 /* Define o2n_vtx_gnum for the current rank and
988 apply back periodic transformation if needed */
989
990 _prepare_update_after_merge(this_join,
991 mesh,
992 local_jmesh,
993 &work_jmesh,
994 &work_join_edges,
995 n_g_new_vertices,
996 init_max_vtx_gnum,
997 &o2n_vtx_gnum); /* Defined for the local rank */
998
999 /* Update cs_mesh_t structure after the vertex merge */
1000
1001 cs_join_update_mesh_after_merge(param,
1002 selection,
1003 o2n_vtx_gnum, /* free inside */
1004 local_jmesh,
1005 mesh);
1006
1007 /* Clean meshes after update (empty edges, degenerate edges, ...) */
1008
1009 cs_join_mesh_clean(work_jmesh, param.verbosity);
1010 cs_join_mesh_clean(local_jmesh, param.verbosity);
1011
1012 /* Define a new cs_join_edges_t structure */
1013
1014 cs_join_mesh_destroy_edges(&work_join_edges);
1015 work_join_edges = cs_join_mesh_define_edges(work_jmesh);
1016
1017 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1018 cs_join_mesh_dump_edges(work_join_edges, work_jmesh);
1019 #endif
1020
1021 t1 = cs_timer_time();
1022 cs_timer_counter_add_diff(&(this_join->stats.t_u_merge_vtx), &t0, &t1);
1023
1024 if (param.verbosity > 0) {
1025 bft_printf(_("\n"
1026 " Merge vertices and mesh update done.\n"));
1027 bft_printf_flush();
1028 }
1029
1030 /* Post if required and level of verbosity is reached */
1031
1032 if (param.visualization > 2)
1033 cs_join_post_dump_mesh("MergeWorkMesh", work_jmesh, param);
1034
1035 /* Set return pointers */
1036
1037 *p_local_jmesh = local_jmesh;
1038 *p_work_jmesh = work_jmesh;
1039 *p_work_join_edges = work_join_edges;
1040 }
1041
1042 /*----------------------------------------------------------------------------
1043 * Prepare mesh update after split operation. Invert face history and update
1044 * local_jmesh in case of periodicity by applying back split operation.
1045 *
1046 * parameters:
1047 * this_join <-- pointer to a cs_join_t structure
1048 * local_jmesh <-> pointer to a local cs_join_mesh_t structure
1049 * mesh <-- pointer to a cs_mesh_t structure
1050 * mesh_builder <-- pointer to a cs_mesh_builder_t structure
1051 * p_history <-> pointer to the history of face splitting
1052 * in: old -> new face links
1053 * out: new -> old face links
1054 *---------------------------------------------------------------------------*/
1055 static void
_prepare_update_after_split(cs_join_t * this_join,cs_join_mesh_t * local_jmesh,cs_mesh_t * mesh,cs_mesh_builder_t * mesh_builder,cs_join_gset_t ** p_history)1056 _prepare_update_after_split(cs_join_t *this_join,
1057 cs_join_mesh_t *local_jmesh,
1058 cs_mesh_t *mesh,
1059 cs_mesh_builder_t *mesh_builder,
1060 cs_join_gset_t **p_history)
1061 {
1062 cs_join_param_t param = this_join->param;
1063 cs_join_gset_t *o2n_hist = *p_history, *n2o_hist = NULL;
1064
1065 const cs_lnum_t n_ranks = cs_glob_n_ranks;
1066
1067 /* Invert face historic */
1068
1069 n2o_hist = cs_join_gset_invert(o2n_hist);
1070
1071 #if defined(HAVE_MPI)
1072 if (n_ranks > 1) {
1073
1074 cs_join_gset_t *n2o_sync_block = NULL;
1075
1076 MPI_Comm comm = cs_glob_mpi_comm;
1077
1078 n2o_sync_block = cs_join_gset_block_sync(local_jmesh->n_g_faces,
1079 n2o_hist,
1080 comm);
1081
1082 cs_join_gset_block_update(local_jmesh->n_g_faces,
1083 n2o_sync_block,
1084 n2o_hist,
1085 comm);
1086
1087 cs_join_gset_destroy(&n2o_sync_block);
1088
1089 }
1090 #endif
1091
1092 /* Build an array keeping relation between old/new global vertex num. */
1093
1094 if (param.perio_type != FVM_PERIODICITY_NULL)
1095 cs_join_perio_split_back(this_join,
1096 local_jmesh,
1097 mesh,
1098 mesh_builder,
1099 o2n_hist,
1100 &n2o_hist);
1101
1102 cs_join_gset_destroy(&o2n_hist);
1103
1104 /* Return pointer */
1105
1106 *p_history = n2o_hist;
1107 }
1108
1109 /*----------------------------------------------------------------------------
1110 * Split faces and update cs_mesh_t structure.
1111 *
1112 * parameters:
1113 * this_join <-- pointer to a cs_join_t structure
1114 * join_type <-> join type as detected so far
1115 * work_join_edges <-- pointer to a cs_join_edges_t structure
1116 * work_face_normal <-- normal based on the original face definition
1117 * rank_face_gnum_index <-- index on face global numering to determine the
1118 * related rank
1119 * p_work_jmesh <-> pointer to a cs_join_mesh_t structure
1120 * local_jmesh <-- pointer to a cs_join_mesh_t structure
1121 * mesh <-> pointer to cs_mesh_t structure
1122 * mesh_builder <-> pointer to cs_mesh_builder structure
1123 *---------------------------------------------------------------------------*/
1124
1125 static void
_split_faces(cs_join_t * this_join,cs_join_type_t * join_type,cs_join_edges_t * work_join_edges,cs_coord_t * work_face_normal,cs_join_mesh_t ** p_work_jmesh,cs_join_mesh_t * local_jmesh,cs_mesh_t * mesh,cs_mesh_builder_t * mesh_builder)1126 _split_faces(cs_join_t *this_join,
1127 cs_join_type_t *join_type,
1128 cs_join_edges_t *work_join_edges,
1129 cs_coord_t *work_face_normal,
1130 cs_join_mesh_t **p_work_jmesh,
1131 cs_join_mesh_t *local_jmesh,
1132 cs_mesh_t *mesh,
1133 cs_mesh_builder_t *mesh_builder)
1134 {
1135 cs_join_gset_t *history = NULL;
1136 cs_join_param_t param = this_join->param;
1137 cs_join_select_t *selection = this_join->selection;
1138 cs_gnum_t *rank_face_gnum_index = selection->compact_rank_index;
1139
1140 cs_timer_t t0 = cs_timer_time();
1141
1142 cs_join_split_faces(param,
1143 work_face_normal,
1144 work_join_edges,
1145 p_work_jmesh,
1146 &history); /* old -> new */
1147
1148 /* Send back to the original rank the new face description */
1149
1150 cs_join_split_update_struct(param,
1151 *p_work_jmesh,
1152 rank_face_gnum_index,
1153 &history, /* old -> new */
1154 &local_jmesh);
1155
1156 _prepare_update_after_split(this_join,
1157 local_jmesh,
1158 mesh,
1159 mesh_builder,
1160 &history); /* in: old->new, out: new -> old */
1161
1162 /* Update cs_mesh_t structure after the face splitting */
1163
1164 cs_join_update_mesh_after_split(param,
1165 selection,
1166 history, /* new -> old */
1167 local_jmesh,
1168 mesh,
1169 mesh_builder);
1170
1171 cs_timer_t t1 = cs_timer_time();
1172
1173 cs_timer_counter_add_diff(&(this_join->stats.t_split_faces), &t0, &t1);
1174
1175 bft_printf_flush();
1176
1177 /* Post if required and level of verbosity is reached */
1178
1179 if (param.visualization > 2)
1180 cs_join_post_dump_mesh("SplitWorkMesh", *p_work_jmesh, param);
1181
1182 /* Check join type if required */
1183
1184 if (cs_glob_mesh->verbosity > 0 && *join_type == CS_JOIN_TYPE_NULL) {
1185
1186 int _join_type = CS_JOIN_TYPE_NULL;
1187 for (cs_lnum_t i = 0; i < history->n_elts; i++) {
1188 if (history->index[i+1] != i+1) {
1189 _join_type = (int)CS_JOIN_TYPE_NON_CONFORMING;
1190 break;
1191 }
1192 else if (local_jmesh->face_gnum[i] != history->g_elts[i]) {
1193 _join_type = (int)CS_JOIN_TYPE_CONFORMING;
1194 break;
1195 }
1196 }
1197 cs_parall_max(1, CS_INT_TYPE, &_join_type);
1198
1199 *join_type = _join_type;
1200
1201 if (*join_type == CS_JOIN_TYPE_NULL)
1202 bft_printf(_("\n Joining operation is null.\n"));
1203 else if (*join_type == CS_JOIN_TYPE_CONFORMING)
1204 bft_printf(_("\n Joining operation is conforming.\n"));
1205 else if (*join_type == CS_JOIN_TYPE_NON_CONFORMING)
1206 bft_printf(_("\n Joining operation is non-conforming.\n"));
1207 bft_printf_flush();
1208
1209 }
1210
1211 /* Free memory */
1212
1213 cs_join_gset_destroy(&history);
1214 }
1215
1216 /*----------------------------------------------------------------------------
1217 * Print information on a mesh structure.
1218 *
1219 * parameters:
1220 * mesh <-- pointer to mesh structure.
1221 * name <-- associated name.
1222 *----------------------------------------------------------------------------*/
1223
1224 static void
_print_mesh_info(const cs_mesh_t * mesh,const char * name)1225 _print_mesh_info(const cs_mesh_t *mesh,
1226 const char *name)
1227 {
1228 bft_printf(_(" %s\n"
1229 " Number of cells: %llu\n"
1230 " Number of interior faces: %llu\n"
1231 " Number of boundary faces: %llu\n"
1232 " Number of vertices: %llu\n"),
1233 name,
1234 (unsigned long long)(mesh->n_g_cells),
1235 (unsigned long long)(mesh->n_g_i_faces),
1236 (unsigned long long)(mesh->n_g_b_faces),
1237 (unsigned long long)(mesh->n_g_vertices));
1238 }
1239
1240 /*----------------------------------------------------------------------------
1241 * Print initial joining info into log file
1242 *---------------------------------------------------------------------------*/
1243
1244 static void
_print_join_info(cs_mesh_t * mesh,cs_join_t * this_join,cs_join_param_t join_param)1245 _print_join_info(cs_mesh_t *mesh,
1246 cs_join_t *this_join,
1247 cs_join_param_t join_param)
1248 {
1249 bft_printf(_("\n -------------------------------------------------------\n"
1250 " Joining number %d:\n\n"), join_param.num);
1251
1252 if (join_param.perio_type != FVM_PERIODICITY_NULL) {
1253 double m[3][4];
1254 memcpy(m, join_param.perio_matrix, sizeof(double)*12);
1255 bft_printf(_(" Periodicity type: %s\n"),
1256 _(fvm_periodicity_type_name[join_param.perio_type]));
1257 bft_printf(_(" Transformation matrix: %12.5g %12.5g %12.5g %12.5g\n"
1258 " %12.5g %12.5g %12.5g %12.5g\n"
1259 " %12.5g %12.5g %12.5g %12.5g\n\n"),
1260 m[0][0], m[0][1], m[0][2], m[0][3],
1261 m[1][0], m[1][1], m[1][2], m[1][3],
1262 m[2][0], m[2][1], m[2][2], m[2][3]);
1263 }
1264
1265 bft_printf(_(" Selection criteria: \"%s\"\n"), this_join->criteria);
1266
1267 if (join_param.verbosity > 0) {
1268 bft_printf(_("\n"
1269 " Parameters for the joining operation:\n"
1270 " Shortest incident edge fraction: %8.5f\n"
1271 " Maximum angle between joined face planes: %8.5f\n\n"),
1272 join_param.fraction, join_param.plane);
1273
1274 bft_printf(_(" Advanced joining parameters:\n"
1275 " Verbosity level: %8d\n"
1276 " Visualization level: %8d\n"
1277 " Deepest level reachable in tree building: %8d\n"
1278 " Max boxes by leaf: %8d\n"
1279 " Max ratio of linked boxes / init. boxes: %8.5f\n"
1280 " Max ratio of boxes for distribution: %8.5f\n"
1281 " Merge step tolerance multiplier: %8.5f\n"
1282 " Pre-merge factor: %8.5f\n"
1283 " Tolerance computation mode: %8d\n"
1284 " Intersection computation mode: %8d\n"
1285 " Max. number of equiv. breaks: %8d\n"
1286 " Max. number of subfaces by face: %8d\n\n"),
1287 join_param.verbosity,
1288 join_param.visualization,
1289 join_param.tree_max_level,
1290 join_param.tree_n_max_boxes,
1291 join_param.tree_max_box_ratio,
1292 join_param.tree_max_box_ratio_distrib,
1293 join_param.merge_tol_coef,
1294 join_param.pre_merge_factor,
1295 join_param.tcm, join_param.icm,
1296 join_param.n_max_equiv_breaks,
1297 join_param.max_sub_faces);
1298
1299 _print_mesh_info(mesh, _(" Before joining"));
1300 bft_printf("\n");
1301 }
1302 }
1303
1304 /*----------------------------------------------------------------------------
1305 * Set advanced parameters to user-defined values.
1306 *
1307 * Out-of range values are silently set to minimum acceptable values
1308 * where those are possible.
1309 *
1310 * parameters:
1311 * join <-> pointer to a cs_join_t structure to update
1312 * mtf <-- merge tolerance coefficient
1313 * pmf <-- pre-merge factor
1314 * tcm <-- tolerance computation mode
1315 * icm <-- intersection computation mode
1316 * maxbrk <-- max number of equivalences to break (merge step)
1317 * max_sub_faces <-- max. possible number of sub-faces by splitting a face
1318 * tml <-- tree max level
1319 * tmb <-- tree max boxes
1320 * tmr <-- tree max ratio
1321 * tmr_distrib <-- tree max ratio for distribution
1322 *---------------------------------------------------------------------------*/
1323
1324 static void
_set_advanced_param(cs_join_t * join,double mtf,double pmf,int tcm,int icm,int maxbrk,int max_sub_faces,int tml,int tmb,double tmr,double tmr_distrib)1325 _set_advanced_param(cs_join_t *join,
1326 double mtf,
1327 double pmf,
1328 int tcm,
1329 int icm,
1330 int maxbrk,
1331 int max_sub_faces,
1332 int tml,
1333 int tmb,
1334 double tmr,
1335 double tmr_distrib)
1336 {
1337 /* Deepest level reachable during tree building */
1338
1339 if (tml < 1)
1340 tml = 1;
1341
1342 join->param.tree_max_level = tml;
1343
1344 /* Max. number of boxes which can be related to a leaf of the tree
1345 if level != tree_max_level */
1346
1347 if (tmb < 1)
1348 tmb = 1;
1349
1350 join->param.tree_n_max_boxes = tmb;
1351
1352 /* Stop tree building if:
1353 n_linked_boxes > tree_max_box_ratio*n_init_boxes */
1354
1355 if (tmr < 1.0 )
1356 tmr = 1.0;
1357
1358 join->param.tree_max_box_ratio = tmr;
1359
1360 if (tmr_distrib < 1.0 )
1361 tmr_distrib = 1.0;
1362
1363 join->param.tree_max_box_ratio_distrib = tmr_distrib;
1364
1365 /* Coef. used to modify the tolerance associated to each vertex BEFORE the
1366 merge operation.
1367 If coef = 0.0 => no vertex merge
1368 If coef < 1.0 => reduce vertex merge
1369 If coef = 1.0 => no change
1370 If coef > 1.0 => increase vertex merge */
1371
1372 if (mtf < 0.0)
1373 mtf = 0.0;
1374
1375 join->param.merge_tol_coef = mtf;
1376
1377 /* Maximum number of equivalence breaks */
1378
1379 if (maxbrk < 0)
1380 maxbrk = 0;
1381
1382 join->param.n_max_equiv_breaks = maxbrk;
1383
1384 /* Pre-merge factor. This parameter is used to define a limit
1385 under which two vertices are merged before the merge step.
1386 Tolerance limit for the pre-merge = pmf * fraction
1387 Default value: 0.10 */
1388
1389 join->param.pre_merge_factor = pmf;
1390
1391 /* Tolerance computation mode */
1392
1393 if ( (tcm)%10 < 1 || (tcm)%10 > 2)
1394 bft_error(__FILE__, __LINE__, 0,
1395 _("Mesh joining:"
1396 " Forbidden value for the tcm parameter.\n"
1397 " It must be 1, 2, 11, or 12 and not: %d\n"), tcm);
1398
1399 join->param.tcm = tcm;
1400
1401 /* Intersection computation mode */
1402
1403 if (icm != 1 && icm != 2)
1404 bft_error(__FILE__, __LINE__, 0,
1405 _("Mesh joining:"
1406 " Forbidden value for icm parameter.\n"
1407 " It must be 1 or 2 and not: %d\n"), icm);
1408
1409 join->param.icm = icm;
1410
1411 /* Maximum number of sub-faces */
1412
1413 if (max_sub_faces < 1)
1414 bft_error(__FILE__, __LINE__, 0,
1415 _("Mesh joining:"
1416 " Forbidden value for the maxsf parameter.\n"
1417 " It must be > 0 and not: %d\n"), max_sub_faces);
1418
1419 join->param.max_sub_faces = max_sub_faces;
1420
1421 }
1422
1423 /*----------------------------------------------------------------------------
1424 * Log statistics and timings for a given joining.
1425 *
1426 * parameters:
1427 * join <-- pointer to a cs_join_t structure
1428 *---------------------------------------------------------------------------*/
1429
1430 static void
_join_performance_log(const cs_join_t * this_join)1431 _join_performance_log(const cs_join_t *this_join)
1432 {
1433 char buf[80];
1434
1435 const cs_join_stats_t *stats = &(this_join->stats);
1436
1437 cs_log_printf(CS_LOG_PERFORMANCE,
1438 _("\nJoining number %d:\n\n"), this_join->param.num);
1439
1440 if (this_join->stats.n_calls > 1)
1441 cs_log_printf(CS_LOG_PERFORMANCE,
1442 _("\n Number of calls (statistics are cumulative): %d:\n\n"),
1443 this_join->stats.n_calls);
1444
1445 cs_log_printf(CS_LOG_PERFORMANCE,
1446 _(" Determination of possible face intersections:\n\n"
1447 " bounding-box tree layout: %dD\n"), stats->bbox_layout);
1448
1449 if (cs_glob_n_ranks > 1 || stats->n_calls > 1) {
1450
1451 cs_gnum_t n = CS_MAX(stats->n_calls, 1);
1452
1453 if (stats->n_calls > 1)
1454 strncpy(buf, _(" rank mean"), 79);
1455 else if (cs_glob_n_ranks <= 1)
1456 strncpy(buf, _(" call mean"), 79);
1457 else
1458 strncpy(buf, _(" rank/call mean"), 79);
1459
1460 cs_log_printf
1461 (CS_LOG_PERFORMANCE,
1462 _("%s minimum maximum\n"
1463 " depth: %10llu | %10llu | %10llu\n"
1464 " number of leaves: %10llu | %10llu | %10llu\n"
1465 " number of boxes: %10llu | %10llu | %10llu\n"
1466 " leaves over threshold: %10llu | %10llu | %10llu\n"
1467 " boxes per leaf: %10llu | %10llu | %10llu\n"
1468 " Memory footprint (kb):\n"
1469 " final search structure: %10llu | %10llu | %10llu\n"
1470 " temporary search structure: %10llu | %10llu | %10llu\n\n"),
1471 buf,
1472 (unsigned long long)(stats->bbox_depth[0] / n),
1473 (unsigned long long)stats->bbox_depth[1],
1474 (unsigned long long)stats->bbox_depth[2],
1475 (unsigned long long)(stats->n_leaves[0] / n),
1476 (unsigned long long)stats->n_leaves[1],
1477 (unsigned long long)stats->n_leaves[2],
1478 (unsigned long long)(stats->n_boxes[0] / n),
1479 (unsigned long long)stats->n_boxes[1],
1480 (unsigned long long)stats->n_boxes[2],
1481 (unsigned long long)(stats->n_th_leaves[0] / n),
1482 (unsigned long long)stats->n_th_leaves[1],
1483 (unsigned long long)stats->n_th_leaves[2],
1484 (unsigned long long)(stats->n_leaf_boxes[0] / n),
1485 (unsigned long long)stats->n_leaf_boxes[1],
1486 (unsigned long long)stats->n_leaf_boxes[2],
1487 (unsigned long long)(stats->box_mem_final[0] / n),
1488 (unsigned long long)stats->box_mem_final[1],
1489 (unsigned long long)stats->box_mem_final[2],
1490 (unsigned long long)(stats->box_mem_required[0] / n),
1491 (unsigned long long)stats->box_mem_required[1],
1492 (unsigned long long)stats->box_mem_required[2]);
1493
1494 }
1495 else
1496 cs_log_printf
1497 (CS_LOG_PERFORMANCE,
1498 _(" depth: %10llu\n"
1499 " number of leaves: %10llu\n"
1500 " number of boxes: %10llu\n"
1501 " leaves over threshold: %10llu\n"
1502 " boxes per leaf: %10llu mean [%llu min, %llu max]\n"
1503 " Memory footprint (kb):\n"
1504 " final search structure: %10llu\n"
1505 " temporary search structure: %10llu\n\n"),
1506 (unsigned long long)stats->bbox_depth[0],
1507 (unsigned long long)stats->n_leaves[0],
1508 (unsigned long long)stats->n_boxes[0],
1509 (unsigned long long)stats->n_th_leaves[0],
1510 (unsigned long long)stats->n_leaf_boxes[0],
1511 (unsigned long long)stats->n_leaf_boxes[1],
1512 (unsigned long long)stats->n_leaf_boxes[2],
1513 (unsigned long long)stats->box_mem_final[0],
1514 (unsigned long long)stats->box_mem_required[0]);
1515
1516 cs_log_printf
1517 (CS_LOG_PERFORMANCE,
1518 _(" Associated times:\n"
1519 " Face bounding boxes tree construction: %10.3g\n"
1520 " Face bounding boxes neighborhood query: %10.3g\n"),
1521 stats->t_box_build.nsec*1.e-9,
1522 stats->t_box_query.nsec*1.e-9);
1523
1524 cs_log_printf
1525 (CS_LOG_PERFORMANCE,
1526 _(" Sorting possible intersections between faces: %10.3g\n"),
1527 stats->t_inter_sort.nsec*1e-9);
1528
1529 cs_log_printf
1530 (CS_LOG_PERFORMANCE,
1531 _(" Definition of local joining mesh: %10.3g\n"),
1532 stats->t_l_join_mesh.nsec*1e-9);
1533
1534 cs_log_printf
1535 (CS_LOG_PERFORMANCE,
1536 _(" Edge intersections: %10.3g\n"),
1537 stats->t_edge_inter.nsec*1e-9);
1538
1539 cs_log_printf
1540 (CS_LOG_PERFORMANCE,
1541 _(" Creation of new vertices: %10.3g\n"),
1542 stats->t_new_vtx.nsec*1e-9);
1543
1544 cs_log_printf
1545 (CS_LOG_PERFORMANCE,
1546 _(" Merging vertices: %10.3g\n"),
1547 stats->t_merge_vtx.nsec*1e-9);
1548
1549 cs_log_printf
1550 (CS_LOG_PERFORMANCE,
1551 _(" Updating structures with vertex merging: %10.3g\n"),
1552 stats->t_u_merge_vtx.nsec*1e-9);
1553
1554 cs_log_printf
1555 (CS_LOG_PERFORMANCE,
1556 _(" Split old faces and reconstruct new faces: %10.3g\n"),
1557 stats->t_split_faces.nsec*1e-9);
1558
1559 cs_log_printf
1560 (CS_LOG_PERFORMANCE,
1561 _("\n"
1562 " Complete treatment for joining %2d:\n"
1563 " wall clock time: %10.3g\n"),
1564 this_join->param.num,
1565 stats->t_total.nsec*1e-9);
1566
1567 cs_log_printf_flush(CS_LOG_PERFORMANCE);
1568 }
1569
1570 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1571
1572 /*============================================================================
1573 * Public function definitions
1574 *===========================================================================*/
1575
1576 /*----------------------------------------------------------------------------*/
1577 /*!
1578 * Add a cs_join_t structure to the list of pending joinings.
1579 *
1580 * \param[in] sel_criteria boundary face selection criteria
1581 * \param[in] fraction value of the fraction parameter
1582 * \param[in] plane value of the plane parameter
1583 * \param[in] verbosity level of verbosity required
1584 * \param[in] visualization level of visualization required
1585 *
1586 * \return number (1 to n) associated with new joining
1587 */
1588 /*----------------------------------------------------------------------------*/
1589
1590 int
cs_join_add(const char * sel_criteria,float fraction,float plane,int verbosity,int visualization)1591 cs_join_add(const char *sel_criteria,
1592 float fraction,
1593 float plane,
1594 int verbosity,
1595 int visualization)
1596 {
1597 /* Allocate and initialize a cs_join_t structure */
1598
1599 BFT_REALLOC(cs_glob_join_array, cs_glob_n_joinings + 1, cs_join_t *);
1600
1601 cs_glob_join_array[cs_glob_n_joinings]
1602 = cs_join_create(cs_glob_n_joinings + 1,
1603 sel_criteria,
1604 fraction,
1605 plane,
1606 FVM_PERIODICITY_NULL,
1607 NULL,
1608 verbosity,
1609 visualization,
1610 true);
1611
1612 cs_glob_join_count++; /* Store number of joining (without periodic ones) */
1613 cs_glob_n_joinings++;
1614
1615 return cs_glob_n_joinings;
1616 }
1617
1618 /*----------------------------------------------------------------------------*/
1619 /*!
1620 * \brief Set advanced parameters for the joining algorithm.
1621 *
1622 * \param[in] join_num joining operation number
1623 * \param[in] mtf merge tolerance factor
1624 * \param[in] pmf pre-merge factor
1625 * \param[in] tcm tolerance computation mode
1626 * \param[in] icm intersection computation mode
1627 * \param[in] max_break max number of equivalences to break
1628 * (merge step)
1629 * \param[in] max_sub_faces max. possible number of sub-faces by
1630 * splitting a face
1631 * \param[in] tml tree max level
1632 * \param[in] tmb tree max boxes
1633 * \param[in] tmr tree max ratio
1634 * \param[in] tmr_distrib tree max ratio for distribution
1635 *
1636 * The first set of parameters allow modifing the intersection tolerance
1637 * detection and interection merge behavior:
1638 *
1639 * * `mtf`: merge tolerance factor, locally modifies the tolerance associated
1640 * to each vertex __before__ the merge step. The following cases apply:
1641 * - if *mtf = 0*, no vertex merge
1642 * - if *mtf < 1*, the vertex merge is stricter. It may increase the number
1643 * of tolerance reductions and therefore define smaller subsets of
1644 * vertices to merge together, possibly leading to more small edges.
1645 * - *mtf = 1* is the default
1646 * - if *mtf > 1*, the vertex merging is less strict. The subset of
1647 * vertices which can be merged is larger.
1648 * * `pmf`: a pre-merge factor. This parameter is used to define a limit under
1649 * which two vertices are merged before the merge step
1650 * (tolerance limit for the pre-merge = pmf * fraction),
1651 * * `tcm`: a tolerance computation mode. If its value is:
1652 * - `1` (default), the tolerance is the minimal edge length related to a
1653 * vertex, multiplied by a fraction.
1654 * - `2`, the tolerance is computed as in `1` with an additional
1655 * multiplication by a coefficient equal to *max(sin(e1), sin(e2)*,
1656 * where *e1* and *e2* are two edges sharing a same vertex *v* for
1657 * the tolerance is computed.
1658 * - `11`: similar to `1` but taking into account only the selected faces.
1659 * - `12`: similar to `2` but taking into account only the selected faces.
1660 * * `icm`: the intersection computation mode. If its value is:
1661 * - `1` (default), the original algorithm is used. Care should be taken
1662 * to clip the intersection to an extremity.
1663 * - `2`, a new intersection algorithm is used. Caution should be used to
1664 * avoid clipping the intersection to an extremity.
1665 * * `maxbrk`: defines the maximum number of equivalence breaks which is
1666 * enabled during the merge step,
1667 * * `max_sub_faces`: defines the maximum number of sub-faces allowed when
1668 * splitting a selected face.
1669 *
1670 * The following parameters used in the search algorithm for face
1671 * intersections between selected faces (octree-like structure). They allow
1672 * modifying the memory/speed tradeoffs, and are useful in case of excess
1673 * memory usage:
1674 *
1675 * * `tml`: the tree maximum level is the deepest level the tree can reach,
1676 * * `tmb`: the tree maximum boxes is the maximum number of bounding boxes (BB)
1677 * which can be linked to a leaf of the tree (not necessary true for
1678 * the deepest level),
1679 * * `tmr`: the tree maximum ratio. The building of the tree structure stops
1680 * when the number of bounding boxes is superior to the product of
1681 * `tmr` with the number of faces to locate. This is efficient
1682 * to reduce memory consumption.
1683 */
1684 /*----------------------------------------------------------------------------*/
1685
1686 void
cs_join_set_advanced_param(int join_num,double mtf,double pmf,int tcm,int icm,int max_break,int max_sub_faces,int tml,int tmb,double tmr,double tmr_distrib)1687 cs_join_set_advanced_param(int join_num,
1688 double mtf,
1689 double pmf,
1690 int tcm,
1691 int icm,
1692 int max_break,
1693 int max_sub_faces,
1694 int tml,
1695 int tmb,
1696 double tmr,
1697 double tmr_distrib)
1698 {
1699 int i, join_id = -1;
1700 cs_join_t *join = NULL;
1701
1702 /* Search for the joining structure related to join_num */
1703
1704 for (i = 0; i < cs_glob_n_joinings; i++) {
1705
1706 join = cs_glob_join_array[i];
1707 if (join_num == join->param.num) {
1708 join_id = i;
1709 break;
1710 }
1711
1712 }
1713
1714 if (join_id < 0)
1715 bft_error(__FILE__, __LINE__, 0,
1716 _(" Joining number %d is not defined.\n"), join_num);
1717
1718 assert(join != NULL);
1719
1720 _set_advanced_param(join,
1721 mtf,
1722 pmf,
1723 tcm,
1724 icm,
1725 max_break,
1726 max_sub_faces,
1727 tml,
1728 tmb,
1729 tmr,
1730 tmr_distrib);
1731 }
1732
1733 /*----------------------------------------------------------------------------
1734 * Apply all the defined joining operations.
1735 *
1736 * parameters:
1737 * preprocess <-- true if we are in the preprocessing stage
1738 *---------------------------------------------------------------------------*/
1739
1740 void
cs_join_all(bool preprocess)1741 cs_join_all(bool preprocess)
1742 {
1743 int join_id;
1744 double full_clock_start, full_clock_end;
1745
1746 cs_join_type_t join_type = CS_JOIN_TYPE_NULL;
1747 cs_mesh_t *mesh = cs_glob_mesh;
1748 cs_mesh_builder_t *mesh_builder = cs_glob_mesh_builder;
1749
1750 if (cs_glob_n_joinings < 1)
1751 return;
1752
1753 /* Sanity checks */
1754
1755 assert(sizeof(cs_lnum_t) == sizeof(cs_lnum_t));
1756 assert(sizeof(double) == sizeof(cs_real_t));
1757
1758 full_clock_start = cs_timer_wtime();
1759
1760 /* Disable all writers until explicitely enabled for this stage */
1761
1762 cs_post_disable_writer(0);
1763
1764 cs_join_post_init();
1765
1766 /* Loop on each defined joining to deal with */
1767
1768 for (join_id = 0; join_id < cs_glob_n_joinings; join_id++) {
1769
1770 if (cs_glob_join_array[join_id] == NULL)
1771 continue;
1772
1773 cs_join_t *this_join = cs_glob_join_array[join_id];
1774
1775 cs_join_param_t join_param = this_join->param;
1776
1777 if (join_param.preprocessing != preprocess)
1778 continue;
1779
1780 cs_timer_t t0 = cs_timer_time(); /* Start timer */
1781
1782 /* Open log file if required */
1783
1784 if (this_join->log_name != NULL) {
1785 cs_glob_join_log = fopen(this_join->log_name, "w");
1786 if (cs_glob_join_log == NULL)
1787 bft_error(__FILE__, __LINE__, errno,
1788 _("Unable to open file: \"%s\" for logging."),
1789 this_join->log_name);
1790 }
1791
1792 /* Print informations into log file */
1793
1794 if (mesh->verbosity > 0)
1795 _print_join_info(mesh, this_join, join_param);
1796
1797 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1798 _dump_mesh(join_param.num, "InitMesh", mesh);
1799 #endif
1800
1801 /* Build arrays and structures required for selection;
1802 will be destroyed after joining and rebuilt for each new join
1803 operation in order to take into account mesh modification */
1804
1805 _select_entities(this_join, mesh);
1806
1807 /* Now execute the joining operation */
1808
1809 if (this_join->selection->n_g_faces > 0) {
1810
1811 cs_lnum_t n_iwm_vertices; /* iwm: initial work mesh */
1812 cs_gnum_t init_max_vtx_gnum, n_g_new_vertices;
1813
1814 cs_real_t *work_face_normal = NULL;
1815 cs_join_gset_t *edge_edge_visibility = NULL;
1816 cs_join_mesh_t *work_jmesh = NULL, *local_jmesh = NULL;
1817 cs_join_edges_t *work_join_edges = NULL;
1818 cs_join_eset_t *vtx_eset = NULL;
1819 cs_join_inter_edges_t *inter_edges = NULL;
1820
1821 if (join_param.perio_type != FVM_PERIODICITY_NULL) {
1822
1823 cs_join_perio_init(this_join, mesh, &mesh_builder);
1824
1825 if (cs_glob_mesh_builder == NULL)
1826 cs_glob_mesh_builder = mesh_builder;
1827 }
1828
1829 _build_join_structures(this_join,
1830 mesh,
1831 &local_jmesh,
1832 &work_jmesh,
1833 &work_join_edges,
1834 &work_face_normal,
1835 &edge_edge_visibility);
1836
1837 if (mesh->verbosity > 0 && join_param.verbosity > 2)
1838 bft_printf(_("\n Number of faces to treat locally: %ld\n"),
1839 (long)work_jmesh->n_faces);
1840
1841 /*
1842
1843 Define new vertices and/or update old vertices from the real
1844 intersection found between edges,
1845 Keep the relation between two intersecting edges through an
1846 equivalence between the vertex of each edge.
1847 Store also the new description of the initial edges through a
1848 cs_join_inter_edges_t structure and synchronize it to get all the
1849 possible equivalences between vertices.
1850 Work mesh structure is not yet fully updated by the new vertices
1851 because the synchronization step has to be done.
1852
1853 */
1854
1855 n_iwm_vertices = work_jmesh->n_vertices;
1856
1857 init_max_vtx_gnum = mesh->n_g_vertices;
1858 if (join_param.perio_type != FVM_PERIODICITY_NULL)
1859 init_max_vtx_gnum += this_join->selection->n_g_vertices;
1860
1861 join_type
1862 = _intersect_edges(this_join,
1863 work_jmesh,
1864 work_join_edges,
1865 &edge_edge_visibility, /* free during this step */
1866 init_max_vtx_gnum,
1867 &n_g_new_vertices,
1868 &vtx_eset,
1869 &inter_edges);
1870
1871 /*
1872 Merge vertices from equivalences found between vertices.
1873 Update work structures after the merge step.
1874 Keep the evolution of the global numbering of initial vertices and
1875 get the sync_block cs_join_inter_edges_t structure to enable the
1876 local structure update after the merge step.
1877 */
1878
1879 if (join_type != CS_JOIN_TYPE_NULL)
1880 _merge_vertices(this_join,
1881 n_iwm_vertices,
1882 init_max_vtx_gnum,
1883 n_g_new_vertices,
1884 &vtx_eset, /* free during this step */
1885 &inter_edges, /* free during this step */
1886 &work_jmesh,
1887 &work_join_edges,
1888 &local_jmesh,
1889 mesh);
1890
1891 else {
1892 cs_join_eset_destroy(&vtx_eset);
1893 cs_join_inter_edges_destroy(&inter_edges);
1894 }
1895
1896 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1897 _dump_mesh(join_param.num, "MeshAfterMerge", mesh);
1898 #endif
1899
1900 /*
1901 Split faces in work_jmesh. Apply modification to the
1902 local_jmesh. Keep a history between old --> new faces.
1903 Update cs_mesh_t structure.
1904 */
1905
1906 _split_faces(this_join,
1907 &join_type,
1908 work_join_edges,
1909 work_face_normal,
1910 &work_jmesh,
1911 local_jmesh,
1912 mesh,
1913 mesh_builder);
1914
1915 /* Free memory */
1916
1917 cs_join_mesh_destroy(&local_jmesh);
1918 cs_join_mesh_destroy(&work_jmesh);
1919 cs_join_mesh_destroy_edges(&work_join_edges);
1920
1921 BFT_FREE(work_face_normal);
1922
1923 /* Clean mesh (delete redundant edge definition) */
1924
1925 cs_join_update_mesh_clean(join_param, mesh);
1926
1927 }
1928 else
1929 bft_printf(_("\nStop joining algorithm: no face selected...\n"));
1930
1931 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1932 _dump_mesh(join_param.num, "FinalMesh", mesh);
1933 #endif
1934
1935 if (mesh->verbosity > 0 && join_param.verbosity > 0) {
1936 bft_printf("\n");
1937 _print_mesh_info(mesh, _(" After joining"));
1938 bft_printf("\n");
1939 }
1940
1941 /* Free temporary structures */
1942
1943 cs_join_select_destroy(this_join->param, &(this_join->selection));
1944
1945 /* Optional synchronization (to be safe) */
1946
1947 #if defined(HAVE_MPI)
1948 if (cs_glob_n_ranks > 1)
1949 MPI_Barrier(cs_glob_mpi_comm);
1950 #endif
1951
1952 /* Close log file if present */
1953
1954 if (cs_glob_join_log != NULL) {
1955 if (fclose(cs_glob_join_log) != 0)
1956 bft_error(__FILE__, __LINE__, errno,
1957 _("Error closing log file for joining: %d."),
1958 this_join->param.num);
1959 cs_glob_join_log = NULL;
1960 }
1961
1962 /* Timing */
1963
1964 cs_timer_t t1 = cs_timer_time();
1965 cs_timer_counter_add_diff(&(this_join->stats.t_total), &t0, &t1);
1966
1967 if (mesh->verbosity > 0) {
1968 bft_printf(_("\n"
1969 " Joining %2d completed (%.3g s).\n"),
1970 join_param.num,
1971 this_join->stats.t_total.nsec*1e-9);
1972 bft_printf_flush();
1973 }
1974
1975 /* Free memory */
1976
1977 this_join->stats.n_calls += 1;
1978
1979 if (join_param.preprocessing) {
1980 _join_performance_log(this_join);
1981 cs_join_destroy(&this_join);
1982 cs_glob_join_array[join_id] = NULL;
1983 }
1984
1985 /* Force joining visualization level to 0 when the joining is not
1986 a part of preprocessing to avoid the increase of postprocessing
1987 meshes during the computation */
1988
1989 if (!join_param.preprocessing) {
1990 this_join->param.visualization = 0;
1991 }
1992
1993 /* Set mesh modification flag */
1994
1995 if (join_type != CS_JOIN_TYPE_NULL)
1996 mesh->modified |= CS_MESH_MODIFIED;
1997
1998 } /* End of loop on joinings */
1999
2000 /* Destroy all remaining structures relative to joining operation
2001 if not needed anymore */
2002
2003 for (join_id = 0; join_id < cs_glob_n_joinings; join_id++) {
2004 if (cs_glob_join_array[join_id] != NULL)
2005 break;
2006 }
2007 if (join_id >= cs_glob_n_joinings) {
2008 BFT_FREE(cs_glob_join_array);
2009 cs_glob_n_joinings = 0;
2010 }
2011
2012 /* Re-enable writers disabled when entering this stage */
2013
2014 cs_post_enable_writer(0);
2015
2016 full_clock_end = cs_timer_wtime();
2017
2018 if (mesh->verbosity > 0) {
2019
2020 bft_printf(_("\n"
2021 " All joining operations successfully finished:\n"
2022 "\n"
2023 " Wall clock time: %10.3g\n\n"),
2024 full_clock_end - full_clock_start);
2025 bft_printf_flush();
2026
2027 cs_log_printf(CS_LOG_PERFORMANCE,
2028 _("\n"
2029 "Joining operations time summary:\n"
2030 " wall clock time: %10.3g\n\n"),
2031 full_clock_end - full_clock_start);
2032
2033 cs_log_separator(CS_LOG_PERFORMANCE);
2034
2035 }
2036 }
2037
2038 /*----------------------------------------------------------------------------
2039 * Clear remaining memory for defined joining operations.
2040 *---------------------------------------------------------------------------*/
2041
2042 void
cs_join_finalize()2043 cs_join_finalize()
2044 {
2045 bool have_log = false;
2046
2047 for (int join_id = 0; join_id < cs_glob_n_joinings; join_id++) {
2048 if (cs_glob_join_array[join_id] != NULL) {
2049 have_log = true;
2050 _join_performance_log(cs_glob_join_array[join_id]);
2051 cs_join_destroy(&(cs_glob_join_array[join_id]));
2052 }
2053 }
2054
2055 BFT_FREE(cs_glob_join_array);
2056 cs_glob_n_joinings = 0;
2057
2058 if (have_log) {
2059 cs_log_printf(CS_LOG_PERFORMANCE, "\n");
2060 cs_log_separator(CS_LOG_PERFORMANCE);
2061 }
2062 }
2063
2064 /*----------------------------------------------------------------------------
2065 * Flag boundary faces that will be selected for joining.
2066 *
2067 * parameters:
2068 * mesh <-- pointer to mesh structure
2069 * preprocess <-- true if we are in the preprocessing stage
2070 * b_select_flag <-> false for boundary faces not selected for joining,
2071 * true for those selected
2072 *---------------------------------------------------------------------------*/
2073
2074 void
cs_join_mark_selected_faces(const cs_mesh_t * mesh,bool preprocess,bool b_select_flag[])2075 cs_join_mark_selected_faces(const cs_mesh_t *mesh,
2076 bool preprocess,
2077 bool b_select_flag[])
2078 {
2079 for (cs_lnum_t face_id = 0; face_id < mesh->n_b_faces; face_id++)
2080 b_select_flag[face_id] = false;
2081
2082 /* Count active joinings */
2083
2084 int n_joinings = 0;
2085
2086 for (int join_id = 0; join_id < cs_glob_n_joinings; join_id++) {
2087 if (cs_glob_join_array[join_id] == NULL)
2088 continue;
2089 cs_join_t *this_join = cs_glob_join_array[join_id];
2090 cs_join_param_t join_param = this_join->param;
2091 if (join_param.preprocessing == preprocess)
2092 n_joinings++;
2093 }
2094
2095 if (n_joinings < 1)
2096 return;
2097
2098 /* Prepare selection structures */
2099
2100 cs_lnum_t *b_face_list;
2101 BFT_MALLOC(b_face_list, mesh->n_b_faces, cs_lnum_t);
2102
2103 cs_real_t *b_face_cog = NULL, *b_face_normal = NULL;
2104
2105 cs_mesh_quantities_b_faces(mesh, &b_face_cog, &b_face_normal);
2106
2107 /* Build temporary selection structures */
2108
2109 const fvm_group_class_set_t *class_defs = mesh->class_defs;
2110 fvm_group_class_set_t *_class_defs = NULL;
2111
2112 if (class_defs == NULL) {
2113 _class_defs = fvm_group_class_set_create();
2114 class_defs = _class_defs;
2115 }
2116
2117 fvm_selector_t *select_b_faces = fvm_selector_create(mesh->dim,
2118 mesh->n_b_faces,
2119 class_defs,
2120 mesh->b_face_family,
2121 1,
2122 b_face_cog,
2123 b_face_normal);
2124
2125 /* Loop on each defined joining to deal with */
2126
2127 for (int join_id = 0; join_id < cs_glob_n_joinings; join_id++) {
2128
2129 if (cs_glob_join_array[join_id] == NULL)
2130 continue;
2131
2132 cs_join_t *this_join = cs_glob_join_array[join_id];
2133
2134 cs_join_param_t join_param = this_join->param;
2135
2136 if (join_param.preprocessing != preprocess)
2137 continue;
2138
2139 /* Extract and mark selected boundary faces */
2140
2141 cs_lnum_t n_b_faces = 0;
2142
2143 fvm_selector_get_list(select_b_faces,
2144 this_join->criteria,
2145 1,
2146 &n_b_faces,
2147 b_face_list);
2148
2149 for (cs_lnum_t i = 0; i < n_b_faces; i++) {
2150 cs_lnum_t face_id = b_face_list[i] - 1;
2151 b_select_flag[face_id] = true;
2152 }
2153
2154 } /* End of loop on joinings */
2155
2156 /* Free arrays and structures needed for selection */
2157
2158 BFT_FREE(b_face_cog);
2159 BFT_FREE(b_face_normal);
2160
2161 select_b_faces = fvm_selector_destroy(select_b_faces);
2162
2163 if (_class_defs != NULL)
2164 _class_defs = fvm_group_class_set_destroy(_class_defs);
2165
2166 BFT_FREE(b_face_list);
2167 }
2168
2169 /*---------------------------------------------------------------------------*/
2170
2171 #if 0 && defined(DEBUG) && !defined(NDEBUG)
2172 cs_debug_glob_mesh_dump("FinalGlobalVertices", mesh);
2173 #endif
2174
2175 END_C_DECLS
2176