1 /*============================================================================
2  * Extract nodal connectivity mesh representations from a native mesh.
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 <stdio.h>
35 #include <stdlib.h>
36 #include <string.h>
37 
38 /*----------------------------------------------------------------------------
39  * Local headers
40  *----------------------------------------------------------------------------*/
41 
42 #include "bft_error.h"
43 #include "bft_mem.h"
44 
45 #include "cs_base.h"
46 #include "cs_mesh.h"
47 #include "cs_sort.h"
48 
49 #include "fvm_defs.h"
50 #include "fvm_nodal.h"
51 #include "fvm_nodal_from_desc.h"
52 #include "fvm_nodal_order.h"
53 
54 /*----------------------------------------------------------------------------
55  *  Header for the current file
56  *----------------------------------------------------------------------------*/
57 
58 #include "cs_mesh_connect.h"
59 
60 /*----------------------------------------------------------------------------*/
61 
62 BEGIN_C_DECLS
63 
64 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
65 
66 /*=============================================================================
67  * Local Macro Definitions
68  *============================================================================*/
69 
70 /*=============================================================================
71  * Local Type Definitions
72  *============================================================================*/
73 
74 /*============================================================================
75  *  Global variables
76  *============================================================================*/
77 
78 /*============================================================================
79  * Private function definitions
80  *============================================================================*/
81 
82 /*----------------------------------------------------------------------------
83  * Add a subset of a mesh's faces to an external mesh.
84  *
85  * The lists of faces to extract are optional (if none is given, boundary
86  * faces are extracted by default); they do not need to be ordered on input,
87  * but they are always ordered on exit (as faces are extracted by increasing
88  * number traversal, the lists are reordered to ensure the coherency of
89  * the extracted mesh's link to its parent faces, built using these lists).
90  *
91  * parameters:
92  *   mesh             <-- base mesh
93  *   extr_mesh        <-- extracted mesh name
94  *   boundary_flag    <-- -1 if unspecified, 0 if faces are not on boundary,
95  *                        1 if faces are on boundary
96  *   include_families <-- include family info if true
97  *   i_face_list_size <-- size of i_face_list[] array
98  *   b_face_list_size <-- size of b_face_list[] array
99  *   i_face_list      <-> list of interior faces (1 to n), or NULL
100  *   b_face_list      <-> list of boundary faces (1 to n), or NULL
101  *----------------------------------------------------------------------------*/
102 
103 static void
_add_faces_to_nodal(const cs_mesh_t * mesh,fvm_nodal_t * extr_mesh,int boundary_flag,bool include_families,cs_lnum_t i_face_list_size,cs_lnum_t b_face_list_size,cs_lnum_t i_face_list[],cs_lnum_t b_face_list[])104 _add_faces_to_nodal(const cs_mesh_t  *mesh,
105                     fvm_nodal_t      *extr_mesh,
106                     int               boundary_flag,
107                     bool              include_families,
108                     cs_lnum_t         i_face_list_size,
109                     cs_lnum_t         b_face_list_size,
110                     cs_lnum_t         i_face_list[],
111                     cs_lnum_t         b_face_list[])
112 {
113   cs_lnum_t   face_id, i;
114 
115   cs_lnum_t   n_max_faces = 0;
116   cs_lnum_t   b_face_count = 0;
117   cs_lnum_t   i_face_count = 0;
118   cs_lnum_t   extr_face_count = 0;
119   cs_lnum_t  *extr_face_idx = NULL;
120   cs_lnum_t  *extr_face_list = NULL;
121 
122   cs_lnum_t   face_num_shift[3];
123   cs_lnum_t  *face_vertices_idx[2];
124   cs_lnum_t  *face_vertices_num[2];
125   const int   *_face_families[2];
126   const int   **face_families = NULL;
127 
128   /* Count the number of faces to convert */
129 
130   n_max_faces = mesh->n_i_faces + mesh->n_b_faces;
131   BFT_MALLOC(extr_face_idx, n_max_faces, cs_lnum_t);
132 
133   /* Initialize index as marker */
134 
135   for (face_id = 0; face_id < n_max_faces; face_id++)
136     extr_face_idx[face_id] = -1;
137 
138   if (b_face_list_size == mesh->n_b_faces) {
139     for (face_id = 0; face_id < mesh->n_b_faces; face_id++)
140       extr_face_idx[face_id] = 1;
141   }
142   else if (b_face_list != NULL) {
143     for (face_id = 0; face_id < b_face_list_size; face_id++)
144       extr_face_idx[b_face_list[face_id] - 1] = 1;
145   }
146 
147   if (i_face_list_size == mesh->n_i_faces) {
148     for (face_id = 0; face_id < mesh->n_i_faces; face_id++)
149       extr_face_idx[face_id + mesh->n_b_faces] = 1;
150   }
151   else if (i_face_list != NULL) {
152     for (face_id = 0; face_id < i_face_list_size; face_id++)
153       extr_face_idx[i_face_list[face_id] - 1 + mesh->n_b_faces] = 1;
154   }
155 
156   /* Convert marked ids to indexes (1 to n) and reconstruct values of
157      b_face_list[] and i_face_list[] to ensure that they are ordered. */
158 
159   b_face_count = 0;
160   i_face_count = 0;
161 
162   if (b_face_list != NULL) {
163     for (face_id = 0; face_id < mesh->n_b_faces; face_id++) {
164       if (extr_face_idx[face_id] == 1) {
165         b_face_list[b_face_count] = face_id + 1;
166         b_face_count++;
167       }
168     }
169   }
170   else
171     b_face_count = CS_MIN(b_face_list_size, mesh->n_b_faces);
172 
173   if (i_face_list != NULL) {
174     for (face_id = 0, i = mesh->n_b_faces;
175          face_id < mesh->n_i_faces;
176          face_id++, i++) {
177       if (extr_face_idx[i] == 1) {
178         i_face_list[i_face_count] = face_id + 1;
179         i_face_count++;
180       }
181     }
182   }
183   else
184     i_face_count = CS_MIN(i_face_list_size, mesh->n_i_faces);
185 
186   BFT_FREE(extr_face_idx);
187 
188   /* Build a contiguous list (boundary faces, interior faces) */
189 
190   extr_face_count = b_face_count + i_face_count;
191 
192   BFT_MALLOC(extr_face_list, extr_face_count, cs_lnum_t);
193 
194   if (b_face_list != NULL) {
195     for (face_id = 0; face_id < b_face_count; face_id++)
196       extr_face_list[face_id] = b_face_list[face_id];
197   }
198   else if (b_face_list == NULL) { /* boundary faces by default if no list */
199     for (face_id = 0; face_id < b_face_count; face_id++)
200       extr_face_list[face_id] = face_id + 1;
201   }
202 
203   if (i_face_list != NULL) {
204     for (face_id = 0, i = b_face_count; face_id < i_face_count; face_id++, i++)
205       extr_face_list[i] = i_face_list[face_id] + mesh->n_b_faces;
206   }
207   else if (i_face_list == NULL) {
208     for (face_id = 0, i = b_face_count; face_id < i_face_count; face_id++, i++)
209       extr_face_list[i] = face_id + mesh->n_b_faces + 1;
210   }
211 
212   if (include_families) {
213     _face_families[0] = mesh->b_face_family;
214     _face_families[1] = mesh->i_face_family;
215     face_families = _face_families;
216   }
217 
218   /* Build the nodal connectivity */
219 
220   face_num_shift[0] = 0;
221   face_num_shift[1] = mesh->n_b_faces + face_num_shift[0];
222   face_num_shift[2] = mesh->n_i_faces + face_num_shift[1];
223 
224   face_vertices_idx[0] = mesh->b_face_vtx_idx;
225   face_vertices_idx[1] = mesh->i_face_vtx_idx;
226   face_vertices_num[0] = mesh->b_face_vtx_lst;
227   face_vertices_num[1] = mesh->i_face_vtx_lst;
228 
229   fvm_nodal_from_desc_add_faces(extr_mesh,
230                                 boundary_flag,
231                                 extr_face_count,
232                                 extr_face_list,
233                                 2,
234                                 face_num_shift,
235                                 (const cs_lnum_t **)face_vertices_idx,
236                                 (const cs_lnum_t **)face_vertices_num,
237                                 face_families,
238                                 NULL);
239 
240   BFT_FREE(extr_face_list);
241 }
242 
243 /*----------------------------------------------------------------------------
244  * Order added faces in an external mesh.
245  *
246  * parameters:
247  *   mesh             <-- base mesh
248  *   extr_mesh        <-- extracted mesh name
249  *----------------------------------------------------------------------------*/
250 
251 static void
_order_nodal_faces(const cs_mesh_t * mesh,fvm_nodal_t * extr_mesh)252 _order_nodal_faces(const cs_mesh_t  *mesh,
253                    fvm_nodal_t      *extr_mesh)
254 {
255   cs_lnum_t   face_id, i;
256 
257   cs_lnum_t   n_max_faces = 0;
258 
259   cs_gnum_t  *num_glob_fac = NULL;
260 
261   /* Count the number of faces to convert */
262 
263   n_max_faces = mesh->n_i_faces + mesh->n_b_faces;
264 
265   /* In case of parallelism or face renumbering, sort faces by
266      increasing global number */
267 
268   if (mesh->global_i_face_num != NULL || mesh->global_b_face_num != NULL) {
269 
270     BFT_MALLOC(num_glob_fac, n_max_faces, cs_gnum_t);
271 
272     if (mesh->global_b_face_num == NULL) {
273       for (face_id = 0; face_id < mesh->n_b_faces; face_id++)
274         num_glob_fac[face_id] = face_id + 1;
275     }
276     else {
277       for (face_id = 0; face_id < mesh->n_b_faces; face_id++)
278         num_glob_fac[face_id] = mesh->global_b_face_num[face_id];
279     }
280 
281     assert(mesh->n_g_b_faces + mesh->n_g_i_faces > 0);
282 
283     if (mesh->global_i_face_num == NULL) {
284       for (face_id = 0, i = mesh->n_b_faces;
285            face_id < mesh->n_i_faces;
286            face_id++, i++)
287         num_glob_fac[i] = face_id + 1 + mesh->n_g_b_faces;
288     }
289     else {
290       for (face_id = 0, i = mesh->n_b_faces;
291            face_id < mesh->n_i_faces;
292            face_id++, i++)
293         num_glob_fac[i] = mesh->global_i_face_num[face_id] + mesh->n_g_b_faces;
294     }
295 
296   }
297 
298   /* Sort faces by increasing global number */
299 
300   fvm_nodal_order_faces(extr_mesh, num_glob_fac);
301   fvm_nodal_init_io_num(extr_mesh, num_glob_fac, 2);
302 
303   if (num_glob_fac != NULL)
304     BFT_FREE(num_glob_fac);
305 }
306 
307 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
308 
309 /*============================================================================
310  * Public function definitions
311  *============================================================================*/
312 
313 /*----------------------------------------------------------------------------
314  * Extract a mesh's "cells -> faces" connectivity.
315  *
316  * We consider a common numbering for internal and boundary faces, in which
317  * boundary faces are defined first. The common id for the i-th boundary
318  * face is thus i, and that of the j-th interior face is n_b_faces + j.
319  *
320  * If ind_cel_extr != NULL, then:
321  * --- ind_cel_extr[cell_id] = id in the list to extract (0 to n-1)
322  *     if cell cell_id should be extracted
323  * --- ind_cel_extr[cell_id] = -1 if cells cell_id should be ignored
324  *
325  * parameters:
326  *   mesh             <-- pointer to mesh structure
327  *   extr_cell_size   <-- size of extr_cell_id[] array
328  *   extr_cell_id     <-- extr_cell_id = ids of extracted cells, or -1
329  *   p_cell_faces_idx --> cells -> faces index
330  *   p_cell_faces_val --> cells -> faces connectivity
331  *----------------------------------------------------------------------------*/
332 
333 void
cs_mesh_connect_get_cell_faces(const cs_mesh_t * mesh,cs_lnum_t extr_cell_size,const cs_lnum_t extr_cell_id[],cs_lnum_t ** const p_cell_faces_idx,cs_lnum_t ** const p_cell_faces_val)334 cs_mesh_connect_get_cell_faces(const cs_mesh_t         *mesh,
335                                cs_lnum_t                extr_cell_size,
336                                const cs_lnum_t          extr_cell_id[],
337                                cs_lnum_t        **const p_cell_faces_idx,
338                                cs_lnum_t        **const p_cell_faces_val)
339 {
340   cs_lnum_t  cell_id, c_id1, c_id2, face_id, n_loc_cells;
341 
342   cs_lnum_t  *cell_face_count = NULL;
343   cs_lnum_t  *cell_faces_idx = NULL;
344   cs_lnum_t  *cell_faces_val = NULL;
345 
346   /* Allocate and initialize cell ->faces index */
347 
348   n_loc_cells = mesh->n_cells;
349   if (extr_cell_id != NULL)
350     n_loc_cells = extr_cell_size;
351 
352   BFT_MALLOC(cell_faces_idx, n_loc_cells + 1, cs_lnum_t);
353 
354   for (cell_id = 0; cell_id < n_loc_cells + 1; cell_id++)
355     cell_faces_idx[cell_id] = 0;
356 
357   /* Count number of faces per cell (we assign the temporary counter
358      corresponding to cell_id to cell_faces_idx[cell_id + 1] rather than
359      cell_faces_idx[cell_id] to simplify the next step) */
360 
361   /* Remark: test if cell_id < mesh->n_cells on internal faces so
362      as to ignore ghost cells */
363 
364   for (face_id = 0; face_id < mesh->n_b_faces; face_id++) {
365     cell_id = mesh->b_face_cells[face_id];
366     if (extr_cell_id != NULL)
367       cell_id = extr_cell_id[cell_id];
368     if (cell_id > -1)
369       cell_faces_idx[cell_id + 1] += 1;
370   }
371 
372   for (face_id = 0; face_id < mesh->n_i_faces; face_id++) {
373     c_id1 = mesh->i_face_cells[face_id][0];
374     c_id2 = mesh->i_face_cells[face_id][1];
375     if (extr_cell_id != NULL) {
376       if (c_id1 < mesh->n_cells)
377         c_id1 = extr_cell_id[c_id1];
378       else
379         c_id1 = -1;
380       if (c_id2 < mesh->n_cells)
381         c_id2 = extr_cell_id[c_id2];
382       else
383         c_id2 = -1;
384     }
385     if (c_id1 > -1 && c_id1 < mesh->n_cells)
386       cell_faces_idx[c_id1 + 1] += 1;
387     if (c_id2 > -1 && c_id2 < mesh->n_cells)
388       cell_faces_idx[c_id2 + 1] += 1;
389   }
390 
391   /* Build cell -> faces index */
392 
393   cell_faces_idx[0] = 1;
394   for (cell_id = 0; cell_id < n_loc_cells; cell_id++)
395     cell_faces_idx[cell_id + 1] += cell_faces_idx[cell_id];
396 
397   /* Build array of values */
398 
399   BFT_MALLOC(cell_faces_val, cell_faces_idx[n_loc_cells] - 1, cs_lnum_t);
400   BFT_MALLOC(cell_face_count, n_loc_cells, cs_lnum_t);
401 
402   for (cell_id = 0; cell_id < n_loc_cells; cell_id++)
403     cell_face_count[cell_id] = 0;
404 
405   for (face_id = 0; face_id < mesh->n_b_faces; face_id++) {
406     cell_id = mesh->b_face_cells[face_id];
407     if (extr_cell_id != NULL)
408       cell_id = extr_cell_id[cell_id];
409     if (cell_id > -1) {
410       cell_faces_val[cell_faces_idx[cell_id] + cell_face_count[cell_id] - 1]
411         = face_id + 1;
412       cell_face_count[cell_id] += 1;
413     }
414   }
415 
416   for (face_id = 0; face_id < mesh->n_i_faces; face_id++) {
417     c_id1 = mesh->i_face_cells[face_id][0];
418     c_id2 = mesh->i_face_cells[face_id][1];
419     if (extr_cell_id != NULL) {
420       if (c_id1 < mesh->n_cells)
421         c_id1 = extr_cell_id[c_id1];
422       else
423         c_id1 = -1;
424       if (c_id2 < mesh->n_cells)
425         c_id2 = extr_cell_id[c_id2];
426       else
427         c_id2 = -1;
428     }
429     if (c_id1 > -1 && c_id1 < mesh->n_cells) {
430       cell_faces_val[cell_faces_idx[c_id1] + cell_face_count[c_id1] - 1]
431         =   face_id + mesh->n_b_faces + 1;
432       cell_face_count[c_id1] += 1;
433     }
434     if (c_id2 > -1 && c_id2 < mesh->n_cells) {
435       cell_faces_val[cell_faces_idx[c_id2] + cell_face_count[c_id2] - 1]
436         = -(face_id + mesh->n_b_faces + 1);
437       cell_face_count[c_id2] += 1;
438     }
439   }
440 
441   BFT_FREE(cell_face_count);
442 
443   /* Return values */
444 
445   *p_cell_faces_idx = cell_faces_idx;
446   *p_cell_faces_val = cell_faces_val;
447 
448 #if 0 && defined(DEBUG) && !defined(NDEBUG)
449  {
450    cs_lnum_t ipos, ival;
451    /* Print arrays */
452    bft_printf("dbg : cs_mesh_ret_cel_fac\n"
453               "nombre de cellules extraites = %d\n", extr_cell_size);
454    for (ipos = 0; ipos < extr_cell_size; ipos++) {
455      bft_printf("  cellule %d\n", ipos);
456      bft_printf("    cell_faces_idx[%d] = %d\n", ipos, cell_faces_idx[ipos]);
457      for (ival = cell_faces_idx[ipos]     - 1;
458           ival < cell_faces_idx[ipos + 1] - 1;
459           ival++)
460        bft_printf("      cell_faces_val[%d] = %d\n",
461                   ival, cell_faces_val[ival]);
462    }
463    bft_printf("  cell_faces_idx[%d] = %d\n", ipos, cell_faces_idx[ipos]);
464  }
465 #endif
466 
467 }
468 
469 /*----------------------------------------------------------------------------
470  * Build a nodal connectivity structure from a subset of a mesh's cells.
471  *
472  * The list of cells to extract is optional (if none is given, all cells
473  * faces are extracted by default); it does not need to be ordered on input,
474  * but is always ordered on exit (as cells are extracted by increasing number
475  * traversal, the list is reordered to ensure the coherency of the extracted
476  * mesh's link to its parent cells, built using this list).
477  *
478  * parameters:
479  *   mesh             <-- base mesh
480  *   name             <-- extracted mesh name
481  *   include_families <-- include family info if true
482  *   cell_list_size   <-- size of cell_list[] array
483  *   cell_list        <-> list of cells (1 to n), or NULL
484  *
485  * returns:
486  *   pointer to extracted nodal mesh
487  *----------------------------------------------------------------------------*/
488 
489 fvm_nodal_t  *
cs_mesh_connect_cells_to_nodal(const cs_mesh_t * mesh,const char * name,bool include_families,cs_lnum_t cell_list_size,cs_lnum_t cell_list[])490 cs_mesh_connect_cells_to_nodal(const cs_mesh_t  *mesh,
491                                const char       *name,
492                                bool              include_families,
493                                cs_lnum_t         cell_list_size,
494                                cs_lnum_t         cell_list[])
495 {
496   cs_lnum_t   face_id, cell_id;
497 
498   int         null_family = 0;
499   cs_lnum_t   extr_cell_count = 0, i_face_count = 0, b_face_count = 0;
500   cs_lnum_t  *extr_cell_idx = NULL;
501 
502   cs_lnum_t  *cell_face_idx = NULL;
503   cs_lnum_t  *cell_face_num = NULL;
504 
505   cs_lnum_t  *i_face_list= NULL, *b_face_list = NULL;
506 
507   cs_lnum_t  face_num_shift[3];
508   cs_lnum_t  *face_vertices_idx[2];
509   cs_lnum_t  *face_vertices_num[2];
510   cs_lnum_t  *polyhedra_faces = NULL;
511   const int  *cell_family = NULL;
512 
513   fvm_nodal_t  *extr_mesh;
514 
515   /* If a family has no attributes, it must be 1st by construction
516      (as families are sorted when merging duplicates) */
517 
518   if (mesh->n_families > 0) {
519     if (mesh->family_item[0] == 0)
520       null_family = 1;
521   }
522 
523   /* Check that the mesh contains face -> vertices connectivity */
524 
525   if (mesh->b_face_vtx_idx == NULL || mesh->i_face_vtx_idx == NULL)
526     bft_error(__FILE__, __LINE__, 0,
527               _("The main mesh does not contain any face -> vertices\n"
528                 "connectivity, necessary for the nodal connectivity\n"
529                 "reconstruction (cs_mesh_connect_cells_to_nodal)."));
530 
531   if (include_families) {
532     BFT_MALLOC(i_face_list, mesh->n_i_faces, cs_lnum_t);
533     BFT_MALLOC(b_face_list, mesh->n_b_faces, cs_lnum_t);
534   }
535 
536   /* Count the number of cells to convert */
537 
538   if (cell_list != NULL) {
539 
540     BFT_MALLOC(extr_cell_idx, mesh->n_cells_with_ghosts, cs_lnum_t);
541 
542     /* Initialize index as marker */
543 
544     for (cell_id = 0; cell_id < mesh->n_cells_with_ghosts; cell_id++)
545       extr_cell_idx[cell_id] = -1;
546     for (cell_id = 0; cell_id < cell_list_size; cell_id++) {
547       if (cell_list[cell_id] <= mesh->n_cells)
548         extr_cell_idx[cell_list[cell_id] - 1] = 1;
549     }
550 
551     /* Also mark faces bearing families */
552 
553     if (include_families) {
554 
555       for (face_id = 0; face_id < mesh->n_i_faces; face_id++) {
556         cs_lnum_t c_id_0 = mesh->i_face_cells[face_id][0];
557         cs_lnum_t c_id_1 = mesh->i_face_cells[face_id][1];
558         if (   (extr_cell_idx[c_id_0] == 1 || extr_cell_idx[c_id_1] == 1)
559             && (mesh->i_face_family[face_id] != null_family)) {
560           i_face_list[i_face_count++] = face_id + 1;
561         }
562       }
563       BFT_REALLOC(i_face_list, i_face_count, cs_lnum_t);
564 
565       for (face_id = 0; face_id < mesh->n_b_faces; face_id++) {
566         cs_lnum_t c_id = mesh->b_face_cells[face_id];
567         if (   (extr_cell_idx[c_id] == 1)
568             && (mesh->b_face_family[face_id] != null_family)) {
569           b_face_list[b_face_count++] = face_id + 1;
570         }
571       }
572       BFT_REALLOC(b_face_list, b_face_count, cs_lnum_t);
573 
574     }
575 
576     /* Convert marked ids to indexes (1 to n) and reconstruct values of
577        cell_list[] to ensure that it is ordered. */
578 
579     extr_cell_count = 0;
580     for (cell_id = 0; cell_id < mesh->n_cells; cell_id++) {
581       if (extr_cell_idx[cell_id] == 1) {
582         cell_list[extr_cell_count] = cell_id + 1;
583         extr_cell_idx[cell_id] = extr_cell_count++;
584       }
585     }
586 
587     assert(extr_cell_count <= cell_list_size);
588 
589   }
590   else {
591     extr_cell_count = CS_MIN(mesh->n_cells, cell_list_size);
592     extr_cell_idx = NULL;
593 
594     if (include_families && extr_cell_count > 0) {
595 
596       for (face_id = 0; face_id < mesh->n_i_faces; face_id++) {
597         cs_lnum_t c_id_0 = mesh->i_face_cells[face_id][0];
598         cs_lnum_t c_id_1 = mesh->i_face_cells[face_id][1];
599         if (   (c_id_0 < extr_cell_count || c_id_1 < extr_cell_count)
600             && (mesh->i_face_family[face_id] != null_family)) {
601           i_face_list[i_face_count++] = face_id + 1;
602         }
603       }
604       BFT_REALLOC(i_face_list, i_face_count, cs_lnum_t);
605 
606       for (face_id = 0; face_id < mesh->n_b_faces; face_id++) {
607         cs_lnum_t c_id = mesh->b_face_cells[face_id];
608         if (   (c_id < extr_cell_count)
609             && (mesh->b_face_family[face_id] != null_family)) {
610           b_face_list[b_face_count++] = face_id + 1;
611         }
612       }
613       BFT_REALLOC(b_face_list, b_face_count, cs_lnum_t);
614 
615     }
616   }
617 
618   /* Extract "cells -> faces" connectivity */
619 
620   cs_mesh_connect_get_cell_faces(mesh,
621                                  extr_cell_count,
622                                  extr_cell_idx,
623                                  &cell_face_idx,
624                                  &cell_face_num);
625 
626   if (extr_cell_idx != NULL)
627     BFT_FREE(extr_cell_idx);
628 
629   /* Build nodal connectivity */
630 
631   face_num_shift[0] = 0;
632   face_num_shift[1] = mesh->n_b_faces + face_num_shift[0];
633   face_num_shift[2] = mesh->n_i_faces + face_num_shift[1];
634 
635   face_vertices_idx[0] = mesh->b_face_vtx_idx;
636   face_vertices_idx[1] = mesh->i_face_vtx_idx;
637   face_vertices_num[0] = mesh->b_face_vtx_lst;
638   face_vertices_num[1] = mesh->i_face_vtx_lst;
639 
640   extr_mesh = fvm_nodal_create(name, 3);
641 
642   fvm_nodal_set_parent(extr_mesh, mesh);
643 
644   if (include_families)
645     cell_family = mesh->cell_family;
646 
647   fvm_nodal_from_desc_add_cells(extr_mesh,
648                                 extr_cell_count,
649                                 2,
650                                 face_num_shift,
651                                 (const cs_lnum_t **)face_vertices_idx,
652                                 (const cs_lnum_t **)face_vertices_num,
653                                 cell_face_idx,
654                                 cell_face_num,
655                                 cell_family,
656                                 cell_list,
657                                 &polyhedra_faces);
658 
659   /* Also add faces bearing families */
660 
661   if (include_families) {
662 
663     _add_faces_to_nodal(mesh,
664                         extr_mesh,
665                         1,
666                         include_families,
667                         0,
668                         b_face_count,
669                         NULL,
670                         b_face_list);
671 
672     _add_faces_to_nodal(mesh,
673                         extr_mesh,
674                         0,
675                         include_families,
676                         i_face_count,
677                         0,
678                         i_face_list,
679                         NULL);
680 
681     _order_nodal_faces(mesh, extr_mesh);
682 
683     BFT_FREE(i_face_list);
684     BFT_FREE(b_face_list);
685   }
686 
687   fvm_nodal_set_shared_vertices(extr_mesh, mesh->vtx_coord);
688   fvm_nodal_set_group_class_set(extr_mesh, mesh->class_defs);
689 
690   /* Free memory */
691 
692   BFT_FREE(polyhedra_faces);
693 
694   BFT_FREE(cell_face_idx);
695   BFT_FREE(cell_face_num);
696 
697   /* Sort cells by increasing global number */
698 
699   fvm_nodal_order_cells(extr_mesh, mesh->global_cell_num);
700   fvm_nodal_init_io_num(extr_mesh, mesh->global_cell_num, 3);
701 
702   /* Sort vertices by increasing global number */
703 
704   fvm_nodal_order_vertices(extr_mesh, mesh->global_vtx_num);
705   fvm_nodal_init_io_num(extr_mesh, mesh->global_vtx_num, 0);
706 
707   /* We are done */
708 
709   return extr_mesh;
710 }
711 
712 /*----------------------------------------------------------------------------
713  * Build a nodal connectivity structure from a subset of a mesh's faces.
714  *
715  * The lists of faces to extract are optional (if none is given, boundary
716  * faces are extracted by default); they do not need to be ordered on input,
717  * but they are always ordered on exit (as faces are extracted by increasing
718  * number traversal, the lists are reordered to ensure the coherency of
719  * the extracted mesh's link to its parent faces, built using these lists).
720  *
721  * parameters:
722  *   mesh             <-- base mesh
723  *   name             <-- extracted mesh name
724  *   include_families <-- include family info if true
725  *   i_face_list_size <-- size of i_face_list[] array
726  *   b_face_list_size <-- size of b_face_list[] array
727  *   i_face_list      <-> list of interior faces (1 to n), or NULL
728  *   b_face_list      <-> list of boundary faces (1 to n), or NULL
729  *
730  * returns:
731  *   pointer to extracted nodal mesh
732  *----------------------------------------------------------------------------*/
733 
734 fvm_nodal_t *
cs_mesh_connect_faces_to_nodal(const cs_mesh_t * mesh,const char * name,bool include_families,cs_lnum_t i_face_list_size,cs_lnum_t b_face_list_size,cs_lnum_t i_face_list[],cs_lnum_t b_face_list[])735 cs_mesh_connect_faces_to_nodal(const cs_mesh_t  *mesh,
736                                const char       *name,
737                                bool              include_families,
738                                cs_lnum_t         i_face_list_size,
739                                cs_lnum_t         b_face_list_size,
740                                cs_lnum_t         i_face_list[],
741                                cs_lnum_t         b_face_list[])
742 {
743   fvm_nodal_t  *extr_mesh = NULL;
744 
745   /* Check that the mesh contains face -> vertices connectivity */
746 
747   if (mesh->b_face_vtx_idx == NULL || mesh->i_face_vtx_idx == NULL)
748     bft_error(__FILE__, __LINE__, 0,
749               _("The main mesh does not contain any face -> vertices\n"
750                 "connectivity, necessary for the nodal connectivity\n"
751                 "reconstruction (cs_mesh_connect_faces_to_nodal)."));
752 
753   extr_mesh = fvm_nodal_create(name, 3);
754 
755   fvm_nodal_set_parent(extr_mesh, mesh);
756 
757   _add_faces_to_nodal(mesh,
758                       extr_mesh,
759                       -1,
760                       include_families,
761                       i_face_list_size,
762                       b_face_list_size,
763                       i_face_list,
764                       b_face_list);
765 
766   _order_nodal_faces(mesh, extr_mesh);
767 
768   fvm_nodal_set_shared_vertices(extr_mesh,
769                                 mesh->vtx_coord);
770 
771   /* Sort vertices by increasing global number */
772 
773   fvm_nodal_order_vertices(extr_mesh, mesh->global_vtx_num);
774   fvm_nodal_init_io_num(extr_mesh, mesh->global_vtx_num, 0);
775 
776   /* Define group classes */
777 
778   if (include_families)
779     fvm_nodal_set_group_class_set(extr_mesh, mesh->class_defs);
780 
781   /* We are done */
782 
783   return extr_mesh;
784 }
785 
786 /*----------------------------------------------------------------------------*/
787 /*!
788  * \brief Build a vertex to cell connectivity for marked vertices only.
789  *
790  * It is the caller's responsibility to free the v2c_idx and v2c arrays,
791  * which are allocated by this function.
792  *
793  * \param[in]    mesh      pointer to mesh structure
794  * \param[in]    v_flag    vertex selection flag (0: not selected, 1: selected)
795  * \param[out]   v2c_idx   vertex to cells index (size: mesh->n_vertices +1)
796  * \param[out]   v2c       vertex to cells
797  */
798 /*----------------------------------------------------------------------------*/
799 
800 void
cs_mesh_connect_vertices_to_cells(cs_mesh_t * mesh,const char v_flag[],cs_lnum_t ** v2c_idx,cs_lnum_t ** v2c)801 cs_mesh_connect_vertices_to_cells(cs_mesh_t    *mesh,
802                                   const char    v_flag[],
803                                   cs_lnum_t   **v2c_idx,
804                                   cs_lnum_t   **v2c)
805 {
806   cs_lnum_t n_vertices = mesh->n_vertices;
807 
808   /* Mark vertices which may be split (vertices lying on new boundary faces) */
809 
810   cs_lnum_t  *_v2c_idx;
811   BFT_MALLOC(_v2c_idx, n_vertices+1, cs_lnum_t);
812 
813   _v2c_idx[0] = 0;
814   for (cs_lnum_t i = 0; i < n_vertices; i++)
815     _v2c_idx[i+1] = 0;
816 
817   /* Now build vertex -> cells index
818      (which will contain duplicate entries at first) */
819 
820   for (cs_lnum_t f_id = 0; f_id < mesh->n_i_faces; f_id++) {
821     cs_lnum_t s_id = mesh->i_face_vtx_idx[f_id];
822     cs_lnum_t e_id = mesh->i_face_vtx_idx[f_id+1];
823     for (cs_lnum_t i = s_id; i < e_id; i++) {
824       cs_lnum_t vtx_id = mesh->i_face_vtx_lst[i];
825       if (v_flag[vtx_id] != 0) {
826         if (mesh->i_face_cells[f_id][0] > -1)
827           _v2c_idx[vtx_id + 1] += 1;
828         if (mesh->i_face_cells[f_id][1] > -1)
829           _v2c_idx[vtx_id + 1] += 1;
830       }
831     }
832   }
833 
834   for (cs_lnum_t f_id = 0; f_id < mesh->n_b_faces; f_id++) {
835     cs_lnum_t s_id = mesh->b_face_vtx_idx[f_id];
836     cs_lnum_t e_id = mesh->b_face_vtx_idx[f_id+1];
837     for (cs_lnum_t i = s_id; i < e_id; i++) {
838       cs_lnum_t vtx_id = mesh->b_face_vtx_lst[i];
839       if (v_flag[vtx_id] != 0)
840         _v2c_idx[vtx_id + 1] += 1;
841     }
842   }
843 
844   /* Transform count to index */
845 
846   for (cs_lnum_t i = 0; i < n_vertices; i++)
847     _v2c_idx[i+1] += _v2c_idx[i];
848 
849   /* Now define selected vertex->cell adjacency */
850 
851   cs_lnum_t *_v2c;
852   BFT_MALLOC(_v2c, _v2c_idx[n_vertices], cs_lnum_t);
853 
854   cs_lnum_t *v2c_count;
855   BFT_MALLOC(v2c_count, n_vertices, cs_lnum_t);
856   for (cs_lnum_t i = 0; i < n_vertices; i++)
857     v2c_count[i] = 0;
858 
859   for (cs_lnum_t f_id = 0; f_id < mesh->n_i_faces; f_id++) {
860     cs_lnum_t s_id = mesh->i_face_vtx_idx[f_id];
861     cs_lnum_t e_id = mesh->i_face_vtx_idx[f_id+1];
862     for (cs_lnum_t i = s_id; i < e_id; i++) {
863       cs_lnum_t vtx_id = mesh->i_face_vtx_lst[i];
864       if (v_flag[vtx_id] != 0) {
865         cs_lnum_t c_id_0 = mesh->i_face_cells[f_id][0];
866         cs_lnum_t c_id_1 = mesh->i_face_cells[f_id][1];
867         cs_lnum_t j = _v2c_idx[vtx_id] + v2c_count[vtx_id];
868         if (c_id_0 > -1) {
869           _v2c[j++] = c_id_0;
870           v2c_count[vtx_id] += 1;
871         }
872         if (c_id_1 > -1) {
873           _v2c[j++] = c_id_1;
874           v2c_count[vtx_id] += 1;
875         }
876       }
877     }
878   }
879 
880   for (cs_lnum_t f_id = 0; f_id < mesh->n_b_faces; f_id++) {
881     cs_lnum_t s_id = mesh->b_face_vtx_idx[f_id];
882     cs_lnum_t e_id = mesh->b_face_vtx_idx[f_id+1];
883     for (cs_lnum_t i = s_id; i < e_id; i++) {
884       cs_lnum_t vtx_id = mesh->b_face_vtx_lst[i];
885       if (v_flag[vtx_id] != 0) {
886         cs_lnum_t c_id_0 = mesh->b_face_cells[f_id];
887         cs_lnum_t j = _v2c_idx[vtx_id] + v2c_count[vtx_id];
888         _v2c[j] = c_id_0;
889         v2c_count[vtx_id] += 1;
890       }
891     }
892   }
893 
894   BFT_FREE(v2c_count);
895 
896   /* Order and compact adjacency array */
897 
898   cs_sort_indexed(n_vertices, _v2c_idx, _v2c);
899 
900   cs_lnum_t *tmp_v2c_idx = NULL;
901 
902   BFT_MALLOC(tmp_v2c_idx, n_vertices+1, cs_lnum_t);
903   memcpy(tmp_v2c_idx, _v2c_idx, (n_vertices+1)*sizeof(cs_lnum_t));
904 
905   cs_lnum_t k = 0;
906 
907   for (cs_lnum_t i = 0; i < n_vertices; i++) {
908     cs_lnum_t js = tmp_v2c_idx[i];
909     cs_lnum_t je = tmp_v2c_idx[i+1];
910     cs_lnum_t v2c_prev = -1;
911     _v2c_idx[i] = k;
912     for (cs_lnum_t j = js; j < je; j++) {
913       if (v2c_prev != _v2c[j]) {
914         _v2c[k++] = _v2c[j];
915         v2c_prev = _v2c[j];
916       }
917     }
918   }
919   _v2c_idx[n_vertices] = k;
920 
921   assert(_v2c_idx[n_vertices] <= tmp_v2c_idx[n_vertices]);
922 
923   BFT_FREE(tmp_v2c_idx);
924   BFT_REALLOC(_v2c, _v2c_idx[n_vertices], cs_lnum_t);
925 
926   *v2c_idx = _v2c_idx;
927   *v2c = _v2c;
928 }
929 
930 /*----------------------------------------------------------------------------*/
931 
932 END_C_DECLS
933