1 /*============================================================================
2 * Manipulation of a cs_join_mesh_t structure
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 <string.h>
35 #include <math.h>
36 #include <float.h>
37
38 /*----------------------------------------------------------------------------
39 * Local headers
40 *---------------------------------------------------------------------------*/
41
42 #include "bft_mem.h"
43 #include "bft_printf.h"
44
45 #include "fvm_io_num.h"
46 #include "fvm_nodal.h"
47 #include "fvm_nodal_from_desc.h"
48 #include "fvm_nodal_order.h"
49
50 #include "cs_all_to_all.h"
51 #include "cs_order.h"
52 #include "cs_search.h"
53 #include "cs_join_post.h"
54 #include "cs_join_set.h"
55 #include "cs_join_util.h"
56 #include "cs_parall.h"
57
58 /*----------------------------------------------------------------------------
59 * Header for the current file
60 *---------------------------------------------------------------------------*/
61
62 #include "cs_join_mesh.h"
63
64 /*---------------------------------------------------------------------------*/
65
66 BEGIN_C_DECLS
67
68 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
69
70 /*============================================================================
71 * Macro and type definitions
72 *===========================================================================*/
73
74 /*============================================================================
75 * Private function definitions
76 *===========================================================================*/
77
78 /*----------------------------------------------------------------------------
79 * Return a string related to a state
80 *
81 * parameters:
82 * state <-- state to display
83 *
84 * returns:
85 * related string
86 *---------------------------------------------------------------------------*/
87
88 static const char *
_print_state(cs_join_state_t state)89 _print_state(cs_join_state_t state)
90 {
91 if (state == CS_JOIN_STATE_UNDEF)
92 return "UDF";
93 else if (state == CS_JOIN_STATE_NEW)
94 return "NEW";
95 else if (state == CS_JOIN_STATE_ORIGIN)
96 return "ORI";
97 else if (state == CS_JOIN_STATE_PERIO)
98 return "PER";
99 else if (state == CS_JOIN_STATE_MERGE)
100 return "MRG";
101 else if (state == CS_JOIN_STATE_PERIO_MERGE)
102 return "PMG";
103 else if (state == CS_JOIN_STATE_SPLIT)
104 return "SPL";
105 else
106 return "ERR";
107 }
108
109 /*----------------------------------------------------------------------------
110 * Compute the cross product of two vectors.
111 *
112 * parameters:
113 * v1 <-- first vector
114 * v2 <-- second vector
115 *
116 * returns:
117 * the resulting cross product (v1 x v2)
118 *---------------------------------------------------------------------------*/
119
120 inline static void
_cross_product(const double v1[],const double v2[],double result[])121 _cross_product(const double v1[],
122 const double v2[],
123 double result[])
124 {
125 result[0] = v1[1] * v2[2] - v2[1] * v1[2];
126 result[1] = v2[0] * v1[2] - v1[0] * v2[2];
127 result[2] = v1[0] * v2[1] - v2[0] * v1[1];
128 }
129
130 /*----------------------------------------------------------------------------
131 * Compute the dot product of two vectors.
132 *
133 * parameters:
134 * v1 <-- first vector
135 * v2 <-- second vector
136 *
137 * returns:
138 * the resulting dot product (v1.v2)
139 *---------------------------------------------------------------------------*/
140
141 inline static double
_dot_product(const double v1[],const double v2[])142 _dot_product(const double v1[],
143 const double v2[])
144 {
145 int i;
146 double result = 0.0;
147
148 for (i = 0; i < 3; i++)
149 result += v1[i] * v2[i];
150
151 return result;
152 }
153
154 /*----------------------------------------------------------------------------
155 * Get the norm of a vector.
156 *
157 * parameters:
158 * v <-> vector to work with.
159 *---------------------------------------------------------------------------*/
160
161 inline static double
_norm(double v[])162 _norm(double v[])
163 {
164 return sqrt(_dot_product(v, v));
165 }
166
167 /*----------------------------------------------------------------------------
168 * Compute the distance between two vertices.
169 *
170 * parameters:
171 * id <-- local id of the vertex to deal with
172 * quantities <-- array keeping edge vector and its length
173 *
174 * returns:
175 * sine in radians between edges sharing vertex id
176 *---------------------------------------------------------------------------*/
177
178 inline static double
_compute_sine(int id,const double quantities[])179 _compute_sine(int id,
180 const double quantities[])
181 {
182 int i;
183
184 double sine;
185 double norm_a, norm_b, a[3], b[3], c[3];
186
187 for (i = 0; i < 3; i++) {
188 a[i] = -quantities[4*id+i];
189 b[i] = quantities[4*(id+1)+i];
190 }
191
192 norm_a = quantities[4*id+3];
193 norm_b = quantities[4*(id+1)+3];
194
195 _cross_product(a, b, c);
196
197 sine = _norm(c) / (norm_a * norm_b);
198
199 return sine;
200 }
201
202 /*----------------------------------------------------------------------------
203 * Compute the distance between two vertices.
204 *
205 * parameters:
206 * a <-- coordinates of the first vertex.
207 * b <-- coordinates of the second vertex.
208 *
209 * returns:
210 * distance between a and b.
211 *---------------------------------------------------------------------------*/
212
213 inline static double
_compute_length(const double a[3],const double b[3])214 _compute_length(const double a[3],
215 const double b[3])
216 {
217 double length;
218
219 length = sqrt( (b[0] - a[0])*(b[0] - a[0])
220 + (b[1] - a[1])*(b[1] - a[1])
221 + (b[2] - a[2])*(b[2] - a[2]));
222
223 return length;
224 }
225
226 /*----------------------------------------------------------------------------
227 * Compute tolerance (mode 2)
228 * tolerance = min[ edge length * sine(v1v2) * fraction]
229 *
230 * parameters:
231 * vertex_coords <-- coordinates of vertices.
232 * vertex_tolerance <-> local tolerance affected to each vertex and
233 * to be updated
234 * n_faces <-- number of selected faces
235 * face_lst <-- list of faces selected to compute the tolerance
236 * f2v_idx <-- "face -> vertex" connect. index
237 * f2v_lst <-- "face -> vertex" connect. list
238 * fraction <-- parameter used to compute the tolerance
239 *---------------------------------------------------------------------------*/
240
241 static void
_compute_tolerance2(const cs_real_t vtx_coords[],double vtx_tolerance[],const cs_lnum_t n_faces,const cs_lnum_t face_lst[],const cs_lnum_t f2v_idx[],const cs_lnum_t f2v_lst[],double fraction)242 _compute_tolerance2(const cs_real_t vtx_coords[],
243 double vtx_tolerance[],
244 const cs_lnum_t n_faces,
245 const cs_lnum_t face_lst[],
246 const cs_lnum_t f2v_idx[],
247 const cs_lnum_t f2v_lst[],
248 double fraction)
249 {
250 int i, j, k, coord;
251 double tolerance, sine;
252 double a[3], b[3];
253
254 int n_max_face_vertices = 0;
255 int *face_connect = NULL;
256 double *edge_quantities = NULL;
257
258 for (i = 0; i < n_faces; i++) {
259 int fid = face_lst[i] - 1;
260 n_max_face_vertices = CS_MAX(n_max_face_vertices,
261 f2v_idx[fid+1] - f2v_idx[fid]);
262 }
263
264 BFT_MALLOC(face_connect, n_max_face_vertices + 1, int);
265 BFT_MALLOC(edge_quantities, 4 * (n_max_face_vertices + 1), double);
266
267 for (i = 0; i < n_faces; i++) {
268
269 int face_id = face_lst[i] - 1;
270 int start = f2v_idx[face_id];
271 int end = f2v_idx[face_id + 1];
272 int n_face_vertices = end - start;
273
274 /* Keep face connect */
275
276 for (k = 0, j = start; j < end; j++, k++)
277 face_connect[k] = f2v_lst[j];
278 face_connect[k] = f2v_lst[start];
279
280 /* Keep edge lengths and edge vectors:
281 - edge_quantities[4*k+0..2] = edge vector
282 - edge_quantities[4*k+3] = edge length */
283
284 for (k = 0; k < n_face_vertices; k++) {
285
286 for (coord = 0; coord < 3; coord++) {
287 a[coord] = vtx_coords[3*face_connect[k] + coord];
288 b[coord] = vtx_coords[3*face_connect[k+1] + coord];
289 edge_quantities[4*(k+1)+coord] = b[coord] - a[coord];
290 }
291 edge_quantities[4*(k+1)+3] = _compute_length(a, b);
292
293 }
294
295 for (coord = 0; coord < 4; coord++)
296 edge_quantities[coord] = edge_quantities[4*k+coord];
297
298 /* Loop on the vertices of the face to update tolerance on
299 each vertex */
300
301 for (k = 0; k < n_face_vertices; k++) {
302
303 int vid = face_connect[k];
304
305 tolerance = fraction * CS_MIN(edge_quantities[4*k+3],
306 edge_quantities[4*(k+1)+3]);
307 sine = _compute_sine(k, edge_quantities);
308
309 vtx_tolerance[vid] = CS_MIN(vtx_tolerance[vid], sine*tolerance);
310
311 }
312
313 } /* End of loop on faces */
314
315 BFT_FREE(face_connect);
316 BFT_FREE(edge_quantities);
317
318 }
319
320 /*----------------------------------------------------------------------------
321 * Compute tolerance (mode 1)
322 * tolerance = shortest edge length * fraction
323 *
324 * parameters:
325 * vertex_coords <-- coordinates of vertices.
326 * vertex_tolerance <-> local tolerance affected to each vertex and
327 * to be updated
328 * n_faces <-- number of selected faces
329 * face_lst <-- list of faces selected to compute the tolerance
330 * face_vtx_idx <-- "face -> vertex" connect. index
331 * face_vtx_lst <-- "face -> vertex" connect. list
332 * fraction <-- parameter used to compute the tolerance
333 *---------------------------------------------------------------------------*/
334
335 static void
_compute_tolerance1(const cs_real_t vtx_coords[],double vtx_tolerance[],const cs_lnum_t n_faces,const cs_lnum_t face_lst[],const cs_lnum_t face_vtx_idx[],const cs_lnum_t face_vtx_lst[],double fraction)336 _compute_tolerance1(const cs_real_t vtx_coords[],
337 double vtx_tolerance[],
338 const cs_lnum_t n_faces,
339 const cs_lnum_t face_lst[],
340 const cs_lnum_t face_vtx_idx[],
341 const cs_lnum_t face_vtx_lst[],
342 double fraction)
343 {
344 cs_lnum_t i, j, k, start, end, face_id, vtx_id1, vtx_id2;
345 double length, tolerance;
346 double a[3], b[3];
347
348 for (i = 0; i < n_faces; i++) {
349
350 face_id = face_lst[i] - 1;
351 start = face_vtx_idx[face_id];
352 end = face_vtx_idx[face_id + 1];
353
354 /* Loop on the vertices of the face */
355
356 for (j = start; j < end - 1; j++) {
357
358 vtx_id1 = face_vtx_lst[j];
359 vtx_id2 = face_vtx_lst[j+1];
360
361 for (k = 0; k < 3; k++) {
362 a[k] = vtx_coords[3*vtx_id1 + k];
363 b[k] = vtx_coords[3*vtx_id2 + k];
364 }
365
366 length = _compute_length(a, b);
367 tolerance = length * fraction;
368 vtx_tolerance[vtx_id1] = CS_MIN(vtx_tolerance[vtx_id1], tolerance);
369 vtx_tolerance[vtx_id2] = CS_MIN(vtx_tolerance[vtx_id2], tolerance);
370
371 }
372
373 /* Case end - start */
374
375 vtx_id1 = face_vtx_lst[end-1];
376 vtx_id2 = face_vtx_lst[start];
377
378 for (k = 0; k < 3; k++) {
379 a[k] = vtx_coords[3*vtx_id1 + k];
380 b[k] = vtx_coords[3*vtx_id2 + k];
381 }
382
383 length = _compute_length(a, b);
384 tolerance = length * fraction;
385 vtx_tolerance[vtx_id1] = CS_MIN(vtx_tolerance[vtx_id1], tolerance);
386 vtx_tolerance[vtx_id2] = CS_MIN(vtx_tolerance[vtx_id2], tolerance);
387
388 } /* End of loop on faces */
389
390 }
391
392 /*----------------------------------------------------------------------------
393 * Define for each vertex a tolerance which is the radius of the
394 * sphere in which the vertex can be fused with another vertex.
395 * This tolerance is computed from the given list of faces (interior or border)
396 *
397 * parameters:
398 * param <-- set of user-defined parameters for the joining
399 * vertex_coords <-- coordinates of vertices.
400 * n_faces <-- number of selected faces
401 * face_lst <-- list of faces selected to compute the tolerance
402 * face_vtx_idx <-- "face -> vertex" connect. index
403 * face_vtx_lst <-- "face -> vertex" connect. list
404 * vertex_tolerance <-> local tolerance affected to each vertex and
405 * to be updated
406 *---------------------------------------------------------------------------*/
407
408 static void
_get_local_tolerance(cs_join_param_t param,const cs_real_t vtx_coords[],const cs_lnum_t n_faces,const cs_lnum_t face_lst[],const cs_lnum_t face_vtx_idx[],const cs_lnum_t face_vtx_lst[],double vtx_tolerance[])409 _get_local_tolerance(cs_join_param_t param,
410 const cs_real_t vtx_coords[],
411 const cs_lnum_t n_faces,
412 const cs_lnum_t face_lst[],
413 const cs_lnum_t face_vtx_idx[],
414 const cs_lnum_t face_vtx_lst[],
415 double vtx_tolerance[])
416 {
417
418 if (param.tcm % 10 == 1) {
419
420 /* tol = min(edge length * fraction) */
421
422 _compute_tolerance1(vtx_coords,
423 vtx_tolerance,
424 n_faces,
425 face_lst,
426 face_vtx_idx,
427 face_vtx_lst,
428 param.fraction);
429
430 }
431 else if (param.tcm % 10 == 2) {
432
433 /* tol = min(edge length * sin(v1v2) * fraction) */
434
435 _compute_tolerance2(vtx_coords,
436 vtx_tolerance,
437 n_faces,
438 face_lst,
439 face_vtx_idx,
440 face_vtx_lst,
441 param.fraction);
442
443 }
444 else
445 bft_error(__FILE__, __LINE__, 0,
446 " Tolerance computation mode (%d) is not defined\n",
447 param.tcm);
448
449 }
450
451 #if defined(HAVE_MPI)
452
453 /*----------------------------------------------------------------------------
454 * Exchange local vertex tolerances to get a global vertex tolerance.
455 *
456 * parameters:
457 * n_vertices <-- number of local selected vertices
458 * select_vtx_io_num <-- fvm_io_num_t structure for the selected vertices
459 * vertex_data <-> data associated to each selected vertex
460 *---------------------------------------------------------------------------*/
461
462 static void
_get_global_tolerance(cs_lnum_t n_vertices,const fvm_io_num_t * select_vtx_io_num,cs_join_vertex_t vtx_data[])463 _get_global_tolerance(cs_lnum_t n_vertices,
464 const fvm_io_num_t *select_vtx_io_num,
465 cs_join_vertex_t vtx_data[])
466 {
467 double *g_vtx_tolerance = NULL, *part_tolerance = NULL;
468 cs_gnum_t n_g_vertices = fvm_io_num_get_global_count(select_vtx_io_num);
469 const cs_gnum_t *part_gnum = fvm_io_num_get_global_num(select_vtx_io_num);
470
471 MPI_Comm mpi_comm = cs_glob_mpi_comm;
472 const int local_rank = CS_MAX(cs_glob_rank_id, 0);
473 const int n_ranks = cs_glob_n_ranks;
474
475 cs_block_dist_info_t
476 bi = cs_block_dist_compute_sizes(local_rank,
477 n_ranks,
478 1,
479 0,
480 n_g_vertices);
481
482 cs_all_to_all_t
483 *d = cs_all_to_all_create_from_block(n_vertices,
484 0, /* flags */
485 part_gnum,
486 bi,
487 mpi_comm);
488
489 /* Send the global numbering for each vertex */
490
491 cs_gnum_t *block_gnum = cs_all_to_all_copy_array(d,
492 CS_GNUM_TYPE,
493 1,
494 false, /* reverse */
495 part_gnum,
496 NULL);
497
498 cs_lnum_t n_recv = cs_all_to_all_n_elts_dest(d);
499
500 /* Send the vertex tolerance for each vertex */
501
502 BFT_MALLOC(part_tolerance, n_vertices, double);
503
504 for (cs_lnum_t i = 0; i < n_vertices; i++)
505 part_tolerance[i] = vtx_data[i].tolerance;
506
507 double *recv_tolerance = cs_all_to_all_copy_array(d,
508 CS_REAL_TYPE,
509 1,
510 false, /* reverse */
511 part_tolerance,
512 NULL);
513
514 /* Define the global tolerance array */
515
516 BFT_MALLOC(g_vtx_tolerance, bi.block_size, double);
517
518 for (cs_lnum_t i = 0; i < bi.block_size; i++)
519 g_vtx_tolerance[i] = DBL_MAX;
520
521 const cs_gnum_t first_vtx_gnum = bi.gnum_range[0];
522
523 for (cs_lnum_t i = 0; i < n_recv; i++) {
524 cs_lnum_t vtx_id = block_gnum[i] - first_vtx_gnum;
525 g_vtx_tolerance[vtx_id] = CS_MIN(g_vtx_tolerance[vtx_id], recv_tolerance[i]);
526 }
527
528 /* Replace local vertex tolerance by the new computed global tolerance */
529
530 for (cs_lnum_t i = 0; i < n_recv; i++) {
531 cs_lnum_t vtx_id = block_gnum[i] - first_vtx_gnum;
532 recv_tolerance[i] = g_vtx_tolerance[vtx_id];
533 }
534
535 BFT_FREE(g_vtx_tolerance);
536
537 cs_all_to_all_copy_array(d,
538 CS_DOUBLE,
539 1,
540 true, /* reverse */
541 recv_tolerance,
542 part_tolerance);
543
544 for (cs_lnum_t i = 0; i < n_vertices; i++)
545 vtx_data[i].tolerance = part_tolerance[i];
546
547 /* Free memory */
548
549 BFT_FREE(part_tolerance);
550 BFT_FREE(recv_tolerance);
551 BFT_FREE(block_gnum);
552
553 cs_all_to_all_destroy(&d);
554 }
555 #endif /* HAVE_MPI */
556
557 /*----------------------------------------------------------------------------
558 * Define for each vertex a tolerance which is the radius of the
559 * sphere in which the vertex can be fused with another vertex.
560 * This tolerance is computed from the given list of faces (interior or border)
561 *
562 * parameters:
563 * param <-- set of user-defined parameters for the joining
564 * selection <-- pointer to a struct. keeping selected entities
565 * b_f2v_idx <-- border "face -> vertex" connectivity index
566 * b_f2v_lst <-- border "face -> vertex" connectivity
567 * i_f2v_idx <-- interior "face -> vertex" connectivity index
568 * i_f2v_lst <-- interior "face -> vertex" connectivity
569 * n_vertices <-- number of vertices in the parent mesh
570 * vtx_coord <-- coordinates of vertices in parent mesh
571 * vtx_gnum <-- global numbering of vertices
572 *
573 * returns:
574 * a pointer to an array of cs_join_vertex_t
575 *---------------------------------------------------------------------------*/
576
577 static cs_join_vertex_t *
_define_vertices(cs_join_param_t param,cs_join_select_t * selection,const cs_lnum_t b_f2v_idx[],const cs_lnum_t b_f2v_lst[],const cs_lnum_t i_f2v_idx[],const cs_lnum_t i_f2v_lst[],const cs_lnum_t n_vertices,const cs_real_t vtx_coord[],const cs_gnum_t vtx_gnum[])578 _define_vertices(cs_join_param_t param,
579 cs_join_select_t *selection,
580 const cs_lnum_t b_f2v_idx[],
581 const cs_lnum_t b_f2v_lst[],
582 const cs_lnum_t i_f2v_idx[],
583 const cs_lnum_t i_f2v_lst[],
584 const cs_lnum_t n_vertices,
585 const cs_real_t vtx_coord[],
586 const cs_gnum_t vtx_gnum[])
587 {
588 int i;
589
590 cs_join_vertex_t *vertices = NULL;
591
592 assert(selection != NULL);
593
594 const int n_ranks = cs_glob_n_ranks;
595
596 /*
597 Define a tolerance around each vertex in the selection list.
598 Tolerance is the radius of the sphere in which the vertex can be merged
599 with another vertex. Radius is the min(fraction * edge_length) on all
600 edges connected to a vertex.
601 Store all data about a vertex in a cs_join_vertex_t structure.
602 */
603
604 if (selection->n_vertices > 0) {
605
606 /* Initialize vertices array */
607
608 BFT_MALLOC(vertices, selection->n_vertices, cs_join_vertex_t);
609
610 #if defined(DEBUG) && !defined(NDEBUG)
611 /* Avoid Valgrind warnings in byte copies due to padding */
612 memset(vertices, 0, selection->n_vertices*sizeof(cs_join_vertex_t));
613 #endif
614
615 for (i = 0; i < selection->n_vertices; i++) {
616
617 cs_lnum_t vtx_id = selection->vertices[i]-1;
618
619 if (n_ranks > 1)
620 vertices[i].gnum = vtx_gnum[vtx_id];
621 else
622 vertices[i].gnum = vtx_id + 1;
623
624 vertices[i].coord[0] = vtx_coord[3*vtx_id];
625 vertices[i].coord[1] = vtx_coord[3*vtx_id+1];
626 vertices[i].coord[2] = vtx_coord[3*vtx_id+2];
627 vertices[i].tolerance = 0.0; /* Default value */
628 vertices[i].state = CS_JOIN_STATE_ORIGIN; /* Default value */
629
630 }
631
632
633 /* Compute the tolerance for each vertex of the mesh */
634
635 if (param.fraction > 0.0) {
636
637 double *vtx_tolerance = NULL;
638
639 BFT_MALLOC(vtx_tolerance, n_vertices, double);
640
641 for (i = 0; i < n_vertices; i++)
642 vtx_tolerance[i] = DBL_MAX;
643
644 /* Define local tolerance */
645
646 _get_local_tolerance(param,
647 vtx_coord,
648 selection->n_faces,
649 selection->faces,
650 b_f2v_idx,
651 b_f2v_lst,
652 vtx_tolerance);
653
654 if (param.tcm / 10 == 0) {
655
656 /* Update local tolerance with adjacent border faces */
657
658 _get_local_tolerance(param,
659 vtx_coord,
660 selection->n_b_adj_faces,
661 selection->b_adj_faces,
662 b_f2v_idx,
663 b_f2v_lst,
664 vtx_tolerance);
665
666 /* Update local tolerance with adjacent interior faces */
667
668 _get_local_tolerance(param,
669 vtx_coord,
670 selection->n_i_adj_faces,
671 selection->i_adj_faces,
672 i_f2v_idx,
673 i_f2v_lst,
674 vtx_tolerance);
675
676 } /* Include adjacent faces in the computation of the vertex tolerance */
677
678 for (i = 0; i < selection->n_vertices; i++)
679 vertices[i].tolerance = vtx_tolerance[selection->vertices[i]-1];
680
681 BFT_FREE(vtx_tolerance);
682
683 } /* End if tolerance > 0.0 */
684
685 #if 0 && defined(DEBUG) && !defined(NDEBUG) /* Sanity check */
686 for (i = 0; i < selection->n_vertices; i++)
687 if (vertices[i].tolerance > (DBL_MAX - 1.))
688 bft_error(__FILE__, __LINE__, 0,
689 "Incompatible value for the \"vertex tolerance\" item\n"
690 "Value must be lower than DBL_MAX and current value is : %f"
691 " (global numbering : %u)\n",
692 vertices[i].tolerance, vertices[i].gnum);
693 #endif
694
695 } /* End if selection->n_vertices > 0 */
696
697 #if defined(HAVE_MPI)
698 if (n_ranks > 1) { /* Parallel treatment : synchro over the ranks */
699
700 /* Global number of selected vertices and associated
701 fvm_io_num_t structure */
702
703 fvm_io_num_t *select_vtx_io_num = fvm_io_num_create(selection->vertices,
704 vtx_gnum,
705 selection->n_vertices,
706 0);
707
708 selection->n_g_vertices = fvm_io_num_get_global_count(select_vtx_io_num);
709
710 _get_global_tolerance(selection->n_vertices,
711 select_vtx_io_num,
712 vertices);
713
714 if (param.verbosity > 1)
715 bft_printf(_(" Global number of selected vertices: %11llu\n\n"),
716 (unsigned long long)(selection->n_g_vertices));
717
718 fvm_io_num_destroy(select_vtx_io_num);
719
720 }
721 #endif /* defined(HAVE_MPI) */
722
723 if (n_ranks == 1) {
724 selection->n_g_vertices = selection->n_vertices;
725
726 if (param.verbosity > 1)
727 bft_printf(_(" Number of selected vertices: %11llu\n\n"),
728 (unsigned long long)(selection->n_g_vertices));
729
730 }
731
732 return vertices;
733 }
734
735 #if defined(HAVE_MPI)
736 /*----------------------------------------------------------------------------
737 * Find for each face of the list its related rank
738 *
739 * parameters:
740 * n_elts <-- number of elements in the glob. list
741 * glob_list <-- global numbering list (must be ordered)
742 * rank_index <-- index defining a range of global numbers related
743 * to a rank
744 *
745 * returns:
746 * an array of size n_elts
747 *---------------------------------------------------------------------------*/
748
749 static int *
_get_rank_from_index(cs_lnum_t n_elts,const cs_gnum_t glob_list[],const cs_gnum_t rank_index[])750 _get_rank_from_index(cs_lnum_t n_elts,
751 const cs_gnum_t glob_list[],
752 const cs_gnum_t rank_index[])
753 {
754 cs_lnum_t i, rank;
755
756 int *rank_list = NULL;
757
758 if (n_elts == 0)
759 return NULL;
760
761 BFT_MALLOC(rank_list, n_elts, int);
762
763 for (i = 0, rank = 0; i < n_elts; i++) {
764
765 for (;rank_index[rank+1] < glob_list[i]; rank++);
766
767 assert(rank < cs_glob_n_ranks);
768 rank_list[i] = rank;
769
770 } /* End of loop on elements */
771
772 return rank_list;
773 }
774
775 /*----------------------------------------------------------------------------
776 * Get the index on ranks and the list of faces to send from a list of global
777 * faces to receive.
778 *
779 * parameters:
780 * gnum_rank_index <-- index on ranks for the global elements
781 * n_elts <-- number of elements to get
782 * glob_list <-- global number of faces to get (ordered)
783 * n_send --> number of face/rank couples to send
784 * send_rank --> list of ranks for the faces to send
785 * send_faces --> list of face ids to send
786 *---------------------------------------------------------------------------*/
787
788 static void
_get_send_faces(const cs_gnum_t gnum_rank_index[],cs_lnum_t n_elts,const cs_gnum_t glob_list[],cs_lnum_t * n_send,int * send_rank[],cs_lnum_t * send_faces[])789 _get_send_faces(const cs_gnum_t gnum_rank_index[],
790 cs_lnum_t n_elts,
791 const cs_gnum_t glob_list[],
792 cs_lnum_t *n_send,
793 int *send_rank[],
794 cs_lnum_t *send_faces[])
795 {
796 cs_lnum_t *_send_faces = NULL;
797
798 MPI_Comm comm = cs_glob_mpi_comm;
799
800 const int local_rank = CS_MAX(cs_glob_rank_id, 0);
801
802 /* Sanity checks */
803
804 assert(gnum_rank_index != NULL);
805
806 /* Find for each element of the list, the rank which owns the element */
807
808 int *gface_ranks = _get_rank_from_index(n_elts, glob_list, gnum_rank_index);
809
810 cs_gnum_t first_gface_id = gnum_rank_index[local_rank];
811
812 int flags = CS_ALL_TO_ALL_NEED_SRC_RANK;
813
814 cs_all_to_all_t
815 *d = cs_all_to_all_create(n_elts,
816 flags,
817 NULL, /* dest_id */
818 gface_ranks,
819 comm);
820
821 cs_all_to_all_transfer_dest_rank(d, &gface_ranks);
822
823 /* Exchange list of global num. to exchange */
824
825 cs_gnum_t *gfaces_to_send = cs_all_to_all_copy_array(d,
826 CS_GNUM_TYPE,
827 1,
828 false, /* reverse */
829 glob_list,
830 NULL);
831
832 cs_lnum_t _n_send = cs_all_to_all_n_elts_dest(d);
833
834 int *_send_rank = cs_all_to_all_get_src_rank(d);
835
836 cs_all_to_all_destroy(&d);
837
838 BFT_MALLOC(_send_faces, _n_send, cs_lnum_t);
839
840 /* Define face ids to send */
841
842 for (cs_lnum_t i = 0; i < _n_send; i++)
843 _send_faces[i] = gfaces_to_send[i] - 1 - first_gface_id;
844
845 /* Free memory */
846
847 BFT_FREE(gfaces_to_send);
848
849 /* Set return pointers */
850
851 *n_send = _n_send;
852 *send_rank = _send_rank;
853 *send_faces = _send_faces;
854 }
855
856 #endif /* HAVE_MPI */
857
858 /*----------------------------------------------------------------------------
859 * Clean the given cs_join_mesh_t structure: remove empty edges.
860 *
861 * parameters:
862 * mesh <-> pointer to the cs_join_mesh_t structure to clean
863 * verbosity <-- level of display
864 *---------------------------------------------------------------------------*/
865
866 static void
_remove_empty_edges(cs_join_mesh_t * mesh,int verbosity)867 _remove_empty_edges(cs_join_mesh_t *mesh,
868 int verbosity)
869 {
870 cs_lnum_t i, j, n_face_vertices;
871
872 cs_lnum_t shift = 0, n_simplified_faces = 0;
873 cs_lnum_t *new_face_vtx_idx = NULL;
874
875 assert(mesh != NULL);
876
877 BFT_MALLOC(new_face_vtx_idx, mesh->n_faces + 1, cs_lnum_t);
878
879 new_face_vtx_idx[0] = 0;
880
881 for (i = 0; i < mesh->n_faces; i++) {
882
883 cs_lnum_t s = mesh->face_vtx_idx[i];
884 cs_lnum_t e = mesh->face_vtx_idx[i+1];
885
886 if (mesh->face_vtx_lst[e-1] != mesh->face_vtx_lst[s])
887 mesh->face_vtx_lst[shift++] = mesh->face_vtx_lst[s];
888
889 /* Loop on face vertices */
890
891 for (j = s; j < e - 1; j++)
892 if (mesh->face_vtx_lst[j] != mesh->face_vtx_lst[j+1])
893 mesh->face_vtx_lst[shift++] = mesh->face_vtx_lst[j+1];
894
895 new_face_vtx_idx[i+1] = shift;
896
897 n_face_vertices = new_face_vtx_idx[i+1] - new_face_vtx_idx[i];
898
899 if (n_face_vertices < e - s) {
900
901 n_simplified_faces++;
902 if (verbosity > 3)
903 bft_printf(" Simplified face %ld (%llu)\n", (long)i+1,
904 (unsigned long long)(mesh->face_gnum[i]));
905
906 if (n_face_vertices < 3)
907 bft_error(__FILE__, __LINE__, 0,
908 _(" The simplified face has less than 3 vertices.\n"
909 " Check your joining parameters.\n"
910 " Face %ld (%llu)\n"), (long)i+1,
911 (unsigned long long)(mesh->face_gnum[i]));
912 }
913
914 } /* End of loop on faces */
915
916 BFT_FREE(mesh->face_vtx_idx);
917 mesh->face_vtx_idx = new_face_vtx_idx;
918
919 BFT_REALLOC(mesh->face_vtx_lst,
920 new_face_vtx_idx[mesh->n_faces], cs_lnum_t);
921
922 if (verbosity > 1) {
923 cs_gnum_t n_g_simplified_faces = n_simplified_faces;
924 cs_parall_counter(&n_g_simplified_faces, 1);
925 bft_printf(_("\n Number of simplified faces: %llu\n"),
926 (unsigned long long)n_simplified_faces);
927 }
928 }
929
930 /*----------------------------------------------------------------------------
931 * Clean the given cs_join_mesh_t structure: remove degenerate edges.
932 *
933 * parameters:
934 * mesh <-> pointer to the cs_join_mesh_t structure to clean
935 * verbosity <-- level of display
936 *---------------------------------------------------------------------------*/
937
938 static void
_remove_degenerate_edges(cs_join_mesh_t * mesh,int verbosity)939 _remove_degenerate_edges(cs_join_mesh_t *mesh,
940 int verbosity)
941 {
942 /*
943 - In the definition of faces based on new edges, a same edge may be
944 traversed twice, in the opposite direction; this is due to merging
945 of edges, as shown below.
946
947 x x
948 |\ |
949 | \ |
950 a2| \a3 A2|
951 | \ |
952 | \ a4 Merge of |
953 ---s1----s2------ vertices x
954 | \ s1 and s2 / \
955 | \ / \
956 a1| \a4 A1/ \A3
957 | \ / \
958 | \ / \
959 x-----------x x-----------x
960 a5 A4
961
962
963 Face: a1 a2 a3 a4 a5 Face: A1 A2 -A2 A3 A4
964
965
966 Caution: the final configuration may be
967 A2 A1 A3 A4 -A2
968 where the references of edges to delete may be at the
969 beginning or end of the face definition
970
971 Remark: several edge pairs may possibly be referenced twice,
972 in the form
973 ... A1 A2 -A2 -A1 ... (where the removal of A2 will
974 make ... A1 -A1 ... appear); We thus run as many passes
975 as necessary on a given face.
976 */
977
978 assert(mesh != NULL);
979
980 cs_lnum_t i, j, k, count, n_face_vertices;
981
982 cs_lnum_t shift = 0;
983 cs_lnum_t n_faces = mesh->n_faces;
984 cs_lnum_t n_modified_faces = 0;
985 cs_gnum_t n_g_modified_faces = 0;
986 cs_join_rset_t *tmp = NULL;
987 cs_join_rset_t *kill = NULL;
988
989 tmp = cs_join_rset_create(8);
990 kill = cs_join_rset_create(8);
991
992 for (i = 0; i < n_faces; i++) {
993
994 cs_lnum_t start_id = mesh->face_vtx_idx[i];
995 cs_lnum_t end_id = mesh->face_vtx_idx[i+1];
996 cs_lnum_t n_init_vertices = end_id - start_id;
997 cs_lnum_t n_elts = n_init_vertices + 2;
998
999 assert(n_init_vertices > 2);
1000
1001 /* Build a temporary list based on the face connectivity */
1002
1003 cs_join_rset_resize(&tmp, n_elts);
1004 cs_join_rset_resize(&kill, n_elts);
1005
1006 for (j = start_id, k = 0; j < end_id; j++, k++) {
1007 tmp->array[k] = mesh->face_vtx_lst[j] + 1;
1008 kill->array[k] = 0;
1009 }
1010
1011 tmp->array[k] = mesh->face_vtx_lst[start_id] + 1;
1012 kill->array[k++] = 0;
1013 tmp->array[k] = mesh->face_vtx_lst[start_id+1] + 1;
1014 kill->array[k++] = 0;
1015
1016 assert(n_elts == k);
1017 tmp->n_elts = n_elts;
1018 kill->n_elts = n_elts;
1019
1020 /* Find degenerate edges */
1021
1022 count = 1;
1023 n_face_vertices = n_init_vertices;
1024
1025 while (count > 0) {
1026
1027 count = 0;
1028 for (j = 0; j < n_face_vertices; j++) {
1029 if (tmp->array[j] == tmp->array[j+2]) {
1030 count++;
1031 kill->array[j] = 1;
1032 kill->array[(j+1)%n_face_vertices] = 1;
1033 }
1034 }
1035
1036 tmp->n_elts = 0;
1037 for (j = 0; j < n_face_vertices; j++) {
1038 if (kill->array[j] == 0)
1039 tmp->array[tmp->n_elts++] = tmp->array[j];
1040 }
1041
1042 n_face_vertices = tmp->n_elts;
1043 tmp->array[tmp->n_elts++] = tmp->array[0];
1044 tmp->array[tmp->n_elts++] = tmp->array[1];
1045
1046 kill->n_elts = tmp->n_elts;
1047 for (j = 0; j < kill->n_elts; j++)
1048 kill->array[j] = 0;
1049
1050 } /* End of while */
1051
1052 if (n_face_vertices != n_init_vertices) {
1053
1054 n_modified_faces += 1;
1055
1056 /* Display the degenerate face */
1057
1058 if (verbosity > 5) {
1059
1060 bft_printf("\n Remove edge for face: %ld [%llu]:",
1061 (long)i+1, (unsigned long long)(mesh->face_gnum[i]));
1062 bft_printf("\n Initial def: ");
1063 for (j = start_id; j < end_id; j++) {
1064 cs_lnum_t v_id = mesh->face_vtx_lst[j];
1065 bft_printf(" %ld (%llu) ", (long)v_id+1,
1066 (unsigned long long)(mesh->vertices[v_id].gnum));
1067 }
1068 bft_printf("\n Final def: ");
1069 for (j = 0; j < n_face_vertices; j++) {
1070 cs_lnum_t v_id = tmp->array[j] - 1;
1071 bft_printf(" %ld (%llu) ", (long)v_id+1,
1072 (unsigned long long)(mesh->vertices[v_id].gnum));
1073 }
1074 bft_printf("\n");
1075 bft_printf_flush();
1076
1077 }
1078
1079 if (n_face_vertices < 3)
1080 bft_error(__FILE__, __LINE__, 0,
1081 _(" The simplified face has less than 3 vertices.\n"
1082 " Check your joining parameters.\n"
1083 " Face %ld (%llu)\n"), (long)i+1,
1084 (unsigned long long)(mesh->face_gnum[i]));
1085
1086 } /* End if n_face_vertices != n_init_vertices */
1087
1088 for (j = 0; j < n_face_vertices; j++)
1089 mesh->face_vtx_lst[shift++] = tmp->array[j] - 1;
1090 mesh->face_vtx_idx[i] = shift;
1091
1092 } /* End of loop on faces */
1093
1094 n_g_modified_faces = n_modified_faces;
1095 cs_parall_counter(&n_g_modified_faces, 1);
1096
1097 if (verbosity > 0)
1098 bft_printf("\n Edge removed for %llu faces (global).\n"
1099 " Join mesh cleaning done.\n",
1100 (unsigned long long)n_g_modified_faces);
1101
1102 for (i = n_faces; i > 0; i--)
1103 mesh->face_vtx_idx[i] = mesh->face_vtx_idx[i-1];
1104 mesh->face_vtx_idx[0] = 0;
1105
1106 BFT_REALLOC(mesh->face_vtx_lst, mesh->face_vtx_idx[n_faces], cs_lnum_t);
1107
1108 /* Free memory */
1109
1110 cs_join_rset_destroy(&tmp);
1111 cs_join_rset_destroy(&kill);
1112 }
1113
1114 /*----------------------------------------------------------------------------
1115 * Count the number of new vertices to add in the new face definition
1116 *
1117 * parameters:
1118 * v1_id <-- first vertex id
1119 * v2_id <-- second vertex id
1120 * old2new <-- indirection array between old and new numbering
1121 * edges <-- cs_join_edges_t structure
1122 * edge_index <-- edge -> new added vertex connectivity index
1123 * edge_new_vtx_lst <-- edge -> new added vertex connectivity list
1124 *
1125 * returns:
1126 * a number of vertices to add
1127 *---------------------------------------------------------------------------*/
1128
1129 static int
_count_new_added_vtx_to_edge(cs_lnum_t v1_id,cs_lnum_t v2_id,const cs_lnum_t old2new[],const cs_join_edges_t * edges,const cs_lnum_t edge_index[],const cs_lnum_t edge_new_vtx_lst[])1130 _count_new_added_vtx_to_edge(cs_lnum_t v1_id,
1131 cs_lnum_t v2_id,
1132 const cs_lnum_t old2new[],
1133 const cs_join_edges_t *edges,
1134 const cs_lnum_t edge_index[],
1135 const cs_lnum_t edge_new_vtx_lst[])
1136 {
1137 cs_lnum_t i, edge_id, edge_num;
1138
1139 cs_lnum_t new_v1_id = old2new[v1_id];
1140 cs_lnum_t new_v2_id = old2new[v2_id];
1141 cs_lnum_t n_adds = 0;
1142
1143 assert(v1_id >= 0);
1144 assert(v2_id >= 0);
1145 assert(new_v1_id >= 0);
1146 assert(new_v2_id >= 0);
1147 assert(edge_index != NULL);
1148 assert(edges != NULL);
1149
1150 /* Find the related edge */
1151
1152 edge_num = cs_join_mesh_get_edge(v1_id+1, v2_id+1, edges);
1153 edge_id = CS_ABS(edge_num) - 1;
1154
1155 if (v1_id == v2_id)
1156 bft_error(__FILE__, __LINE__, 0,
1157 _("\n Problem in mesh connectivity.\n"
1158 " Detected when updating connectivity.\n"
1159 " Edge number: %ld (%llu) - (%ld, %ld) in old numbering.\n"),
1160 (long)edge_num, (unsigned long long)(edges->gnum[edge_id]),
1161 (long)v1_id, (long)v2_id);
1162
1163 /* Add the first vertex (new_v1_id) */
1164
1165 n_adds = 1;
1166
1167 /* Add another vertices if needed */
1168
1169 for (i = edge_index[edge_id]; i < edge_index[edge_id+1]; i++) {
1170
1171 cs_lnum_t new_vtx_id = edge_new_vtx_lst[i] - 1;
1172
1173 if (new_vtx_id != new_v1_id && new_vtx_id != new_v2_id)
1174 n_adds++;
1175
1176 }
1177
1178 return n_adds;
1179 }
1180
1181 /*----------------------------------------------------------------------------
1182 * Add new vertex to the face -> vertex connectivity
1183 *
1184 * parameters:
1185 * v1_id <-- first vertex id
1186 * v2_id <-- second vertex id
1187 * old2new <-- indirection array between old and new numbering
1188 * edges <-- cs_join_edges_t structure
1189 * edge_index <-- edge -> new added vertex connectivity index
1190 * edge_new_vtx_lst <-- edge -> new added vertex connectivity list
1191 * new_face_vtx_lst <-> new face -> vertex connectivity list
1192 * p_shift <-> pointer to the shift in the connectivity list
1193 *---------------------------------------------------------------------------*/
1194
1195 static void
_add_new_vtx_to_edge(cs_lnum_t v1_id,cs_lnum_t v2_id,const cs_lnum_t old2new[],const cs_join_edges_t * edges,const cs_lnum_t edge_index[],const cs_lnum_t edge_new_vtx_lst[],cs_lnum_t new_face_vtx_lst[],cs_lnum_t * p_shift)1196 _add_new_vtx_to_edge(cs_lnum_t v1_id,
1197 cs_lnum_t v2_id,
1198 const cs_lnum_t old2new[],
1199 const cs_join_edges_t *edges,
1200 const cs_lnum_t edge_index[],
1201 const cs_lnum_t edge_new_vtx_lst[],
1202 cs_lnum_t new_face_vtx_lst[],
1203 cs_lnum_t *p_shift)
1204 {
1205 cs_lnum_t new_v1_id = old2new[v1_id];
1206 cs_lnum_t shift = *p_shift;
1207
1208 assert(edges != NULL);
1209
1210 /* Add first vertex num to the connectivity list */
1211
1212 new_face_vtx_lst[shift++] = new_v1_id;
1213
1214 if (edge_new_vtx_lst != NULL) {
1215
1216 cs_lnum_t i, edge_id, edge_num, e_start, e_end;
1217
1218 cs_lnum_t new_v2_id = old2new[v2_id];
1219
1220 /* Find the related edge */
1221
1222 edge_num = cs_join_mesh_get_edge(v1_id+1, v2_id+1, edges);
1223 edge_id = CS_ABS(edge_num) - 1;
1224 e_start = edge_index[edge_id];
1225 e_end = edge_index[edge_id+1];
1226
1227 /* Add a vertex if needed */
1228
1229 if (edge_num > 0) {
1230
1231 for (i = e_start; i < e_end; i++) {
1232
1233 cs_lnum_t new_vtx_id = edge_new_vtx_lst[i] - 1;
1234
1235 if (new_vtx_id != new_v1_id && new_vtx_id != new_v2_id)
1236 new_face_vtx_lst[shift++] = new_vtx_id;
1237
1238 }
1239 }
1240 else { /* edge_num < 0 */
1241
1242 for (i = e_end - 1; i > e_start - 1; i--) {
1243
1244 cs_lnum_t new_vtx_id = edge_new_vtx_lst[i] - 1;
1245
1246 if (new_vtx_id != new_v1_id && new_vtx_id != new_v2_id)
1247 new_face_vtx_lst[shift++] = new_vtx_id;
1248
1249 }
1250
1251 } /* End if edge_num < 0 */
1252
1253 } /* End if edge_new_vtx_lst != NULL */
1254
1255 /* Return pointer */
1256
1257 *p_shift = shift;
1258 }
1259
1260 #if defined(HAVE_MPI)
1261
1262 /*----------------------------------------------------------------------------
1263 * Dump a cs_join_vertex_t structure into a file.
1264 *
1265 * parameters:
1266 * f <-- handle to output file
1267 * vertex <-- cs_join_vertex_t structure to dump
1268 *---------------------------------------------------------------------------*/
1269
1270 static void
_log_vertex(cs_join_vertex_t vertex)1271 _log_vertex(cs_join_vertex_t vertex)
1272 {
1273 assert(vertex.gnum > 0);
1274 assert(vertex.tolerance >= 0.0);
1275
1276 bft_printf(" %10llu | %11.6f | % 12.10e % 12.10e % 12.10e | %s\n",
1277 (unsigned long long)vertex.gnum, vertex.tolerance,
1278 vertex.coord[0], vertex.coord[1], vertex.coord[2],
1279 _print_state(vertex.state));
1280 }
1281
1282 /*----------------------------------------------------------------------------
1283 * Create a MPI_Datatype for the cs_join_vertex_t structure.
1284 *
1285 * returns:
1286 * a MPI_Datatype associated to the cs_join_vertex_t structure
1287 *---------------------------------------------------------------------------*/
1288
1289 static MPI_Datatype
_create_vtx_datatype(void)1290 _create_vtx_datatype(void)
1291 {
1292 int j;
1293 cs_join_vertex_t v_data;
1294 MPI_Datatype new_type;
1295
1296 int blocklengths[4] = {1, 1, 1, 3};
1297 MPI_Aint displacements[4] = {0 , 0, 0, 0};
1298 MPI_Datatype types[4] = {MPI_INT, CS_MPI_GNUM, MPI_DOUBLE, CS_MPI_COORD};
1299
1300 /* Initialize vertex structure */
1301
1302 v_data.state = CS_JOIN_STATE_UNDEF;
1303 v_data.gnum = 1;
1304 v_data.tolerance = 0.0;
1305 for (j = 0; j < 3; j++)
1306 v_data.coord[j] = 0.0;
1307
1308 /* Define array of displacements */
1309
1310 #if defined(MPI_VERSION) && (MPI_VERSION >= 2)
1311 MPI_Get_address(&v_data, displacements);
1312 MPI_Get_address(&v_data.gnum, displacements + 1);
1313 MPI_Get_address(&v_data.tolerance, displacements + 2);
1314 MPI_Get_address(&v_data.coord, displacements + 3);
1315 #else
1316 MPI_Address(&v_data, displacements);
1317 MPI_Address(&v_data.gnum, displacements + 1);
1318 MPI_Address(&v_data.tolerance, displacements + 2);
1319 MPI_Address(&v_data.coord, displacements + 3);
1320 #endif
1321
1322 displacements[3] -= displacements[0];
1323 displacements[2] -= displacements[0];
1324 displacements[1] -= displacements[0];
1325 displacements[0] = 0;
1326
1327 /* Create new datatype */
1328
1329 #if (MPI_VERSION >= 2)
1330 MPI_Type_create_struct(4, blocklengths, displacements, types, &new_type);
1331 #else
1332 MPI_Type_struct(4, blocklengths, displacements, types, &new_type);
1333 #endif
1334
1335 MPI_Type_commit(&new_type);
1336
1337 return new_type;
1338 }
1339
1340 /*----------------------------------------------------------------------------
1341 * Create a function to define an operator for MPI reduction operation
1342 *
1343 * parameters:
1344 * in <-- input vertices
1345 * inout <-> in/out vertices (vertex with the min. tolerance)
1346 * len <-- size of input array
1347 * datatype <-- MPI_datatype associated to cs_join_vertex_t
1348 *---------------------------------------------------------------------------*/
1349
1350 static void
_mpi_vertex_min(cs_join_vertex_t * in,cs_join_vertex_t * inout,int * len,MPI_Datatype * datatype)1351 _mpi_vertex_min(cs_join_vertex_t *in,
1352 cs_join_vertex_t *inout,
1353 int *len,
1354 MPI_Datatype *datatype)
1355 {
1356 CS_UNUSED(datatype);
1357
1358 int i, j;
1359
1360 assert(in != NULL && inout != NULL);
1361
1362 for (i = 0; i < *len; i++) {
1363
1364 if (in->tolerance <= inout->tolerance) {
1365
1366 if (in->tolerance < inout->tolerance) {
1367
1368 inout->gnum = in->gnum;
1369 for (j = 0; j < 3; j++)
1370 inout->coord[j] = in->coord[j];
1371 inout->tolerance = in->tolerance;
1372 inout->state = in->state;
1373
1374 }
1375 else {
1376
1377 if (in->gnum < inout->gnum) {
1378 inout->gnum = in->gnum;
1379 for (j = 0; j < 3; j++)
1380 inout->coord[j] = in->coord[j];
1381 inout->tolerance = in->tolerance;
1382 inout->state = in->state;
1383 }
1384
1385 }
1386
1387 } /* in.tol <= inout.tol */
1388
1389 } /* End of loop on array */
1390
1391 }
1392
1393 /*----------------------------------------------------------------------------
1394 * Create a function to define an operator for MPI reduction operation
1395 *
1396 * parameters:
1397 * in <-- input vertices
1398 * inout <-> in/out vertices (vertex with the max. toelrance)
1399 * len <-- size of input array
1400 * datatype <-- MPI_datatype associated to cs_join_vertex_t
1401 *---------------------------------------------------------------------------*/
1402
1403 static void
_mpi_vertex_max(cs_join_vertex_t * in,cs_join_vertex_t * inout,int * len,MPI_Datatype * datatype)1404 _mpi_vertex_max(cs_join_vertex_t *in,
1405 cs_join_vertex_t *inout,
1406 int *len,
1407 MPI_Datatype *datatype)
1408 {
1409 CS_UNUSED(datatype);
1410
1411 int i, j;
1412
1413 assert(in != NULL && inout != NULL);
1414
1415 for (i = 0; i < *len; i++) {
1416
1417 if (in->tolerance >= inout->tolerance) {
1418
1419 if (in->tolerance > inout->tolerance) {
1420
1421 inout->gnum = in->gnum;
1422 for (j = 0; j < 3; j++)
1423 inout->coord[j] = in->coord[j];
1424 inout->tolerance = in->tolerance;
1425 inout->state = in->state;
1426
1427 }
1428 else {
1429
1430 if (in->gnum < inout->gnum) {
1431 inout->gnum = in->gnum;
1432 for (j = 0; j < 3; j++)
1433 inout->coord[j] = in->coord[j];
1434 inout->tolerance = in->tolerance;
1435 inout->state = in->state;
1436
1437 }
1438
1439 }
1440
1441 } /* in.tol >= inout.tol */
1442
1443 } /* End of loop on array */
1444
1445 }
1446
1447 #endif /* HAVE_MPI */
1448
1449 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1450
1451 /*============================================================================
1452 * Public function definitions
1453 *===========================================================================*/
1454
1455 /*----------------------------------------------------------------------------
1456 * Allocate and initialize a new cs_join_mesh_t structure.
1457 *
1458 * parameters:
1459 * name <-- name of the mesh
1460 *
1461 * returns:
1462 * a pointer to a cs_join_mesh_t structure.
1463 *---------------------------------------------------------------------------*/
1464
1465 cs_join_mesh_t *
cs_join_mesh_create(const char * name)1466 cs_join_mesh_create(const char *name)
1467 {
1468 cs_join_mesh_t *new_mesh = NULL;
1469
1470 BFT_MALLOC(new_mesh, 1, cs_join_mesh_t);
1471
1472 if (name != NULL) {
1473
1474 int len = strlen(name);
1475
1476 BFT_MALLOC(new_mesh->name, len + 1, char);
1477 strncpy(new_mesh->name, name, len);
1478 new_mesh->name[len] = '\0';
1479
1480 }
1481 else
1482 new_mesh->name = NULL;
1483
1484 new_mesh->n_faces = 0;
1485 new_mesh->n_g_faces = 0;
1486 new_mesh->face_gnum = NULL;
1487 new_mesh->face_vtx_idx = NULL;
1488 new_mesh->face_vtx_lst = NULL;
1489 new_mesh->n_vertices = 0;
1490 new_mesh->n_g_vertices = 0;
1491 new_mesh->vertices = NULL;
1492
1493 return new_mesh;
1494 }
1495
1496 /*----------------------------------------------------------------------------
1497 * Get a cs_join_mesh_t structure with the given list of global faces inside.
1498 *
1499 * Exchange between ranks to get the connectivity associated to each
1500 * face of the global numbering list.
1501 *
1502 * parameters:
1503 * mesh_name <-- name of the created mesh
1504 * n_elts <-- number of elements in the global list
1505 * glob_sel <-- list of global elements (ordered)
1506 * gnum_rank_index <-- index on ranks for the global elements
1507 * local_mesh <-- pointer to the local part of the distributed
1508 * cs_join_mesh_t structure on selected elements
1509 *
1510 * returns:
1511 * a pointer to a new allocated cs_join_mesh_t structure
1512 *---------------------------------------------------------------------------*/
1513
1514 cs_join_mesh_t *
cs_join_mesh_create_from_glob_sel(const char * mesh_name,cs_lnum_t n_elts,const cs_gnum_t glob_sel[],const cs_gnum_t gnum_rank_index[],const cs_join_mesh_t * local_mesh)1515 cs_join_mesh_create_from_glob_sel(const char *mesh_name,
1516 cs_lnum_t n_elts,
1517 const cs_gnum_t glob_sel[],
1518 const cs_gnum_t gnum_rank_index[],
1519 const cs_join_mesh_t *local_mesh)
1520 {
1521 cs_join_mesh_t *new_mesh = NULL;
1522
1523 const int n_ranks = cs_glob_n_ranks;
1524
1525 if (n_ranks == 1) {
1526
1527 cs_lnum_t i;
1528 cs_lnum_t *loc_sel = NULL;
1529
1530 BFT_MALLOC(loc_sel, n_elts, cs_lnum_t);
1531
1532 for (i = 0; i < n_elts; i++)
1533 loc_sel[i] = glob_sel[i];
1534
1535 new_mesh = cs_join_mesh_create_from_subset(mesh_name,
1536 n_elts,
1537 loc_sel,
1538 local_mesh);
1539
1540 BFT_FREE(loc_sel);
1541 }
1542
1543 #if defined(HAVE_MPI)
1544
1545 else { /* Parallel mode */
1546
1547 int *send_rank = NULL;
1548 cs_lnum_t n_send_faces = 0;
1549 cs_lnum_t *send_faces = NULL;
1550
1551 new_mesh = cs_join_mesh_create(mesh_name);
1552
1553 /* Define a send list (face ids to send) from the global list of faces
1554 to receive. */
1555
1556 _get_send_faces(gnum_rank_index,
1557 n_elts,
1558 glob_sel,
1559 &n_send_faces,
1560 &send_rank,
1561 &send_faces);
1562
1563 /* Get useful connectivity on ranks for computing local intersections */
1564
1565 cs_join_mesh_exchange(n_send_faces,
1566 send_rank,
1567 send_faces,
1568 local_mesh,
1569 new_mesh,
1570 cs_glob_mpi_comm);
1571
1572 BFT_FREE(send_faces);
1573 BFT_FREE(send_rank);
1574
1575 cs_join_mesh_face_order(new_mesh);
1576 }
1577
1578 #endif
1579
1580 return new_mesh;
1581 }
1582
1583 /*----------------------------------------------------------------------------
1584 * Allocate and define a cs_join_mesh_t structure relative to an extraction
1585 * of selected faces.
1586 *
1587 * The selection must be ordered.
1588 *
1589 * parameters:
1590 * mesh_name <-- name of the name to create
1591 * subset_size <-- number of selected faces in the subset
1592 * selection <-- list of selected faces. Numbering in parent mesh
1593 * parent_mesh <-- parent cs_join_mesh_t structure
1594 *
1595 * returns:
1596 * a pointer to a cs_join_mesh_t structure
1597 *---------------------------------------------------------------------------*/
1598
1599 cs_join_mesh_t *
cs_join_mesh_create_from_subset(const char * mesh_name,cs_lnum_t subset_size,const cs_lnum_t selection[],const cs_join_mesh_t * parent_mesh)1600 cs_join_mesh_create_from_subset(const char *mesh_name,
1601 cs_lnum_t subset_size,
1602 const cs_lnum_t selection[],
1603 const cs_join_mesh_t *parent_mesh)
1604 {
1605 cs_lnum_t i, j, shift, parent_id, start, end;
1606
1607 cs_lnum_t n_select_vertices = 0;
1608 cs_lnum_t *select_vtx_id = NULL;
1609
1610 cs_join_mesh_t *mesh = NULL;
1611
1612 assert(parent_mesh != NULL);
1613
1614 /* Get the selected vertices relative to the subset selection */
1615
1616 BFT_MALLOC(select_vtx_id, parent_mesh->n_vertices, cs_lnum_t);
1617
1618 for (i = 0; i < parent_mesh->n_vertices; i++)
1619 select_vtx_id[i] = -1;
1620
1621 for (i = 0; i < subset_size; i++) {
1622
1623 parent_id = selection[i] - 1;
1624
1625 for (j = parent_mesh->face_vtx_idx[parent_id];
1626 j < parent_mesh->face_vtx_idx[parent_id+1];
1627 j++) {
1628 select_vtx_id[parent_mesh->face_vtx_lst[j]] = 0;
1629 }
1630
1631 }
1632
1633 n_select_vertices = 0;
1634 for (i = 0; i < parent_mesh->n_vertices; i++) {
1635 if (select_vtx_id[i] > -1)
1636 select_vtx_id[i] = n_select_vertices++;
1637 }
1638
1639 /* Create a new cs_join_mesh_t structure */
1640
1641 mesh = cs_join_mesh_create(mesh_name);
1642
1643 mesh->n_faces = subset_size;
1644
1645 /* Build face_vtx_idx, and face global numbering */
1646
1647 BFT_MALLOC(mesh->face_vtx_idx, mesh->n_faces + 1, cs_lnum_t);
1648 BFT_MALLOC(mesh->face_gnum, mesh->n_faces, cs_gnum_t);
1649
1650 for (i = 0; i < mesh->n_faces; i++) {
1651
1652 parent_id = selection[i] - 1;
1653
1654 mesh->face_vtx_idx[i+1] = parent_mesh->face_vtx_idx[parent_id+1]
1655 - parent_mesh->face_vtx_idx[parent_id];
1656 mesh->face_gnum[i] = parent_mesh->face_gnum[parent_id];
1657
1658 }
1659
1660 mesh->face_vtx_idx[0] = 0;
1661 for (i = 0; i < mesh->n_faces; i++)
1662 mesh->face_vtx_idx[i+1] += mesh->face_vtx_idx[i];
1663
1664 BFT_MALLOC(mesh->face_vtx_lst,
1665 mesh->face_vtx_idx[mesh->n_faces], cs_lnum_t);
1666
1667 /* Build face_vtx_lst */
1668
1669 for (i = 0; i < mesh->n_faces; i++) {
1670
1671 parent_id = selection[i] - 1;
1672 start = parent_mesh->face_vtx_idx[parent_id];
1673 end = parent_mesh->face_vtx_idx[parent_id+1];
1674 shift = mesh->face_vtx_idx[i];
1675
1676 for (j = start; j < end; j++) {
1677
1678 mesh->face_vtx_lst[shift + j - start]
1679 = select_vtx_id[parent_mesh->face_vtx_lst[j]];
1680
1681 }
1682
1683 } /* End of loop on selected faces */
1684
1685 /* Define vertices */
1686
1687 mesh->n_vertices = n_select_vertices;
1688
1689 BFT_MALLOC(mesh->vertices, n_select_vertices, cs_join_vertex_t);
1690
1691 n_select_vertices = 0;
1692 for (i = 0; i < parent_mesh->n_vertices; i++) {
1693 if (select_vtx_id[i] > -1)
1694 mesh->vertices[n_select_vertices++] = parent_mesh->vertices[i];
1695 }
1696
1697 /* Define global numbering */
1698
1699 if (cs_glob_n_ranks == 1) {
1700
1701 mesh->n_g_faces = mesh->n_faces;
1702 mesh->n_g_vertices = mesh->n_vertices;
1703
1704 }
1705 else {
1706
1707 fvm_io_num_t *io_num = NULL;
1708 cs_gnum_t *vtx_gnum = NULL;
1709
1710 const cs_gnum_t *io_gnum = NULL;
1711
1712 /* Get the global number of faces in the subset */
1713
1714 io_num = fvm_io_num_create(NULL, mesh->face_gnum, subset_size, 0);
1715
1716 mesh->n_g_faces = fvm_io_num_get_global_count(io_num);
1717
1718 io_num = fvm_io_num_destroy(io_num);
1719
1720 /* Get the global number of vertices in the subset */
1721
1722 BFT_MALLOC(vtx_gnum, mesh->n_vertices, cs_gnum_t);
1723
1724 for (i = 0; i < mesh->n_vertices; i++)
1725 vtx_gnum[i] = mesh->vertices[i].gnum;
1726
1727 io_num = fvm_io_num_create(NULL, vtx_gnum, mesh->n_vertices, 0);
1728
1729 mesh->n_g_vertices = fvm_io_num_get_global_count(io_num);
1730
1731 io_gnum = fvm_io_num_get_global_num(io_num);
1732
1733 for (i = 0; i < mesh->n_vertices; i++)
1734 mesh->vertices[i].gnum = io_gnum[i];
1735
1736 io_num = fvm_io_num_destroy(io_num);
1737
1738 BFT_FREE(vtx_gnum);
1739 }
1740
1741 /* Free memory */
1742
1743 BFT_FREE(select_vtx_id);
1744
1745 /* Order faces by increasing global number */
1746
1747 cs_join_mesh_face_order(mesh);
1748
1749 return mesh;
1750 }
1751
1752 /*----------------------------------------------------------------------------
1753 * Define a cs_join_mesh_t structure from a selection of faces and its
1754 * related vertices.
1755 *
1756 * parameters:
1757 * name <-- mesh name of the resulting cs_join_mesh_t structure
1758 * param <-- set of user-defined parameters for the joining
1759 * selection <-> selected entities
1760 * b_f2v_idx <-- border "face -> vertex" connectivity index
1761 * b_f2v_lst <-- border "face -> vertex" connectivity
1762 * i_f2v_idx <-- interior "face -> vertex" connectivity index
1763 * i_f2v_lst <-- interior "face -> vertex" connectivity
1764 * n_vertices <-- number of vertices in the parent mesh
1765 * vtx_coord <-- coordinates of vertices in parent mesh
1766 * vtx_gnum <-- global numbering of vertices
1767 *
1768 * returns:
1769 * a pointer to a cs_join_mesh_t structure
1770 *---------------------------------------------------------------------------*/
1771
1772 cs_join_mesh_t *
cs_join_mesh_create_from_select(const char * name,const cs_join_param_t param,cs_join_select_t * selection,const cs_lnum_t b_f2v_idx[],const cs_lnum_t b_f2v_lst[],const cs_lnum_t i_f2v_idx[],const cs_lnum_t i_f2v_lst[],const cs_lnum_t n_vertices,const cs_real_t vtx_coord[],const cs_gnum_t vtx_gnum[])1773 cs_join_mesh_create_from_select(const char *name,
1774 const cs_join_param_t param,
1775 cs_join_select_t *selection,
1776 const cs_lnum_t b_f2v_idx[],
1777 const cs_lnum_t b_f2v_lst[],
1778 const cs_lnum_t i_f2v_idx[],
1779 const cs_lnum_t i_f2v_lst[],
1780 const cs_lnum_t n_vertices,
1781 const cs_real_t vtx_coord[],
1782 const cs_gnum_t vtx_gnum[])
1783 {
1784 int i, j, id, face_id, start, end, shift;
1785
1786 cs_join_vertex_t *vertices = NULL;
1787 cs_join_mesh_t *mesh = NULL;
1788
1789 assert(selection != NULL);
1790
1791 mesh = cs_join_mesh_create(name);
1792
1793 /* Define face connectivity */
1794
1795 mesh->n_faces = selection->n_faces;
1796 mesh->n_g_faces = selection->n_g_faces;
1797
1798 /* Define face_vtx_idx */
1799
1800 BFT_MALLOC(mesh->face_vtx_idx, selection->n_faces + 1, cs_lnum_t);
1801
1802 for (i = 0; i < selection->n_faces; i++) {
1803 face_id = selection->faces[i] - 1;
1804 mesh->face_vtx_idx[i+1] = b_f2v_idx[face_id+1] - b_f2v_idx[face_id];
1805 }
1806
1807 mesh->face_vtx_idx[0] = 0;
1808 for (i = 0; i < selection->n_faces; i++)
1809 mesh->face_vtx_idx[i+1] += mesh->face_vtx_idx[i];
1810
1811 BFT_MALLOC(mesh->face_vtx_lst,
1812 mesh->face_vtx_idx[mesh->n_faces], cs_lnum_t);
1813
1814 /* Define face_vtx_lst */
1815
1816 for (i = 0; i < selection->n_faces; i++) {
1817
1818 shift = mesh->face_vtx_idx[i];
1819 face_id = selection->faces[i] - 1;
1820 start = b_f2v_idx[face_id];
1821 end = b_f2v_idx[face_id+1];
1822
1823 for (j = 0; j < end - start; j++) {
1824
1825 id = cs_search_binary(selection->n_vertices,
1826 b_f2v_lst[start + j] + 1,
1827 selection->vertices);
1828
1829 mesh->face_vtx_lst[shift + j] = id;
1830
1831 }
1832
1833 } /* End of loop on selected faces */
1834
1835 /* Define global face numbering */
1836
1837 BFT_MALLOC(mesh->face_gnum, mesh->n_faces, cs_gnum_t);
1838
1839 if (selection->compact_face_gnum == NULL)
1840 for (i = 0; i < selection->n_faces; i++)
1841 mesh->face_gnum[i] = selection->faces[i];
1842 else
1843 for (i = 0; i < selection->n_faces; i++)
1844 mesh->face_gnum[i] = selection->compact_face_gnum[i];
1845
1846 /* Define vertices */
1847
1848 vertices = _define_vertices(param,
1849 selection,
1850 b_f2v_idx,
1851 b_f2v_lst,
1852 i_f2v_idx,
1853 i_f2v_lst,
1854 n_vertices,
1855 vtx_coord,
1856 vtx_gnum);
1857
1858 mesh->n_vertices = selection->n_vertices;
1859 mesh->n_g_vertices = selection->n_g_vertices;
1860 mesh->vertices = vertices;
1861
1862 return mesh;
1863 }
1864
1865 /*----------------------------------------------------------------------------
1866 * Destroy a cs_join_mesh_t structure.
1867 *
1868 * parameters:
1869 * mesh <-> pointer to pointer to cs_join_mesh_t structure to destroy
1870 *---------------------------------------------------------------------------*/
1871
1872 void
cs_join_mesh_destroy(cs_join_mesh_t ** mesh)1873 cs_join_mesh_destroy(cs_join_mesh_t **mesh)
1874 {
1875 if (*mesh != NULL) {
1876
1877 cs_join_mesh_t *m = *mesh;
1878
1879 BFT_FREE(m->name);
1880 BFT_FREE(m->face_vtx_idx);
1881 BFT_FREE(m->face_vtx_lst);
1882 BFT_FREE(m->face_gnum);
1883 BFT_FREE(m->vertices);
1884 BFT_FREE(*mesh);
1885
1886 }
1887 }
1888
1889 /*----------------------------------------------------------------------------
1890 * Re-initialize an existing cs_join_mesh_t structure.
1891 *
1892 * parameters:
1893 * mesh <-> pointer to a cs_join_mesh_t structure
1894 *---------------------------------------------------------------------------*/
1895
1896 void
cs_join_mesh_reset(cs_join_mesh_t * mesh)1897 cs_join_mesh_reset(cs_join_mesh_t *mesh)
1898 {
1899 if (mesh == NULL)
1900 return;
1901
1902 mesh->n_faces = 0;
1903 mesh->n_g_faces = 0;
1904
1905 BFT_FREE(mesh->face_gnum);
1906 BFT_FREE(mesh->face_vtx_lst);
1907 BFT_FREE(mesh->face_vtx_idx);
1908
1909 mesh->n_vertices = 0;
1910 mesh->n_g_vertices = 0;
1911
1912 BFT_FREE(mesh->vertices);
1913 }
1914
1915 /*----------------------------------------------------------------------------
1916 * Copy a cs_join_mesh_t structure into another.
1917 *
1918 * parameters:
1919 * mesh <-> pointer to a cs_join_mesh_t structure to fill
1920 * ref_mesh <-- pointer to the reference
1921 *---------------------------------------------------------------------------*/
1922
1923 void
cs_join_mesh_copy(cs_join_mesh_t ** mesh,const cs_join_mesh_t * ref_mesh)1924 cs_join_mesh_copy(cs_join_mesh_t **mesh,
1925 const cs_join_mesh_t *ref_mesh)
1926 {
1927 cs_lnum_t i;
1928 cs_join_mesh_t *_mesh = *mesh;
1929
1930 if (ref_mesh == NULL) {
1931 cs_join_mesh_destroy(mesh);
1932 return;
1933 }
1934
1935 if (_mesh == NULL)
1936 _mesh = cs_join_mesh_create(ref_mesh->name);
1937
1938 _mesh->n_faces = ref_mesh->n_faces;
1939 _mesh->n_g_faces = ref_mesh->n_g_faces;
1940
1941 BFT_REALLOC(_mesh->face_gnum, _mesh->n_faces, cs_gnum_t);
1942 BFT_REALLOC(_mesh->face_vtx_idx, _mesh->n_faces + 1, cs_lnum_t);
1943
1944 _mesh->face_vtx_idx[0] = 0;
1945
1946 for (i = 0; i < _mesh->n_faces; i++) {
1947 _mesh->face_gnum[i] = ref_mesh->face_gnum[i];
1948 _mesh->face_vtx_idx[i+1] = ref_mesh->face_vtx_idx[i+1];
1949 }
1950
1951 BFT_REALLOC(_mesh->face_vtx_lst,
1952 _mesh->face_vtx_idx[_mesh->n_faces],
1953 cs_lnum_t);
1954
1955 for (i = 0; i < _mesh->face_vtx_idx[_mesh->n_faces]; i++)
1956 _mesh->face_vtx_lst[i] = ref_mesh->face_vtx_lst[i];
1957
1958 _mesh->n_vertices = ref_mesh->n_vertices;
1959 _mesh->n_g_vertices = ref_mesh->n_g_vertices;
1960
1961 BFT_REALLOC(_mesh->vertices, _mesh->n_vertices, cs_join_vertex_t);
1962
1963 memcpy(_mesh->vertices,
1964 ref_mesh->vertices,
1965 _mesh->n_vertices*sizeof(cs_join_vertex_t));
1966
1967 /* Set return pointer */
1968
1969 *mesh = _mesh;
1970 }
1971
1972 /*----------------------------------------------------------------------------
1973 * Compute the global min/max tolerance defined on vertices and display it
1974 *
1975 * parameters:
1976 * param <-- user-defined parameters for the joining algorithm
1977 * mesh <-- pointer to a cs_join_mesh_t structure
1978 *---------------------------------------------------------------------------*/
1979
1980 void
cs_join_mesh_minmax_tol(cs_join_param_t param,cs_join_mesh_t * mesh)1981 cs_join_mesh_minmax_tol(cs_join_param_t param,
1982 cs_join_mesh_t *mesh)
1983 {
1984 cs_lnum_t i;
1985 cs_join_vertex_t _min, _max, g_min, g_max;
1986
1987 assert(mesh != NULL);
1988
1989 const int n_ranks = cs_glob_n_ranks;
1990
1991 _min.state = CS_JOIN_STATE_UNDEF;
1992 _max.state = CS_JOIN_STATE_UNDEF;
1993 _min.gnum = 0;
1994 _max.gnum = 0;
1995 _min.tolerance = DBL_MAX;
1996 _max.tolerance = -DBL_MAX;
1997 for (i = 0; i < 3; i++) {
1998 _min.coord[i] = DBL_MAX;
1999 _max.coord[i] = DBL_MAX;
2000 }
2001
2002 g_min = _min;
2003 g_max = _max;
2004
2005 /* Compute local min/max */
2006
2007 if (mesh->n_vertices > 0) {
2008
2009 for (i = 0; i < mesh->n_vertices; i++) {
2010
2011 if (_min.tolerance > mesh->vertices[i].tolerance)
2012 _min = mesh->vertices[i];
2013 if (_max.tolerance < mesh->vertices[i].tolerance)
2014 _max = mesh->vertices[i];
2015
2016 }
2017
2018 if (param.verbosity > 3) {
2019 fprintf(cs_glob_join_log,
2020 "\n Local min/max. tolerance:\n\n"
2021 " Glob. Num. | Tolerance | Coordinates\n");
2022 cs_join_mesh_dump_vertex(cs_glob_join_log, _min);
2023 cs_join_mesh_dump_vertex(cs_glob_join_log, _max);
2024 }
2025
2026 }
2027
2028 #if defined(HAVE_MPI)
2029 #if !defined(_WIN32)
2030 if (n_ranks > 1) {
2031
2032 MPI_Datatype MPI_JOIN_VERTEX = _create_vtx_datatype();
2033 MPI_Op MPI_Vertex_min, MPI_Vertex_max;
2034
2035 MPI_Op_create((MPI_User_function *)_mpi_vertex_min,
2036 true, &MPI_Vertex_min);
2037 MPI_Op_create((MPI_User_function *)_mpi_vertex_max,
2038 false, &MPI_Vertex_max);
2039
2040 MPI_Allreduce(&_min, &g_min, 1, MPI_JOIN_VERTEX, MPI_Vertex_min,
2041 cs_glob_mpi_comm);
2042 MPI_Allreduce(&_max, &g_max, 1, MPI_JOIN_VERTEX, MPI_Vertex_max,
2043 cs_glob_mpi_comm);
2044
2045 bft_printf(_(" Global min/max. tolerance:\n\n"
2046 " Glob. Num. | Tolerance | Coordinates\n\n"));
2047 _log_vertex(g_min);
2048 _log_vertex(g_max);
2049
2050 MPI_Op_free(&MPI_Vertex_min);
2051 MPI_Op_free(&MPI_Vertex_max);
2052 MPI_Type_free(&MPI_JOIN_VERTEX);
2053
2054 }
2055 #endif
2056 #endif
2057
2058 }
2059
2060
2061 #if defined(HAVE_MPI)
2062
2063 /*----------------------------------------------------------------------------
2064 * Get the connectivity of a list of global elements distributed over the
2065 * ranks.
2066 *
2067 * parameters:
2068 * n_send <-- number of face/rank couples to send
2069 * send_rank <-- index on ranks for the face distribution
2070 * send_faces <-- list of face ids to send
2071 * send_mesh <-- pointer to the sending cs_join_mesh_t structure
2072 * recv_mesh <-> pointer to the receiving cs_join_mesh_t structure
2073 * comm <-- mpi communicator on which take places comm.
2074 *---------------------------------------------------------------------------*/
2075
2076 void
cs_join_mesh_exchange(cs_lnum_t n_send,const int send_rank[],const cs_lnum_t send_faces[],const cs_join_mesh_t * send_mesh,cs_join_mesh_t * recv_mesh,MPI_Comm comm)2077 cs_join_mesh_exchange(cs_lnum_t n_send,
2078 const int send_rank[],
2079 const cs_lnum_t send_faces[],
2080 const cs_join_mesh_t *send_mesh,
2081 cs_join_mesh_t *recv_mesh,
2082 MPI_Comm comm)
2083 {
2084 assert(send_mesh != NULL);
2085 assert(recv_mesh != NULL);
2086
2087 /* Try to use datatype larger than char if for join vertex type
2088 to reduce risk of reaching maximum index size with cs_lnum_t
2089 (with current size of 48 bytes, it would take more than
2090 40 million vertices in the connectvity to approach
2091 the 2 Gb limit, but a safety factor of 4 or 8 is still better).
2092 */
2093
2094 size_t vtx_t_size = sizeof(cs_join_vertex_t);
2095 cs_datatype_t vtx_type = CS_CHAR;
2096
2097 if (vtx_t_size % 8 == 0) {
2098 vtx_t_size /= 8;
2099 vtx_type = CS_UINT64;
2100 }
2101 else if (vtx_t_size % 4 == 0) {
2102 vtx_t_size /= 4;
2103 vtx_type = CS_UINT32;
2104 }
2105
2106 /* Sanity checks */
2107
2108 assert(send_mesh != NULL);
2109 assert(recv_mesh != NULL);
2110 assert(send_rank != NULL || n_send == 0);
2111
2112 cs_all_to_all_t *d
2113 = cs_all_to_all_create(n_send,
2114 CS_ALL_TO_ALL_ORDER_BY_SRC_RANK, /* needed ? */
2115 NULL, /* dest_id */
2116 send_rank,
2117 comm);
2118
2119 /* The mesh doesn't change from a global point of view.
2120 It's only a redistribution of the elements according to the send_faces
2121 list. */
2122
2123 recv_mesh->n_g_faces = send_mesh->n_g_faces;
2124 recv_mesh->n_g_vertices = send_mesh->n_g_vertices;
2125
2126 /* Exchange face connect. count */
2127
2128 cs_lnum_t *send_index;
2129 BFT_MALLOC(send_index, n_send+1, cs_lnum_t);
2130 send_index[0] = 0;
2131
2132 cs_gnum_t *send_gbuf;
2133 BFT_MALLOC(send_gbuf, n_send*2, cs_gnum_t);
2134
2135 for (cs_lnum_t i = 0; i < n_send; i++) {
2136
2137 cs_lnum_t face_id = send_faces[i];
2138 cs_lnum_t n_f_vtx = send_mesh->face_vtx_idx[face_id+1]
2139 - send_mesh->face_vtx_idx[face_id];
2140
2141 send_gbuf[i*2] = send_mesh->face_gnum[face_id];
2142 send_gbuf[i*2+1] = n_f_vtx;
2143
2144 send_index[i+1] = send_index[i] + n_f_vtx;
2145
2146 } /* End of loop on elements to send */
2147
2148 cs_gnum_t *recv_gbuf
2149 = cs_all_to_all_copy_array(d,
2150 CS_GNUM_TYPE,
2151 2,
2152 false, /* reverse */
2153 send_gbuf,
2154 NULL);
2155
2156 BFT_FREE(send_gbuf);
2157
2158 /* Update cs_join_mesh_t structure */
2159
2160 cs_lnum_t n_recv = cs_all_to_all_n_elts_dest(d);
2161
2162 recv_mesh->n_faces = n_recv;
2163
2164 BFT_MALLOC(recv_mesh->face_gnum, n_recv, cs_gnum_t);
2165 BFT_MALLOC(recv_mesh->face_vtx_idx, n_recv + 1, cs_lnum_t);
2166
2167 /* Scan recv_gbuf to build face->vertex connect. index */
2168
2169 recv_mesh->face_vtx_idx[0] = 0;
2170
2171 for (cs_lnum_t i = 0; i < n_recv; i++) {
2172
2173 recv_mesh->face_gnum[i] = recv_gbuf[i*2];
2174 recv_mesh->face_vtx_idx[i+1] = recv_mesh->face_vtx_idx[i]
2175 + recv_gbuf[i*2+1];
2176
2177 } /* End of loop on received elements */
2178
2179 BFT_FREE(recv_gbuf);
2180
2181 /* Store vtx_shift data into send_count in order to re-use it */
2182
2183 cs_join_vertex_t *send_vbuf;
2184 BFT_MALLOC(send_vbuf, send_index[n_send], cs_join_vertex_t);
2185
2186 /* Exchange vertex buffers */
2187
2188 for (cs_lnum_t i = 0; i < n_send; i++) {
2189
2190 cs_lnum_t face_id = send_faces[i];
2191
2192 cs_lnum_t s_id = send_mesh->face_vtx_idx[face_id];
2193 cs_lnum_t e_id = send_mesh->face_vtx_idx[face_id+1];
2194
2195 size_t shift = send_index[i];
2196
2197 for (cs_lnum_t j = s_id; j < e_id; j++) {
2198 cs_lnum_t vtx_id = send_mesh->face_vtx_lst[j];
2199 send_vbuf[shift++] = send_mesh->vertices[vtx_id];
2200 }
2201
2202 }
2203
2204 for (cs_lnum_t i = 0; i < n_send; i++)
2205 send_index[i+1] *= vtx_t_size;
2206 for (cs_lnum_t i = 0; i < n_recv; i++)
2207 recv_mesh->face_vtx_idx[i+1] *= vtx_t_size;
2208
2209 recv_mesh->vertices
2210 = cs_all_to_all_copy_indexed(d,
2211 vtx_type,
2212 false, /* reverse */
2213 send_index,
2214 send_vbuf,
2215 recv_mesh->face_vtx_idx,
2216 NULL);
2217
2218 for (cs_lnum_t i = 0; i < n_recv; i++)
2219 recv_mesh->face_vtx_idx[i+1] /= vtx_t_size;
2220
2221 BFT_FREE(send_vbuf);
2222 BFT_FREE(send_index);
2223
2224 /* Update cs_join_mesh_t structure */
2225
2226 recv_mesh->n_vertices = recv_mesh->face_vtx_idx[n_recv];
2227
2228 cs_lnum_t recv_lst_size = recv_mesh->face_vtx_idx[n_recv];
2229 BFT_MALLOC(recv_mesh->face_vtx_lst, recv_lst_size, cs_lnum_t);
2230
2231 for (cs_lnum_t i = 0; i < recv_lst_size; i++)
2232 recv_mesh->face_vtx_lst[i] = i;
2233
2234 /* Delete vertices which appear several times */
2235
2236 cs_join_mesh_vertex_clean(recv_mesh);
2237
2238 /* Free memory */
2239
2240 cs_all_to_all_destroy(&d);
2241 }
2242
2243 #endif /* HAVE_MPI */
2244
2245 /*----------------------------------------------------------------------------
2246 * Destroy a cs_join_edges_t structure.
2247 *
2248 * parameters:
2249 * edges <-> pointer to pointer to cs_join_edges_t structure to destroy
2250 *---------------------------------------------------------------------------*/
2251
2252 void
cs_join_mesh_destroy_edges(cs_join_edges_t ** edges)2253 cs_join_mesh_destroy_edges(cs_join_edges_t **edges)
2254 {
2255 if (*edges != NULL) {
2256
2257 cs_join_edges_t *e = *edges;
2258
2259 if (e->n_edges > 0) {
2260
2261 BFT_FREE(e->gnum);
2262 BFT_FREE(e->def);
2263 }
2264
2265 BFT_FREE(e->vtx_idx);
2266 BFT_FREE(e->adj_vtx_lst);
2267 BFT_FREE(e->edge_lst);
2268
2269 BFT_FREE(*edges);
2270 }
2271 }
2272
2273 /*----------------------------------------------------------------------------
2274 * Order a cs_join_mesh_t structure according to its global face numbering
2275 * Delete redundancies.
2276 *
2277 * parameters:
2278 * mesh <-> pointer to a cs_join_mesh_t structure to order
2279 *---------------------------------------------------------------------------*/
2280
2281 void
cs_join_mesh_face_order(cs_join_mesh_t * mesh)2282 cs_join_mesh_face_order(cs_join_mesh_t *mesh)
2283 {
2284 int i, j, o_id;
2285 cs_lnum_t shift, start, end, n_new_faces;
2286 cs_gnum_t prev, cur;
2287
2288 cs_lnum_t n_faces = mesh->n_faces;
2289 cs_lnum_t *num_buf = NULL, *selection = NULL;
2290 cs_lnum_t *order = NULL;
2291 cs_gnum_t *gnum_buf = NULL;
2292
2293 assert(mesh != NULL);
2294
2295 if (n_faces == 0)
2296 return;
2297
2298 /* Order faces according to their global numbering */
2299
2300 BFT_MALLOC(order, n_faces, cs_lnum_t);
2301
2302 cs_order_gnum_allocated(NULL, mesh->face_gnum, order, n_faces);
2303
2304 /* Order global face numbering */
2305
2306 BFT_MALLOC(gnum_buf, n_faces, cs_gnum_t);
2307 BFT_MALLOC(selection, n_faces, cs_lnum_t);
2308
2309 for (i = 0; i < n_faces; i++)
2310 gnum_buf[i] = mesh->face_gnum[i];
2311
2312 prev = 0;
2313 n_new_faces = 0;
2314
2315 for (i = 0; i < n_faces; i++) {
2316
2317 o_id = order[i];
2318 cur = gnum_buf[o_id];
2319
2320 if (prev != cur) {
2321 prev = cur;
2322 selection[n_new_faces] = o_id;
2323 mesh->face_gnum[n_new_faces] = cur;
2324 n_new_faces++;
2325 }
2326
2327 }
2328
2329 mesh->n_faces = n_new_faces;
2330
2331 BFT_FREE(gnum_buf);
2332 BFT_FREE(order);
2333
2334 BFT_REALLOC(mesh->face_gnum, n_new_faces, cs_gnum_t);
2335 BFT_REALLOC(selection, n_new_faces, cs_lnum_t);
2336
2337 /* Order face -> vertex connectivity list */
2338
2339 BFT_MALLOC(num_buf, mesh->face_vtx_idx[n_faces], cs_lnum_t);
2340
2341 for (i = 0; i < mesh->face_vtx_idx[n_faces]; i++)
2342 num_buf[i] = mesh->face_vtx_lst[i];
2343
2344 shift = 0;
2345
2346 for (i = 0; i < n_new_faces; i++) {
2347
2348 o_id = selection[i];
2349 start = mesh->face_vtx_idx[o_id];
2350 end = mesh->face_vtx_idx[o_id+1];
2351
2352 for (j = start; j < end; j++)
2353 mesh->face_vtx_lst[shift++] = num_buf[j];
2354
2355 } /* End of loop on faces */
2356
2357 BFT_REALLOC(num_buf, n_faces, cs_lnum_t);
2358
2359 for (i = 0; i < n_faces; i++)
2360 num_buf[i] = mesh->face_vtx_idx[i+1] - mesh->face_vtx_idx[i];
2361
2362 for (i = 0; i < n_new_faces; i++) {
2363 o_id = selection[i];
2364 mesh->face_vtx_idx[i+1] = mesh->face_vtx_idx[i] + num_buf[o_id];
2365 }
2366
2367 /* Memory management */
2368
2369 BFT_FREE(selection);
2370 BFT_FREE(num_buf);
2371 BFT_REALLOC(mesh->face_vtx_idx, n_new_faces+1, cs_lnum_t);
2372 BFT_REALLOC(mesh->face_vtx_lst, mesh->face_vtx_idx[n_new_faces], cs_lnum_t);
2373 }
2374
2375 #if defined(HAVE_MPI)
2376
2377 /*----------------------------------------------------------------------------
2378 * Synchronize vertices definition over the rank. For a vertex with the same
2379 * global number but a different tolerance, we keep the minimal tolerance.
2380 *
2381 * parameters:
2382 * mesh <-> pointer to the cs_join_mesh_t structure to synchronize
2383 *---------------------------------------------------------------------------*/
2384
2385 void
cs_join_mesh_sync_vertices(cs_join_mesh_t * mesh)2386 cs_join_mesh_sync_vertices(cs_join_mesh_t *mesh)
2387 {
2388 cs_lnum_t *order = NULL;
2389 cs_gnum_t *recv_gnum = NULL;
2390
2391 MPI_Comm mpi_comm = cs_glob_mpi_comm;
2392
2393 const int n_ranks = cs_glob_n_ranks;
2394 const int local_rank = CS_MAX(cs_glob_rank_id, 0);
2395
2396 assert(n_ranks > 1);
2397 assert(mesh != NULL);
2398
2399 /* Get the max global number */
2400
2401 cs_gnum_t l_max_gnum = 0, g_max_gnum = 0;
2402 for (cs_lnum_t i = 0; i < mesh->n_vertices; i++)
2403 l_max_gnum = CS_MAX(l_max_gnum, mesh->vertices[i].gnum);
2404
2405 MPI_Allreduce(&l_max_gnum, &g_max_gnum, 1, CS_MPI_GNUM, MPI_MAX, mpi_comm);
2406
2407 cs_block_dist_info_t bi = cs_block_dist_compute_sizes(local_rank,
2408 n_ranks,
2409 1,
2410 0,
2411 g_max_gnum);
2412
2413 int *dest_rank = NULL;
2414 BFT_MALLOC(dest_rank, mesh->n_vertices, int);
2415
2416 for (cs_lnum_t i = 0; i < mesh->n_vertices; i++)
2417 dest_rank[i] = ((mesh->vertices[i].gnum - 1)/bi.block_size) * bi.rank_step;
2418
2419 cs_all_to_all_t
2420 *d = cs_all_to_all_create(mesh->n_vertices,
2421 0, /* flags */
2422 NULL,
2423 dest_rank,
2424 mpi_comm);
2425
2426 cs_all_to_all_transfer_dest_rank(d, &dest_rank);
2427
2428 /* Send vertices to sync and receive its part of work */
2429
2430 int stride = sizeof(cs_join_vertex_t);
2431
2432 cs_join_vertex_t *recv_vertices = cs_all_to_all_copy_array(d,
2433 CS_CHAR,
2434 stride,
2435 false, /* reverse */
2436 mesh->vertices,
2437 NULL);
2438
2439 cs_lnum_t n_recv = cs_all_to_all_n_elts_dest(d);
2440
2441 /* Order vertices by increasing global number */
2442
2443 BFT_MALLOC(recv_gnum, n_recv, cs_gnum_t);
2444 BFT_MALLOC(order, n_recv, cs_lnum_t);
2445
2446 for (cs_lnum_t i = 0; i < n_recv; i++)
2447 recv_gnum[i] = recv_vertices[i].gnum;
2448
2449 cs_order_gnum_allocated(NULL, recv_gnum, order, n_recv);
2450
2451 /* Sync. vertices sharing the same global number */
2452
2453 cs_lnum_t start = 0, end = 0;
2454
2455 while (start < n_recv) {
2456
2457 cs_lnum_t i;
2458 double min_tol = recv_vertices[order[start]].tolerance;
2459 cs_gnum_t ref_gnum = recv_vertices[order[start]].gnum;
2460 for (i = start;
2461 i < n_recv && ref_gnum == recv_vertices[order[i]].gnum;
2462 i++);
2463 end = i;
2464
2465 /* Get min tolerance */
2466 for (i = start; i < end; i++)
2467 min_tol = CS_MIN(min_tol, recv_vertices[order[i]].tolerance);
2468
2469 /* Set min tolerance to all vertices sharing the same global number */
2470 for (i = start; i < end; i++)
2471 recv_vertices[order[i]].tolerance = min_tol;
2472
2473 start = end;
2474 }
2475
2476 /* Send back vertices after synchronization */
2477
2478 cs_all_to_all_copy_array(d,
2479 CS_CHAR,
2480 stride,
2481 true, /* reverse */
2482 recv_vertices,
2483 mesh->vertices);
2484
2485 /* Free buffers */
2486
2487 BFT_FREE(recv_gnum);
2488 BFT_FREE(order);
2489 BFT_FREE(recv_vertices);
2490
2491 cs_all_to_all_destroy(&d);
2492 }
2493
2494 #endif /* HAVE_MPI */
2495
2496 /*----------------------------------------------------------------------------
2497 * Delete vertices which appear several times (same global number) and
2498 * vertices which are not used in face definition.
2499 *
2500 * parameters:
2501 * mesh <-> pointer to cs_join_mesh_t structure to clean
2502 *---------------------------------------------------------------------------*/
2503
2504 void
cs_join_mesh_vertex_clean(cs_join_mesh_t * mesh)2505 cs_join_mesh_vertex_clean(cs_join_mesh_t *mesh)
2506 {
2507 cs_lnum_t i, j, shift, n_init_vertices, n_final_vertices;
2508 cs_gnum_t prev, cur;
2509
2510 cs_lnum_t *order = NULL;
2511 cs_lnum_t *init2final = NULL, *tag = NULL;
2512 cs_gnum_t *gnum_buf = NULL;
2513 cs_join_vertex_t *final_vertices = NULL;
2514
2515 assert(mesh != NULL);
2516
2517 n_init_vertices = mesh->n_vertices;
2518
2519 if (n_init_vertices < 2)
2520 return;
2521
2522 /* Count the final number of vertices */
2523
2524 BFT_MALLOC(order, n_init_vertices, cs_lnum_t);
2525 BFT_MALLOC(tag, n_init_vertices, cs_lnum_t);
2526 BFT_MALLOC(gnum_buf, n_init_vertices, cs_gnum_t);
2527
2528 for (i = 0; i < n_init_vertices; i++) {
2529 gnum_buf[i] = mesh->vertices[i].gnum;
2530 tag[i] = 0;
2531 }
2532
2533 /* Tag vertices really used in the mesh definition */
2534
2535 for (i = 0; i < mesh->n_faces; i++)
2536 for (j = mesh->face_vtx_idx[i]; j < mesh->face_vtx_idx[i+1]; j++)
2537 tag[mesh->face_vtx_lst[j]] = 1;
2538
2539 /* Order vertices by increasing global number */
2540
2541 cs_order_gnum_allocated(NULL, gnum_buf, order, n_init_vertices);
2542
2543 n_final_vertices = 0;
2544 prev = 0;
2545
2546 for (i = 0; i < n_init_vertices; i++) {
2547
2548 shift = order[i];
2549 cur = gnum_buf[shift];
2550
2551 if (prev != cur && tag[i] > 0) {
2552 n_final_vertices++;
2553 prev = cur;
2554 }
2555
2556 }
2557
2558 /* Define the final vertices structure and indirection between
2559 initial numbering and final numbering */
2560
2561 BFT_MALLOC(final_vertices, n_final_vertices, cs_join_vertex_t);
2562 BFT_MALLOC(init2final, n_init_vertices, cs_lnum_t);
2563
2564 n_final_vertices = 0;
2565 prev = 0;
2566
2567 for (i = 0; i < n_init_vertices; i++) {
2568
2569 shift = order[i];
2570 cur = gnum_buf[shift];
2571
2572 if (prev != cur && tag[i] > 0) {
2573
2574 final_vertices[n_final_vertices++] = mesh->vertices[shift];
2575 prev = cur;
2576
2577 }
2578
2579 init2final[shift] = n_final_vertices - 1;
2580
2581 }
2582
2583 BFT_FREE(mesh->vertices);
2584
2585 mesh->vertices = final_vertices;
2586 mesh->n_vertices = n_final_vertices;
2587
2588 /* Update face->vertex connectivity list */
2589
2590 for (i = 0; i < mesh->n_faces; i++) {
2591
2592 for (j = mesh->face_vtx_idx[i]; j < mesh->face_vtx_idx[i+1]; j++)
2593 mesh->face_vtx_lst[j] = init2final[mesh->face_vtx_lst[j]];
2594
2595 } /* end of loop on faces */
2596
2597 BFT_FREE(init2final);
2598 BFT_FREE(gnum_buf);
2599 BFT_FREE(tag);
2600 BFT_FREE(order);
2601 }
2602
2603 /*----------------------------------------------------------------------------
2604 * Clean the given cs_join_mesh_t structure, removing degenerate edges.
2605 *
2606 * parameters:
2607 * mesh <-> pointer to the cs_join_mesh_t structure to clean
2608 * verbosity <-- level of display
2609 *---------------------------------------------------------------------------*/
2610
2611 void
cs_join_mesh_clean(cs_join_mesh_t * mesh,int verbosity)2612 cs_join_mesh_clean(cs_join_mesh_t *mesh,
2613 int verbosity)
2614 {
2615 assert(mesh != NULL);
2616
2617 /* Delete empty edge:
2618 These edges are generated during the merge step. If two vertices
2619 sharing the same edge are fused, we have to delete the resulting
2620 edge.
2621 */
2622
2623 _remove_empty_edges(mesh, verbosity);
2624
2625 /* Delete degenerate edge:
2626 These edges are generated during the merge step.
2627
2628
2629 x x
2630 |\ |
2631 | \ |
2632 a2| \a3 A2|
2633 | \ Merge |
2634 | \ a4 of |
2635 ---s1----s2------ vertices x
2636 | \ s1 and s2 / \
2637 | \ / \
2638 a1| \a4 A1/ \A3
2639 | \ / \
2640 | \ / \
2641 x-----------x x-----------x
2642 a5 A4
2643
2644 */
2645
2646 _remove_degenerate_edges(mesh, verbosity);
2647 }
2648
2649 /*----------------------------------------------------------------------------
2650 * Define a list of edges associated to a cs_join_mesh_t structure.
2651 *
2652 * parameters:
2653 * mesh <-- pointer to a cs_join_mesh_t structure
2654 *
2655 * returns:
2656 * a pointer to the new defined cs_join_edges_t structure.
2657 *---------------------------------------------------------------------------*/
2658
2659 cs_join_edges_t *
cs_join_mesh_define_edges(const cs_join_mesh_t * mesh)2660 cs_join_mesh_define_edges(const cs_join_mesh_t *mesh)
2661 {
2662 int i, j;
2663 cs_lnum_t v1_num, v2_num, o_id1, o_id2;
2664 cs_lnum_t edge_shift, shift, n_init_edges;
2665 cs_gnum_t v1_gnum, v2_gnum;
2666
2667 cs_lnum_t *order = NULL;
2668 cs_lnum_t *vtx_counter = NULL, *vtx_lst = NULL;
2669 cs_gnum_t *adjacency = NULL;
2670 cs_join_edges_t *edges = NULL;
2671
2672 if (mesh == NULL)
2673 return edges;
2674
2675 /* Initialization and structure allocation */
2676
2677 BFT_MALLOC(edges, 1, cs_join_edges_t);
2678
2679 edges->n_edges = 0;
2680 edges->def = NULL;
2681 edges->gnum = NULL;
2682 edges->n_vertices = mesh->n_vertices;
2683 edges->vtx_idx = NULL;
2684 edges->adj_vtx_lst = NULL;
2685 edges->edge_lst = NULL;
2686
2687 /* Define edges */
2688
2689 n_init_edges = mesh->face_vtx_idx[mesh->n_faces];
2690
2691 BFT_MALLOC(edges->def, 2*n_init_edges, cs_lnum_t);
2692 BFT_MALLOC(edges->vtx_idx, mesh->n_vertices + 1, cs_lnum_t);
2693
2694 for (i = 0; i < mesh->n_vertices + 1; i++)
2695 edges->vtx_idx[i] = 0;
2696
2697 /* Loop on faces to initialize edge list */
2698
2699 BFT_MALLOC(vtx_lst, 2*n_init_edges, cs_lnum_t);
2700 BFT_MALLOC(adjacency, 2*n_init_edges, cs_gnum_t);
2701
2702 for (shift = 0, i = 0; i < mesh->n_faces; i++) {
2703
2704 cs_lnum_t start = mesh->face_vtx_idx[i];
2705 cs_lnum_t end = mesh->face_vtx_idx[i+1];
2706
2707 assert(end-start > 0);
2708
2709 for (j = start; j < end - 1; j++) {
2710
2711 v1_num = mesh->face_vtx_lst[j] + 1;
2712 v1_gnum = (mesh->vertices[v1_num-1]).gnum;
2713 v2_num = mesh->face_vtx_lst[j+1] + 1;
2714 v2_gnum = (mesh->vertices[v2_num-1]).gnum;
2715
2716 if (v1_gnum > v2_gnum) {
2717
2718 vtx_lst[2*shift] = v2_num;
2719 adjacency[2*shift] = v2_gnum;
2720 vtx_lst[2*shift+1] = v1_num;
2721 adjacency[2*shift+1] = v1_gnum;
2722
2723 }
2724 else {
2725
2726 vtx_lst[2*shift] = v1_num;
2727 adjacency[2*shift] = v1_gnum;
2728 vtx_lst[2*shift+1] = v2_num;
2729 adjacency[2*shift+1] = v2_gnum;
2730
2731 }
2732
2733 shift++;
2734
2735 } /* End of loop on n-1 first vertices */
2736
2737 v1_num = mesh->face_vtx_lst[end-1] + 1;
2738 v1_gnum = (mesh->vertices[v1_num-1]).gnum;
2739 v2_num = mesh->face_vtx_lst[start] + 1;
2740 v2_gnum = (mesh->vertices[v2_num-1]).gnum;
2741
2742 if (v1_gnum > v2_gnum) {
2743
2744 vtx_lst[2*shift] = v2_num;
2745 adjacency[2*shift] = v2_gnum;
2746 vtx_lst[2*shift+1] = v1_num;
2747 adjacency[2*shift+1] = v1_gnum;
2748
2749 }
2750 else {
2751
2752 vtx_lst[2*shift] = v1_num;
2753 adjacency[2*shift] = v1_gnum;
2754 vtx_lst[2*shift+1] = v2_num;
2755 adjacency[2*shift+1] = v2_gnum;
2756
2757 }
2758
2759 shift++;
2760
2761 } /* End of loop on faces */
2762
2763 assert(shift == n_init_edges);
2764
2765 BFT_MALLOC(order, n_init_edges, cs_lnum_t);
2766
2767 cs_order_gnum_allocated_s(NULL, adjacency, 2, order, n_init_edges);
2768
2769 if (n_init_edges > 0) {
2770
2771 /* Fill cs_join_edges_t structure */
2772
2773 o_id1 = order[0];
2774 edges->def[0] = vtx_lst[2*o_id1];
2775 edges->def[1] = vtx_lst[2*o_id1+1];
2776 edges->vtx_idx[vtx_lst[2*o_id1]] += 1;
2777 edges->vtx_idx[vtx_lst[2*o_id1+1]] += 1;
2778 edge_shift = 1;
2779
2780 for (i = 1; i < n_init_edges; i++) {
2781
2782 o_id1 = order[i-1];
2783 o_id2 = order[i];
2784
2785 if ( vtx_lst[2*o_id1] != vtx_lst[2*o_id2]
2786 || vtx_lst[2*o_id1+1] != vtx_lst[2*o_id2+1]) {
2787
2788 edges->vtx_idx[vtx_lst[2*o_id2]] += 1;
2789 edges->vtx_idx[vtx_lst[2*o_id2+1]] += 1;
2790 edges->def[2*edge_shift] = vtx_lst[2*o_id2];
2791 edges->def[2*edge_shift+1] = vtx_lst[2*o_id2+1];
2792 edge_shift++;
2793
2794 }
2795
2796 } /* End of loop on edges */
2797
2798 edges->n_edges = edge_shift;
2799 BFT_REALLOC(edges->def, 2*edges->n_edges, cs_lnum_t);
2800
2801 } /* If n_init_edges > 0 */
2802
2803 /* Build adj_vtx_lst and edge_lst */
2804
2805 BFT_MALLOC(vtx_counter, mesh->n_vertices, cs_lnum_t);
2806
2807 for (i = 0; i < mesh->n_vertices; i++) {
2808 edges->vtx_idx[i+1] += edges->vtx_idx[i];
2809 vtx_counter[i] = 0;
2810 }
2811
2812 BFT_MALLOC(edges->adj_vtx_lst, edges->vtx_idx[mesh->n_vertices], cs_lnum_t);
2813 BFT_MALLOC(edges->edge_lst, edges->vtx_idx[mesh->n_vertices], cs_lnum_t);
2814
2815 if (n_init_edges > 0) {
2816
2817 cs_lnum_t vtx_id_a, vtx_id_b, shift_a, shift_b;
2818 cs_gnum_t vtx_gnum_a, vtx_gnum_b;
2819
2820 cs_lnum_t cur_edge_num = 1;
2821
2822 /* Initiate edge_lst and adj_vtx_lst building */
2823
2824 o_id1 = order[0];
2825
2826 vtx_id_a = vtx_lst[2*o_id1]-1;
2827 vtx_id_b = vtx_lst[2*o_id1+1]-1;
2828
2829 vtx_gnum_a = (mesh->vertices[vtx_id_a]).gnum;
2830 vtx_gnum_b = (mesh->vertices[vtx_id_a]).gnum;
2831
2832 shift_a = edges->vtx_idx[vtx_id_a];
2833 shift_b = edges->vtx_idx[vtx_id_b];
2834
2835 vtx_counter[vtx_id_a] += 1;
2836 vtx_counter[vtx_id_b] += 1;
2837
2838 edges->adj_vtx_lst[shift_a] = vtx_id_b;
2839 edges->adj_vtx_lst[shift_b] = vtx_id_a;
2840
2841 if (vtx_gnum_a > vtx_gnum_b) {
2842 edges->edge_lst[shift_a] = -cur_edge_num;
2843 edges->edge_lst[shift_b] = cur_edge_num;
2844 }
2845 else {
2846 edges->edge_lst[shift_a] = cur_edge_num;
2847 edges->edge_lst[shift_b] = -cur_edge_num;
2848 }
2849
2850 cur_edge_num++;
2851
2852 for (i = 1; i < n_init_edges; i++) {
2853
2854 o_id1 = order[i-1];
2855 o_id2 = order[i];
2856
2857 if ( vtx_lst[2*o_id1] != vtx_lst[2*o_id2]
2858 || vtx_lst[2*o_id1+1] != vtx_lst[2*o_id2+1]) {
2859
2860 vtx_id_a = vtx_lst[2*o_id2]-1;
2861 vtx_id_b = vtx_lst[2*o_id2+1]-1;
2862
2863 vtx_gnum_a = (mesh->vertices[vtx_id_a]).gnum;
2864 vtx_gnum_b = (mesh->vertices[vtx_id_a]).gnum;
2865
2866 shift_a = edges->vtx_idx[vtx_id_a] + vtx_counter[vtx_id_a];
2867 shift_b = edges->vtx_idx[vtx_id_b] + vtx_counter[vtx_id_b];
2868
2869 vtx_counter[vtx_id_a] += 1;
2870 vtx_counter[vtx_id_b] += 1;
2871
2872 edges->adj_vtx_lst[shift_a] = vtx_id_b;
2873 edges->adj_vtx_lst[shift_b] = vtx_id_a;
2874
2875 if (vtx_gnum_a > vtx_gnum_b) {
2876 edges->edge_lst[shift_a] = -cur_edge_num;
2877 edges->edge_lst[shift_b] = cur_edge_num;
2878 }
2879 else {
2880 edges->edge_lst[shift_a] = cur_edge_num;
2881 edges->edge_lst[shift_b] = -cur_edge_num;
2882 }
2883
2884 cur_edge_num++;
2885
2886 }
2887
2888 } /* End of loop on edges */
2889
2890 assert(cur_edge_num - 1 == edges->n_edges);
2891
2892 } /* End of adj_vtx_lst and edge_lst building if n_init_edges > 0 */
2893
2894 /* Partial clean-up */
2895
2896 BFT_FREE(vtx_lst);
2897 BFT_FREE(vtx_counter);
2898
2899 /* Define a global numbering on edges */
2900
2901 BFT_MALLOC(edges->gnum, edges->n_edges, cs_gnum_t);
2902 BFT_REALLOC(adjacency, 2*edges->n_edges, cs_gnum_t);
2903
2904 for (i = 0; i < edges->n_edges; i++) {
2905
2906 cs_lnum_t v1_id = edges->def[2*i] - 1;
2907 cs_lnum_t v2_id = edges->def[2*i+1] - 1;
2908
2909 v1_gnum = (mesh->vertices[v1_id]).gnum;
2910 v2_gnum = (mesh->vertices[v2_id]).gnum;
2911
2912 if (v1_gnum > v2_gnum) {
2913 adjacency[2*i] = v2_gnum;
2914 adjacency[2*i+1] = v1_gnum;
2915 }
2916 else {
2917 adjacency[2*i] = v1_gnum;
2918 adjacency[2*i+1] = v2_gnum;
2919 }
2920
2921 } /* End of loop on edges */
2922
2923 /* Order vtx_lst and build an order list into adjacency */
2924
2925 cs_order_gnum_allocated_s(NULL, adjacency, 2, order, edges->n_edges);
2926
2927 if (cs_glob_n_ranks == 1) { /* Serial treatment */
2928
2929 edges->n_g_edges = edges->n_edges;
2930
2931 for (i = 0; i < edges->n_edges; i++) {
2932
2933 cs_lnum_t o_id = order[i];
2934
2935 edges->gnum[i] = o_id+1;
2936
2937 }
2938
2939 }
2940 else { /* Parallel treatment */
2941
2942 cs_gnum_t *order_couples = NULL;
2943 fvm_io_num_t *edge_io_num = NULL;
2944 const cs_gnum_t *edges_gnum = NULL;
2945
2946 assert(cs_glob_n_ranks > 1);
2947
2948 BFT_MALLOC(order_couples, 2*edges->n_edges, cs_gnum_t);
2949
2950 for (i = 0; i < edges->n_edges; i++) {
2951
2952 cs_lnum_t o_id = order[i];
2953
2954 order_couples[2*i] = adjacency[2*o_id];
2955 order_couples[2*i+1] = adjacency[2*o_id+1];
2956
2957 } /* End of loop on edges */
2958
2959 edge_io_num = fvm_io_num_create_from_adj_s(NULL,
2960 order_couples,
2961 edges->n_edges,
2962 2);
2963
2964 edges->n_g_edges = fvm_io_num_get_global_count(edge_io_num);
2965 edges_gnum = fvm_io_num_get_global_num(edge_io_num);
2966
2967 for (i = 0; i < edges->n_edges; i++)
2968 edges->gnum[i] = edges_gnum[i];
2969
2970 /* Partial Clean-up */
2971
2972 BFT_FREE(order_couples);
2973 fvm_io_num_destroy(edge_io_num);
2974
2975 }
2976
2977 /* Memory management */
2978
2979 BFT_FREE(adjacency);
2980 BFT_FREE(order);
2981
2982 /* Return pointers */
2983
2984 return edges;
2985 }
2986
2987 /*----------------------------------------------------------------------------
2988 * Get the edge number relative to a couple of vertex numbers.
2989 *
2990 * edge_num > 0 if couple is in the same order as the edge->def
2991 * edge_num < 0 otherwise
2992 *
2993 * parameters:
2994 * v1_num <-- vertex number for the first vertex
2995 * v2_num <-- vertex number for the second vertex
2996 * edges <-- pointer to a cs_join_edges_t structure
2997 *
2998 * returns:
2999 * an edge number relative to the couple of vertices
3000 *---------------------------------------------------------------------------*/
3001
3002 cs_lnum_t
cs_join_mesh_get_edge(cs_lnum_t v1_num,cs_lnum_t v2_num,const cs_join_edges_t * edges)3003 cs_join_mesh_get_edge(cs_lnum_t v1_num,
3004 cs_lnum_t v2_num,
3005 const cs_join_edges_t *edges)
3006 {
3007 cs_lnum_t i;
3008 cs_lnum_t edge_num = 0;
3009
3010 assert(edges != NULL);
3011 assert(v1_num > 0);
3012 assert(v2_num > 0);
3013
3014 if (edges->vtx_idx[v1_num] - edges->vtx_idx[v1_num-1] == 0)
3015 bft_error(__FILE__, __LINE__, 0,
3016 _(" The given vertex number: %ld is not defined"
3017 " in the edge structure (edges->vtx_idx).\n"), (long)v1_num);
3018
3019 for (i = edges->vtx_idx[v1_num-1]; i < edges->vtx_idx[v1_num]; i++) {
3020 if (edges->adj_vtx_lst[i] == v2_num - 1) {
3021 edge_num = edges->edge_lst[i];
3022 break;
3023 }
3024 }
3025
3026 if (edge_num == 0)
3027 bft_error(__FILE__, __LINE__, 0,
3028 _(" The given couple of vertex numbers :\n"
3029 " vertex 1 : %ld\n"
3030 " vertex 2 : %ld\n"
3031 " is not defined in the edge structure.\n"),
3032 (long)v1_num, (long)v2_num);
3033
3034 assert(edge_num != 0);
3035
3036 return edge_num;
3037 }
3038
3039 /*----------------------------------------------------------------------------
3040 * Re-organize the cs_join_mesh_t structure after a renumbering of
3041 * the vertices following the merge operation + a new description of each
3042 * face.
3043 *
3044 * parameters:
3045 * mesh <-> pointer to the cs_join_mesh_t structure to update
3046 * edges <-- pointer to a cs_join_edges_t structure
3047 * edge_index <-- index on edges for the new vertices
3048 * edge_new_vtx_lst <-- list of new vertices for each edge
3049 * n_new_vertices <-- new local number of vertices after merge
3050 * old2new <-- array storing the relation between old/new vertex id
3051 *---------------------------------------------------------------------------*/
3052
3053 void
cs_join_mesh_update(cs_join_mesh_t * mesh,const cs_join_edges_t * edges,const cs_lnum_t edge_index[],const cs_lnum_t edge_new_vtx_lst[],cs_lnum_t n_new_vertices,const cs_lnum_t old2new[])3054 cs_join_mesh_update(cs_join_mesh_t *mesh,
3055 const cs_join_edges_t *edges,
3056 const cs_lnum_t edge_index[],
3057 const cs_lnum_t edge_new_vtx_lst[],
3058 cs_lnum_t n_new_vertices,
3059 const cs_lnum_t old2new[])
3060 {
3061 cs_lnum_t i, j, n_adds;
3062
3063 cs_join_vertex_t *new_vertices = NULL;
3064 cs_lnum_t *new_face_vtx_idx = NULL, *new_face_vtx_lst = NULL;
3065
3066 /* Sanity checks */
3067
3068 assert(mesh != NULL);
3069 assert(edges != NULL);
3070
3071 /* Update description and numbering for face -> vertex connectivity */
3072
3073 if (edge_new_vtx_lst != NULL) {
3074
3075 BFT_MALLOC(new_face_vtx_idx, mesh->n_faces + 1, cs_lnum_t);
3076
3077 for (i = 0; i < mesh->n_faces + 1; i++)
3078 new_face_vtx_idx[i] = 0;
3079
3080 /* Update face -> vertex connectivity.
3081 add new vertices between the existing one */
3082
3083 for (i = 0; i < mesh->n_faces; i++) {
3084
3085 cs_lnum_t start_id = mesh->face_vtx_idx[i];
3086 cs_lnum_t end_id = mesh->face_vtx_idx[i+1];
3087
3088 for (j = start_id; j < end_id - 1; j++) {
3089
3090 n_adds = _count_new_added_vtx_to_edge(mesh->face_vtx_lst[j],
3091 mesh->face_vtx_lst[j+1],
3092 old2new,
3093 edges,
3094 edge_index,
3095 edge_new_vtx_lst);
3096
3097 new_face_vtx_idx[i+1] += n_adds;
3098
3099 }
3100
3101 /* Case end - start */
3102
3103 n_adds = _count_new_added_vtx_to_edge(mesh->face_vtx_lst[end_id-1],
3104 mesh->face_vtx_lst[start_id],
3105 old2new,
3106 edges,
3107 edge_index,
3108 edge_new_vtx_lst);
3109
3110 new_face_vtx_idx[i+1] += n_adds;
3111
3112 } /* End of loop on faces */
3113
3114 /* Build new face_vtx_idx */
3115
3116 new_face_vtx_idx[0] = 0;
3117 for (i = 0; i < mesh->n_faces; i++) {
3118
3119 new_face_vtx_idx[i+1] += new_face_vtx_idx[i];
3120
3121 if (new_face_vtx_idx[i+1] < 3)
3122 bft_error(__FILE__, __LINE__, 0,
3123 _(" Problem in mesh connectivity."
3124 " Face: %llu\n"
3125 " Problem detected during connectivity update:\n"
3126 " The face is defined by less than 3 points"
3127 " (excessive merging has occured).\n\n"
3128 " Modify joining parameters to reduce merging"
3129 " (fraction & merge).\n"),
3130 (unsigned long long)(mesh->face_gnum[i]));
3131
3132 }
3133
3134 /* Build new_face_vtx_lst */
3135
3136 BFT_MALLOC(new_face_vtx_lst, new_face_vtx_idx[mesh->n_faces], cs_lnum_t);
3137
3138 } /* End if edge_new_vtx_lst != NULL */
3139
3140 else { /* edge_new_vtx_lst == NULL */
3141
3142 new_face_vtx_idx = mesh->face_vtx_idx;
3143 new_face_vtx_lst = mesh->face_vtx_lst;
3144
3145 }
3146
3147 for (i = 0; i < mesh->n_faces; i++) {
3148
3149 cs_lnum_t start_id = mesh->face_vtx_idx[i];
3150 cs_lnum_t end_id = mesh->face_vtx_idx[i+1];
3151 cs_lnum_t shift = new_face_vtx_idx[i];
3152
3153 for (j = start_id; j < end_id-1; j++)
3154 _add_new_vtx_to_edge(mesh->face_vtx_lst[j],
3155 mesh->face_vtx_lst[j+1],
3156 old2new,
3157 edges,
3158 edge_index,
3159 edge_new_vtx_lst,
3160 new_face_vtx_lst,
3161 &shift);
3162
3163 /* Case end - start */
3164
3165 _add_new_vtx_to_edge(mesh->face_vtx_lst[end_id-1],
3166 mesh->face_vtx_lst[start_id],
3167 old2new,
3168 edges,
3169 edge_index,
3170 edge_new_vtx_lst,
3171 new_face_vtx_lst,
3172 &shift);
3173
3174 } /* End of loop on faces */
3175
3176 if (edge_new_vtx_lst != NULL) {
3177
3178 BFT_FREE(mesh->face_vtx_idx);
3179 BFT_FREE(mesh->face_vtx_lst);
3180
3181 mesh->face_vtx_idx = new_face_vtx_idx;
3182 mesh->face_vtx_lst = new_face_vtx_lst;
3183
3184 }
3185
3186 /* Define the new_vertices structure */
3187
3188 BFT_MALLOC(new_vertices, n_new_vertices, cs_join_vertex_t);
3189
3190 for (i = 0; i < mesh->n_vertices; i++)
3191 new_vertices[old2new[i]] = mesh->vertices[i];
3192
3193 #if 0 && defined(DEBUG) && !defined(NDEBUG)
3194 bft_printf("\n\n Dump Old2New array : "
3195 "n_old_vertices = %d - n_new_vertices = %d\n",
3196 mesh->n_vertices, n_new_vertices);
3197 for (i = 0; i < mesh->n_vertices; i++)
3198 bft_printf("Old num : %7d (%9u) => New num : %7d (%9u)\n",
3199 i+1 , mesh->vertices[i].gnum,
3200 old2new[i]+1, new_vertices[old2new[i]].gnum);
3201 bft_printf_flush();
3202 #endif
3203
3204 /* Update mesh structure */
3205
3206 BFT_FREE(mesh->vertices);
3207
3208 mesh->n_vertices = n_new_vertices;
3209 mesh->n_g_vertices = n_new_vertices;
3210 mesh->vertices = new_vertices;
3211
3212 #if defined(HAVE_MPI)
3213
3214 if (cs_glob_n_ranks > 1) {
3215
3216 cs_gnum_t *vtx_gnum = NULL;
3217 fvm_io_num_t *io_num = NULL;
3218
3219 /* Global number of selected vertices and associated
3220 fvm_io_num_t structure */
3221
3222 BFT_MALLOC(vtx_gnum, n_new_vertices, cs_gnum_t);
3223
3224 for (i = 0; i < n_new_vertices; i++)
3225 vtx_gnum[i] = (mesh->vertices[i]).gnum;
3226
3227 io_num = fvm_io_num_create(NULL, vtx_gnum, n_new_vertices, 0);
3228
3229 mesh->n_g_vertices = fvm_io_num_get_global_count(io_num);
3230
3231 fvm_io_num_destroy(io_num);
3232
3233 BFT_FREE(vtx_gnum);
3234
3235 }
3236 #endif
3237 }
3238
3239 /*----------------------------------------------------------------------------
3240 * Compute for each face of the cs_join_mesh_t structure the face normal.
3241 * || face_normal || = 1 (divided by the area of the face)
3242 *
3243 * The caller is responsible for freeing the returned array.
3244 *
3245 * parameters:
3246 * mesh <-- pointer to a cs_join_mesh_t structure
3247 *
3248 * Pi+1
3249 * *---------* B : barycenter of the polygon
3250 * / . . \
3251 * / . . \ Pi : vertices of the polygon
3252 * / . . \
3253 * / . . Ti \ Ti : triangle
3254 * *.........B.........* Pi
3255 * Pn-1 \ . . /
3256 * \ . . /
3257 * \ . . /
3258 * \ . T0 . /
3259 * *---------*
3260 * P0
3261 *
3262 *
3263 * returns:
3264 * an array with the face normal for each face of the mesh
3265 *---------------------------------------------------------------------------*/
3266
3267 cs_real_t *
cs_join_mesh_get_face_normal(const cs_join_mesh_t * mesh)3268 cs_join_mesh_get_face_normal(const cs_join_mesh_t *mesh)
3269 {
3270 cs_lnum_t i, j, k, vid;
3271 double inv_norm;
3272
3273 cs_lnum_t n_max_vertices = 0;
3274 cs_real_t *face_vtx_coord = NULL;
3275 cs_real_t *face_normal = NULL;
3276
3277 if (mesh == NULL)
3278 return face_normal;
3279
3280 if (mesh->n_faces == 0)
3281 return face_normal;
3282
3283 BFT_MALLOC(face_normal, 3*mesh->n_faces, cs_real_t);
3284
3285 for (i = 0; i < 3*mesh->n_faces; i++)
3286 face_normal[i] = 0.0;
3287
3288 /* Compute n_max_vertices */
3289
3290 for (i = 0; i < mesh->n_faces; i++)
3291 n_max_vertices = CS_MAX(n_max_vertices,
3292 mesh->face_vtx_idx[i+1] - mesh->face_vtx_idx[i]);
3293
3294 BFT_MALLOC(face_vtx_coord, 3*(n_max_vertices+1), cs_real_t);
3295
3296 for (i = 0; i < mesh->n_faces; i++) {
3297
3298 cs_real_t v1[3], v2[3], tri_normal[3];
3299
3300 cs_lnum_t shift = 0;
3301 cs_lnum_t s = mesh->face_vtx_idx[i];
3302 cs_lnum_t e = mesh->face_vtx_idx[i+1];
3303 cs_lnum_t n_face_vertices = e - s;
3304 double inv_n_face_vertices = 1/(double)n_face_vertices;
3305
3306 cs_real_t bary[3] = { 0.0, 0.0, 0.0};
3307 cs_real_t fnorm[3] = { 0.0, 0.0, 0.0};
3308
3309 /* Fill face_vtx_coord */
3310
3311 for (j = s; j < e; j++) {
3312 vid = mesh->face_vtx_lst[j];
3313 for (k = 0; k < 3; k++)
3314 face_vtx_coord[shift++] = mesh->vertices[vid].coord[k];
3315 }
3316
3317 vid = mesh->face_vtx_lst[s];
3318 for (k = 0; k < 3; k++)
3319 face_vtx_coord[shift++] = mesh->vertices[vid].coord[k];
3320
3321 /* Compute the barycenter of the face */
3322
3323 for (j = 0; j < n_face_vertices; j++)
3324 for (k = 0; k < 3; k++)
3325 bary[k] += face_vtx_coord[3*j+k];
3326
3327 for (k = 0; k < 3; k++)
3328 bary[k] *= inv_n_face_vertices;
3329
3330 /* Loop on the triangles of the face defined by an edge of the face
3331 and the barycenter */
3332
3333 for (j = 0; j < n_face_vertices; j++) {
3334
3335 /* Computation of the normal of each triangle Ti:
3336 -> --> -->
3337 N(Ti) = 1/2 ( BPi X BPi+1 )
3338 */
3339
3340 for (k = 0; k < 3; k++) {
3341 v1[k] = face_vtx_coord[3*j + k] - bary[k];
3342 v2[k] = face_vtx_coord[3*(j+1)+ k] - bary[k];
3343 }
3344
3345 _cross_product(v1, v2, tri_normal);
3346
3347 /* Computation of the normal of the polygon
3348 => vectorial sum of normals of each triangle
3349
3350 -> n-1 ->
3351 N(P) = Sum ( N(Ti) )
3352 i=0
3353 */
3354
3355 for (k = 0; k < 3; k++)
3356 fnorm[k] += 0.5 * tri_normal[k];
3357
3358 } /* End of loop on vertices of the face */
3359
3360 inv_norm = 1/sqrt(_dot_product(fnorm, fnorm));
3361
3362 for (k = 0; k < 3; k++)
3363 face_normal[3*i+k] = inv_norm * fnorm[k];
3364
3365 #if 0 && defined(DEBUG) && !defined(NDEBUG)
3366 bft_printf(" Face_num: %5d (%u)- face_normal [%8.4e, %8.4e, %8.4e]\n",
3367 i+1, mesh->face_gnum[i],
3368 face_normal[3*i], face_normal[3*i+1], face_normal[3*i+2]);
3369 bft_printf_flush();
3370 #endif
3371
3372 } /* End of loop on faces */
3373
3374 /* Free memory */
3375
3376 BFT_FREE(face_vtx_coord);
3377
3378 return face_normal;
3379 }
3380
3381 /*----------------------------------------------------------------------------
3382 * Allocate and define an "edge -> face" connectivity
3383 *
3384 * parameters:
3385 * mesh <-- pointer to a cs_join_mesh_t structure
3386 * edges <-- pointer to a cs_join_edges_t structure
3387 * edge_face_idx --> pointer to the edge -> face connect. index
3388 * edge_face_lst --> pointer to the edge -> face connect. list
3389 *---------------------------------------------------------------------------*/
3390
3391 void
cs_join_mesh_get_edge_face_adj(const cs_join_mesh_t * mesh,const cs_join_edges_t * edges,cs_lnum_t * edge_face_idx[],cs_lnum_t * edge_face_lst[])3392 cs_join_mesh_get_edge_face_adj(const cs_join_mesh_t *mesh,
3393 const cs_join_edges_t *edges,
3394 cs_lnum_t *edge_face_idx[],
3395 cs_lnum_t *edge_face_lst[])
3396 {
3397 cs_lnum_t i, j, k, edge_id, shift;
3398 cs_lnum_t n_edges, n_faces;
3399
3400 cs_lnum_t n_max_vertices = 0;
3401 cs_lnum_t *counter = NULL, *face_connect = NULL;
3402 cs_lnum_t *_edge_face_idx = NULL, *_edge_face_lst = NULL;
3403
3404 if (mesh == NULL || edges == NULL)
3405 return;
3406
3407 n_edges = edges->n_edges;
3408 n_faces = mesh->n_faces;
3409
3410 /* Compute n_max_vertices */
3411
3412 for (i = 0; i < n_faces; i++)
3413 n_max_vertices = CS_MAX(n_max_vertices,
3414 mesh->face_vtx_idx[i+1]-mesh->face_vtx_idx[i]);
3415
3416 BFT_MALLOC(face_connect, n_max_vertices + 1, cs_lnum_t);
3417 BFT_MALLOC(counter, n_edges, cs_lnum_t);
3418
3419 /* Build an edge -> face connectivity */
3420
3421 BFT_MALLOC(_edge_face_idx, n_edges+1, cs_lnum_t);
3422
3423 for (i = 0; i < n_edges+1; i++)
3424 _edge_face_idx[i] = 0;
3425
3426 for (i = 0; i < n_edges; i++)
3427 counter[i] = 0;
3428
3429 /* Build index */
3430
3431 for (i = 0; i < n_faces; i++) {
3432
3433 cs_lnum_t start_id = mesh->face_vtx_idx[i];
3434 cs_lnum_t end_id = mesh->face_vtx_idx[i+1];
3435 cs_lnum_t n_face_vertices = end_id - start_id;
3436
3437 for (j = start_id, k = 0; j < end_id; j++, k++)
3438 face_connect[k] = mesh->face_vtx_lst[j];
3439 face_connect[n_face_vertices] = mesh->face_vtx_lst[start_id];
3440
3441 assert(n_face_vertices == k);
3442
3443 for (j = 0; j < n_face_vertices; j++) {
3444
3445 cs_lnum_t vtx_id1 = face_connect[j];
3446
3447 for (k = edges->vtx_idx[vtx_id1]; k < edges->vtx_idx[vtx_id1+1]; k++)
3448 if (edges->adj_vtx_lst[k] == face_connect[j+1])
3449 break;
3450
3451 assert(k != edges->vtx_idx[vtx_id1+1]);
3452
3453 _edge_face_idx[CS_ABS(edges->edge_lst[k])] += 1;
3454
3455 } /* End of loop on vertices of the face */
3456
3457 } /* End of loop on faces */
3458
3459 for (i = 0; i < n_edges; i++)
3460 _edge_face_idx[i+1] += _edge_face_idx[i];
3461
3462 BFT_MALLOC(_edge_face_lst, _edge_face_idx[n_edges], cs_lnum_t);
3463
3464 /* Fill "edge -> face" connectivity list */
3465
3466 for (i = 0; i < n_faces; i++) {
3467
3468 cs_lnum_t start_id = mesh->face_vtx_idx[i];
3469 cs_lnum_t end_id = mesh->face_vtx_idx[i+1];
3470 cs_lnum_t n_face_vertices = end_id - start_id;
3471
3472 for (j = start_id, k = 0; j < end_id; j++, k++)
3473 face_connect[k] = mesh->face_vtx_lst[j];
3474 face_connect[n_face_vertices] = mesh->face_vtx_lst[start_id];
3475
3476 for (j = 0; j < n_face_vertices; j++) {
3477
3478 cs_lnum_t vtx_id1 = face_connect[j];
3479
3480 for (k = edges->vtx_idx[vtx_id1];
3481 k < edges->vtx_idx[vtx_id1+1]; k++)
3482 if (edges->adj_vtx_lst[k] == face_connect[j+1])
3483 break;
3484
3485 edge_id = CS_ABS(edges->edge_lst[k]) - 1;
3486 shift = _edge_face_idx[edge_id] + counter[edge_id];
3487 _edge_face_lst[shift] = i+1;
3488 counter[edge_id] += 1;
3489
3490 } /* End of loop on vertices of the face */
3491
3492 } /* End of loop on faces */
3493
3494 #if 0 && defined(DEBUG) && !defined(NDEBUG)
3495 bft_printf("\n DUMP EDGE -> FACE CONNECTIVITY:\n\n");
3496
3497 for (i = 0; i < n_edges; i++) {
3498
3499 cs_lnum_t start = _edge_face_idx[i];
3500 cs_lnum_t end = _edge_face_idx[i+1];
3501 cs_lnum_t v1_id = edges->def[2*i] - 1;
3502 cs_lnum_t v2_id = edges->def[2*i+1] - 1;
3503
3504 bft_printf(" edge_num: %6d (%9u) [%9u - %9u]: size: %4d, faces: ",
3505 i+1, edges->gnum[i],
3506 mesh->vertices[v1_id].gnum, mesh->vertices[v2_id].gnum,
3507 end-start);
3508 for (j = start; j < end; j++)
3509 bft_printf(" %u ", mesh->face_gnum[_edge_face_lst[j]-1]);
3510 bft_printf("\n");
3511 bft_printf_flush();
3512 }
3513 bft_printf("\n");
3514 #endif
3515
3516 /* Set return pointers */
3517
3518 *edge_face_idx = _edge_face_idx;
3519 *edge_face_lst = _edge_face_lst;
3520
3521 /* Free memory*/
3522
3523 BFT_FREE(counter);
3524 BFT_FREE(face_connect);
3525 }
3526
3527 /*----------------------------------------------------------------------------
3528 * Dump a cs_join_vertex_t structure into a file.
3529 *
3530 * parameters:
3531 * f <-- handle to output file
3532 * vertex <-- cs_join_vertex_t structure to dump
3533 *---------------------------------------------------------------------------*/
3534
3535 void
cs_join_mesh_dump_vertex(FILE * f,const cs_join_vertex_t vertex)3536 cs_join_mesh_dump_vertex(FILE *f,
3537 const cs_join_vertex_t vertex)
3538 {
3539 assert(vertex.gnum > 0);
3540 assert(vertex.tolerance >= 0.0);
3541
3542 fprintf(f," %10llu | %11.6f | % 12.10e % 12.10e % 12.10e | %s\n",
3543 (unsigned long long)vertex.gnum, vertex.tolerance,
3544 vertex.coord[0], vertex.coord[1], vertex.coord[2],
3545 _print_state( vertex.state));
3546 }
3547
3548 /*----------------------------------------------------------------------------
3549 * Dump a cs_join_mesh_t structure into a file.
3550 *
3551 * parameters:
3552 * f <-- handle to output file
3553 * mesh <-- pointer to cs_join_mesh_t structure to dump
3554 *---------------------------------------------------------------------------*/
3555
3556 void
cs_join_mesh_dump(FILE * f,const cs_join_mesh_t * mesh)3557 cs_join_mesh_dump(FILE *f,
3558 const cs_join_mesh_t *mesh)
3559 {
3560 int i, j;
3561
3562 if (mesh == NULL) {
3563 fprintf(f,
3564 "\n\n -- Dump a cs_join_mesh_t structure: (%p) --\n",
3565 (const void *)mesh);
3566 return;
3567 }
3568
3569 fprintf(f, "\n\n -- Dump a cs_join_mesh_t structure: %s (%p) --\n",
3570 mesh->name, (const void *)mesh);
3571 fprintf(f, "\n mesh->n_faces: %11ld\n", (long)mesh->n_faces);
3572 fprintf(f, " mesh->n_g_faces: %11llu\n\n",
3573 (unsigned long long)mesh->n_g_faces);
3574
3575 if (mesh->face_vtx_idx != NULL) {
3576
3577 for (i = 0; i < mesh->n_faces; i++) {
3578
3579 cs_lnum_t start = mesh->face_vtx_idx[i];
3580 cs_lnum_t end = mesh->face_vtx_idx[i+1];
3581
3582 fprintf(f, "\n face_id: %9ld gnum: %10llu n_vertices : %4ld\n",
3583 (long)i, (unsigned long long)mesh->face_gnum[i],
3584 (long)(end-start));
3585
3586 for (j = start; j < end; j++) {
3587
3588 cs_lnum_t vtx_id = mesh->face_vtx_lst[j];
3589 cs_join_vertex_t v_data = mesh->vertices[vtx_id];
3590
3591 fprintf(f," %8ld - %10llu - [ % 7.5e % 7.5e % 7.5e] - %s\n",
3592 (long)vtx_id+1, (unsigned long long)v_data.gnum,
3593 v_data.coord[0], v_data.coord[1], v_data.coord[2],
3594 _print_state(v_data.state));
3595
3596 }
3597 fprintf(f,"\n");
3598
3599 /* Check if there is no incoherency in the mesh definition */
3600
3601 for (j = start; j < end - 1; j++) {
3602
3603 cs_lnum_t vtx_id1 = mesh->face_vtx_lst[j];
3604 cs_lnum_t vtx_id2 = mesh->face_vtx_lst[j+1];
3605
3606 if (vtx_id1 == vtx_id2) {
3607 fprintf(f,
3608 " Incoherency found in the current mesh definition\n"
3609 " Face number: %ld (global: %llu)\n"
3610 " Vertices: local (%ld, %ld), global (%llu, %llu)"
3611 " are defined twice\n",
3612 (long)i+1, (unsigned long long)mesh->face_gnum[i],
3613 (long)vtx_id1+1, (long)vtx_id2+1,
3614 (unsigned long long)(mesh->vertices[vtx_id1]).gnum,
3615 (unsigned long long)(mesh->vertices[vtx_id2]).gnum);
3616 fflush(f);
3617 assert(0);
3618 }
3619
3620 }
3621
3622 {
3623 cs_lnum_t vtx_id1 = mesh->face_vtx_lst[end-1];
3624 cs_lnum_t vtx_id2 = mesh->face_vtx_lst[start];
3625
3626 if (vtx_id1 == vtx_id2) {
3627 fprintf(f,
3628 " Incoherency found in the current mesh definition\n"
3629 " Face number: %ld (global: %llu)\n"
3630 " Vertices: local (%ld, %ld), global (%llu, %llu)"
3631 " are defined twice\n",
3632 (long)i+1, (unsigned long long)(mesh->face_gnum[i]),
3633 (long)vtx_id1+1, (long)vtx_id2+1,
3634 (unsigned long long)((mesh->vertices[vtx_id1]).gnum),
3635 (unsigned long long)((mesh->vertices[vtx_id2]).gnum));
3636 fflush(f);
3637 assert(0);
3638 }
3639 }
3640
3641 } /* End of loop on faces */
3642
3643 } /* End if face_vtx_idx != NULL */
3644
3645 fprintf(f,
3646 "\n Dump vertex data\n"
3647 " mesh->vertices : %p\n"
3648 " mesh->n_vertices : %11ld\n"
3649 " mesh->n_g_vertices : %11llu\n\n",
3650 (const void *)mesh->vertices, (long)mesh->n_vertices,
3651 (unsigned long long)mesh->n_g_vertices);
3652
3653 if (mesh->n_vertices > 0) {
3654
3655 fprintf(f,
3656 " Local Num | Global Num | Tolerance | Coordinates\n\n");
3657
3658 for (i = 0; i < mesh->n_vertices; i++) {
3659 fprintf(f," %9d |", i+1);
3660 cs_join_mesh_dump_vertex(f, mesh->vertices[i]);
3661 }
3662
3663 }
3664 fprintf(f,"\n");
3665 fflush(f);
3666 }
3667
3668 /*----------------------------------------------------------------------------
3669 * Dump a list of cs_join_edge_t structures.
3670 *
3671 * parameters:
3672 * f <-- handle to output file
3673 * edges <-- cs_join_edges_t structure to dump
3674 * mesh <-- associated cs_join_mesh_t structure
3675 *---------------------------------------------------------------------------*/
3676
3677 void
cs_join_mesh_dump_edges(FILE * f,const cs_join_edges_t * edges,const cs_join_mesh_t * mesh)3678 cs_join_mesh_dump_edges(FILE *f,
3679 const cs_join_edges_t *edges,
3680 const cs_join_mesh_t *mesh)
3681 {
3682 cs_lnum_t i, j;
3683
3684 if (edges == NULL)
3685 return;
3686
3687 fprintf(f, "\n Edge connectivity used in the joining operation:\n");
3688 fprintf(f, " Number of edges: %8ld\n", (long)edges->n_edges);
3689 fprintf(f, " Number of vertices: %8ld\n", (long)edges->n_vertices);
3690
3691 for (i = 0; i < edges->n_edges; i++) { /* Dump edge connectivity */
3692
3693 cs_lnum_t v1_id = edges->def[2*i] - 1;
3694 cs_lnum_t v2_id = edges->def[2*i+1] - 1;
3695 cs_gnum_t v1_gnum = (mesh->vertices[v1_id]).gnum;
3696 cs_gnum_t v2_gnum = (mesh->vertices[v2_id]).gnum;
3697
3698 fprintf(f, " Edge %6ld (%8llu) <Vertex> [%8llu %8llu]\n",
3699 (long)i+1, (unsigned long long)edges->gnum[i],
3700 (unsigned long long)v1_gnum,
3701 (unsigned long long)v2_gnum);
3702
3703 /* Check coherency */
3704
3705 if (v1_id == v2_id) {
3706 fprintf(f, " Incoherency found in the current edge definition\n"
3707 " Edge number: %ld\n"
3708 " Vertices: local (%ld, %ld), global (%llu, %llu)"
3709 " are defined twice\n",
3710 (long)i+1, (long)v1_id+1, (long)v2_id+1,
3711 (unsigned long long)v1_gnum, (unsigned long long)v2_gnum);
3712 fflush(f);
3713 assert(0);
3714 }
3715
3716 if (v1_gnum == v2_gnum) {
3717 fprintf(f, " Incoherency found in the current edge definition\n"
3718 " Edge number: %ld\n"
3719 " Vertices: local (%ld, %ld), global (%llu, %llu)"
3720 " are defined twice\n",
3721 (long)i+1, (long)v1_id+1, (long)v2_id+1,
3722 (unsigned long long)v1_gnum, (unsigned long long)v2_gnum);
3723 fflush(f);
3724 assert(0);
3725 }
3726
3727 } /* End of loop on edges */
3728
3729 fprintf(f, "\n Vertex -> Vertex connectivity :\n\n");
3730
3731 for (i = 0; i < mesh->n_vertices; i++) {
3732
3733 cs_lnum_t start = edges->vtx_idx[i];
3734 cs_lnum_t end = edges->vtx_idx[i+1];
3735
3736 fprintf(f, " Vertex %6ld (%7llu) - %3d - ",
3737 (long)i+1, (unsigned long long)(mesh->vertices[i]).gnum,
3738 (int)(end - start));
3739
3740 for (j = start; j < end; j++) {
3741 if (edges->edge_lst[j] > 0)
3742 fprintf
3743 (f, " [ v: %7llu, e: %7llu] ",
3744 (unsigned long long)(mesh->vertices[edges->adj_vtx_lst[j]]).gnum,
3745 (unsigned long long)edges->gnum[edges->edge_lst[j] - 1]);
3746 else
3747 fprintf
3748 (f, " [ v: %7llu, e: %7llu] ",
3749 (unsigned long long)(mesh->vertices[edges->adj_vtx_lst[j]]).gnum,
3750 (unsigned long long)edges->gnum[- edges->edge_lst[j] - 1]);
3751 }
3752 fprintf(f, "\n");
3753
3754 }
3755
3756 fflush(f);
3757 }
3758
3759 /*---------------------------------------------------------------------------*/
3760
3761 END_C_DECLS
3762