1 /*============================================================================
2 * Fortran interfaces of functions needing a synchronization of the extended
3 * neighborhood.
4 *============================================================================*/
5
6 /*
7 This file is part of Code_Saturne, a general-purpose CFD tool.
8
9 Copyright (C) 1998-2021 EDF S.A.
10
11 This program is free software; you can redistribute it and/or modify it under
12 the terms of the GNU General Public License as published by the Free Software
13 Foundation; either version 2 of the License, or (at your option) any later
14 version.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
19 details.
20
21 You should have received a copy of the GNU General Public License along with
22 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
23 Street, Fifth Floor, Boston, MA 02110-1301, USA.
24 */
25
26 /*----------------------------------------------------------------------------*/
27
28 #include "cs_defs.h"
29
30 /*----------------------------------------------------------------------------
31 * Standard C library headers
32 *----------------------------------------------------------------------------*/
33
34 #include <assert.h>
35 #include <errno.h>
36 #include <stdio.h>
37 #include <stdarg.h>
38 #include <string.h>
39 #include <math.h>
40 #include <float.h>
41
42 #if defined(HAVE_MPI)
43 #include <mpi.h>
44 #endif
45
46 /*----------------------------------------------------------------------------
47 * Local headers
48 *----------------------------------------------------------------------------*/
49
50 #include "bft_error.h"
51 #include "bft_mem.h"
52 #include "bft_printf.h"
53
54 #include "cs_halo.h"
55 #include "cs_halo_perio.h"
56 #include "cs_math.h"
57 #include "cs_mesh.h"
58 #include "cs_mesh_adjacencies.h"
59 #include "cs_mesh_quantities.h"
60 #include "cs_parall.h"
61 #include "cs_sort.h"
62
63 /*----------------------------------------------------------------------------
64 * Header for the current file
65 *----------------------------------------------------------------------------*/
66
67 #include "cs_ext_neighborhood.h"
68
69 /*----------------------------------------------------------------------------*/
70
71 BEGIN_C_DECLS
72
73 /*=============================================================================
74 * Additional doxygen documentation
75 *============================================================================*/
76
77 /*!
78 \file cs_ext_neighborhood.c
79 Extended cell neighborhood.
80
81 \enum cs_ext_neighborhood_type_t
82
83 \brief Type of extended neighborhood associated with the mesh
84
85 Extended cell neighborhoods are based on cells adjacent through a shared
86 vertex but with no shared face.
87
88 \var CS_EXT_NEIGHBORHOOD_NONE
89 No extended neighborhood
90
91 \image html ext_neighborhood_none.svg "No extended neighborhood"
92
93 \var CS_EXT_NEIGHBORHOOD_COMPLETE
94 Full extended neighborhood
95
96 \image html ext_neighborhood_complete.svg "Full extended neighborhood"
97
98 This option should lead to the smoothest gradient, as it uses information
99 from all neighbors, but is quite costly. On average, a hexahedral mesh
100 has about 21 vertex neighbors per cell, a tetrahedral mesh
101 around 150 vertex neighbors.
102
103 \var CS_EXT_NEIGHBORHOOD_CELL_CENTER_OPPOSITE
104 Cells with centers best aligned opposite to face-adjacent cell centers
105
106 \image html ext_neighborhood_cell_center_opposite.svg "Opposite cell centers"
107
108 Add cells whose centers are closest to the half-line prolonging
109 the [face-adjacent cell-center, cell center segment]. The number of
110 additional cells in the extended neighborhood is at most equal to the
111 number of cell faces.
112
113 \var CS_EXT_NEIGHBORHOOD_NON_ORTHO_MAX
114 Cells adjacent to faces whose non-orthogonality exceeds
115 a given threshold (45 degrees by default)
116
117 \image html ext_neighborhood_non_ortho_max.svg "Maximum non-orthogonality"
118
119 Add cells adjacent to vertices of faces whose non-orthogonality exceeds
120 a given threshold.
121
122 This is the legacy reduced extended neighborhood option.
123
124 Depending on the configuration, information may be heavily weighted to
125 some sides of a cell and quite limited on other sides, which may explain
126 why user feedback seems to dismiss this option.
127 */
128
129 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
130
131 /*=============================================================================
132 * Local Macro definitions
133 *============================================================================*/
134
135 /* Cache line multiple, in cs_real_t units */
136
137 #define CS_CL (CS_CL_SIZE/8)
138
139 /* Vector size, in cs_real_t units (large enough for AVX512) */
140
141 #define CS_VS 8
142
143 /*============================================================================
144 * Type definition
145 *============================================================================*/
146
147 /*============================================================================
148 * Global variables
149 *============================================================================*/
150
151 /* Short names for extended neighborhood types */
152
153 const char *cs_ext_neighborhood_type_name[]
154 = {N_("none"),
155 N_("complete"),
156 N_("opposite cell centers"),
157 N_("non-orthogonality threshold")};
158
159 static cs_ext_neighborhood_type_t _ext_nbh_type
160 = CS_EXT_NEIGHBORHOOD_CELL_CENTER_OPPOSITE;
161 static cs_real_t _non_ortho_max = 45;
162 static bool _full_nb_boundary = false;
163
164 /*============================================================================
165 * Private function definitions
166 *============================================================================*/
167
168 /*----------------------------------------------------------------------------*/
169 /*!
170 * \brief Compute the square norm of the distance from a point to a line
171 * defined by a vector.
172 *
173 * \param[in] u vector of 3 real values aligned with line
174 * \param[in] v vector of 3 real values from one pont one line
175 * to given point
176 *
177 * \return square of distance, given by: |u x v| / |v|
178 */
179 /*----------------------------------------------------------------------------*/
180
181 static inline cs_real_t
_cross_product_sq_norm(const cs_real_t u[3],const cs_real_t v[3])182 _cross_product_sq_norm(const cs_real_t u[3],
183 const cs_real_t v[3])
184 {
185 cs_real_t uv[3];
186
187 uv[0] = u[1]*v[2] - u[2]*v[1];
188 uv[1] = u[2]*v[0] - u[0]*v[2];
189 uv[2] = u[0]*v[1] - u[1]*v[0];
190
191 cs_real_t d2n = uv[0]*uv[0] + uv[1]*uv[1] + uv[2]*uv[2];
192 cs_real_t d2d = u[0]*u[0] + u[1]*u[1] + u[2]*u[2];
193
194 return d2n / d2d;
195 }
196
197 /*----------------------------------------------------------------------------
198 * Extract a mesh's "cell -> internal faces" connectivity.
199 *
200 * parameters:
201 * mesh <-- pointer to a cs_mesh_t structure
202 * p_cell_i_faces_idx --> pointer to the "cell -> faces" connectivity index
203 * (0 to n-1 numbering)
204 * p_cell_i_faces_lst --> pointer to the "cell -> faces" connectivity list
205 * (0 to n-1, no orientation)
206 *----------------------------------------------------------------------------*/
207
208 static void
_get_cell_i_faces_connectivity(const cs_mesh_t * mesh,cs_lnum_t ** const p_cell_i_faces_idx,cs_lnum_t ** const p_cell_i_faces_lst)209 _get_cell_i_faces_connectivity(const cs_mesh_t *mesh,
210 cs_lnum_t **const p_cell_i_faces_idx,
211 cs_lnum_t **const p_cell_i_faces_lst)
212 {
213
214 cs_lnum_t i, j, j1, j2;
215
216 cs_lnum_t *cell_faces_count = NULL;
217 cs_lnum_t *cell_faces_idx = NULL;
218 cs_lnum_t *cell_faces_lst = NULL;
219
220 /* Allocate and initialize index */
221
222 BFT_MALLOC(cell_faces_idx, mesh->n_cells + 1, cs_lnum_t);
223
224 for (i = 0; i < mesh->n_cells + 1; i++)
225 cell_faces_idx[i] = 0;
226
227 /* Count number of faces per cell (we assign the temporary counter
228 * to i and cell_faces_idx[i + 1] instead of cell_faces_idx[i]
229 * to simplify the next stage) */
230
231 /* Note: test if j < mesh->n_cells on internal faces to ignore
232 parallel and/or periodic ghost cells */
233
234 for (i = 0; i < mesh->n_i_faces; i++) {
235 j1 = mesh->i_face_cells[i][0];
236 j2 = mesh->i_face_cells[i][1];
237 if (j1 < mesh->n_cells)
238 cell_faces_idx[j1 + 1] += 1;
239 if (j2 < mesh->n_cells)
240 cell_faces_idx[j2 + 1] += 1;
241 }
242
243 /* Build position index */
244
245 cell_faces_idx[0] = 0;
246 for (j = 0; j < mesh->n_cells; j++)
247 cell_faces_idx[j + 1] += cell_faces_idx[j];
248
249 /* Build array of values */
250
251 BFT_MALLOC(cell_faces_lst, cell_faces_idx[mesh->n_cells], cs_lnum_t);
252 BFT_MALLOC(cell_faces_count, mesh->n_cells, cs_lnum_t);
253
254 for (i = 0; i < mesh->n_cells; i++)
255 cell_faces_count[i] = 0;
256
257 for (i = 0; i < mesh->n_i_faces; i++) {
258 j1 = mesh->i_face_cells[i][0];
259 j2 = mesh->i_face_cells[i][1];
260 if (j1 < mesh->n_cells) {
261 cell_faces_lst[cell_faces_idx[j1] + cell_faces_count[j1]] = i;
262 cell_faces_count[j1] += 1;
263 }
264 if (j2 < mesh->n_cells) {
265 cell_faces_lst[cell_faces_idx[j2] + cell_faces_count[j2]] = i;
266 cell_faces_count[j2] += 1;
267 }
268 }
269
270 BFT_FREE(cell_faces_count);
271
272 /* Set return values */
273
274 *p_cell_i_faces_idx = cell_faces_idx;
275 *p_cell_i_faces_lst = cell_faces_lst;
276 }
277
278 /*----------------------------------------------------------------------------
279 * Extract a mesh's "cell -> boundary faces" connectivity.
280 *
281 * parameters:
282 * mesh <-- pointer to a cs_mesh_t structure
283 * p_cell_b_faces_idx --> pointer to the "cell -> faces" connectivity index
284 * (0 to n-1 numbering)
285 * p_cell_b_faces_lst --> pointer to the "cell -> faces" connectivity list
286 * (0 to n-1, no orientation)
287 *----------------------------------------------------------------------------*/
288
289 static void
_get_cell_b_faces_connectivity(const cs_mesh_t * mesh,cs_lnum_t ** const p_cell_b_faces_idx,cs_lnum_t ** const p_cell_b_faces_lst)290 _get_cell_b_faces_connectivity(const cs_mesh_t *mesh,
291 cs_lnum_t **const p_cell_b_faces_idx,
292 cs_lnum_t **const p_cell_b_faces_lst)
293 {
294
295 cs_lnum_t i, j, j1;
296
297 cs_lnum_t *cell_faces_count = NULL;
298 cs_lnum_t *cell_faces_idx = NULL;
299 cs_lnum_t *cell_faces_lst = NULL;
300
301 /* Allocate and initialize index */
302
303 BFT_MALLOC(cell_faces_idx, mesh->n_cells + 1, cs_lnum_t);
304
305 for (i = 0; i < mesh->n_cells + 1; i++)
306 cell_faces_idx[i] = 0;
307
308 /* Count number of faces per cell (we assign the temporary counter
309 * to i and cell_faces_idx[i + 1] instead of cell_faces_idx[i]
310 * to simplify the next stage) */
311
312 /* Note: test if j < mesh->n_cells on internal faces to ignore
313 parallel and/or periodic ghost cells */
314
315 for (i = 0; i < mesh->n_b_faces; i++) {
316 j1 = mesh->b_face_cells[i];
317 if (j1 < mesh->n_cells)
318 cell_faces_idx[j1 + 1] += 1;
319 }
320
321 /* Build position index */
322
323 cell_faces_idx[0] = 0;
324 for (j = 0; j < mesh->n_cells; j++)
325 cell_faces_idx[j + 1] += cell_faces_idx[j];
326
327 /* Build array of values */
328
329 BFT_MALLOC(cell_faces_lst, cell_faces_idx[mesh->n_cells], cs_lnum_t);
330 BFT_MALLOC(cell_faces_count, mesh->n_cells, cs_lnum_t);
331
332 for (i = 0; i < mesh->n_cells; i++)
333 cell_faces_count[i] = 0;
334
335 for (i = 0; i < mesh->n_b_faces; i++) {
336 j1 = mesh->b_face_cells[i];
337 if (j1 < mesh->n_cells) {
338 cell_faces_lst[cell_faces_idx[j1] + cell_faces_count[j1]] = i;
339 cell_faces_count[j1] += 1;
340 }
341 }
342
343 BFT_FREE(cell_faces_count);
344
345 /* Set return values */
346
347 *p_cell_b_faces_idx = cell_faces_idx;
348 *p_cell_b_faces_lst = cell_faces_lst;
349 }
350
351 /*----------------------------------------------------------------------------
352 * Create a "vertex -> cells" connectivity.
353 *
354 * parameters:
355 * mesh <-- pointer to a cs_mesh_t structure
356 * p_vtx_cells_idx --> pointer to the "vtx -> cells" connectivity index
357 * p_vtx_cells_lst --> pointer to the "vtx -> cells" connectivity list
358 *----------------------------------------------------------------------------*/
359
360 static void
_create_vtx_cells_connect(cs_mesh_t * mesh,cs_lnum_t * p_vtx_cells_idx[],cs_lnum_t * p_vtx_cells_lst[])361 _create_vtx_cells_connect(cs_mesh_t *mesh,
362 cs_lnum_t *p_vtx_cells_idx[],
363 cs_lnum_t *p_vtx_cells_lst[])
364 {
365 cs_lnum_t i, j, idx;
366 cs_lnum_t vtx_id, face_id, cell_id;
367
368 bool already_seen;
369
370 cs_lnum_t vtx_cells_connect_size = 0;
371
372 cs_lnum_t *vtx_faces_idx = NULL, *vtx_faces_lst = NULL;
373 cs_lnum_t *vtx_cells_idx = NULL, *vtx_cells_lst = NULL;
374
375 const cs_lnum_t n_vertices = mesh->n_vertices;
376 const cs_lnum_t n_faces = mesh->n_i_faces;
377 const cs_lnum_t *face_vtx_idx = mesh->i_face_vtx_idx;
378 const cs_lnum_t *face_vtx_lst = mesh->i_face_vtx_lst;
379 const cs_lnum_2_t *face_cells = (const cs_lnum_2_t *)(mesh->i_face_cells);
380
381 cs_lnum_t vtx_cells_estimated_connect_size = 3 * n_vertices;
382
383 BFT_MALLOC(vtx_cells_idx, n_vertices + 1, cs_lnum_t);
384 BFT_MALLOC(vtx_faces_idx, n_vertices + 1, cs_lnum_t);
385
386 for (vtx_id = 0; vtx_id < n_vertices + 1; vtx_id++) {
387 vtx_cells_idx[vtx_id] = 0;
388 vtx_faces_idx[vtx_id] = 0;
389 }
390
391 /* Define vtx -> faces connectivity index */
392
393 for (face_id = 0; face_id < n_faces; face_id++) {
394
395 for (i = face_vtx_idx[face_id]; i < face_vtx_idx[face_id+1]; i++) {
396 vtx_id = face_vtx_lst[i];
397 vtx_faces_idx[vtx_id + 1] += 1;
398 }
399
400 } /* End of loop on faces */
401
402 vtx_faces_idx[0] = 0;
403 for (vtx_id = 0; vtx_id < n_vertices; vtx_id++)
404 vtx_faces_idx[vtx_id + 1] += vtx_faces_idx[vtx_id];
405
406 /* Allocation and definiton of "vtx -> faces" connectivity list */
407
408 BFT_MALLOC(vtx_faces_lst, vtx_faces_idx[n_vertices], cs_lnum_t);
409
410 for (face_id = 0; face_id < n_faces; face_id++) {
411
412 for (i = face_vtx_idx[face_id]; i < face_vtx_idx[face_id+1]; i++) {
413
414 vtx_id = face_vtx_lst[i];
415 vtx_faces_lst[vtx_faces_idx[vtx_id]] = face_id;
416 vtx_faces_idx[vtx_id] += 1;
417
418 }
419
420 } /* End of loop on faces */
421
422 for (vtx_id = n_vertices; vtx_id > 0; vtx_id--)
423 vtx_faces_idx[vtx_id] = vtx_faces_idx[vtx_id-1];
424 vtx_faces_idx[0] = 0;
425
426 /* Define "vertex -> cells" connectivity.
427 Use "vertex -> faces" connectivity and "face -> cells" connectivity */
428
429 BFT_MALLOC(vtx_cells_lst, vtx_cells_estimated_connect_size, cs_lnum_t);
430
431 vtx_cells_idx[0] = 0;
432
433 for (vtx_id = 0; vtx_id < n_vertices; vtx_id++) {
434
435 for (i = vtx_faces_idx[vtx_id]; i < vtx_faces_idx[vtx_id+1]; i++) {
436
437 face_id = vtx_faces_lst[i];
438
439 for (j = 0; j < 2; j++) { /* For the cells sharing this face */
440
441 cell_id = face_cells[face_id][j];
442
443 already_seen = false;
444 idx = vtx_cells_idx[vtx_id];
445
446 while ((already_seen == false) && (idx < vtx_cells_connect_size)) {
447 if (cell_id == vtx_cells_lst[idx])
448 already_seen = true;
449 idx++;
450 }
451
452 if (already_seen == false) {
453
454 if (vtx_cells_connect_size + 1 > vtx_cells_estimated_connect_size) {
455 vtx_cells_estimated_connect_size *= 2;
456 BFT_REALLOC(vtx_cells_lst,
457 vtx_cells_estimated_connect_size, cs_lnum_t);
458 }
459
460 vtx_cells_lst[vtx_cells_connect_size] = cell_id;
461 vtx_cells_connect_size += 1;
462
463 }
464
465 } /* End of loop on cells sharing this face */
466
467 } /* End of loop on faces sharing this vertex */
468
469 vtx_cells_idx[vtx_id+1] = vtx_cells_connect_size;
470
471 } /* End of loop on vertices */
472
473 BFT_REALLOC(vtx_cells_lst, vtx_cells_connect_size, cs_lnum_t);
474
475 /* Free memory */
476
477 BFT_FREE(vtx_faces_idx);
478 BFT_FREE(vtx_faces_lst);
479
480 *p_vtx_cells_idx = vtx_cells_idx;
481 *p_vtx_cells_lst = vtx_cells_lst;
482
483 }
484
485 /*----------------------------------------------------------------------------
486 * Create a "vertex -> cells" connectivity.
487 *
488 * parameters:
489 * face_id <-- identification number for the face
490 * cell_id <-- identification number for the cell sharing this face
491 * mesh <-- pointer to a cs_mesh_t structure
492 * cell_cells_tag <-- pointer to a tag array
493 * vtx_cells_idx --> pointer to the "vtx -> cells" connectivity index
494 * vtx_cells_lst --> pointer to the "vtx -> cells" connectivity list
495 * vtx_gcells_idx --> pointer to the "vtx -> gcells" connectivity index
496 * vtx_gcells_lst --> pointer to the "vtx -> gcells" connectivity list
497 *----------------------------------------------------------------------------*/
498
499 static void
_tag_cells(cs_lnum_t face_id,cs_lnum_t cell_id,const cs_mesh_t * mesh,char cell_cells_tag[],cs_lnum_t vtx_cells_idx[],cs_lnum_t vtx_cells_lst[],cs_lnum_t vtx_gcells_idx[],cs_lnum_t vtx_gcells_lst[])500 _tag_cells(cs_lnum_t face_id,
501 cs_lnum_t cell_id,
502 const cs_mesh_t *mesh,
503 char cell_cells_tag[],
504 cs_lnum_t vtx_cells_idx[],
505 cs_lnum_t vtx_cells_lst[],
506 cs_lnum_t vtx_gcells_idx[],
507 cs_lnum_t vtx_gcells_lst[])
508 {
509 const cs_lnum_t *cell_cells_lst = mesh->cell_cells_lst;
510 const cs_lnum_t *cell_cells_idx = mesh->cell_cells_idx;
511
512 const cs_lnum_t n_cells = mesh->n_cells;
513 const cs_lnum_t *face_vtx_idx = mesh->i_face_vtx_idx;
514 const cs_lnum_t *face_vtx_lst = mesh->i_face_vtx_lst;
515
516 if (cell_id < n_cells) {
517
518 for (cs_lnum_t i = cell_cells_idx[cell_id];
519 i < cell_cells_idx[cell_id+1];
520 i++) {
521
522 /* For cell not tagged yet */
523
524 if (cell_cells_tag[i] == 0) {
525
526 cs_lnum_t ext_cell_id = cell_cells_lst[i];
527
528 /* Cells sharing a vertex with the face */
529
530 for (cs_lnum_t j = face_vtx_idx[face_id];
531 j < face_vtx_idx[face_id+1];
532 j++) {
533
534 cs_lnum_t vtx_id = face_vtx_lst[j];
535
536 if (ext_cell_id < n_cells) { /* true cell */
537
538 for (cs_lnum_t k = vtx_cells_idx[vtx_id];
539 k < vtx_cells_idx[vtx_id+1];
540 k++) {
541
542 cs_lnum_t loc_cell_id = vtx_cells_lst[k];
543
544 /* Comparison and selection */
545
546 if (loc_cell_id == ext_cell_id && cell_cells_tag[i] == 0)
547 cell_cells_tag[i] = 1;
548
549 } /* End of loop on cells sharing this vertex */
550
551 }
552 else { /* ext_cell_num >= n_cells; i.e. ghost cell */
553
554 for (cs_lnum_t k = vtx_gcells_idx[vtx_id];
555 k < vtx_gcells_idx[vtx_id+1];
556 k++) {
557
558 cs_lnum_t loc_cell_id = vtx_gcells_lst[k] + n_cells;
559
560 /* Comparison and selection */
561
562 if (loc_cell_id == ext_cell_id && cell_cells_tag[i] == 0)
563 cell_cells_tag[i] = 1;
564
565 }
566
567 }
568
569 } /* End of loop on vertices of the face */
570
571 }
572
573 } /* End of loop on cells in the extended neighborhood */
574
575 } /* End if cell_id < n_cells */
576
577 }
578
579 /*---------------------------------------------------------------------------
580 * Reverse "ghost cell -> vertex" connectivity into "vertex -> ghost cells"
581 * connectivity for halo elements.
582 * Build the connectivity index.
583 *
584 * parameters:
585 * halo <-- pointer to a cs_halo_t structure
586 * rank_id <-- rank number to work with
587 * checker <-> temporary array to check vertices
588 * gcell_vtx_idx <-- "ghost cell -> vertices" connectivity index
589 * gcell_vtx_lst <-- "ghost cell -> vertices" connectivity list
590 * vtx_gcells_idx <-> "vertex -> ghost cells" connectivity index
591 *---------------------------------------------------------------------------*/
592
593 static void
_reverse_connectivity_idx(cs_halo_t * halo,cs_lnum_t n_vertices,cs_lnum_t rank_id,cs_lnum_t * checker,cs_lnum_t * gcell_vtx_idx,cs_lnum_t * gcell_vtx_lst,cs_lnum_t * vtx_gcells_idx)594 _reverse_connectivity_idx(cs_halo_t *halo,
595 cs_lnum_t n_vertices,
596 cs_lnum_t rank_id,
597 cs_lnum_t *checker,
598 cs_lnum_t *gcell_vtx_idx,
599 cs_lnum_t *gcell_vtx_lst,
600 cs_lnum_t *vtx_gcells_idx)
601 {
602 cs_lnum_t i, j, id, vtx_id, start_idx, end_idx;
603
604 /* Initialize index */
605
606 vtx_gcells_idx[0] = 0;
607 for (i = 0; i < n_vertices; i++) {
608 vtx_gcells_idx[i+1] = 0;
609 checker[i] = -1;
610 }
611
612 if (rank_id == -1) {
613 start_idx = 0;
614 end_idx = halo->n_elts[CS_HALO_EXTENDED];
615 }
616 else { /* Call with rank_id > 1 for standard halo */
617 start_idx = halo->index[2*rank_id];
618 end_idx = halo->index[2*rank_id+1];
619 }
620
621 /* Define index */
622
623 for (id = start_idx; id < end_idx; id++) {
624
625 for (j = gcell_vtx_idx[id]; j < gcell_vtx_idx[id+1]; j++) {
626
627 vtx_id = gcell_vtx_lst[j];
628
629 if (checker[vtx_id] != id) {
630 checker[vtx_id] = id;
631 vtx_gcells_idx[vtx_id+1] += 1;
632 }
633
634 }
635
636 } /* End of loop of ghost cells */
637
638 for (i = 0; i < n_vertices; i++)
639 vtx_gcells_idx[i+1] += vtx_gcells_idx[i];
640 }
641
642 /*---------------------------------------------------------------------------
643 * Reverse "ghost cells -> vertex" connectivity into "vertex -> ghost cells"
644 * connectivity for halo elements.
645 * Build the connectivity list.
646 *
647 * parameters:
648 * halo <-- pointer to a cs_halo_t structure
649 * n_vertices <-- number of vertices
650 * rank_id <-- rank number to work with
651 * counter <-> temporary array to count vertices
652 * checker <-> temporary array to check vertices
653 * gcell_vtx_idx <-- "ghost cell -> vertices" connectivity index
654 * gcell_vtx_lst <-- "ghost cell -> vertices" connectivity list
655 * vtx_gcells_idx <-- "vertex -> ghost cells" connectivity index
656 * vtx_gcells_lst <-> "vertex -> ghost cells" connectivity list
657 *---------------------------------------------------------------------------*/
658
659 static void
_reverse_connectivity_lst(cs_halo_t * halo,cs_lnum_t n_vertices,cs_lnum_t rank_id,cs_lnum_t * counter,cs_lnum_t * checker,cs_lnum_t * gcell_vtx_idx,cs_lnum_t * gcell_vtx_lst,cs_lnum_t * vtx_gcells_idx,cs_lnum_t * vtx_gcells_lst)660 _reverse_connectivity_lst(cs_halo_t *halo,
661 cs_lnum_t n_vertices,
662 cs_lnum_t rank_id,
663 cs_lnum_t *counter,
664 cs_lnum_t *checker,
665 cs_lnum_t *gcell_vtx_idx,
666 cs_lnum_t *gcell_vtx_lst,
667 cs_lnum_t *vtx_gcells_idx,
668 cs_lnum_t *vtx_gcells_lst)
669 {
670 cs_lnum_t i, j, id, shift, vtx_id, start_idx, end_idx;
671
672 /* Initialize buffers */
673
674 for (i = 0; i < n_vertices; i++) {
675 counter[i] = 0;
676 checker[i] = -1;
677 }
678
679 if (rank_id == -1) {
680 start_idx = 0;
681 end_idx = halo->n_elts[CS_HALO_EXTENDED];
682 }
683 else {
684 start_idx = halo->index[2*rank_id];
685 end_idx = halo->index[2*rank_id+1];
686 }
687
688 /* Fill the connectivity list */
689
690 for (id = start_idx; id < end_idx; id++) {
691
692 for (j = gcell_vtx_idx[id]; j < gcell_vtx_idx[id+1]; j++) {
693
694 vtx_id = gcell_vtx_lst[j];
695
696 if (checker[vtx_id] != id) {
697
698 checker[vtx_id] = id;
699 shift = vtx_gcells_idx[vtx_id] + counter[vtx_id];
700 vtx_gcells_lst[shift] = id;
701 counter[vtx_id] += 1;
702
703 }
704
705 }
706
707 } /* End of loop of ghost cells */
708
709 }
710
711 /*---------------------------------------------------------------------------
712 * Create a "vertex -> ghost cells" connectivity.
713 * Add mesh->n_cells to all the elements of vtx_gcells_lst to obtain
714 * the local cell numbering.
715 *
716 * parameters:
717 * halo <-- pointer to a cs_halo_t structure
718 * n_vertices <-- number of vertices
719 * gcell_vtx_idx <-- "ghost cell -> vertices" connectivity index
720 * gcell_vtx_lst <-- "ghost cell -> vertices" connectivity list
721 * p_vtx_gcells_idx --> pointer to "vertex -> ghost cells" index
722 * p_vtx_gcells_lst --> pointer to "vertex -> ghost cells" list
723 *---------------------------------------------------------------------------*/
724
725 static void
_create_vtx_gcells_connect(cs_halo_t * halo,cs_lnum_t n_vertices,cs_lnum_t * gcells_vtx_idx,cs_lnum_t * gcells_vtx_lst,cs_lnum_t * p_vtx_gcells_idx[],cs_lnum_t * p_vtx_gcells_lst[])726 _create_vtx_gcells_connect(cs_halo_t *halo,
727 cs_lnum_t n_vertices,
728 cs_lnum_t *gcells_vtx_idx,
729 cs_lnum_t *gcells_vtx_lst,
730 cs_lnum_t *p_vtx_gcells_idx[],
731 cs_lnum_t *p_vtx_gcells_lst[])
732 {
733 cs_lnum_t *vtx_buffer = NULL, *vtx_counter = NULL, *vtx_checker = NULL;
734 cs_lnum_t *vtx_gcells_idx = NULL, *vtx_gcells_lst = NULL;
735
736 BFT_MALLOC(vtx_buffer, 2*n_vertices, cs_lnum_t);
737 vtx_counter = &(vtx_buffer[0]);
738 vtx_checker = &(vtx_buffer[n_vertices]);
739
740 BFT_MALLOC(vtx_gcells_idx, n_vertices + 1, cs_lnum_t);
741
742 /* Create a vertex -> ghost cells connectivity */
743
744 _reverse_connectivity_idx(halo,
745 n_vertices,
746 -1,
747 vtx_checker,
748 gcells_vtx_idx,
749 gcells_vtx_lst,
750 vtx_gcells_idx);
751
752 BFT_MALLOC(vtx_gcells_lst, vtx_gcells_idx[n_vertices], cs_lnum_t);
753
754 _reverse_connectivity_lst(halo,
755 n_vertices,
756 -1,
757 vtx_counter,
758 vtx_checker,
759 gcells_vtx_idx,
760 gcells_vtx_lst,
761 vtx_gcells_idx,
762 vtx_gcells_lst);
763
764 *p_vtx_gcells_idx = vtx_gcells_idx;
765 *p_vtx_gcells_lst = vtx_gcells_lst;
766
767 /* Free memory */
768
769 BFT_FREE(vtx_buffer);
770
771 }
772
773 /*---------------------------------------------------------------------------
774 * Create a "vertex -> cells" connectivity.
775 *
776 * parameters:
777 * mesh <-- pointer to cs_mesh_t structure
778 * cell_i_faces_idx <-- "cell -> internal faces" connectivity index
779 * cell_i_faces_lst <-- "cell -> internal faces" connectivity list
780 * p_vtx_cells_idx --> pointer to "vertex -> cells" connectivity index
781 * p_vtx_cells_lst --> pointer to "vertex -> cells" connectivity list
782 *---------------------------------------------------------------------------*/
783
784 static void
_create_vtx_cells_connect2(cs_mesh_t * mesh,cs_lnum_t * cell_i_faces_idx,cs_lnum_t * cell_i_faces_lst,cs_lnum_t * p_vtx_cells_idx[],cs_lnum_t * p_vtx_cells_lst[])785 _create_vtx_cells_connect2(cs_mesh_t *mesh,
786 cs_lnum_t *cell_i_faces_idx,
787 cs_lnum_t *cell_i_faces_lst,
788 cs_lnum_t *p_vtx_cells_idx[],
789 cs_lnum_t *p_vtx_cells_lst[])
790 {
791 cs_lnum_t i, cell_id, fac_id, i_vtx;
792 cs_lnum_t shift, vtx_id;
793
794 cs_lnum_t *vtx_buffer = NULL, *vtx_count = NULL, *vtx_tag = NULL;
795 cs_lnum_t *vtx_cells_idx = NULL, *vtx_cells_lst = NULL;
796
797 const cs_lnum_t n_cells = mesh->n_cells;
798 const cs_lnum_t n_vertices = mesh->n_vertices;
799 const cs_lnum_t *fac_vtx_idx = mesh->i_face_vtx_idx;
800 const cs_lnum_t *fac_vtx_lst = mesh->i_face_vtx_lst;
801
802 /* Initialize buffers */
803
804 BFT_MALLOC(vtx_buffer, 2*n_vertices, cs_lnum_t);
805 BFT_MALLOC(vtx_cells_idx, n_vertices + 1, cs_lnum_t);
806
807 vtx_count = &(vtx_buffer[0]);
808 vtx_tag = &(vtx_buffer[n_vertices]);
809
810 vtx_cells_idx[0] = 0;
811 for (i = 0; i < n_vertices; i++) {
812
813 vtx_cells_idx[i + 1] = 0;
814 vtx_tag[i] = -1;
815 vtx_count[i] = 0;
816
817 }
818
819 /* Define index */
820
821 for (cell_id = 0; cell_id < n_cells; cell_id++) {
822
823 for (i = cell_i_faces_idx[cell_id]; i < cell_i_faces_idx[cell_id+1]; i++) {
824
825 fac_id = cell_i_faces_lst[i];
826
827 for (i_vtx = fac_vtx_idx[fac_id];
828 i_vtx < fac_vtx_idx[fac_id+1];
829 i_vtx++) {
830
831 vtx_id = fac_vtx_lst[i_vtx];
832
833 if (vtx_tag[vtx_id] != cell_id) {
834
835 vtx_cells_idx[vtx_id+1] += 1;
836 vtx_tag[vtx_id] = cell_id;
837
838 } /* Add this cell to the connectivity of the vertex */
839
840 } /* End of loop on vertices */
841
842 } /* End of loop on cell's faces */
843
844 } /* End of loop on cells */
845
846 for (i = 0; i < n_vertices; i++) {
847
848 vtx_cells_idx[i+1] += vtx_cells_idx[i];
849 vtx_tag[i] = -1;
850
851 }
852
853 BFT_MALLOC(vtx_cells_lst, vtx_cells_idx[n_vertices], cs_lnum_t);
854
855 /* Fill list */
856
857 for (cell_id = 0; cell_id < n_cells; cell_id++) {
858
859 for (i = cell_i_faces_idx[cell_id]; i < cell_i_faces_idx[cell_id+1]; i++) {
860
861 fac_id = cell_i_faces_lst[i];
862
863 for (i_vtx = fac_vtx_idx[fac_id];
864 i_vtx < fac_vtx_idx[fac_id+1];
865 i_vtx++) {
866
867 vtx_id = fac_vtx_lst[i_vtx];
868
869 if (vtx_tag[vtx_id] != cell_id) {
870
871 shift = vtx_cells_idx[vtx_id] + vtx_count[vtx_id];
872 vtx_tag[vtx_id] = cell_id;
873 vtx_cells_lst[shift] = cell_id;
874 vtx_count[vtx_id] += 1;
875
876 } /* Add this cell to the connectivity of the vertex */
877
878 } /* End of loop on vertices */
879
880 } /* End of loop on cell's faces */
881
882 } /* End of loop on cells */
883
884 BFT_FREE(vtx_buffer);
885
886 *p_vtx_cells_idx = vtx_cells_idx;
887 *p_vtx_cells_lst = vtx_cells_lst;
888
889 }
890
891 /*---------------------------------------------------------------------------
892 * Create a "cell -> cells" connectivity.
893 *
894 * parameters:
895 * mesh <-- pointer to cs_mesh_t structure
896 * cell_i_faces_idx <-- "cell -> faces" connectivity index
897 * cell_i_faces_lst <-- "cell -> faces" connectivity list
898 * vtx_gcells_idx --> "vertex -> ghost cells" connectivity index
899 * vtx_gcells_lst --> "vertex -> ghost cells" connectivity list
900 * vtx_cells_idx --> "vertex -> cells" connectivity index
901 * vtx_cells_lst --> "vertex -> cells" connectivity list
902 * p_cell_cells_idx --> pointer to "cell -> cells" connectivity index
903 * p_cell_cells_lst --> pointer to "cell -> cells" connectivity list
904 *---------------------------------------------------------------------------*/
905
906 static void
_create_cell_cells_connect(cs_mesh_t * mesh,cs_lnum_t * cell_i_faces_idx,cs_lnum_t * cell_i_faces_lst,cs_lnum_t * vtx_gcells_idx,cs_lnum_t * vtx_gcells_lst,cs_lnum_t * vtx_cells_idx,cs_lnum_t * vtx_cells_lst,cs_lnum_t * p_cell_cells_idx[],cs_lnum_t * p_cell_cells_lst[])907 _create_cell_cells_connect(cs_mesh_t *mesh,
908 cs_lnum_t *cell_i_faces_idx,
909 cs_lnum_t *cell_i_faces_lst,
910 cs_lnum_t *vtx_gcells_idx,
911 cs_lnum_t *vtx_gcells_lst,
912 cs_lnum_t *vtx_cells_idx,
913 cs_lnum_t *vtx_cells_lst,
914 cs_lnum_t *p_cell_cells_idx[],
915 cs_lnum_t *p_cell_cells_lst[])
916 {
917 cs_lnum_t i, j, i_cel, fac_id, i_vtx;
918 cs_lnum_t cell_id, vtx_id, shift;
919
920 cs_lnum_t *cell_buffer = NULL, *cell_tag = NULL, *cell_count = NULL;
921 cs_lnum_t *cell_cells_idx = NULL, *cell_cells_lst = NULL;
922
923 const cs_lnum_t n_cells = mesh->n_cells;
924 const cs_lnum_t n_cells_wghosts = mesh->n_cells_with_ghosts;
925 const cs_lnum_2_t *face_cells = (const cs_lnum_2_t *)(mesh->i_face_cells);
926 const cs_lnum_t *fac_vtx_idx = mesh->i_face_vtx_idx;
927 const cs_lnum_t *fac_vtx_lst = mesh->i_face_vtx_lst;
928
929 /* Allocate and initialize buffers */
930
931 BFT_MALLOC(cell_cells_idx, n_cells + 1, cs_lnum_t);
932 BFT_MALLOC(cell_buffer, n_cells_wghosts + n_cells, cs_lnum_t);
933
934 cell_tag = &(cell_buffer[0]);
935 cell_count = &(cell_buffer[n_cells_wghosts]);
936
937 cell_cells_idx[0] = 0;
938 for (i = 0; i < n_cells; i++) {
939 cell_cells_idx[i+1] = 0;
940 cell_count[i] = 0;
941 }
942
943 for (i = 0; i < n_cells_wghosts; i++)
944 cell_tag[i] = -1;
945
946 /* Define index */
947
948 for (i_cel = 0; i_cel < n_cells; i_cel++) {
949
950 /* First loop on faces to tag cells sharing a face */
951
952 for (i = cell_i_faces_idx[i_cel]; i < cell_i_faces_idx[i_cel+1]; i++) {
953
954 fac_id = cell_i_faces_lst[i];
955
956 cell_tag[face_cells[fac_id][0]] = i_cel;
957 cell_tag[face_cells[fac_id][1]] = i_cel;
958
959 } /* End of loop on cell's faces */
960
961 /* Second loop on faces to update index */
962
963 for (i = cell_i_faces_idx[i_cel]; i < cell_i_faces_idx[i_cel+1]; i++) {
964
965 fac_id = cell_i_faces_lst[i];
966
967 for (i_vtx = fac_vtx_idx[fac_id];
968 i_vtx < fac_vtx_idx[fac_id+1];
969 i_vtx++) {
970
971 vtx_id = fac_vtx_lst[i_vtx];
972
973 /* For cells belonging to this rank, get vertex -> cells connect. */
974
975 for (j = vtx_cells_idx[vtx_id]; j < vtx_cells_idx[vtx_id+1]; j++) {
976
977 cell_id = vtx_cells_lst[j];
978
979 if (cell_tag[cell_id] != i_cel) {
980 cell_cells_idx[i_cel+1] += 1;
981 cell_tag[cell_id] = i_cel;
982 }
983
984 }
985
986 if (n_cells_wghosts - n_cells > 0) { /* If there are ghost cells */
987
988 /* For ghost cells, get vertex -> ghost cells connect. */
989
990 for (j = vtx_gcells_idx[vtx_id]; j < vtx_gcells_idx[vtx_id+1]; j++) {
991
992 cell_id = vtx_gcells_lst[j] + n_cells;
993
994 if (cell_tag[cell_id] != i_cel) {
995 cell_cells_idx[i_cel+1] += 1;
996 cell_tag[cell_id] = i_cel;
997 }
998
999 }
1000
1001 } /* If there are ghost cells */
1002
1003 } /* End of loop on vertices */
1004
1005 } /* End of loop on cell's faces */
1006
1007 } /* End of loop on cells */
1008
1009 /* Create index */
1010
1011 for (i = 0; i < n_cells; i++)
1012 cell_cells_idx[i+1] += cell_cells_idx[i];
1013
1014 for (i = 0; i < n_cells_wghosts; i++)
1015 cell_tag[i] = -1;
1016
1017 BFT_MALLOC(cell_cells_lst, cell_cells_idx[n_cells], cs_lnum_t);
1018
1019 /* Fill list */
1020
1021 for (i_cel = 0; i_cel < n_cells; i_cel++) {
1022
1023 /* First loop on faces to tag cells sharing a face */
1024
1025 for (i = cell_i_faces_idx[i_cel]; i < cell_i_faces_idx[i_cel+1]; i++) {
1026
1027 fac_id = cell_i_faces_lst[i];
1028
1029 cell_tag[face_cells[fac_id][0]] = i_cel;
1030 cell_tag[face_cells[fac_id][1]] = i_cel;
1031
1032 } /* End of loop on cell's faces */
1033
1034 /* Second loop on faces to update list */
1035
1036 for (i = cell_i_faces_idx[i_cel]; i < cell_i_faces_idx[i_cel+1]; i++) {
1037
1038 fac_id = cell_i_faces_lst[i];
1039
1040 for (i_vtx = fac_vtx_idx[fac_id];
1041 i_vtx < fac_vtx_idx[fac_id+1];
1042 i_vtx++) {
1043
1044 vtx_id = fac_vtx_lst[i_vtx];
1045
1046 /* For cells belonging to this rank, get vertex -> cells connect. */
1047
1048 for (j = vtx_cells_idx[vtx_id]; j < vtx_cells_idx[vtx_id+1]; j++) {
1049
1050 cell_id = vtx_cells_lst[j];
1051
1052 if (cell_tag[cell_id] != i_cel) {
1053
1054 shift = cell_cells_idx[i_cel] + cell_count[i_cel];
1055 cell_cells_lst[shift] = cell_id;
1056 cell_tag[cell_id] = i_cel;
1057 cell_count[i_cel] += 1;
1058
1059 } /* Add this cell_id */
1060
1061 }
1062
1063 if (n_cells_wghosts - n_cells > 0) { /* If there are ghost cells */
1064
1065 /* For ghost cells, get vertex -> ghost cells connect. */
1066
1067 for (j = vtx_gcells_idx[vtx_id]; j < vtx_gcells_idx[vtx_id+1]; j++){
1068
1069 cell_id = vtx_gcells_lst[j] + n_cells;
1070
1071 if (cell_tag[cell_id] != i_cel) {
1072
1073 shift = cell_cells_idx[i_cel] + cell_count[i_cel];
1074 cell_cells_lst[shift] = cell_id;
1075 cell_tag[cell_id] = i_cel;
1076 cell_count[i_cel] += 1;
1077
1078 }
1079
1080 }
1081
1082 } /* If there are ghost cells */
1083
1084 } /* End of loop on vertices */
1085
1086 } /* End of loop on cell's faces */
1087
1088 } /* End of loop on cells */
1089
1090 /* Sort line elements by column id (for better access patterns) */
1091
1092 bool unique = cs_sort_indexed(n_cells, cell_cells_idx, cell_cells_lst);
1093
1094 assert(unique == true);
1095
1096 *p_cell_cells_idx = cell_cells_idx;
1097 *p_cell_cells_lst = cell_cells_lst;
1098
1099 /* Free memory */
1100
1101 BFT_FREE(cell_buffer);
1102 }
1103
1104 /*----------------------------------------------------------------------------*/
1105 /*!
1106 * \brief Reduce the "cell -> cells" connectivity for the
1107 * extended neighborhood using a non-orthogonality criterion.
1108 *
1109 * Note: Only cells sharing only a vertex or vertices (not a face)
1110 * belong to the "cell -> cells" connectivity.
1111 *
1112 * \param[in] mesh pointer to mesh structure
1113 * \param[in] mesh_quantities associated mesh quantities
1114 * \param[in, out] cell_cells_tag tag adjacencies to retain
1115 */
1116 /*----------------------------------------------------------------------------*/
1117
1118 static void
_neighborhood_reduce_anomax(cs_mesh_t * mesh,cs_mesh_quantities_t * mesh_quantities,char cell_cells_tag[])1119 _neighborhood_reduce_anomax(cs_mesh_t *mesh,
1120 cs_mesh_quantities_t *mesh_quantities,
1121 char cell_cells_tag[])
1122 {
1123 cs_lnum_t i, face_id, cell_i, cell_j;
1124
1125 cs_real_t v_ij[3];
1126 cs_real_t face_normal[3];
1127 cs_real_t norm_ij, face_norm, cos_ij_fn;
1128 cs_real_t dprod;
1129
1130 cs_lnum_t *vtx_cells_idx = NULL, *vtx_cells_lst = NULL;
1131 cs_lnum_t *vtx_gcells_idx = NULL, *vtx_gcells_lst = NULL;
1132
1133 /* Convert degrees to radians */
1134 const cs_real_t non_ortho_max = _non_ortho_max * cs_math_pi / 180;
1135
1136 const cs_lnum_t n_faces = mesh->n_i_faces;
1137 const cs_lnum_2_t *face_cells = (const cs_lnum_2_t *)(mesh->i_face_cells);
1138
1139 const cs_real_t cos_ij_fn_min = cos(non_ortho_max);
1140 const cs_real_t *cell_cen = mesh_quantities->cell_cen;
1141
1142 enum {X, Y, Z};
1143
1144 assert(mesh->dim == 3);
1145
1146 if (non_ortho_max <= 0)
1147 return;
1148
1149 /*
1150 For each internal face, we select in the extended neighborhood
1151 of the two cells sharing this face all the cells sharing a
1152 vertex of this face if the non-orthogonality angle is above a
1153 criterion.
1154 */
1155
1156 /*
1157 First: re-build a "vertex -> cells" connectivity
1158 ------------------------------------------------
1159 We have to invert the "face -> vertices" connectivity and then
1160 we will use the "face -> cells" conectivity.
1161 */
1162
1163 _create_vtx_cells_connect(mesh,
1164 &vtx_cells_idx,
1165 &vtx_cells_lst);
1166
1167 if (cs_mesh_n_g_ghost_cells(mesh) > 0)
1168 _create_vtx_gcells_connect(mesh->halo,
1169 mesh->n_vertices,
1170 mesh->gcell_vtx_idx,
1171 mesh->gcell_vtx_lst,
1172 &vtx_gcells_idx,
1173 &vtx_gcells_lst);
1174
1175 for (face_id = 0; face_id < n_faces; face_id++) {
1176
1177 /* We compute the cosine of the non-orthogonality angle
1178 of internal faces (angle between the normal of the face
1179 and the line between I (center of the cell I) and J (center
1180 of the cell J) */
1181
1182 /* Vector IJ and normal of the face */
1183
1184 cell_i = face_cells[face_id][0];
1185 cell_j = face_cells[face_id][1];
1186 dprod = 0;
1187
1188 for (i = 0; i < 3; i++) {
1189 v_ij[i] = cell_cen[3*cell_j + i] - cell_cen[3*cell_i + i];
1190 face_normal[i] = mesh_quantities->i_face_normal[3*face_id + i];
1191 dprod += v_ij[i]*face_normal[i];
1192 }
1193
1194 norm_ij = cs_math_3_norm(v_ij);
1195 face_norm = cs_math_3_norm(face_normal);
1196
1197 assert(norm_ij > 0.);
1198 assert(face_norm > 0.);
1199
1200 /* Dot product : norm_ij . face_norm */
1201
1202 cos_ij_fn = dprod / (norm_ij * face_norm);
1203
1204 /* Comparison to a predefined limit.
1205 This is non-orthogonal if we are below the limit and so we keep
1206 the cell in the extended neighborhood of the two cells sharing
1207 the face. (The cell is tagged (<0) then we will change the sign
1208 and eliminate all cells < 0 */
1209
1210 if (cos_ij_fn <= cos_ij_fn_min) {
1211
1212 /* For each cell sharing the face : intersection between
1213 cells in the extended neighborhood and cells sharing a
1214 vertex of the face. */
1215
1216 _tag_cells(face_id, cell_i, mesh,
1217 cell_cells_tag,
1218 vtx_cells_idx, vtx_cells_lst,
1219 vtx_gcells_idx, vtx_gcells_lst);
1220 _tag_cells(face_id, cell_j, mesh,
1221 cell_cells_tag,
1222 vtx_cells_idx, vtx_cells_lst,
1223 vtx_gcells_idx, vtx_gcells_lst);
1224
1225 }
1226
1227 } /* End of loop on faces */
1228
1229 /* Free "vertex -> cells" connectivity */
1230
1231 BFT_FREE(vtx_cells_idx);
1232 BFT_FREE(vtx_cells_lst);
1233
1234 if (cs_mesh_n_g_ghost_cells(mesh) > 0) {
1235 BFT_FREE(vtx_gcells_idx);
1236 BFT_FREE(vtx_gcells_lst);
1237 }
1238 }
1239
1240 /*----------------------------------------------------------------------------*/
1241 /*!
1242 * \brief Reduce extended neighborhood based on cell center opposite.
1243 *
1244 * \param[in, out] mesh pointer to a mesh structure
1245 * \param[in] mq associated mesh quantities
1246 * \param[in, out] cell_cells_tag tag adjacencies to retain
1247 */
1248 /*----------------------------------------------------------------------------*/
1249
1250 static void
_neighborhood_reduce_cell_center_opposite(cs_mesh_t * mesh,cs_mesh_quantities_t * mq,char cell_cells_tag[])1251 _neighborhood_reduce_cell_center_opposite(cs_mesh_t *mesh,
1252 cs_mesh_quantities_t *mq,
1253 char cell_cells_tag[])
1254 {
1255 cs_lnum_t *cell_i_faces_idx = NULL, *cell_i_faces_lst = NULL;
1256 cs_lnum_t *cell_b_faces_idx = NULL, *cell_b_faces_lst = NULL;
1257
1258 const cs_lnum_t *cell_cells_lst = mesh->cell_cells_lst;
1259 const cs_lnum_t *cell_cells_idx = mesh->cell_cells_idx;
1260
1261 const cs_lnum_t n_cells = mesh->n_cells;
1262
1263 const cs_lnum_2_t *i_face_cells = (const cs_lnum_2_t *)mesh->i_face_cells;
1264 const cs_real_3_t *cell_cen = (const cs_real_3_t *)mq->cell_cen;
1265 const cs_real_3_t *b_face_cog = (const cs_real_3_t *)mq->b_face_cog;
1266
1267 /* Get "cell -> faces" connectivity for the local mesh */
1268
1269 _get_cell_i_faces_connectivity(mesh,
1270 &cell_i_faces_idx,
1271 &cell_i_faces_lst);
1272
1273 _get_cell_b_faces_connectivity(mesh,
1274 &cell_b_faces_idx,
1275 &cell_b_faces_lst);
1276
1277 # pragma omp parallel if (n_cells > CS_THR_MIN)
1278 {
1279 cs_lnum_t t_s_id, t_e_id;
1280 cs_parall_thread_range(n_cells, sizeof(cs_real_t), &t_s_id, &t_e_id);
1281
1282 cs_real_3_t *n_c_s = NULL;
1283
1284 cs_lnum_t n_max_c = 0;
1285
1286 /* Loop on cells */
1287
1288 for (cs_lnum_t c_id = t_s_id; c_id < t_e_id; c_id++) {
1289
1290 /* Build cache for extended neighbors */
1291
1292 cs_lnum_t c_s_id = cell_cells_idx[c_id];
1293 cs_lnum_t c_e_id = cell_cells_idx[c_id+1];
1294 cs_lnum_t n_c = c_e_id - c_s_id;
1295
1296 if (n_c > n_max_c) {
1297 n_max_c = n_c*2;
1298 BFT_REALLOC(n_c_s, n_max_c, cs_real_3_t);
1299 }
1300
1301 for (cs_lnum_t i = 0; i < n_c; i++) {
1302 cs_lnum_t c_id_n = cell_cells_lst[c_s_id + i];
1303 for (cs_lnum_t j = 0; j < 3; j++)
1304 n_c_s[i][j] = cell_cen[c_id_n][j] - cell_cen[c_id][j];
1305 }
1306
1307 /* Loop on interior cell faces */
1308
1309 cs_lnum_t i_f_s_id = cell_i_faces_idx[c_id];
1310 cs_lnum_t i_f_e_id = cell_i_faces_idx[c_id+1];
1311
1312 cs_lnum_t b_f_s_id = cell_b_faces_idx[c_id];
1313 cs_lnum_t b_f_e_id = cell_b_faces_idx[c_id+1];
1314
1315 for (cs_lnum_t i = i_f_s_id; i < i_f_e_id; i++) {
1316
1317 cs_real_t i_dist_align_min = HUGE_VAL;
1318
1319 cs_lnum_t f_id_0 = cell_i_faces_lst[i];
1320 cs_lnum_t c_id_0 = i_face_cells[f_id_0][0];
1321 if (c_id_0 == c_id)
1322 c_id_0 = i_face_cells[f_id_0][1];
1323
1324 cs_real_t seg_0[3];
1325 for (cs_lnum_t k = 0; k < 3; k++)
1326 seg_0[k] = cell_cen[c_id][k] - cell_cen[c_id_0][k];
1327
1328 /* Initialize using base (face->cells) connectivity) */
1329
1330 for (cs_lnum_t j = i_f_s_id; j < i_f_e_id; j++) {
1331
1332 if (i == j) continue;
1333
1334 cs_real_t seg_1[3];
1335
1336 cs_lnum_t f_id_1 = cell_i_faces_lst[j];
1337 cs_lnum_t c_id_1 = i_face_cells[f_id_1][0];
1338 if (c_id_1 == c_id)
1339 c_id_1 = i_face_cells[f_id_1][1];
1340 for (cs_lnum_t k = 0; k < 3; k++)
1341 seg_1[k] = cell_cen[c_id_1][k] - cell_cen[c_id][k];
1342
1343 if (cs_math_3_dot_product(seg_0, seg_1) <= 0)
1344 continue;
1345
1346 cs_real_t d2 = _cross_product_sq_norm(seg_0, seg_1);
1347
1348 if (d2 < i_dist_align_min)
1349 i_dist_align_min = d2;
1350
1351 }
1352
1353 for (cs_lnum_t j = b_f_s_id; j < b_f_e_id; j++) {
1354
1355 cs_real_t seg_1[3];
1356
1357 cs_lnum_t f_id_1 = cell_b_faces_lst[j];
1358 for (cs_lnum_t k = 0; k < 3; k++)
1359 seg_1[k] = b_face_cog[f_id_1][k] - cell_cen[c_id][k];
1360
1361 if (cs_math_3_dot_product(seg_0, seg_1) <= 0)
1362 continue;
1363
1364 cs_real_t d2 = _cross_product_sq_norm(seg_0, seg_1);
1365
1366 if (d2 < i_dist_align_min)
1367 i_dist_align_min = d2;
1368
1369 }
1370
1371 /* Compare with cells in extended neighborhood */
1372
1373 cs_lnum_t i_cell_align_min = -1;
1374 for (cs_lnum_t j = 0; j < n_c; j++) {
1375
1376 if (cs_math_3_dot_product(seg_0, n_c_s[j]) <= 0)
1377 continue;
1378
1379 cs_real_t d2 = _cross_product_sq_norm(seg_0, n_c_s[j]);
1380
1381 if (d2 < i_dist_align_min) {
1382 i_cell_align_min = j;
1383 i_dist_align_min = d2;
1384 }
1385
1386 }
1387
1388 if (i_cell_align_min > -1)
1389 cell_cells_tag[c_s_id + i_cell_align_min] = 1;
1390
1391 } /* End of loop on interior cell faces */
1392
1393 /* Loop on boundary cell faces */
1394
1395 for (cs_lnum_t i = b_f_s_id; i < b_f_e_id; i++) {
1396
1397 cs_real_t i_dist_align_min = HUGE_VAL;
1398
1399 cs_lnum_t f_id_0 = cell_b_faces_lst[i];
1400
1401 cs_real_t seg_0[3];
1402 for (cs_lnum_t k = 0; k < 3; k++)
1403 seg_0[k] = cell_cen[c_id][k] - b_face_cog[f_id_0][k];
1404
1405 /* Initialize using base (face->cells) connectivity) */
1406
1407 for (cs_lnum_t j = i_f_s_id; j < i_f_e_id; j++) {
1408
1409 cs_real_t seg_1[3];
1410
1411 cs_lnum_t f_id_1 = cell_i_faces_lst[j];
1412 cs_lnum_t c_id_1 = i_face_cells[f_id_1][0];
1413 if (c_id_1 == c_id)
1414 c_id_1 = i_face_cells[f_id_1][1];
1415 for (cs_lnum_t k = 0; k < 3; k++)
1416 seg_1[k] = cell_cen[c_id_1][k] - cell_cen[c_id][k];
1417
1418 if (cs_math_3_dot_product(seg_0, seg_1) <= 0)
1419 continue;
1420
1421 cs_real_t d2 = _cross_product_sq_norm(seg_0, seg_1);
1422
1423 if (d2 < i_dist_align_min)
1424 i_dist_align_min = d2;
1425
1426 }
1427
1428 for (cs_lnum_t j = b_f_s_id; j < b_f_e_id; j++) {
1429
1430 if (i == j) continue;
1431
1432 cs_real_t seg_1[3];
1433
1434 cs_lnum_t f_id_1 = cell_b_faces_lst[j];
1435 for (cs_lnum_t k = 0; k < 3; k++)
1436 seg_1[k] = b_face_cog[f_id_1][k] - cell_cen[c_id][k];
1437
1438 if (cs_math_3_dot_product(seg_0, seg_1) <= 0)
1439 continue;
1440
1441 cs_real_t d2 = _cross_product_sq_norm(seg_0, seg_1);
1442
1443 if (d2 < i_dist_align_min)
1444 i_dist_align_min = d2;
1445
1446 }
1447
1448 /* Compare with cells in extended neighborhood */
1449
1450 cs_lnum_t i_cell_align_min = -1;
1451 for (cs_lnum_t j = 0; j < n_c; j++) {
1452
1453 if (cs_math_3_dot_product(seg_0, n_c_s[j]) <= 0)
1454 continue;
1455
1456 cs_real_t d2 = _cross_product_sq_norm(seg_0, n_c_s[j]);
1457
1458 if (d2 < i_dist_align_min) {
1459 i_cell_align_min = j;
1460 i_dist_align_min = d2;
1461 }
1462
1463 }
1464
1465 if (i_cell_align_min > -1)
1466 cell_cells_tag[c_s_id + i_cell_align_min] = 1;
1467
1468 } /* End of loop on boundary cell faces */
1469
1470 } /* End of loop on cells */
1471
1472 BFT_FREE(n_c_s);
1473
1474 } /* End of Open MP block */
1475
1476 BFT_FREE(cell_i_faces_idx);
1477 BFT_FREE(cell_i_faces_lst);
1478 BFT_FREE(cell_b_faces_idx);
1479 BFT_FREE(cell_b_faces_lst);
1480 }
1481
1482 /*----------------------------------------------------------------------------*/
1483 /*!
1484 * \brief Reduce extended neighborhood based on proximity to boundary
1485 *
1486 * \param[in, out] mesh pointer to a mesh structure
1487 * \param[in, out] cell_cells_tag tag adjacencies to retain
1488 */
1489 /*----------------------------------------------------------------------------*/
1490
1491 static void
_neighborhood_reduce_full_boundary(cs_mesh_t * mesh,char cell_cells_tag[])1492 _neighborhood_reduce_full_boundary(cs_mesh_t *mesh,
1493 char cell_cells_tag[])
1494 {
1495 const cs_lnum_t *cell_cells_idx = mesh->cell_cells_idx;
1496
1497 const cs_lnum_t n_b_cells = mesh->n_b_cells;
1498
1499 # pragma omp parallel if (n_b_cells > CS_THR_MIN)
1500 {
1501 cs_lnum_t t_s_id, t_e_id;
1502 cs_parall_thread_range(n_b_cells, sizeof(cs_real_t), &t_s_id, &t_e_id);
1503
1504 /* Loop on boundary cells */
1505
1506 for (cs_lnum_t b_c_id = t_s_id; b_c_id < t_e_id; b_c_id++) {
1507
1508 cs_lnum_t c_id = mesh->b_cells[b_c_id];
1509
1510 /* Build cache for extended neighbors */
1511
1512 cs_lnum_t c_s_id = cell_cells_idx[c_id];
1513 cs_lnum_t c_e_id = cell_cells_idx[c_id+1];
1514
1515 for (cs_lnum_t i = c_s_id; i < c_e_id; i++)
1516 cell_cells_tag[i] = 1;
1517
1518 } /* End of loop on boundary cells */
1519
1520 } /* End of Open MP block */
1521 }
1522
1523 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
1524
1525 /*============================================================================
1526 * Public function definitions
1527 *============================================================================*/
1528
1529 /*----------------------------------------------------------------------------*/
1530 /*!
1531 * \brief Get the extended neighborhood type.
1532 *
1533 * \return extended neighborhood type
1534 */
1535 /*----------------------------------------------------------------------------*/
1536
1537 cs_ext_neighborhood_type_t
cs_ext_neighborhood_get_type(void)1538 cs_ext_neighborhood_get_type(void)
1539 {
1540 return _ext_nbh_type;
1541 }
1542
1543 /*----------------------------------------------------------------------------*/
1544 /*!
1545 * \brief Set the extended neighborhood type.
1546 *
1547 * \param[in] enh_type extended neighborhood type
1548 */
1549 /*----------------------------------------------------------------------------*/
1550
1551 void
cs_ext_neighborhood_set_type(cs_ext_neighborhood_type_t enh_type)1552 cs_ext_neighborhood_set_type(cs_ext_neighborhood_type_t enh_type)
1553 {
1554 if ( enh_type < CS_EXT_NEIGHBORHOOD_NONE
1555 || enh_type > CS_EXT_NEIGHBORHOOD_NON_ORTHO_MAX)
1556 enh_type = CS_EXT_NEIGHBORHOOD_NONE;
1557
1558 _ext_nbh_type = enh_type;
1559 }
1560
1561 /*----------------------------------------------------------------------------*/
1562 /*!
1563 * \brief Get the non_orthogonality threshold (in degrees) associated with the
1564 * CS_EXT_NEIGHBORHOOD_NON_ORTHO_MAX neighborhood type.
1565 *
1566 * \return non-orthogonality threshold
1567 */
1568 /*----------------------------------------------------------------------------*/
1569
1570 cs_real_t
cs_ext_neighborhood_get_non_ortho_max(void)1571 cs_ext_neighborhood_get_non_ortho_max(void)
1572 {
1573 return _non_ortho_max;
1574 }
1575
1576 /*----------------------------------------------------------------------------*/
1577 /*!
1578 * \brief Set the non_orthogonality threshold (in degrees) associated with the
1579 * CS_EXT_NEIGHBORHOOD_NON_ORTHO_MAX neighborhood type.
1580 *
1581 * \param[in] non_ortho_max non-orthogonality threshold
1582 */
1583 /*----------------------------------------------------------------------------*/
1584
1585 void
cs_ext_neighborhood_set_non_ortho_max(cs_real_t non_ortho_max)1586 cs_ext_neighborhood_set_non_ortho_max(cs_real_t non_ortho_max)
1587 {
1588 _non_ortho_max = non_ortho_max;
1589 }
1590
1591 /*----------------------------------------------------------------------------*/
1592 /*!
1593 * \brief Reduce the "cell -> cells" connectivity for the
1594 * extended neighborhood using a non-orthogonality criterion.
1595 *
1596 * Note: Only cells sharing only a vertex or vertices (not a face)
1597 * belong to the "cell -> cells" connectivity.
1598 *
1599 * \param[in] mesh pointer to mesh structure
1600 * \param[in] mesh_quantities associated mesh quantities
1601 */
1602 /*----------------------------------------------------------------------------*/
1603
1604 void
cs_ext_neighborhood_reduce(cs_mesh_t * mesh,cs_mesh_quantities_t * mesh_quantities)1605 cs_ext_neighborhood_reduce(cs_mesh_t *mesh,
1606 cs_mesh_quantities_t *mesh_quantities)
1607 {
1608 /* Return if there is no extended neighborhood */
1609
1610 if ( mesh->cell_cells_idx == NULL
1611 || (mesh->halo_type == CS_HALO_STANDARD && _full_nb_boundary == false)
1612 || _ext_nbh_type == CS_EXT_NEIGHBORHOOD_COMPLETE)
1613 return;
1614
1615 cs_lnum_t i, cell_id;
1616 cs_gnum_t init_cell_cells_connect_size;
1617
1618 cs_gnum_t n_deleted_cells = 0;
1619 cs_lnum_t previous_idx = 0, new_idx = -1;
1620
1621 cs_lnum_t *cell_cells_idx = mesh->cell_cells_idx;
1622 cs_lnum_t *cell_cells_lst = mesh->cell_cells_lst;
1623
1624 const cs_lnum_t n_cells = mesh->n_cells;
1625
1626 /* Tag cells to eliminate (0 to eliminate, 1 to keep) */
1627
1628 char *cell_cells_tag;
1629 BFT_MALLOC(cell_cells_tag, mesh->cell_cells_idx[n_cells], char);
1630 for (i = 0; i < mesh->cell_cells_idx[n_cells]; i++)
1631 cell_cells_tag[i] = 0;
1632
1633 switch(_ext_nbh_type) {
1634 case CS_EXT_NEIGHBORHOOD_NON_ORTHO_MAX:
1635 _neighborhood_reduce_anomax(mesh,
1636 mesh_quantities,
1637 cell_cells_tag);
1638 break;
1639 case CS_EXT_NEIGHBORHOOD_CELL_CENTER_OPPOSITE:
1640 _neighborhood_reduce_cell_center_opposite(mesh,
1641 mesh_quantities,
1642 cell_cells_tag);
1643 break;
1644 default:
1645 break;
1646 }
1647
1648 if (_full_nb_boundary)
1649 _neighborhood_reduce_full_boundary(mesh, cell_cells_tag);
1650
1651 /* Delete cells with tag == 0 */
1652
1653 init_cell_cells_connect_size = cell_cells_idx[n_cells];
1654
1655 for (cell_id = 0; cell_id < n_cells; cell_id++) {
1656
1657 for (i = previous_idx; i < cell_cells_idx[cell_id+1]; i++) {
1658
1659 if (cell_cells_tag[i] != 0) {
1660 new_idx++;
1661 cell_cells_lst[new_idx] = cell_cells_lst[i];
1662 }
1663 else
1664 n_deleted_cells++;
1665
1666 } /* End of loop on cells in the extended neighborhood of cell_id+1 */
1667
1668 previous_idx = cell_cells_idx[cell_id+1];
1669 cell_cells_idx[cell_id+1] -= n_deleted_cells;
1670
1671 } /* End of loop on cells */
1672
1673 BFT_FREE(cell_cells_tag);
1674
1675 /* Reallocation of cell_cells_lst */
1676
1677 BFT_REALLOC(mesh->cell_cells_lst, cell_cells_idx[n_cells], cs_lnum_t);
1678
1679 /* Output for log */
1680
1681 if (mesh->verbosity > 0) {
1682
1683 #if defined(HAVE_MPI)
1684
1685 if (cs_glob_n_ranks > 1) {
1686
1687 cs_gnum_t count_g[2], count_l[2];
1688
1689 count_l[0] = init_cell_cells_connect_size;
1690 count_l[1] = n_deleted_cells;
1691
1692 MPI_Allreduce(count_l, count_g, 2, CS_MPI_GNUM,
1693 MPI_SUM, cs_glob_mpi_comm);
1694
1695 init_cell_cells_connect_size = count_g[0];
1696 n_deleted_cells = count_g[1];
1697 }
1698
1699 #endif
1700
1701 double ratio = 100;
1702 if (init_cell_cells_connect_size)
1703 ratio *= (init_cell_cells_connect_size - n_deleted_cells)
1704 / init_cell_cells_connect_size;
1705
1706 bft_printf
1707 (_("\n"
1708 " Extended neighborhood reduction: %s\n"
1709 " --------------------------------\n"
1710 "\n"
1711 " Size of complete cell-cell connectivity: %12llu\n"
1712 " Size of filtered cell-cell connectivity: %12llu\n"
1713 " (%2.2g %% of extended connectivity retained)\n"),
1714 _(cs_ext_neighborhood_type_name[_ext_nbh_type]),
1715 (unsigned long long)init_cell_cells_connect_size,
1716 (unsigned long long)(init_cell_cells_connect_size - n_deleted_cells),
1717 ratio);
1718
1719 }
1720
1721 #if 0 /* For debugging purposes */
1722 for (i = 0; i < mesh->n_cells; i++) {
1723 cs_lnum_t j;
1724 bft_printf(" cell %d :: ", i+1);
1725 for (j = mesh->cell_cells_idx[i];
1726 j < mesh->cell_cells_idx[i+1];
1727 j++)
1728 bft_printf(" %d ", mesh->cell_cells_lst[j]);
1729 bft_printf("\n");
1730 }
1731 #endif
1732
1733 cs_sort_indexed(n_cells, mesh->cell_cells_idx, mesh->cell_cells_lst);
1734
1735 cs_mesh_adjacencies_update_cell_cells_e();
1736 }
1737
1738 /*----------------------------------------------------------------------------*/
1739 /*!
1740 * \brief Create the "cell -> cells" connectivity.
1741 *
1742 * \param[in, out] mesh pointer to a mesh structure
1743 */
1744 /*----------------------------------------------------------------------------*/
1745
1746 void
cs_ext_neighborhood_define(cs_mesh_t * mesh)1747 cs_ext_neighborhood_define(cs_mesh_t *mesh)
1748 {
1749 cs_lnum_t *vtx_gcells_idx = NULL, *vtx_gcells_lst = NULL;
1750 cs_lnum_t *vtx_cells_idx = NULL, *vtx_cells_lst = NULL;
1751 cs_lnum_t *cell_i_faces_idx = NULL, *cell_i_faces_lst = NULL;
1752 cs_lnum_t *cell_cells_idx = NULL, *cell_cells_lst = NULL;
1753
1754 cs_halo_t *halo = mesh->halo;
1755
1756 /* Get "cell -> faces" connectivity for the local mesh */
1757
1758 _get_cell_i_faces_connectivity(mesh,
1759 &cell_i_faces_idx,
1760 &cell_i_faces_lst);
1761
1762 /* Create a "vertex -> cell" connectivity */
1763
1764 _create_vtx_cells_connect2(mesh,
1765 cell_i_faces_idx,
1766 cell_i_faces_lst,
1767 &vtx_cells_idx,
1768 &vtx_cells_lst);
1769
1770 if (cs_mesh_n_g_ghost_cells(mesh) > 0) {
1771
1772 /* Create a "vertex -> ghost cells" connectivity */
1773
1774 _create_vtx_gcells_connect(halo,
1775 mesh->n_vertices,
1776 mesh->gcell_vtx_idx,
1777 mesh->gcell_vtx_lst,
1778 &vtx_gcells_idx,
1779 &vtx_gcells_lst);
1780
1781 }
1782
1783 /* Create the "cell -> cells" connectivity for the extended halo */
1784
1785 _create_cell_cells_connect(mesh,
1786 cell_i_faces_idx,
1787 cell_i_faces_lst,
1788 vtx_gcells_idx,
1789 vtx_gcells_lst,
1790 vtx_cells_idx,
1791 vtx_cells_lst,
1792 &cell_cells_idx,
1793 &cell_cells_lst);
1794
1795 mesh->cell_cells_idx = cell_cells_idx;
1796 mesh->cell_cells_lst = cell_cells_lst;
1797
1798 /* Free memory */
1799
1800 BFT_FREE(vtx_gcells_idx);
1801 BFT_FREE(vtx_gcells_lst);
1802
1803 BFT_FREE(cell_i_faces_idx);
1804 BFT_FREE(cell_i_faces_lst);
1805 BFT_FREE(vtx_cells_idx);
1806 BFT_FREE(vtx_cells_lst);
1807 }
1808
1809 /*----------------------------------------------------------------------------*/
1810
1811 END_C_DECLS
1812