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