1 /*============================================================================
2  * Boundary zones handling.
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 <float.h>
35 #include <math.h>
36 #include <stdarg.h>
37 #include <stdio.h>
38 #include <stdlib.h>
39 #include <string.h>
40 
41 #if defined(HAVE_MPI)
42 #include <mpi.h>
43 #endif
44 
45 /*----------------------------------------------------------------------------
46  * Local headers
47  *----------------------------------------------------------------------------*/
48 
49 #include <ple_locator.h>
50 
51 #include "bft_mem.h"
52 #include "bft_error.h"
53 #include "bft_printf.h"
54 
55 #include "cs_base.h"
56 #include "cs_flag_check.h"
57 #include "cs_log.h"
58 #include "cs_map.h"
59 #include "cs_mesh.h"
60 #include "cs_mesh_location.h"
61 #include "cs_mesh_quantities.h"
62 #include "cs_parall.h"
63 #include "cs_parameters.h"
64 
65 /*----------------------------------------------------------------------------
66  *  Header for the current file
67  *----------------------------------------------------------------------------*/
68 
69 #include "cs_boundary_zone.h"
70 
71 /*----------------------------------------------------------------------------*/
72 
73 BEGIN_C_DECLS
74 
75 /*=============================================================================
76  * Additional doxygen documentation
77  *============================================================================*/
78 
79 /*!
80   \file cs_boundary_zone.c
81         Boundary zone handling.
82 */
83 
84 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
85 
86 /*=============================================================================
87  * Local Macro Definitions
88  *============================================================================*/
89 
90 /* Zone descriptor allocation block size */
91 
92 #define _CS_ZONE_S_ALLOC_SIZE       16
93 
94 /*=============================================================================
95  * Local Type Definitions
96  *============================================================================*/
97 
98 /*! Boundary zone ids */
99 
100 typedef struct {
101   cs_lnum_t         n_elts;   /*!< local number of associated elements */
102   const cs_lnum_t  *max_id;   /*!< local face ids, or NULL if trivial */
103 } cs_boundary_zone_id_t;
104 
105 /*============================================================================
106  *  Global variables
107  *============================================================================*/
108 
109 /* Zone definitions */
110 
111 static int  _n_zones = 0;
112 static int  _n_zones_max = 0;
113 static cs_zone_t   **_zones = NULL;
114 static cs_map_name_to_id_t   *_zone_map = NULL;
115 
116 /* Boundary zone id associated with boundary faces */
117 
118 static int *_zone_id = NULL;
119 
120 /* Optional, separate zone classification id for data extraction */
121 
122 static int  _max_zone_class_id = -1;
123 static int *_zone_class_id = NULL;
124 
125 /*============================================================================
126  * Prototypes for functions intended for use only by Fortran wrappers.
127  * (descriptions follow, with function bodies).
128  *============================================================================*/
129 
130 /*============================================================================
131  * Private function definitions
132  *============================================================================*/
133 
134 /*----------------------------------------------------------------------------*/
135 /*!
136  * \brief Return a pointer to a boundary zone based on its name if present.
137  *
138  * If no boundary zone of the given name is defined, NULL is returned.
139  *
140  * \param[in]  name  boundary zone name
141  *
142  * \return  pointer to (read only) zone structure, or NULL
143  */
144 /*----------------------------------------------------------------------------*/
145 
146 static cs_zone_t  *
_zone_by_name_try(const char * name)147 _zone_by_name_try(const char  *name)
148 {
149   cs_zone_t *z = NULL;
150   int id = cs_map_name_to_id_try(_zone_map, name);
151 
152   if (id > -1)
153     z = _zones[id];
154 
155   return z;
156 }
157 
158 /*----------------------------------------------------------------------------*/
159 /*!
160  * \brief Create a boundary zone.
161  *
162  * \param[in]  name  zone name
163  *
164  * \return   pointer to new zone.
165  */
166 /*----------------------------------------------------------------------------*/
167 
168 static cs_zone_t *
_zone_define(const char * name)169 _zone_define(const char  *name)
170 {
171   int zone_id = -1;
172   const char *addr_0 = NULL, *addr_1 = NULL;
173 
174   cs_zone_t *z = _zone_by_name_try(name);
175 
176   /* Check this name was not already used */
177 
178   if (z != NULL)
179     return z;
180 
181   /* Initialize if necessary */
182 
183   if (_zone_map == NULL)
184     _zone_map = cs_map_name_to_id_create();
185 
186   else
187     addr_0 = cs_map_name_to_id_reverse(_zone_map, 0);
188 
189   size_t l = 0;
190   if (name != NULL)
191     l = strlen(name);
192   if (l == 0)
193     bft_error(__FILE__, __LINE__, 0, _("Defining a zone requires a name."));
194 
195   /* Insert entry in map */
196 
197   zone_id = cs_map_name_to_id(_zone_map, name);
198 
199   /* Move name pointers of previous zones if necessary
200      (i.e. reallocation of map names array) */
201 
202   addr_1 = cs_map_name_to_id_reverse(_zone_map, 0);
203 
204   if (addr_1 != addr_0) {
205     ptrdiff_t addr_shift = addr_1 - addr_0;
206     for (int i = 0; i < zone_id; i++)
207       _zones[i]->name += addr_shift;
208   }
209 
210   if (zone_id == _n_zones)
211     _n_zones = zone_id + 1;
212 
213   /* Reallocate zones pointer if necessary */
214 
215   if (_n_zones > _n_zones_max) {
216     if (_n_zones_max == 0)
217       _n_zones_max = 8;
218     else
219       _n_zones_max *= 2;
220     BFT_REALLOC(_zones, _n_zones_max, cs_zone_t *);
221   }
222 
223   /* Allocate zones descriptor block if necessary
224      (to reduce fragmentation and improve locality of zone
225      descriptors, they are allocated in blocks) */
226 
227   int shift_in_alloc_block = zone_id % _CS_ZONE_S_ALLOC_SIZE;
228   if (shift_in_alloc_block == 0)
229     BFT_MALLOC(_zones[zone_id], _CS_ZONE_S_ALLOC_SIZE, cs_zone_t);
230   else
231     _zones[zone_id] =   _zones[zone_id - shift_in_alloc_block]
232                       + shift_in_alloc_block;
233 
234   /* Assign zone */
235 
236   z = _zones[zone_id];
237 
238   z->name = cs_map_name_to_id_reverse(_zone_map, zone_id);
239 
240   z->id = zone_id;
241   z->type = 0;
242 
243   z->location_id = 0;
244 
245   z->n_elts = 0;
246   z->elt_ids = NULL;
247 
248   z->time_varying = false;
249   z->allow_overlay = false;
250 
251   z->n_g_elts = 0;
252 
253   z->measure = -1.;
254   z->f_measure = -1.;
255   z->boundary_measure = -1.;
256   z->f_boundary_measure = -1.;
257 
258   for (int idim = 0; idim < 3; idim++)
259     z->cog[idim] = 0.;
260 
261   return z;
262 }
263 
264 /*----------------------------------------------------------------------------
265  * Add type flag info to the current position in the setup log.
266  *
267  * parameters:
268  *   type <-- type flag
269  *----------------------------------------------------------------------------*/
270 
271 static inline void
_log_type(int type)272 _log_type(int type)
273 {
274   if (type == 0)
275     return;
276 
277   cs_log_printf(CS_LOG_SETUP,
278                 _("    type:                       %d"), type);
279   cs_log_printf(CS_LOG_SETUP, "\n");
280 }
281 
282 /*----------------------------------------------------------------------------*/
283 /*!
284  * \brief Build zone class id array, based on zone class id.
285  */
286 /*----------------------------------------------------------------------------*/
287 
288 static void
_build_zone_class_id(void)289 _build_zone_class_id(void)
290 {
291   cs_mesh_t  *m = cs_glob_mesh;
292 
293   cs_lnum_t n_faces = m->n_b_faces;
294 
295   BFT_REALLOC(_zone_class_id, n_faces, int);
296 # pragma omp parallel for if (n_faces > CS_THR_MIN)
297   for (cs_lnum_t i = 0; i < n_faces; i++)
298     _zone_class_id[i] = _zone_id[i];
299 }
300 
301 /*----------------------------------------------------------------------------*/
302 /*!
303  * \brief Compute geometrical measure of a boundary zone (surface and perimeter)
304  *
305  * For time-varying zones, the associated mesh location is updated.
306  *
307  * \param[in]  mesh_modified indicate if mesh has been modified
308  * \param[in]  z             zone for which measures need to be computed
309  */
310 /*----------------------------------------------------------------------------*/
311 
312 static void
_boundary_zone_compute_metadata(bool mesh_modified,cs_zone_t * z)313 _boundary_zone_compute_metadata(bool       mesh_modified,
314                                 cs_zone_t *z)
315 {
316   /* We recompute values only if mesh is modified or zone is time varying.
317    * FIXME: For the moment, the boundary measure is not computed, but set to -1.
318    * to be improved in the future.
319    * Not that boundary measure of a (2D) surface should be the length of the
320    * perimeter
321    */
322   if (z->time_varying || mesh_modified) {
323 
324     cs_real_t *b_face_surf   = cs_glob_mesh_quantities->b_face_surf;
325     cs_real_t *b_f_face_surf = cs_glob_mesh_quantities->b_f_face_surf;
326     cs_real_3_t *face_cen = (cs_real_3_t *)cs_glob_mesh_quantities->b_face_cog;
327 
328     z->measure = 0.;
329     z->f_measure = 0.;
330     z->boundary_measure = -1.;
331     z->f_boundary_measure = -1.;
332 
333     for (int idim = 0; idim < 3; idim++)
334       z->cog[idim] = 0.;
335 
336     for (cs_lnum_t e_id = 0; e_id < z->n_elts; e_id++) {
337       cs_lnum_t f_id = z->elt_ids[e_id];
338       z->measure   += b_face_surf[f_id];
339       z->f_measure += b_f_face_surf[f_id];
340       for (int idim = 0; idim < 3; idim++)
341         z->cog[idim] += face_cen[f_id][idim] * b_face_surf[f_id];
342     }
343 
344     cs_real_t measures[7] = {z->measure, z->f_measure,
345                              z->boundary_measure, z->f_boundary_measure,
346                              z->cog[0], z->cog[1], z->cog[2]};
347 
348     cs_gnum_t n_g_elts = z->n_elts;
349     cs_parall_sum(1, CS_GNUM_TYPE, &n_g_elts);
350 
351     cs_parall_sum(7, CS_REAL_TYPE, measures);
352 
353     z->n_g_elts = n_g_elts;
354 
355     z->measure = measures[0];
356     z->f_measure = measures[1];
357     z->boundary_measure = measures[2];
358     z->f_boundary_measure = measures[3];
359 
360     /* Avoid a SIGFPE error */
361     if (fabs(measures[0]) > DBL_MIN)
362       for (int idim = 0; idim < 3; idim++)
363         z->cog[idim] = measures[4+idim] / measures[0];
364 
365   } /* Need to compute metadata */
366 
367 }
368 
369 /*============================================================================
370  * Fortran wrapper function definitions
371  *============================================================================*/
372 
373 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
374 
375 /*============================================================================
376  * Public function definitions
377  *============================================================================*/
378 
379 /*----------------------------------------------------------------------------*/
380 /*!
381  * \brief Initialize boundary zone structures.
382  *
383  * This defines a default boundary zone. This is the first function of
384  * the boundary zone handling functions which should be called, and it should
385  * only be called after \ref cs_mesh_location_initialize.
386  */
387 /*----------------------------------------------------------------------------*/
388 
389 void
cs_boundary_zone_initialize(void)390 cs_boundary_zone_initialize(void)
391 {
392   assert(_n_zones == 0);
393 
394   cs_mesh_location_set_explicit_ids(CS_MESH_LOCATION_BOUNDARY_FACES, true);
395 
396   const char *name = cs_mesh_location_get_name(CS_MESH_LOCATION_BOUNDARY_FACES);
397 
398   cs_zone_t *z = _zone_define(name);
399 
400   z->location_id = CS_MESH_LOCATION_BOUNDARY_FACES;
401 
402   z->type = 0;
403 
404   z->allow_overlay = true;
405 
406   assert(z->id == 0);
407 }
408 
409 /*----------------------------------------------------------------------------*/
410 /*!
411  * \brief Free all boundary zone structures.
412  */
413 /*----------------------------------------------------------------------------*/
414 
415 void
cs_boundary_zone_finalize(void)416 cs_boundary_zone_finalize(void)
417 {
418   BFT_FREE(_zone_class_id);
419   BFT_FREE(_zone_id);
420 
421   for (int i = 0; i < _n_zones; i++) {
422     if (i % _CS_ZONE_S_ALLOC_SIZE == 0)
423       BFT_FREE(_zones[i]);
424   }
425 
426   BFT_FREE(_zones);
427 
428   cs_map_name_to_id_destroy(&_zone_map);
429 
430   _n_zones = 0;
431   _n_zones_max = 0;
432 }
433 
434 /*----------------------------------------------------------------------------*/
435 /*!
436  * \brief Return number of boundary zones defined.
437  */
438 /*----------------------------------------------------------------------------*/
439 
440 int
cs_boundary_zone_n_zones(void)441 cs_boundary_zone_n_zones(void)
442 {
443   return _n_zones;
444 }
445 
446 /*----------------------------------------------------------------------------*/
447 /*!
448  * \brief Return number of boundary zones which may vary in time.
449  *
450  * \return  number of zones which may vary in time
451  */
452 /*----------------------------------------------------------------------------*/
453 
454 int
cs_boundary_zone_n_zones_time_varying(void)455 cs_boundary_zone_n_zones_time_varying(void)
456 {
457   int count = 0;
458 
459   for (int i = 0; i < _n_zones; i++) {
460     if (_zones[i]->time_varying)
461       count += 1;
462   }
463 
464   return count;
465 }
466 
467 /*----------------------------------------------------------------------------*/
468 /*!
469  * \brief Update association of a given private boundary zone with a mesh.
470  *
471  * For time-varying zones, the associated mesh location is updated.
472  *
473  * \param[in]  id  zone id
474  */
475 /*----------------------------------------------------------------------------*/
476 
477 void
cs_boundary_zone_build_private(int id)478 cs_boundary_zone_build_private(int  id)
479 {
480   if (id < 0 || id >= _n_zones)
481     bft_error(__FILE__, __LINE__, 0,
482               _("Boundary zone with id %d is not defined."), id);
483 
484   cs_zone_t *z = _zones[id];
485   if (! (z->type & CS_BOUNDARY_ZONE_PRIVATE))
486     return;
487 
488   cs_mesh_t  *m = cs_glob_mesh;
489   cs_mesh_location_build(m, z->location_id);
490 
491   z->n_elts = cs_mesh_location_get_n_elts(z->location_id)[0];
492   z->elt_ids = cs_mesh_location_get_elt_ids(z->location_id);
493 }
494 
495 /*----------------------------------------------------------------------------*/
496 /*!
497  * \brief Update association of boundary zones with a mesh.
498  *
499  * For time-varying zones, the associated mesh location is updated.
500  *
501  * \param[in]  mesh_modified  indicate if mesh has been modified
502  */
503 /*----------------------------------------------------------------------------*/
504 
505 void
cs_boundary_zone_build_all(bool mesh_modified)506 cs_boundary_zone_build_all(bool  mesh_modified)
507 {
508   cs_mesh_t  *m = cs_glob_mesh;
509   bool has_time_varying = false;
510 
511   /* update zone lists */
512 
513   for (int i = 0; i < _n_zones; i++) {
514     cs_zone_t *z = _zones[i];
515     if (z->time_varying) {
516       cs_mesh_location_build(m, z->location_id);
517       if (! (z->type & CS_BOUNDARY_ZONE_PRIVATE))
518         has_time_varying = true;
519     }
520     z->n_elts = cs_mesh_location_get_n_elts(z->location_id)[0];
521     z->elt_ids = cs_mesh_location_get_elt_ids(z->location_id);
522   }
523 
524   /* Assign maximum zone id and check for overlap errors
525      (start with zone 1, as 0 is default) */
526 
527   if (mesh_modified)
528     BFT_REALLOC(_zone_id, m->n_b_faces, int);
529 
530   if (mesh_modified || has_time_varying) {
531 
532     cs_lnum_t n_faces = m->n_b_faces;
533 
534 #   pragma omp parallel for if (n_faces > CS_THR_MIN)
535     for (cs_lnum_t i = 0; i <n_faces; i++)
536       _zone_id[i] = 0;
537 
538     int overlap_error[2] = {_n_zones, _n_zones};
539 
540     for (int i = 1; i < _n_zones; i++) {
541       cs_zone_t *z = _zones[i];
542       if (z->type & CS_BOUNDARY_ZONE_PRIVATE)
543         continue;
544       for (cs_lnum_t j = 0; j < z->n_elts; j++) {
545         cs_lnum_t f_id = z->elt_ids[j];
546         int z_id_prev = _zone_id[f_id];
547         if (z_id_prev == 0)
548           _zone_id[f_id] = z->id;
549         else if (_zones[z_id_prev]->allow_overlay)
550           _zone_id[f_id] = z->id;
551         else if (overlap_error[0] == _n_zones) {
552           overlap_error[0] = z_id_prev;
553           overlap_error[1] = z->id;
554           break;
555         }
556       }
557     }
558 
559     cs_parall_min(2, CS_INT_TYPE, overlap_error);
560 
561     if (overlap_error[0] < _n_zones) {
562 
563       for (int i = 1; i < _n_zones; i++) {
564         cs_zone_t *z = _zones[i];
565         if (z->type & CS_BOUNDARY_ZONE_PRIVATE)
566           continue;
567         for (cs_lnum_t j = 0; j < z->n_elts; j++) {
568           cs_lnum_t f_id = z->elt_ids[j];
569           int z_id_prev = CS_ABS(_zone_id[f_id]);
570           if (z_id_prev == 0)
571             _zone_id[f_id] = z->id;
572           else if (   _zones[z_id_prev]->allow_overlay
573                    && _zone_id[f_id] > 0)
574             _zone_id[f_id] = z->id;
575           else
576             _zone_id[f_id] = -z->id;
577         }
578       }
579 
580       cs_flag_check_error_info(_("face with forbidden zone overlap"),
581                                _("zone id"),
582                                _("zone_id"),
583                                _("Faces with zone error"),
584                                _("Faces with valid zones"),
585                                CS_MESH_LOCATION_BOUNDARY_FACES,
586                                0, /* min_flag */
587                                _zone_id);
588 
589       int i0 = overlap_error[0], i1 = overlap_error[1];
590 
591       bft_error(__FILE__, __LINE__, 0,
592                 _("Boundary zone %i (\"%s\") contains at least\n"
593                   "one face already marked with zone id %d (\"%s\").\n\n"
594                   "Check definitions or allow overlays for this zone."),
595                 i1, _zones[i1]->name, i0, _zones[i0]->name);
596 
597     }
598 
599     if (_max_zone_class_id > -1)
600       _build_zone_class_id();
601 
602     /* Compute or update zone geometrical measures */
603     for (int i = 0; i < _n_zones; i++) {
604       cs_zone_t *z = _zones[i];
605       _boundary_zone_compute_metadata(mesh_modified, z);
606     }
607   }
608 }
609 
610 /*----------------------------------------------------------------------------*/
611 /*!
612  * \brief Define a new boundary zone using a selection criteria string.
613  *
614  * \param[in]  name       name of location to define
615  * \param[in]  criteria   selection criteria for associated elements
616  * \param[in]  type_flag  mask of zone category values
617  *
618  * \return  id of newly defined boundary zone
619  */
620 /*----------------------------------------------------------------------------*/
621 
622 int
cs_boundary_zone_define(const char * name,const char * criteria,int type_flag)623 cs_boundary_zone_define(const char  *name,
624                         const char  *criteria,
625                         int          type_flag)
626 {
627   if (criteria == NULL)
628     bft_error(__FILE__, __LINE__, 0,
629               _("%s: selection criteria string must be non-null."),
630               __func__);
631 
632   cs_zone_t *z = _zone_define(name);
633 
634   if (strcmp(criteria, "all[]"))
635     z->location_id = cs_mesh_location_add(name,
636                                           CS_MESH_LOCATION_BOUNDARY_FACES,
637                                           criteria);
638   else
639     z->location_id = CS_MESH_LOCATION_BOUNDARY_FACES;
640 
641   z->type = type_flag;
642 
643   return z->id;
644 }
645 
646 /*----------------------------------------------------------------------------*/
647 /*!
648  * \brief Define a new mesh location with an associated selection function.
649  *
650  * So as to define a subset of mesh entities of a given type, a pointer
651  * to a selection function may be given.
652  *
653  * This requires more programming but allows finer control than selection
654  * criteria, as the function has access to the complete mesh structure.
655  *
656  * \param[in]  name  name of location to define
657  * \param[in]  func  pointer to selection function for associated elements
658  * \param[in, out]  input  pointer to optional (untyped) value
659  *                         or structure.
660  * \param[in]  type_flag  mask of zone category values
661  *
662  * \return  id of newly defined created mesh location
663  */
664 /*----------------------------------------------------------------------------*/
665 
666 int
cs_boundary_zone_define_by_func(const char * name,cs_mesh_location_select_t * func,void * input,int type_flag)667 cs_boundary_zone_define_by_func(const char                 *name,
668                                 cs_mesh_location_select_t  *func,
669                                 void                       *input,
670                                 int                         type_flag)
671 {
672   if (func == NULL)
673     bft_error(__FILE__, __LINE__, 0,
674               _("%s: selection function pointer must be non-null."),
675               __func__);
676 
677   cs_zone_t *z = _zone_define(name);
678 
679   z->location_id = cs_mesh_location_add_by_func(name,
680                                                 CS_MESH_LOCATION_BOUNDARY_FACES,
681                                                 func,
682                                                 input);
683 
684   z->type = type_flag;
685 
686   return z->id;
687 }
688 
689 /*----------------------------------------------------------------------------*/
690 /*!
691  * \brief Return a pointer to a boundary zone based on its id.
692  *
693  * This function requires that a boundary zone of the given id is defined.
694  *
695  * \param[in]  id   zone id
696  *
697  * \return  pointer to the boundary zone structure
698  */
699 /*----------------------------------------------------------------------------*/
700 
701 const cs_zone_t  *
cs_boundary_zone_by_id(int id)702 cs_boundary_zone_by_id(int  id)
703 {
704   if (id > -1 && id < _n_zones)
705     return _zones[id];
706   else {
707     bft_error(__FILE__, __LINE__, 0,
708               _("Boundary zone with id %d is not defined."), id);
709     return NULL;
710   }
711 }
712 
713 /*----------------------------------------------------------------------------*/
714 /*!
715  * \brief Return a pointer to a boundary zone based on its name if present.
716  *
717  * This function requires that a boundary zone of the given name is defined.
718  *
719  * \param[in]  name  boundary zone name
720  *
721  * \return  pointer to (read-only) zone structure
722  */
723 /*----------------------------------------------------------------------------*/
724 
725 const cs_zone_t  *
cs_boundary_zone_by_name(const char * name)726 cs_boundary_zone_by_name(const char  *name)
727 {
728   int id = cs_map_name_to_id_try(_zone_map, name);
729 
730   if (id > -1)
731     return _zones[id];
732   else {
733     bft_error(__FILE__, __LINE__, 0,
734               _("Boundary zone \"%s\" is not defined."), name);
735     return NULL;
736   }
737 }
738 
739 /*----------------------------------------------------------------------------*/
740 /*!
741  * \brief Return a pointer to a boundary zone based on its name if present.
742  *
743  * If no boundary zone of the given name is defined, NULL is returned.
744  *
745  * \param[in]  name  boundary zone name
746  *
747  * \return  pointer to (read only) zone structure, or NULL
748  */
749 /*----------------------------------------------------------------------------*/
750 
751 const cs_zone_t  *
cs_boundary_zone_by_name_try(const char * name)752 cs_boundary_zone_by_name_try(const char  *name)
753 {
754   const cs_zone_t *z = NULL;
755   int id = cs_map_name_to_id_try(_zone_map, name);
756 
757   if (id > -1)
758     z = _zones[id];
759 
760   return z;
761 }
762 
763 /*----------------------------------------------------------------------------*/
764 /*!
765  * \brief Set type flag for a given boundary zone.
766  *
767  * \param[in]  id         boundary zone id
768  * \param[in]  type_flag  volume zone type flag
769  */
770 /*----------------------------------------------------------------------------*/
771 
772 void
cs_boundary_zone_set_type(int id,int type_flag)773 cs_boundary_zone_set_type(int   id,
774                           int   type_flag)
775 {
776   const cs_zone_t *z0 = cs_boundary_zone_by_id(id);
777 
778   _zones[z0->id]->type |= type_flag;
779 }
780 
781 /*----------------------------------------------------------------------------*/
782 /*!
783  * \brief Set time varying behavior for a given boundary zone.
784  *
785  * \param[in]  id            boundary zone id
786  * \param[in]  time_varying  true if the zone's definition varies in time
787  */
788 /*----------------------------------------------------------------------------*/
789 
790 void
cs_boundary_zone_set_time_varying(int id,bool time_varying)791 cs_boundary_zone_set_time_varying(int   id,
792                                   bool  time_varying)
793 {
794   const cs_zone_t *z0 = cs_boundary_zone_by_id(id);
795 
796   _zones[z0->id]->time_varying = time_varying;
797 }
798 
799 /*----------------------------------------------------------------------------*/
800 /*!
801  * \brief Set overlay behavior for a given boundary zone.
802  *
803  * \param[in]  id             boundary zone id
804  * \param[in]  allow_overlay  true if the zone may be overlayed by another
805  */
806 /*----------------------------------------------------------------------------*/
807 
808 void
cs_boundary_zone_set_overlay(int id,bool allow_overlay)809 cs_boundary_zone_set_overlay(int   id,
810                              bool  allow_overlay)
811 {
812   const cs_zone_t *z0 = cs_boundary_zone_by_id(id);
813 
814   _zones[z0->id]->allow_overlay = allow_overlay;
815 }
816 
817 /*----------------------------------------------------------------------------*/
818 /*!
819  * \brief Return pointer to zone id associated with each boundary face.
820  *
821  * In case of overlayed zones, the highest zone id associated with
822  * a given face is given. Private (automatically defined) zones
823  * are excluded from this definition.
824  */
825 /*----------------------------------------------------------------------------*/
826 
827 const int *
cs_boundary_zone_face_zone_id(void)828 cs_boundary_zone_face_zone_id(void)
829 {
830   return (const int *)_zone_id;
831 }
832 
833 /*----------------------------------------------------------------------------*/
834 /*!
835  * \brief Print info relative to a given boundary zone to log file.
836  *
837  * \param[in]  z   pointer to boundary zone structure
838  */
839 /*----------------------------------------------------------------------------*/
840 
841 void
cs_boundary_zone_log_info(const cs_zone_t * z)842 cs_boundary_zone_log_info(const cs_zone_t  *z)
843 {
844   if (z == NULL)
845     return;
846 
847   /* Global indicators */
848   /*-------------------*/
849 
850   cs_log_printf(CS_LOG_SETUP,
851                 _("\n"
852                   "  Zone: \"%s\"\n"
853                   "    id:                         %d\n"),
854                 z->name, z->id);
855 
856 
857   _log_type(z->type);
858 
859   cs_log_printf(CS_LOG_SETUP,
860                 _("    location_id:                %d\n"),
861                 z->location_id);
862 
863   if (z->time_varying)
864     cs_log_printf(CS_LOG_SETUP, _("    time varying\n"));
865   if (z->type & CS_BOUNDARY_ZONE_PRIVATE)
866     cs_log_printf(CS_LOG_SETUP, _("    private (automatic)\n"));
867   else if (z->allow_overlay)
868     cs_log_printf(CS_LOG_SETUP, _("    allow overlay\n"));
869 
870   cs_mesh_location_def_t ml_def =
871     cs_mesh_location_get_definition_method(z->location_id);
872 
873   if (ml_def == CS_MESH_LOCATION_DEF_SELECTION_STR) {
874     const char *sel_str = cs_mesh_location_get_selection_string(z->location_id);
875 
876     cs_log_printf(CS_LOG_SETUP,
877                   _("    selection criteria:         \"%s\"\n"),
878                   sel_str);
879   }
880   else if (ml_def == CS_MESH_LOCATION_DEF_SELECTION_FUNC) {
881     cs_mesh_location_select_t *select_fp
882       = cs_mesh_location_get_selection_function(z->location_id);
883 
884     cs_log_printf(CS_LOG_SETUP,
885                   _("    selection function:         %p\n"),
886                   (void *)select_fp);
887   }
888   else if (ml_def == CS_MESH_LOCATION_DEF_UNION) {
889     int n_sub_ids = cs_mesh_location_get_n_sub_ids(z->location_id);
890     int *sub_ids  = cs_mesh_location_get_sub_ids(z->location_id);
891 
892     bool is_complement = cs_mesh_location_is_complement(z->location_id);
893     if (!is_complement)
894       cs_log_printf(CS_LOG_SETUP,
895                     _("    Union of %d mesh locations:\n"),
896                     n_sub_ids);
897     else
898       cs_log_printf(CS_LOG_SETUP,
899                     _("    Complement of %d mesh locations:\n"),
900                     n_sub_ids);
901     for (int iid = 0; iid < n_sub_ids; iid++) {
902       cs_log_printf(CS_LOG_SETUP,
903                     _("      sub-location %d/%d\n"),
904                     iid+1,
905                     n_sub_ids);
906 
907       int sub_ml_id = sub_ids[iid];
908 
909       cs_log_printf(CS_LOG_SETUP,
910                     _("        location_id:            %d\n"),
911                     sub_ml_id);
912 
913       cs_mesh_location_def_t sub_ml_def =
914         cs_mesh_location_get_definition_method(sub_ml_id);
915 
916       if (sub_ml_def == CS_MESH_LOCATION_DEF_SELECTION_STR) {
917         const char *sub_sel_str =
918           cs_mesh_location_get_selection_string(sub_ml_id);
919 
920         cs_log_printf(CS_LOG_SETUP,
921                       _("        selection criteria:     \"%s\"\n"),
922                       sub_sel_str);
923       }
924       else if (sub_ml_def == CS_MESH_LOCATION_DEF_SELECTION_FUNC) {
925         cs_mesh_location_select_t *sub_select_fp
926           = cs_mesh_location_get_selection_function(sub_ml_id);
927 
928         cs_log_printf(CS_LOG_SETUP,
929                       _("        selection function:     %p\n"),
930                       (void *)sub_select_fp);
931       }
932     }
933   }
934 }
935 
936 /*----------------------------------------------------------------------------*/
937 /*!
938  * \brief Log setup information relative to defined boundary zones.
939  */
940 /*----------------------------------------------------------------------------*/
941 
942 void
cs_boundary_zone_log_setup(void)943 cs_boundary_zone_log_setup(void)
944 {
945   if (_n_zones == 0)
946     return;
947 
948   cs_log_printf(CS_LOG_SETUP,
949                 _("\n"
950                   "Boundary zones\n"
951                   "--------------\n"));
952 
953   for (int i = 0; i < _n_zones; i++)
954     cs_boundary_zone_log_info(_zones[i]);
955 }
956 
957 /*----------------------------------------------------------------------------*/
958 /*!
959  * \brief Return number of boundary zones associated with a
960  *        given zone flag.
961  *
962  * Private (automatic) zones are excluded from this count.
963  *
964  * \param[in]  type_flag  flag to compare to zone type
965  *
966  * \return  number of zones matching the given type flag
967  */
968 /*----------------------------------------------------------------------------*/
969 
970 int
cs_boundary_zone_n_type_zones(int type_flag)971 cs_boundary_zone_n_type_zones(int  type_flag)
972 {
973   int count = 0;
974 
975   for (int i = 0; i < _n_zones; i++) {
976     if (    (_zones[i]->type & type_flag)
977         && !(_zones[i]->type & CS_BOUNDARY_ZONE_PRIVATE))
978       count += 1;
979   }
980 
981   return count;
982 }
983 
984 /*----------------------------------------------------------------------------*/
985 /*!
986  * \brief Get pointer to optional boundary face class ids.
987  *
988  * For each boundary face, a specific output (logging and postprocessing)
989  * class id may be assigned. This allows realizing logging, postprocessing,
990  * or otherwise extracting data based on this class.
991  *
992  * Using this function at a given point indicates that user-defined class
993  * ids will be used. The face class ids are initially equal to the
994  * face zone ids, but may be modified by the user.
995  *
996  * In the presence of a time-varying mesh or boundary zones, the face
997  * class ids will be reset to the zone ids, so it may be necessary to
998  * update the user definitions.
999  *
1000  * The class id values are arbitrarily chosen by the user, but must be
1001  * positive integers; numbers do not need to be contiguous, but very high
1002  * numbers may also lead to higher memory consumption.
1003  *
1004  * \return  pointer to array of boundary face output zone ids;
1005  */
1006 /*----------------------------------------------------------------------------*/
1007 
1008 int *
cs_boundary_zone_face_class_id(void)1009 cs_boundary_zone_face_class_id(void)
1010 {
1011   if (_max_zone_class_id < 0)
1012     _build_zone_class_id();
1013 
1014   _max_zone_class_id = 0;
1015 
1016   return _zone_class_id;
1017 }
1018 
1019 /*----------------------------------------------------------------------------*/
1020 /*!
1021  * \brief Update boundary face output class ids if present.
1022  *
1023  * Face class ids lower than 0 are replaced by the matching face zone id.
1024  */
1025 /*----------------------------------------------------------------------------*/
1026 
1027 void
cs_boundary_zone_update_face_class_id(void)1028 cs_boundary_zone_update_face_class_id(void)
1029 {
1030   int max_class = -1;
1031 
1032   if (_max_zone_class_id >= 0) {
1033 
1034     const cs_lnum_t n_b_faces = cs_glob_mesh->n_b_faces;
1035 
1036     for (cs_lnum_t i = 0; i < n_b_faces; i++) {
1037       int o_id = _zone_class_id[i];
1038       if (o_id < 0) {
1039         o_id = _zone_id[i];
1040         _zone_class_id[i] = o_id;
1041       }
1042       if (o_id >= max_class)
1043         max_class = o_id;
1044     }
1045 
1046   }
1047 
1048   cs_parall_max(1, CS_INT_TYPE, &max_class);
1049 
1050   _max_zone_class_id = max_class;
1051 }
1052 
1053 /*----------------------------------------------------------------------------*/
1054 /*!
1055  * \brief Get read pointer to optional boundary face class or zone ids.
1056  *
1057  * If no face classes have been defined by \ref cs_boundary_zone_face_class_id
1058  * the boundary face zone id is returned instead.
1059  *
1060  * \return  pointer to array of boundary face output zone ids;
1061  */
1062 /*----------------------------------------------------------------------------*/
1063 
1064 const int *
cs_boundary_zone_face_class_or_zone_id(void)1065 cs_boundary_zone_face_class_or_zone_id(void)
1066 {
1067   const int *retval = _zone_class_id;
1068 
1069   if (retval == NULL)
1070     retval = _zone_id;
1071 
1072   return retval;
1073 }
1074 
1075 /*----------------------------------------------------------------------------*/
1076 /*!
1077  * \brief Return the maximum defined face class or zone id.
1078  *
1079  * This value is valid only if \ref cs_boundary_zone_update_face_class_id
1080  * \return  maximum face class or zone id;
1081  */
1082 /*----------------------------------------------------------------------------*/
1083 
1084 int
cs_boundary_zone_max_class_or_zone_id(void)1085 cs_boundary_zone_max_class_or_zone_id(void)
1086 {
1087   int retval = _n_zones-1;
1088 
1089   if (_max_zone_class_id > retval)
1090     retval = _max_zone_class_id;
1091 
1092   return retval;
1093 }
1094 
1095 /*----------------------------------------------------------------------------*/
1096 /*!
1097  * \brief Print boundary zones information to listing file
1098  */
1099 /*----------------------------------------------------------------------------*/
1100 
1101 void
cs_boundary_zone_print_info(void)1102 cs_boundary_zone_print_info(void)
1103 {
1104 
1105   bft_printf("\n");
1106   bft_printf(" --- Information on boundary zones\n");
1107   cs_real_t *b_face_surf   = cs_glob_mesh_quantities->b_face_surf;
1108   cs_real_t *b_f_face_surf = cs_glob_mesh_quantities->b_f_face_surf;
1109 
1110   for (int i = 0; i < _n_zones; i++) {
1111     cs_zone_t *z = _zones[i];
1112     bft_printf(_("  Boundary zone \"%s\"\n"
1113                  "    id              = %d\n"
1114                  "    Number of faces = %llu\n"
1115                  "    Surface         = %1.5g\n"
1116                  "    Center of gravity = (%1.5g, %1.5g, %1.5g)\n"),
1117                z->name,
1118                z->id,
1119                (unsigned long long)z->n_g_elts,
1120                z->measure,
1121                z->cog[0], z->cog[1], z->cog[2]);
1122     /* Only log fluid fluid when different to surface */
1123     if (b_f_face_surf != b_face_surf && b_f_face_surf != NULL)
1124       bft_printf(_("    Fluid surface   = %1.5g\n"),
1125                  z->f_measure);
1126 
1127     if (z->boundary_measure < 0.) {
1128       bft_printf(_("    Perimeter       = -1 (not computed)\n"));
1129       /* Only log fluid fluid when different to surface */
1130       if (b_f_face_surf != b_face_surf && b_f_face_surf != NULL)
1131         bft_printf(_("    Fluid perimeter = -1 (not computed)\n"));
1132     }
1133     else {
1134       bft_printf(_("    Perimeter       = %1.5g\n"),
1135                  z->boundary_measure);
1136       /* Only log fluid fluid when different to surface */
1137       if (b_f_face_surf != b_face_surf && b_f_face_surf != NULL)
1138         bft_printf(_("    Fluid perimeter = %1.5g\n"),
1139                    z->f_boundary_measure);
1140     }
1141   }
1142 
1143   bft_printf_flush();
1144 }
1145 
1146 /*----------------------------------------------------------------------------*/
1147 
1148 END_C_DECLS
1149