1 /*============================================================================
2 * Manage the definition/setting of properties
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 <ctype.h>
35 #include <float.h>
36 #include <stdio.h>
37 #include <stdlib.h>
38 #include <string.h>
39 #include <math.h>
40
41 /*----------------------------------------------------------------------------
42 * Local headers
43 *----------------------------------------------------------------------------*/
44
45 #include <bft_mem.h>
46
47 #include "cs_defs.h"
48 #include "cs_log.h"
49 #include "cs_param_cdo.h"
50 #include "cs_reco.h"
51 #include "cs_volume_zone.h"
52 #include "cs_xdef_eval.h"
53
54 /*----------------------------------------------------------------------------
55 * Header for the current file
56 *----------------------------------------------------------------------------*/
57
58 #include "cs_property.h"
59
60 /*----------------------------------------------------------------------------*/
61
62 BEGIN_C_DECLS
63
64 /*=============================================================================
65 * Local Macro definitions and structure definitions
66 *============================================================================*/
67
68 #define CS_PROPERTY_DBG 0
69
70 /*============================================================================
71 * Private variables
72 *============================================================================*/
73
74 static const char _err_empty_pty[] =
75 " Stop setting an empty cs_property_t structure.\n"
76 " Please check your settings.\n";
77
78 /* Pointer to shared structures (owned by a cs_domain_t structure) */
79 static const cs_cdo_quantities_t *cs_cdo_quant;
80 static const cs_cdo_connect_t *cs_cdo_connect;
81
82 static int _n_properties = 0;
83 static int _n_max_properties = 0;
84 static cs_property_t **_properties = NULL;
85
86 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
87
88 /*============================================================================
89 * Private function prototypes
90 *============================================================================*/
91
92 /*----------------------------------------------------------------------------*/
93 /*!
94 * \brief Check if the settings are valid
95 *
96 * \param[in] tens values of the tensor
97 */
98 /*----------------------------------------------------------------------------*/
99
100 static inline bool
_is_tensor_symmetric(const cs_real_3_t * tens)101 _is_tensor_symmetric(const cs_real_3_t *tens)
102 {
103 if ((tens[0][1] - tens[1][0]) > cs_math_zero_threshold ||
104 (tens[0][2] - tens[2][0]) > cs_math_zero_threshold ||
105 (tens[1][2] - tens[2][1]) > cs_math_zero_threshold)
106 return false;
107 else
108 return true;
109 }
110
111 /*----------------------------------------------------------------------------*/
112 /*!
113 * \brief Print a tensor
114 *
115 * \param[in] tensor values of the tensor
116 */
117 /*----------------------------------------------------------------------------*/
118
119 static inline void
_print_tensor(const cs_real_3_t * tensor)120 _print_tensor(const cs_real_3_t *tensor)
121 {
122 cs_log_printf(CS_LOG_DEFAULT,
123 "\n Tensor property: | % 8.4e % 8.4e % 8.4e |\n"
124 " | % 8.4e % 8.4e % 8.4e |\n"
125 " | % 8.4e % 8.4e % 8.4e |\n",
126 tensor[0][0], tensor[0][1], tensor[0][2],
127 tensor[1][0], tensor[1][1], tensor[1][2],
128 tensor[2][0], tensor[2][1], tensor[2][2]);
129 }
130
131 /*----------------------------------------------------------------------------*/
132 /*!
133 * \brief Add a new definition to a cs_property_t structure
134 * Sanity checks on the settings related to this definition.
135
136 * \param[in, out] pty pointer to a cs_property_t structure
137 *
138 * \return the new definition id
139 */
140 /*----------------------------------------------------------------------------*/
141
142 inline static int
_add_new_def(cs_property_t * pty)143 _add_new_def(cs_property_t *pty)
144 {
145 int new_id = pty->n_definitions;
146
147 pty->n_definitions += 1;
148 BFT_REALLOC(pty->defs, pty->n_definitions, cs_xdef_t *);
149 BFT_REALLOC(pty->get_eval_at_cell, pty->n_definitions,
150 cs_xdef_eval_t *);
151 BFT_REALLOC(pty->get_eval_at_cell_cw, pty->n_definitions,
152 cs_xdef_cw_eval_t *);
153
154 return new_id;
155 }
156
157 /*----------------------------------------------------------------------------*/
158 /*!
159 * \brief Check if the settings are valid and then invert a tensor
160 *
161 * \param[in, out] tens values of the tensor
162 * \param[in] type type of property to deal with
163 */
164 /*----------------------------------------------------------------------------*/
165
166 static void
_invert_tensor(cs_real_3_t * tens,cs_property_type_t type)167 _invert_tensor(cs_real_3_t *tens,
168 cs_property_type_t type)
169 {
170 #if defined(DEBUG) && !defined(NDEBUG) && CS_PROPERTY_DBG > 0 /* Sanity check */
171 bool is_ok =true;
172 for (int k = 0; k < 3; k++)
173 if (fabs(tensor[k][k]) < cs_math_zero_threshold)
174 is_ok = false;
175
176 if (is_ok)
177 _is_tensor_symmetric(tens);
178
179 if (!is_ok) {
180 _print_tensor((const cs_real_t (*)[3])tens);
181 bft_error(__FILE__, __LINE__, 0,
182 " %s: A problem has been detected during the definition of the"
183 " property %s in the cell %d.\n", __func__, pty->name, c_id);
184 }
185 #endif
186
187 if (type & CS_PROPERTY_ISO || type & CS_PROPERTY_ORTHO)
188 for (int k = 0; k < 3; k++)
189 tens[k][k] = 1.0/tens[k][k];
190
191 else { /* anisotropic (sym. or not) */
192
193 cs_real_33_t invmat;
194
195 cs_math_33_inv_cramer((const cs_real_3_t (*))tens, invmat);
196 for (int k = 0; k < 3; k++)
197 for (int l = 0; l < 3; l++)
198 tens[k][l] = invmat[k][l];
199
200 }
201 }
202
203 /*----------------------------------------------------------------------------*/
204 /*!
205 * \brief Compute the value of a property at the cell center
206 *
207 * \param[in] c_id id of the current cell
208 * \param[in] t_eval physical time at which one evaluates the term
209 * \param[in] pty pointer to a cs_property_t structure
210 *
211 * \return the value of the property for the given cell
212 */
213 /*----------------------------------------------------------------------------*/
214
215 static inline cs_real_t
_get_cell_value(cs_lnum_t c_id,cs_real_t t_eval,const cs_property_t * pty)216 _get_cell_value(cs_lnum_t c_id,
217 cs_real_t t_eval,
218 const cs_property_t *pty)
219 {
220 int def_id = 0;
221 if (pty->n_definitions > 1) {
222 assert(pty->def_ids != NULL);
223 def_id = pty->def_ids[c_id];
224 assert(def_id > -1);
225 }
226
227 assert(pty->get_eval_at_cell[def_id] != NULL);
228
229 cs_xdef_t *def = pty->defs[def_id];
230 cs_real_t result = 0;
231
232 pty->get_eval_at_cell[def_id](1, &c_id, true, /* dense output */
233 cs_glob_mesh,
234 cs_cdo_connect,
235 cs_cdo_quant,
236 t_eval,
237 def->context,
238 &result);
239
240 return result;
241 }
242
243 /*----------------------------------------------------------------------------*/
244 /*!
245 * \brief Compute the value of a property at the cell center
246 * Version using a cs_cell_mesh_t structure
247 *
248 * \param[in] cm pointer to a cs_cell_mesh_t structure
249 * \param[in] pty pointer to a cs_property_t structure
250 * \param[in] t_eval physical time at which one evaluates the term
251 *
252 * \return the value of the property for the given cell
253 */
254 /*----------------------------------------------------------------------------*/
255
256 static inline cs_real_t
_value_in_cell(const cs_cell_mesh_t * cm,const cs_property_t * pty,cs_real_t t_eval)257 _value_in_cell(const cs_cell_mesh_t *cm,
258 const cs_property_t *pty,
259 cs_real_t t_eval)
260 {
261 cs_real_t result = 0;
262 int def_id = 0;
263 if (pty->n_definitions > 1) {
264 def_id = pty->def_ids[cm->c_id];
265 assert(def_id > -1);
266 }
267 cs_xdef_t *def = pty->defs[def_id];
268
269 assert(pty->get_eval_at_cell_cw[def_id] != NULL);
270 pty->get_eval_at_cell_cw[def_id](cm, t_eval, def->context, &result);
271
272 return result;
273 }
274
275 /*----------------------------------------------------------------------------*/
276 /*!
277 * \brief Compute the value of the tensor attached to a property at the cell
278 * center
279 *
280 * \param[in] c_id id of the current cell
281 * \param[in] t_eval physical time at which one evaluates the term
282 * \param[in] pty pointer to a cs_property_t structure
283 * \param[in, out] tensor 3x3 matrix
284 */
285 /*----------------------------------------------------------------------------*/
286
287 static void
_get_cell_tensor(cs_lnum_t c_id,cs_real_t t_eval,const cs_property_t * pty,cs_real_t tensor[3][3])288 _get_cell_tensor(cs_lnum_t c_id,
289 cs_real_t t_eval,
290 const cs_property_t *pty,
291 cs_real_t tensor[3][3])
292 {
293 int def_id = 0;
294 if (pty->n_definitions > 1) {
295 def_id = pty->def_ids[c_id];
296 assert(def_id > -1);
297 }
298 assert(pty->get_eval_at_cell[def_id] != NULL);
299
300 cs_xdef_t *def = pty->defs[def_id];
301
302 if (pty->type & CS_PROPERTY_ISO) {
303
304 double eval;
305 pty->get_eval_at_cell[def_id](1, &c_id, true, /* dense output */
306 cs_glob_mesh,
307 cs_cdo_connect,
308 cs_cdo_quant,
309 t_eval,
310 def->context,
311 &eval);
312
313 tensor[0][0] = tensor[1][1] = tensor[2][2] = eval;
314
315 }
316 else if (pty->type & CS_PROPERTY_ORTHO) {
317
318 double eval[3];
319 pty->get_eval_at_cell[def_id](1, &c_id, true, /* dense output */
320 cs_glob_mesh,
321 cs_cdo_connect,
322 cs_cdo_quant,
323 t_eval,
324 def->context,
325 eval);
326
327 for (int k = 0; k < 3; k++)
328 tensor[k][k] = eval[k];
329
330 }
331 else if (pty->type & CS_PROPERTY_ANISO_SYM) {
332
333 double eval[6];
334 pty->get_eval_at_cell[def_id](1, &c_id, true, /* dense output */
335 cs_glob_mesh,
336 cs_cdo_connect,
337 cs_cdo_quant,
338 t_eval,
339 def->context,
340 eval);
341
342 /* Diag. values */
343 tensor[0][0] = eval[0];
344 tensor[1][1] = eval[1];
345 tensor[2][2] = eval[2];
346
347 /* Extra-diag. values */
348 tensor[0][1] = tensor[1][0] = eval[3];
349 tensor[0][2] = tensor[2][0] = eval[4];
350 tensor[1][2] = tensor[2][1] = eval[5];
351
352 }
353 else {
354
355 assert(pty->type & CS_PROPERTY_ANISO);
356 pty->get_eval_at_cell[def_id](1, &c_id, true, /* dense output */
357 cs_glob_mesh,
358 cs_cdo_connect,
359 cs_cdo_quant,
360 t_eval,
361 def->context,
362 (cs_real_t *)tensor);
363
364 }
365 }
366
367 /*----------------------------------------------------------------------------*/
368 /*!
369 * \brief Compute the value of the tensor attached to a property at the cell
370 * center. Case of a property as the product of two other properties.
371 *
372 * \param[in] c_id id of the current cell
373 * \param[in] t_eval physical time at which one evaluates the term
374 * \param[in] pty pointer to a cs_property_t structure
375 * \param[in, out] tensor 3x3 matrix
376 */
377 /*----------------------------------------------------------------------------*/
378
379 static void
_get_cell_tensor_by_property_product(cs_lnum_t c_id,cs_real_t t_eval,const cs_property_t * pty,cs_real_t tensor[3][3])380 _get_cell_tensor_by_property_product(cs_lnum_t c_id,
381 cs_real_t t_eval,
382 const cs_property_t *pty,
383 cs_real_t tensor[3][3])
384 {
385 assert(pty->related_properties != NULL);
386 const cs_property_t *a = pty->related_properties[0];
387 const cs_property_t *b = pty->related_properties[1];
388
389 cs_real_t tensor_a[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
390 cs_real_t tensor_b[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
391
392 /* Evaluates each property */
393 _get_cell_tensor(c_id, t_eval, a, tensor_a);
394 _get_cell_tensor(c_id, t_eval, b, tensor_b);
395
396 /* Compute the product */
397 if (pty->type & CS_PROPERTY_ISO) {
398 /* a and b are isotropic */
399 tensor[0][0] = tensor[1][1] = tensor[2][2] = tensor_a[0][0]*tensor_b[0][0];
400 }
401 else if (pty->type & CS_PROPERTY_ORTHO) {
402 for (int k = 0; k < 3; k++)
403 tensor[k][k] = tensor_a[k][k]*tensor_b[k][k];
404 }
405 else {
406 assert(pty->type & CS_PROPERTY_ANISO);
407 /* tensor has been initialized by the calling function */
408 cs_math_33_product_add(tensor_a, tensor_b, tensor);
409 }
410 }
411
412 /*----------------------------------------------------------------------------*/
413 /*!
414 * \brief Compute the value of the tensor attached a property at the cell
415 * center
416 * Version using a cs_cell_mesh_t structure
417 *
418 * \param[in] cm pointer to a cs_cell_mesh_t structure
419 * \param[in] pty pointer to a cs_property_t structure
420 * \param[in] t_eval physical time at which one evaluates the term
421 * \param[in, out] tensor 3x3 matrix
422 */
423 /*----------------------------------------------------------------------------*/
424
425 static void
_tensor_in_cell(const cs_cell_mesh_t * cm,const cs_property_t * pty,cs_real_t t_eval,cs_real_t tensor[3][3])426 _tensor_in_cell(const cs_cell_mesh_t *cm,
427 const cs_property_t *pty,
428 cs_real_t t_eval,
429 cs_real_t tensor[3][3])
430 {
431 int def_id = 0;
432 if (pty->n_definitions > 1) {
433 def_id = pty->def_ids[cm->c_id];
434 assert(def_id > -1);
435 }
436 assert(pty->get_eval_at_cell_cw[def_id] != NULL);
437
438 cs_xdef_t *def = pty->defs[def_id];
439
440 if (pty->type & CS_PROPERTY_ISO) {
441
442 double eval;
443 pty->get_eval_at_cell_cw[def_id](cm, t_eval, def->context, &eval);
444 tensor[0][0] = tensor[1][1] = tensor[2][2] = eval;
445
446 }
447 else if (pty->type & CS_PROPERTY_ORTHO) {
448
449 double eval[3];
450 pty->get_eval_at_cell_cw[def_id](cm, t_eval, def->context, eval);
451 for (int k = 0; k < 3; k++)
452 tensor[k][k] = eval[k];
453
454 }
455 else if (pty->type & CS_PROPERTY_ANISO_SYM) {
456
457 double eval[6];
458 pty->get_eval_at_cell_cw[def_id](cm, t_eval, def->context, eval);
459
460 /* Diag. values */
461 tensor[0][0] = eval[0];
462 tensor[1][1] = eval[1];
463 tensor[2][2] = eval[2];
464
465 /* Extra-diag. values */
466 tensor[0][1] = tensor[1][0] = eval[3];
467 tensor[0][2] = tensor[2][0] = eval[4];
468 tensor[1][2] = tensor[2][1] = eval[5];
469
470 }
471 else {
472
473 assert(pty->type & CS_PROPERTY_ANISO);
474 pty->get_eval_at_cell_cw[def_id](cm, t_eval, def->context,
475 (cs_real_t *)tensor);
476
477 }
478 }
479
480 /*----------------------------------------------------------------------------*/
481 /*!
482 * \brief Compute the value of the tensor attached a property at the cell
483 * center.
484 * Version using a cs_cell_mesh_t structure and with a property
485 * defined as the product of two existing ones.
486 *
487 * \param[in] cm pointer to a cs_cell_mesh_t structure
488 * \param[in] pty pointer to a cs_property_t structure
489 * \param[in] t_eval physical time at which one evaluates the term
490 * \param[in, out] tensor 3x3 matrix
491 */
492 /*----------------------------------------------------------------------------*/
493
494 static void
_tensor_in_cell_by_property_product(const cs_cell_mesh_t * cm,const cs_property_t * pty,cs_real_t t_eval,cs_real_t tensor[3][3])495 _tensor_in_cell_by_property_product(const cs_cell_mesh_t *cm,
496 const cs_property_t *pty,
497 cs_real_t t_eval,
498 cs_real_t tensor[3][3])
499 {
500 assert(pty->related_properties != NULL);
501 const cs_property_t *a = pty->related_properties[0];
502 const cs_property_t *b = pty->related_properties[1];
503
504 cs_real_t tensor_a[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
505 cs_real_t tensor_b[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
506
507 /* Evaluates each property */
508 _tensor_in_cell(cm, a, t_eval, tensor_a);
509 _tensor_in_cell(cm, b, t_eval, tensor_b);
510
511 /* Compute the product */
512 if (pty->type & CS_PROPERTY_ISO) {
513 /* a and b are isotropic */
514 tensor[0][0] = tensor[1][1] = tensor[2][2] = tensor_a[0][0]*tensor_b[0][0];
515 }
516 else if (pty->type & CS_PROPERTY_ORTHO) {
517 for (int k = 0; k < 3; k++)
518 tensor[k][k] = tensor_a[k][k]*tensor_b[k][k];
519 }
520 else {
521 assert(pty->type & CS_PROPERTY_ANISO);
522 tensor[0][0] = tensor[1][1] = tensor[2][2] = 0;
523 cs_math_33_product_add(tensor_a, tensor_b, tensor);
524 }
525 }
526
527 /*----------------------------------------------------------------------------*/
528 /*!
529 * \brief Define a property as the product of two existing properties
530 *
531 * \param[in, out] pty resulting property
532 */
533 /*----------------------------------------------------------------------------*/
534
535 static void
_define_pty_by_product(cs_property_t * pty)536 _define_pty_by_product(cs_property_t *pty)
537 {
538 /* Only one definition is added in this case specifying that the definition
539 * relies on other definitions to be defined. The exact way to specify values
540 * is managed by the calling code (with a call to each sub-definition using
541 * the standard algorithm)
542 */
543
544 int id = _add_new_def(pty);
545 assert(id == 0);
546
547 int dim = 1;
548 if (pty->type & CS_PROPERTY_ORTHO)
549 dim = 3;
550 else if (pty->type & CS_PROPERTY_ANISO_SYM)
551 dim = 6;
552 else if (pty->type & CS_PROPERTY_ANISO)
553 dim = 9;
554
555 cs_flag_t state_flag = 0;
556 cs_flag_t meta_flag = 0;
557
558 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_SUB_DEFINITIONS,
559 dim,
560 0, /* zone_id = all cells */
561 state_flag,
562 meta_flag,
563 NULL); /* no input */
564
565 /* Set pointers */
566 pty->defs[id] = d;
567 pty->get_eval_at_cell[id] = NULL;
568 pty->get_eval_at_cell_cw[id] = NULL;
569 }
570
571 /*----------------------------------------------------------------------------*/
572 /*!
573 * \brief Create and initialize a new property structure
574 *
575 * \param[in] name name of the property
576 * \param[in] id id of the property to create
577 * \param[in] type type of property
578 *
579 * \return a pointer to a new allocated cs_property_t structure
580 */
581 /*----------------------------------------------------------------------------*/
582
583 static cs_property_t *
_create_property(const char * name,int id,cs_property_type_t type)584 _create_property(const char *name,
585 int id,
586 cs_property_type_t type)
587 {
588 /* Check the sanity of type */
589
590 int n_types = 0;
591 const int flags[] = {CS_PROPERTY_ISO,
592 CS_PROPERTY_ORTHO,
593 CS_PROPERTY_ANISO_SYM,
594 CS_PROPERTY_ANISO};
595
596 for (int i = 0; i < 4; i++) {
597 if (type & flags[i])
598 n_types += 1;
599 }
600
601 if (n_types > 1) {
602
603 const char *names[] = {"CS_PROPERTY_ISO",
604 "CS_PROPERTY_ORTHO",
605 "CS_PROPERTY_ANISO_SYM",
606 "CS_PROPERTY_ANISO"};
607 int l = 0;
608 char prop_list[256] = "";
609
610 for (int i = 0; i < 4 && l > 0; i++) {
611 if (type & flags[i]) {
612 snprintf(prop_list+l, 256-l, " %s\n", names[i]);
613 prop_list[255] = '\0';
614 l = strlen(prop_list);
615 }
616 }
617
618 }
619 else if (n_types < 1)
620 if ((type & CS_PROPERTY_ANISO) == 0)
621 bft_error(__FILE__, __LINE__, 0,
622 "%s: No known type specified for property %s\n"
623 " Set one among\n"
624 " CS_PROPERTY_ISO,\n"
625 " CS_PROPERTY_ORTHO,\n"
626 " CS_PROPERTY_ANISO_SYM,\n"
627 " CS_PROPERTY_ANISO.\n", __func__, name);
628
629 cs_property_t *pty = NULL;
630 BFT_MALLOC(pty, 1, cs_property_t);
631
632 /* Copy name */
633
634 size_t len = strlen(name);
635 BFT_MALLOC(pty->name, len + 1, char);
636 strncpy(pty->name, name, len);
637 pty->name[len] = '\0';
638
639 pty->id = id;
640 pty->type = type;
641 pty->state_flag = 0;
642 pty->process_flag = 0;
643
644 pty->ref_value = 1.0; /* default setting */
645
646 pty->n_definitions = 0;
647 pty->defs = NULL;
648 pty->def_ids = NULL;
649
650 pty->get_eval_at_cell = NULL;
651 pty->get_eval_at_cell_cw = NULL;
652
653 pty->n_related_properties = 0;
654 pty->related_properties = NULL;
655
656 return pty;
657 }
658
659 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
660
661 /*============================================================================
662 * Public function prototypes
663 *============================================================================*/
664
665 /*----------------------------------------------------------------------------*/
666 /*!
667 * \brief Set shared pointers to main domain members
668 *
669 * \param[in] quant additional mesh quantities struct.
670 * \param[in] connect pointer to a cs_cdo_connect_t struct.
671 */
672 /*----------------------------------------------------------------------------*/
673
674 void
cs_property_set_shared_pointers(const cs_cdo_quantities_t * quant,const cs_cdo_connect_t * connect)675 cs_property_set_shared_pointers(const cs_cdo_quantities_t *quant,
676 const cs_cdo_connect_t *connect)
677 {
678 /* Assign static const pointers */
679 cs_cdo_quant = quant;
680 cs_cdo_connect = connect;
681 }
682
683 /*----------------------------------------------------------------------------*/
684 /*!
685 * \brief Retrieve the number of properties
686 *
687 * \return the number of properties
688 */
689 /*----------------------------------------------------------------------------*/
690
691 int
cs_property_get_n_properties(void)692 cs_property_get_n_properties(void)
693 {
694 return _n_properties;
695 }
696
697 /*----------------------------------------------------------------------------*/
698 /*!
699 * \brief Create and initialize a new property structure
700 *
701 * \param[in] name name of the property
702 * \param[in] type type of property
703 *
704 * \return a pointer to a new allocated cs_property_t structure
705 */
706 /*----------------------------------------------------------------------------*/
707
708 cs_property_t *
cs_property_add(const char * name,cs_property_type_t type)709 cs_property_add(const char *name,
710 cs_property_type_t type)
711 {
712 cs_property_t *pty = cs_property_by_name(name);
713
714 if (pty != NULL) {
715 cs_base_warn(__FILE__, __LINE__);
716 cs_log_printf(CS_LOG_DEFAULT,
717 _(" %s: An existing property has already the name %s.\n"
718 " Stop adding this property.\n"), __func__, name);
719 return pty;
720 }
721
722 int pty_id = _n_properties;
723
724 if (pty_id == 0) {
725 _n_max_properties = 3;
726 BFT_MALLOC(_properties, _n_max_properties, cs_property_t *);
727 }
728
729 _n_properties += 1;
730
731 if (_n_properties > _n_max_properties) {
732 _n_max_properties *= 2;
733 BFT_REALLOC(_properties, _n_max_properties, cs_property_t *);
734 }
735
736 _properties[pty_id] = _create_property(name, pty_id, type);
737
738 return _properties[pty_id];
739 }
740
741 /*----------------------------------------------------------------------------*/
742 /*!
743 * \brief Define a cs_property_t structure thanks to the product of two
744 * properties
745 * The type is infered from that of the related properties
746 * The value of the property is given as:
747 * value_ab = value_a * value_b
748 *
749 * \param[in] name name of the property
750 * \param[in] pty_a pointer to a cs_property_t structure
751 * \param[in] pty_b pointer to a cs_property_t structure
752 *
753 * \return a pointer to a new allocated cs_property_t structure
754 */
755 /*----------------------------------------------------------------------------*/
756
757 cs_property_t *
cs_property_add_as_product(const char * name,const cs_property_t * pty_a,const cs_property_t * pty_b)758 cs_property_add_as_product(const char *name,
759 const cs_property_t *pty_a,
760 const cs_property_t *pty_b)
761 {
762 if (pty_a == NULL || pty_b == NULL)
763 return NULL;
764
765 /* Determine the type of the new property */
766 cs_property_type_t type = CS_PROPERTY_BY_PRODUCT;
767
768 /* pty_b |pty_a -> iso | ortho | aniso
769 * iso | iso | ortho | aniso
770 * ortho | ortho | ortho | aniso
771 * aniso | aniso | aniso | aniso */
772
773 if (pty_a->type & CS_PROPERTY_ISO) {
774 if (pty_b->type & CS_PROPERTY_ISO)
775 type |= CS_PROPERTY_ISO;
776 else if (pty_b->type & CS_PROPERTY_ORTHO)
777 type |= CS_PROPERTY_ORTHO;
778 else if (pty_b->type & CS_PROPERTY_ANISO)
779 type |= CS_PROPERTY_ANISO;
780 else
781 bft_error(__FILE__, __LINE__, 0,
782 " %s: Invalid type of property.", __func__);
783 }
784 else if (pty_a->type & CS_PROPERTY_ANISO) {
785 type |= CS_PROPERTY_ANISO;
786 }
787 else if (pty_a->type & CS_PROPERTY_ORTHO) {
788 if (pty_b->type & CS_PROPERTY_ANISO)
789 type |= CS_PROPERTY_ANISO;
790 else
791 type |= CS_PROPERTY_ORTHO;
792 }
793 else
794 bft_error(__FILE__, __LINE__, 0,
795 " %s: Invalid type of property.", __func__);
796
797 cs_property_t *pty_ab = cs_property_add(name, type);
798
799 pty_ab->n_related_properties = 2;
800 BFT_MALLOC(pty_ab->related_properties, 2, const cs_property_t *);
801
802 pty_ab->related_properties[0] = pty_a;
803 pty_ab->related_properties[1] = pty_b;
804
805 return pty_ab;
806 }
807
808 /*----------------------------------------------------------------------------*/
809 /*!
810 * \brief Find the related property definition from its name
811 *
812 * \param[in] name name of the property to find
813 *
814 * \return NULL if not found otherwise the associated pointer
815 */
816 /*----------------------------------------------------------------------------*/
817
818 cs_property_t *
cs_property_by_name(const char * name)819 cs_property_by_name(const char *name)
820 {
821 if (_n_properties < 0)
822 return NULL;
823 assert(name != NULL);
824
825 for (int i = 0; i < _n_properties; i++) {
826 cs_property_t *pty = _properties[i];
827 assert(pty->name != NULL);
828 if (strcmp(pty->name, name) == 0)
829 return pty;
830 }
831
832 return NULL;
833 }
834
835 /*----------------------------------------------------------------------------*/
836 /*!
837 * \brief Find the related property definition from its id
838 *
839 * \param[in] id id of the property to find
840 *
841 * \return NULL if not found otherwise the associated pointer
842 */
843 /*----------------------------------------------------------------------------*/
844
845 cs_property_t *
cs_property_by_id(int id)846 cs_property_by_id(int id)
847 {
848 if (_n_properties < 0)
849 return NULL;
850 if (id < 0 || id >= _n_max_properties)
851 return NULL;
852
853 return _properties[id];
854 }
855
856 /*----------------------------------------------------------------------------*/
857 /*!
858 * \brief Set optional parameters related to a cs_property_t structure
859 *
860 * \param[in, out] pty pointer to a cs_property_t structure
861 * \param[in] key key related to a setting option
862 */
863 /*----------------------------------------------------------------------------*/
864
865 void
cs_property_set_option(cs_property_t * pty,cs_property_key_t key)866 cs_property_set_option(cs_property_t *pty,
867 cs_property_key_t key)
868 {
869 if (pty == NULL)
870 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
871
872 switch(key) {
873
874 case CS_PTYKEY_POST_FOURIER:
875 pty->process_flag |= CS_PROPERTY_POST_FOURIER;
876 break;
877
878 default:
879 bft_error(__FILE__, __LINE__, 0,
880 _(" Key not implemented for setting a property."));
881 break;
882
883 } /* Switch on keys */
884
885 }
886
887
888 /*----------------------------------------------------------------------------*/
889 /*!
890 * \brief Set the reference value associated to a \ref cs_property_t structure
891 * This is a real number even whatever the type of property is.
892 *
893 * \param[in, out] pty pointer to a cs_property_t structure
894 * \param[in] refval value to set
895 */
896 /*----------------------------------------------------------------------------*/
897
898 void
cs_property_set_reference_value(cs_property_t * pty,double refval)899 cs_property_set_reference_value(cs_property_t *pty,
900 double refval)
901 {
902 if (pty == NULL)
903 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
904
905 pty->ref_value = refval;
906 }
907
908 /*----------------------------------------------------------------------------*/
909 /*!
910 * \brief Free all cs_property_t structures and the array storing all the
911 * structures
912 */
913 /*----------------------------------------------------------------------------*/
914
915 void
cs_property_destroy_all(void)916 cs_property_destroy_all(void)
917 {
918 if (_n_properties == 0)
919 return;
920
921 for (int i = 0; i < _n_properties; i++) {
922
923 cs_property_t *pty = _properties[i];
924
925 if (pty == NULL)
926 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
927
928 BFT_FREE(pty->name);
929 BFT_FREE(pty->def_ids);
930
931 for (int j = 0; j < pty->n_definitions; j++)
932 pty->defs[j] = cs_xdef_free(pty->defs[j]);
933
934 BFT_FREE(pty->defs);
935 BFT_FREE(pty->get_eval_at_cell);
936 BFT_FREE(pty->get_eval_at_cell_cw);
937
938 if (pty->n_related_properties > 0)
939 BFT_FREE(pty->related_properties);
940
941 BFT_FREE(pty);
942
943 } /* Loop on properties */
944
945 BFT_FREE(_properties);
946 _n_properties = 0;
947 _n_max_properties = 0;
948 }
949
950 /*----------------------------------------------------------------------------*/
951 /*!
952 * \brief Last stage of the definition of a property based on several
953 * definitions (i.e. definition by subdomains)
954 */
955 /*----------------------------------------------------------------------------*/
956
957 void
cs_property_finalize_setup(void)958 cs_property_finalize_setup(void)
959 {
960 if (_n_properties == 0)
961 return;
962
963 for (int i = 0; i < _n_properties; i++) {
964
965 cs_property_t *pty = _properties[i];
966
967 if (pty == NULL)
968 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
969
970 if (pty->type & CS_PROPERTY_BY_PRODUCT)
971 continue;
972
973 if (pty->n_definitions > 1) { /* Initialization of def_ids */
974
975 const cs_lnum_t n_cells = cs_cdo_quant->n_cells;
976
977 BFT_MALLOC(pty->def_ids, n_cells, short int);
978
979 # pragma omp parallel for if (n_cells > CS_THR_MIN)
980 for (cs_lnum_t j = 0; j < n_cells; j++)
981 pty->def_ids[j] = -1; /* Unset by default */
982
983 for (int id = 0; id < pty->n_definitions; id++) {
984
985 const cs_xdef_t *def = pty->defs[id];
986
987 assert(def->z_id > 0);
988 assert(def->support == CS_XDEF_SUPPORT_VOLUME);
989
990 const cs_zone_t *z = cs_volume_zone_by_id(def->z_id);
991
992 assert(z != NULL);
993
994 # pragma omp parallel for if (z->n_elts > CS_THR_MIN)
995 for (cs_lnum_t j = 0; j < z->n_elts; j++)
996 pty->def_ids[z->elt_ids[j]] = id;
997
998 } /* Loop on definitions */
999
1000 /* Check if the property is defined everywhere */
1001 for (cs_lnum_t j = 0; j < n_cells; j++)
1002 if (pty->def_ids[j] == -1)
1003 bft_error(__FILE__, __LINE__, 0,
1004 " %s: cell %ld is unset for property %s\n",
1005 __func__, (long)j, pty->name);
1006
1007 }
1008 else if (pty->n_definitions == 0) {
1009
1010 /* Default initialization */
1011 if (pty->type & CS_PROPERTY_ISO)
1012 cs_property_def_iso_by_value(pty, NULL, pty->ref_value);
1013 else if (pty->type & CS_PROPERTY_ORTHO) {
1014 cs_real_t ref[3] = {pty->ref_value, pty->ref_value, pty->ref_value};
1015 cs_property_def_ortho_by_value(pty, NULL, ref);
1016 }
1017 else if (pty->type & CS_PROPERTY_ANISO) {
1018 cs_real_t ref[3][3] = { {pty->ref_value, 0, 0},
1019 {0, pty->ref_value, 0},
1020 {0, 0, pty->ref_value} };
1021 cs_property_def_aniso_by_value(pty, NULL, ref);
1022 }
1023 else
1024 bft_error(__FILE__, __LINE__, 0, "%s: Incompatible property type.",
1025 __func__);
1026
1027 cs_log_printf(CS_LOG_DEFAULT,
1028 "\n Property \"%s\" will be defined using its reference"
1029 " value.\n", pty->name);
1030
1031 }
1032
1033 } /* Loop on properties */
1034
1035 for (int i = 0; i < _n_properties; i++) {
1036
1037 cs_property_t *pty = _properties[i];
1038
1039 if (pty->type & CS_PROPERTY_BY_PRODUCT) {
1040
1041 assert(pty->n_related_properties == 2);
1042
1043 const cs_property_t *pty_a = pty->related_properties[0];
1044 const cs_property_t *pty_b = pty->related_properties[1];
1045
1046 pty->ref_value = pty_a->ref_value * pty_b->ref_value;
1047
1048 _define_pty_by_product(pty);
1049
1050 } /* Only properties defined as a product */
1051
1052 } /* Loop on properties */
1053
1054 }
1055
1056 /*----------------------------------------------------------------------------*/
1057 /*!
1058 * \brief Initialize a \ref cs_property_data_t structure. If property is NULL
1059 * then one considers that this is a unitary property
1060 *
1061 * \param[in] need_tensor true if one needs a tensor-valued evaluation
1062 * \param[in] need_eigen true if one needs an evaluation of eigen values
1063 * \param[in] property pointer to the \ref cs_property_t structure
1064 * \param[in, out] data structure to initialize (already allocated)
1065 */
1066 /*----------------------------------------------------------------------------*/
1067
1068 void
cs_property_data_init(bool need_tensor,bool need_eigen,const cs_property_t * property,cs_property_data_t * data)1069 cs_property_data_init(bool need_tensor,
1070 bool need_eigen,
1071 const cs_property_t *property,
1072 cs_property_data_t *data)
1073 {
1074 if (data == NULL)
1075 return;
1076
1077 data->property = property;
1078
1079 data->is_unity = false;
1080 data->is_iso = false;
1081
1082 if (property == NULL) {
1083 data->is_iso = true;
1084 data->is_unity = true;
1085 }
1086 else {
1087
1088 if (property->type & CS_PROPERTY_ISO) {
1089 data->is_iso = true;
1090
1091 if (property->n_definitions == 1) {
1092 cs_xdef_t *d = property->defs[0];
1093 if (d->type == CS_XDEF_BY_VALUE) {
1094 double *dval = (double *)d->context;
1095 if (fabs(dval[0] - 1) < FLT_MIN)
1096 data->is_iso = true;
1097 }
1098 }
1099 }
1100
1101 }
1102
1103 const cs_real_t ref_val = (property == NULL) ? 1. : property->ref_value;
1104
1105 data->need_eigen = need_eigen;
1106 data->eigen_max = ref_val;
1107 data->eigen_ratio = 1.0;
1108
1109 data->need_tensor = need_tensor;
1110
1111 data->value = ref_val;
1112 data->tensor[0][0] = ref_val, data->tensor[0][1] = 0, data->tensor[0][2] = 0;
1113 data->tensor[1][0] = 0, data->tensor[1][1] = ref_val, data->tensor[1][2] = 0;
1114 data->tensor[2][0] = 0, data->tensor[2][1] = 0, data->tensor[2][2] = ref_val;
1115 }
1116
1117 /*----------------------------------------------------------------------------*/
1118 /*!
1119 * \brief Define a single uniform and steady isotropic definition for the
1120 * given cs_property_t structure.
1121 * This is a specialized variant of \ref cs_property_def_iso_by_value
1122 * since several assumptions are satisfied.
1123 *
1124 * \param[in, out] pty pointer to a cs_property_t structure
1125 * \param[in] val value to set
1126 *
1127 * \return a pointer to the resulting cs_xdef_t structure
1128 */
1129 /*----------------------------------------------------------------------------*/
1130
1131 cs_xdef_t *
cs_property_def_constant_value(cs_property_t * pty,double val)1132 cs_property_def_constant_value(cs_property_t *pty,
1133 double val)
1134 {
1135 if (pty == NULL)
1136 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
1137 if ((pty->type & CS_PROPERTY_ISO) == 0)
1138 bft_error(__FILE__, __LINE__, 0,
1139 " Invalid setting: property %s is not isotropic.\n"
1140 " Please check your settings.", pty->name);
1141
1142 int new_id = _add_new_def(pty);
1143
1144 if (new_id > 0)
1145 bft_error(__FILE__, __LINE__, 0,
1146 " %s: Invalid setting: property %s is assumed to be constant.\n"
1147 " Several definitions have been added.\n"
1148 " Please check your settings.", __func__, pty->name);
1149
1150 cs_flag_t state_flag = CS_FLAG_STATE_UNIFORM | CS_FLAG_STATE_CELLWISE |
1151 CS_FLAG_STATE_STEADY;
1152 cs_flag_t meta_flag = 0; /* metadata */
1153 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_VALUE,
1154 1, /* dim */
1155 0, /* all cells */
1156 state_flag,
1157 meta_flag,
1158 &val); /* context */
1159
1160 pty->defs[new_id] = d;
1161 pty->get_eval_at_cell[new_id] = cs_xdef_eval_scalar_by_val;
1162 pty->get_eval_at_cell_cw[new_id] = cs_xdef_cw_eval_scalar_by_val;
1163
1164 /* Set the state flag */
1165 pty->state_flag |= CS_FLAG_STATE_CELLWISE | CS_FLAG_STATE_STEADY;
1166 pty->state_flag |= CS_FLAG_STATE_UNIFORM;
1167
1168 /* Set automatically the reference value if all cells are selected */
1169 cs_property_set_reference_value(pty, val);
1170
1171 return d;
1172 }
1173
1174 /*----------------------------------------------------------------------------*/
1175 /*!
1176 * \brief Define an isotropic cs_property_t structure by value for entities
1177 * related to a volume zone
1178 *
1179 * \param[in, out] pty pointer to a cs_property_t structure
1180 * \param[in] zname name of the associated zone (if NULL or "" all
1181 * cells are considered)
1182 * \param[in] val value to set
1183 *
1184 * \return a pointer to the resulting cs_xdef_t structure
1185 */
1186 /*----------------------------------------------------------------------------*/
1187
1188 cs_xdef_t *
cs_property_def_iso_by_value(cs_property_t * pty,const char * zname,double val)1189 cs_property_def_iso_by_value(cs_property_t *pty,
1190 const char *zname,
1191 double val)
1192 {
1193 if (pty == NULL)
1194 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
1195 if ((pty->type & CS_PROPERTY_ISO) == 0)
1196 bft_error(__FILE__, __LINE__, 0,
1197 " Invalid setting: property %s is not isotropic.\n"
1198 " Please check your settings.", pty->name);
1199
1200 int new_id = _add_new_def(pty);
1201 int z_id = cs_get_vol_zone_id(zname);
1202 cs_flag_t state_flag = CS_FLAG_STATE_UNIFORM | CS_FLAG_STATE_CELLWISE |
1203 CS_FLAG_STATE_STEADY;
1204 cs_flag_t meta_flag = 0; /* metadata */
1205 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_VALUE,
1206 1, /* dim */
1207 z_id,
1208 state_flag,
1209 meta_flag,
1210 &val); /* context */
1211
1212 pty->defs[new_id] = d;
1213 pty->get_eval_at_cell[new_id] = cs_xdef_eval_scalar_by_val;
1214 pty->get_eval_at_cell_cw[new_id] = cs_xdef_cw_eval_scalar_by_val;
1215
1216 /* Set the state flag */
1217 pty->state_flag |= CS_FLAG_STATE_CELLWISE | CS_FLAG_STATE_STEADY;
1218 if (z_id == 0)
1219 pty->state_flag |= CS_FLAG_STATE_UNIFORM;
1220
1221 /* Set automatically the reference value if all cells are selected */
1222 if (z_id == 0)
1223 cs_property_set_reference_value(pty, val);
1224
1225 return d;
1226 }
1227
1228 /*----------------------------------------------------------------------------*/
1229 /*!
1230 * \brief Define an orthotropic cs_property_t structure by value for entities
1231 * related to a volume zone
1232 *
1233 * \param[in, out] pty pointer to a cs_property_t structure
1234 * \param[in] zname name of the associated zone (if NULL or "" all
1235 * cells are considered)
1236 * \param[in] val values to set (vector of size 3)
1237 *
1238 * \return a pointer to the resulting cs_xdef_t structure
1239 */
1240 /*----------------------------------------------------------------------------*/
1241
1242 cs_xdef_t *
cs_property_def_ortho_by_value(cs_property_t * pty,const char * zname,double val[])1243 cs_property_def_ortho_by_value(cs_property_t *pty,
1244 const char *zname,
1245 double val[])
1246 {
1247 if (pty == NULL)
1248 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
1249 if ((pty->type & CS_PROPERTY_ORTHO) == 0)
1250 bft_error(__FILE__, __LINE__, 0,
1251 " Invalid setting: property %s is not orthotropic.\n"
1252 " Please check your settings.", pty->name);
1253
1254 int new_id = _add_new_def(pty);
1255 int z_id = cs_get_vol_zone_id(zname);
1256 cs_flag_t state_flag = CS_FLAG_STATE_UNIFORM | CS_FLAG_STATE_CELLWISE |
1257 CS_FLAG_STATE_STEADY;
1258 cs_flag_t meta_flag = 0; /* metadata */
1259 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_VALUE,
1260 3, /* dim */
1261 z_id,
1262 state_flag,
1263 meta_flag,
1264 (void *)val);
1265
1266 pty->defs[new_id] = d;
1267 pty->get_eval_at_cell[new_id] = cs_xdef_eval_vector_by_val;
1268 pty->get_eval_at_cell_cw[new_id] = cs_xdef_cw_eval_vector_by_val;
1269
1270 /* Set the state flag */
1271 pty->state_flag |= CS_FLAG_STATE_CELLWISE | CS_FLAG_STATE_STEADY;
1272 if (z_id == 0)
1273 pty->state_flag |= CS_FLAG_STATE_UNIFORM;
1274
1275 return d;
1276 }
1277
1278 /*----------------------------------------------------------------------------*/
1279 /*!
1280 * \brief Define an anisotropic cs_property_t structure by value for entities
1281 * related to a volume zone
1282 *
1283 * \param[in, out] pty pointer to a cs_property_t structure
1284 * \param[in] zname name of the associated zone (if NULL or "" all
1285 * cells are considered)
1286 * \param[in] tens values to set (3x3 tensor)
1287 *
1288 * \return a pointer to the resulting cs_xdef_t structure
1289 */
1290 /*----------------------------------------------------------------------------*/
1291
1292 cs_xdef_t *
cs_property_def_aniso_by_value(cs_property_t * pty,const char * zname,cs_real_t tens[3][3])1293 cs_property_def_aniso_by_value(cs_property_t *pty,
1294 const char *zname,
1295 cs_real_t tens[3][3])
1296 {
1297 if (pty == NULL)
1298 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
1299 if ((pty->type & CS_PROPERTY_ANISO) == 0)
1300 bft_error(__FILE__, __LINE__, 0,
1301 " Invalid setting: property %s is not anisotropic.\n"
1302 " Please check your settings.", pty->name);
1303
1304 /* Check the symmetry */
1305 if (!_is_tensor_symmetric((const cs_real_t (*)[3])tens))
1306 bft_error(__FILE__, __LINE__, 0,
1307 _(" The definition of the tensor related to the"
1308 " property %s is not symmetric.\n"
1309 " This case is not handled. Please check your settings.\n"),
1310 pty->name);
1311
1312 int new_id = _add_new_def(pty);
1313 int z_id = cs_get_vol_zone_id(zname);
1314 cs_flag_t state_flag = CS_FLAG_STATE_UNIFORM | CS_FLAG_STATE_CELLWISE |
1315 CS_FLAG_STATE_STEADY;
1316 cs_flag_t meta_flag = 0; /* metadata */
1317 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_VALUE,
1318 9, /* dim */
1319 z_id,
1320 state_flag,
1321 meta_flag,
1322 tens);
1323
1324 pty->defs[new_id] = d;
1325 pty->get_eval_at_cell[new_id] = cs_xdef_eval_tensor_by_val;
1326 pty->get_eval_at_cell_cw[new_id] = cs_xdef_cw_eval_tensor_by_val;
1327
1328 /* Set the state flag */
1329 pty->state_flag |= CS_FLAG_STATE_CELLWISE | CS_FLAG_STATE_STEADY;
1330 if (z_id == 0)
1331 pty->state_flag |= CS_FLAG_STATE_UNIFORM;
1332
1333 return d;
1334 }
1335
1336 /*----------------------------------------------------------------------------*/
1337 /*!
1338 * \brief Define an anisotropic cs_property_t structure by value for entities
1339 * related to a volume zone
1340 *
1341 * \param[in, out] pty pointer to a cs_property_t structure
1342 * \param[in] zname name of the associated zone (if NULL or "" all
1343 * cells are considered)
1344 * \param[in] symtens values to set (6 values -- symmetric storage)
1345 *
1346 * \return a pointer to the resulting cs_xdef_t structure
1347 */
1348 /*----------------------------------------------------------------------------*/
1349
1350 cs_xdef_t *
cs_property_def_aniso_sym_by_value(cs_property_t * pty,const char * zname,cs_real_t symtens[6])1351 cs_property_def_aniso_sym_by_value(cs_property_t *pty,
1352 const char *zname,
1353 cs_real_t symtens[6])
1354 {
1355 if (pty == NULL)
1356 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
1357 if ((pty->type & CS_PROPERTY_ANISO_SYM) == 0)
1358 bft_error(__FILE__, __LINE__, 0,
1359 " Invalid setting: property %s is not anisotropic"
1360 " with a symmetric storage.\n"
1361 " Please check your settings.", pty->name);
1362
1363 int new_id = _add_new_def(pty);
1364 int z_id = cs_get_vol_zone_id(zname);
1365 cs_flag_t state_flag = CS_FLAG_STATE_UNIFORM | CS_FLAG_STATE_CELLWISE |
1366 CS_FLAG_STATE_STEADY;
1367 cs_flag_t meta_flag = 0; /* metadata */
1368 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_VALUE,
1369 6, /* dim */
1370 z_id,
1371 state_flag,
1372 meta_flag,
1373 symtens);
1374
1375 pty->defs[new_id] = d;
1376 pty->get_eval_at_cell[new_id] = cs_xdef_eval_symtens_by_val;
1377 pty->get_eval_at_cell_cw[new_id] = cs_xdef_cw_eval_symtens_by_val;
1378
1379 /* Set the state flag */
1380 pty->state_flag |= CS_FLAG_STATE_CELLWISE | CS_FLAG_STATE_STEADY;
1381 if (z_id == 0)
1382 pty->state_flag |= CS_FLAG_STATE_UNIFORM;
1383
1384 return d;
1385 }
1386
1387 /*----------------------------------------------------------------------------*/
1388 /*!
1389 * \brief Define a cs_property_t structure thanks to an analytic function in
1390 * a subdomain attached to the mesh location named ml_name
1391 *
1392 * \param[in, out] pty pointer to a cs_property_t structure
1393 * \param[in] zname name of the associated zone (if NULL or "" all
1394 * cells are considered)
1395 * \param[in] func pointer to a cs_analytic_func_t function
1396 * \param[in] input NULL or pointer to a structure cast on-the-fly
1397 *
1398 * \return a pointer to the resulting cs_xdef_t structure
1399 */
1400 /*----------------------------------------------------------------------------*/
1401
1402 cs_xdef_t *
cs_property_def_by_time_func(cs_property_t * pty,const char * zname,cs_time_func_t * func,void * input)1403 cs_property_def_by_time_func(cs_property_t *pty,
1404 const char *zname,
1405 cs_time_func_t *func,
1406 void *input)
1407 {
1408 if (pty == NULL)
1409 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
1410
1411 int new_id = _add_new_def(pty);
1412 int z_id = cs_get_vol_zone_id(zname);
1413 cs_flag_t state_flag = CS_FLAG_STATE_UNIFORM | CS_FLAG_STATE_CELLWISE;
1414 cs_flag_t meta_flag = 0; /* metadata */
1415 cs_xdef_time_func_context_t tfc = { .func = func,
1416 .input = input,
1417 .free_input = NULL };
1418
1419 /* Default initialization */
1420 pty->get_eval_at_cell[new_id] = NULL;
1421 pty->get_eval_at_cell_cw[new_id] = cs_xdef_cw_eval_by_time_func;
1422
1423 int dim = 0;
1424 if (pty->type & CS_PROPERTY_ISO) {
1425 dim = 1;
1426 pty->get_eval_at_cell[new_id] = cs_xdef_eval_scalar_at_cells_by_time_func;
1427 }
1428 else if (pty->type & CS_PROPERTY_ORTHO) {
1429 dim = 3;
1430 pty->get_eval_at_cell[new_id] = cs_xdef_eval_vector_at_cells_by_time_func;
1431 }
1432 else if (pty->type & CS_PROPERTY_ANISO_SYM) {
1433 dim = 6;
1434 pty->get_eval_at_cell[new_id] = cs_xdef_eval_symtens_at_cells_by_time_func;
1435 }
1436 else if (pty->type & CS_PROPERTY_ANISO) {
1437 dim = 9;
1438 pty->get_eval_at_cell[new_id] = cs_xdef_eval_tensor_at_cells_by_time_func;
1439 }
1440 else
1441 bft_error(__FILE__, __LINE__, 0, "%s: Incompatible property type.",
1442 __func__);
1443
1444 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_TIME_FUNCTION,
1445 dim,
1446 z_id,
1447 state_flag,
1448 meta_flag,
1449 &tfc);
1450
1451 pty->defs[new_id] = d;
1452
1453 /* Set the state flag for the property */
1454
1455 pty->state_flag |= CS_FLAG_STATE_CELLWISE;
1456 if (z_id == 0)
1457 pty->state_flag |= CS_FLAG_STATE_UNIFORM;
1458
1459 return d;
1460 }
1461
1462 /*----------------------------------------------------------------------------*/
1463 /*!
1464 * \brief Define a cs_property_t structure thanks to an analytic function in
1465 * a subdomain attached to the mesh location named ml_name
1466 *
1467 * \param[in, out] pty pointer to a cs_property_t structure
1468 * \param[in] zname name of the associated zone (if NULL or "" all
1469 * cells are considered)
1470 * \param[in] func pointer to a cs_analytic_func_t function
1471 * \param[in] input NULL or pointer to a structure cast on-the-fly
1472 *
1473 * \return a pointer to the resulting cs_xdef_t structure
1474 */
1475 /*----------------------------------------------------------------------------*/
1476
1477 cs_xdef_t *
cs_property_def_by_analytic(cs_property_t * pty,const char * zname,cs_analytic_func_t * func,void * input)1478 cs_property_def_by_analytic(cs_property_t *pty,
1479 const char *zname,
1480 cs_analytic_func_t *func,
1481 void *input)
1482 {
1483 if (pty == NULL)
1484 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
1485
1486 cs_flag_t state_flag = 0, meta_flag = 0; /* metadata */
1487 int z_id = cs_get_vol_zone_id(zname);
1488 cs_xdef_analytic_context_t ac = { .z_id = z_id,
1489 .func = func,
1490 .input = input,
1491 .free_input = NULL };
1492
1493 int dim = cs_property_get_dim(pty);
1494
1495 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_ANALYTIC_FUNCTION,
1496 dim,
1497 z_id,
1498 state_flag,
1499 meta_flag,
1500 &ac);
1501
1502 int new_id = _add_new_def(pty);
1503 pty->defs[new_id] = d;
1504 pty->get_eval_at_cell[new_id] = cs_xdef_eval_at_cells_by_analytic;
1505 pty->get_eval_at_cell_cw[new_id] = cs_xdef_cw_eval_by_analytic;
1506
1507 return d;
1508 }
1509
1510 /*----------------------------------------------------------------------------*/
1511 /*!
1512 * \brief Define a cs_property_t structure thanks to law depending on one
1513 * scalar variable in a subdomain attached to the mesh location named
1514 * ml_name
1515 *
1516 * \param[in, out] pty pointer to a cs_property_t structure
1517 * \param[in] zname name of the associated zone (if NULL or "" all
1518 * cells are considered)
1519 * \param[in] context pointer to a structure (may be NULL)
1520 * \param[in] get_eval_at_cell pointer to a function
1521 * \param[in] get_eval_at_cell_cw pointer to a function
1522 *
1523 * \return a pointer to the resulting cs_xdef_t structure
1524 */
1525 /*----------------------------------------------------------------------------*/
1526
1527 cs_xdef_t *
cs_property_def_by_func(cs_property_t * pty,const char * zname,void * context,cs_xdef_eval_t * get_eval_at_cell,cs_xdef_cw_eval_t * get_eval_at_cell_cw)1528 cs_property_def_by_func(cs_property_t *pty,
1529 const char *zname,
1530 void *context,
1531 cs_xdef_eval_t *get_eval_at_cell,
1532 cs_xdef_cw_eval_t *get_eval_at_cell_cw)
1533 {
1534 if (pty == NULL)
1535 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
1536
1537 int def_id = _add_new_def(pty);
1538 int z_id = cs_get_vol_zone_id(zname);
1539 cs_flag_t state_flag = 0;
1540 cs_flag_t meta_flag = 0; /* metadata */
1541
1542 int dim = cs_property_get_dim(pty);
1543
1544 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_FUNCTION,
1545 dim,
1546 z_id,
1547 state_flag,
1548 meta_flag,
1549 context);
1550
1551 pty->defs[def_id] = d;
1552 pty->get_eval_at_cell[def_id] = get_eval_at_cell;
1553 pty->get_eval_at_cell_cw[def_id] = get_eval_at_cell_cw;
1554
1555 return d;
1556 }
1557
1558 /*----------------------------------------------------------------------------*/
1559 /*!
1560 * \brief Define a cs_property_t structure thanks to an array of values
1561 *
1562 * \param[in, out] pty pointer to a cs_property_t structure
1563 * \param[in] loc information to know where are located values
1564 * \param[in] array pointer to an array
1565 * \param[in] is_owner transfer the lifecycle to the cs_xdef_t structure
1566 * (true or false)
1567 * \param[in] index optional pointer to the array index
1568 *
1569 * \return a pointer to the resulting cs_xdef_t structure
1570 */
1571 /*----------------------------------------------------------------------------*/
1572
1573 cs_xdef_t *
cs_property_def_by_array(cs_property_t * pty,cs_flag_t loc,cs_real_t * array,bool is_owner,cs_lnum_t * index)1574 cs_property_def_by_array(cs_property_t *pty,
1575 cs_flag_t loc,
1576 cs_real_t *array,
1577 bool is_owner,
1578 cs_lnum_t *index)
1579 {
1580 if (pty == NULL)
1581 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
1582
1583 int id = _add_new_def(pty);
1584
1585 if (pty->n_definitions > 1)
1586 bft_error(__FILE__, __LINE__, 0,
1587 " When a definition by array is requested, the max. number"
1588 " of subdomains to consider should be equal to 1.\n"
1589 " Current value is %d for property %s.\n"
1590 " Please modify your settings.",
1591 pty->n_definitions, pty->name);
1592
1593 int dim = cs_property_get_dim(pty);
1594 cs_flag_t state_flag = 0; /* Will be updated during the creation */
1595 cs_flag_t meta_flag = 0; /* metadata */
1596
1597 /* z_id = 0 since all the support is selected in this case */
1598 cs_xdef_array_context_t input = { .z_id = 0,
1599 .stride = dim,
1600 .loc = loc,
1601 .values = array,
1602 .is_owner = is_owner,
1603 .index = index };
1604
1605 cs_xdef_t *d = cs_xdef_volume_create(CS_XDEF_BY_ARRAY,
1606 dim,
1607 0, /* zone_id */
1608 state_flag,
1609 meta_flag,
1610 &input);
1611
1612 /* Set pointers */
1613 pty->defs[id] = d;
1614
1615 if (dim == 1)
1616 pty->get_eval_at_cell[id] = cs_xdef_eval_scalar_at_cells_by_array;
1617 else
1618 pty->get_eval_at_cell[id] = cs_xdef_eval_nd_at_cells_by_array;
1619 pty->get_eval_at_cell_cw[id] = cs_xdef_cw_eval_by_array;
1620
1621 if (cs_flag_test(loc, cs_flag_primal_cell) == false &&
1622 cs_flag_test(loc, cs_flag_primal_vtx) == false &&
1623 cs_flag_test(loc, cs_flag_dual_face_byc) == false)
1624 bft_error(__FILE__, __LINE__, 0,
1625 " %s: case not available.\n", __func__);
1626
1627 /* Set the state flag */
1628 if (cs_flag_test(loc, cs_flag_primal_cell))
1629 pty->state_flag |= CS_FLAG_STATE_CELLWISE;
1630
1631 return d;
1632 }
1633
1634 /*----------------------------------------------------------------------------*/
1635 /*!
1636 * \brief Define a cs_property_t structure thanks to a field structure
1637 *
1638 * \param[in, out] pty pointer to a cs_property_t structure
1639 * \param[in] field pointer to a cs_field_t structure
1640 */
1641 /*----------------------------------------------------------------------------*/
1642
1643 void
cs_property_def_by_field(cs_property_t * pty,cs_field_t * field)1644 cs_property_def_by_field(cs_property_t *pty,
1645 cs_field_t *field)
1646 {
1647 if (pty == NULL)
1648 bft_error(__FILE__, __LINE__, 0, _(_err_empty_pty));
1649
1650 if (field == NULL)
1651 return;
1652
1653 int id = _add_new_def(pty);
1654 int dim = cs_property_get_dim(pty);
1655
1656 /* Sanity checks */
1657 assert(dim == field->dim);
1658 assert(id == 0);
1659 /* z_id = 0 since all the support is selected in this case */
1660
1661 const cs_zone_t *z = cs_volume_zone_by_id(0);
1662 if (field->location_id != z->location_id)
1663 bft_error(__FILE__, __LINE__, 0,
1664 " Property defined by field requests that the field location"
1665 " is supported by cells\n"
1666 " Property %s\n", pty->name);
1667 if (pty->n_definitions > 1)
1668 bft_error(__FILE__, __LINE__, 0,
1669 " When a definition by field is requested, the max. number"
1670 " of subdomains to consider should be equal to 1.\n"
1671 " Current value is %d for property %s.\n"
1672 " Please modify your settings.",
1673 pty->n_definitions, pty->name);
1674
1675 cs_flag_t state_flag = CS_FLAG_STATE_CELLWISE;
1676 cs_flag_t meta_flag = 0; /* metadata */
1677
1678 pty->defs[id] = cs_xdef_volume_create(CS_XDEF_BY_FIELD,
1679 dim,
1680 0, /* zone_id */
1681 state_flag,
1682 meta_flag,
1683 field);
1684
1685 pty->get_eval_at_cell[id] = cs_xdef_eval_cell_by_field;
1686 pty->get_eval_at_cell_cw[id] = cs_xdef_cw_eval_by_field;
1687
1688 /* Set the state flag */
1689 pty->state_flag |= CS_FLAG_STATE_CELLWISE;
1690 }
1691
1692 /*----------------------------------------------------------------------------*/
1693 /*!
1694 * \brief Evaluate the value of the property at each cell. Store the
1695 * evaluation in the given array.
1696 *
1697 * \param[in] t_eval physical time at which one evaluates the term
1698 * \param[in] pty pointer to a cs_property_t structure
1699 * \param[out] pty_stride = 0 if uniform, =1 otherwise
1700 * \param[in, out] pty_vals pointer to an array of values. Allocated if not
1701 * The size of the allocation depends on the value
1702 * of the pty_stride
1703 */
1704 /*----------------------------------------------------------------------------*/
1705
1706 void
cs_property_iso_get_cell_values(cs_real_t t_eval,const cs_property_t * pty,int * pty_stride,cs_real_t ** p_pty_vals)1707 cs_property_iso_get_cell_values(cs_real_t t_eval,
1708 const cs_property_t *pty,
1709 int *pty_stride,
1710 cs_real_t **p_pty_vals)
1711 {
1712 if (pty == NULL)
1713 return;
1714 assert(pty->type & CS_PROPERTY_ISO);
1715
1716 bool allocate = (*p_pty_vals == NULL) ? true : false;
1717 cs_real_t *values = *p_pty_vals;
1718
1719 if (cs_property_is_uniform(pty)) {
1720
1721 *pty_stride = 0;
1722 if (allocate)
1723 BFT_MALLOC(values, 1, cs_real_t);
1724 /* Evaluation at c_id = 0. One assumes that there is at least one cell per
1725 MPI rank */
1726 values[0] = cs_property_get_cell_value(0, t_eval, pty);
1727
1728 }
1729 else {
1730
1731 *pty_stride = 1;
1732 if (allocate)
1733 BFT_MALLOC(values, cs_cdo_quant->n_cells, cs_real_t);
1734 cs_property_eval_at_cells(t_eval, pty, values);
1735
1736 }
1737
1738 /* Return the pointer to values */
1739 *p_pty_vals = values;
1740 }
1741
1742 /*----------------------------------------------------------------------------*/
1743 /*!
1744 * \brief Evaluate the value of the property at each cell. Store the
1745 * evaluation in the given array.
1746 *
1747 * \param[in] t_eval physical time at which one evaluates the term
1748 * \param[in] pty pointer to a cs_property_t structure
1749 * \param[in, out] array pointer to an array of values (must be allocated)
1750 */
1751 /*----------------------------------------------------------------------------*/
1752
1753 void
cs_property_eval_at_cells(cs_real_t t_eval,const cs_property_t * pty,cs_real_t * array)1754 cs_property_eval_at_cells(cs_real_t t_eval,
1755 const cs_property_t *pty,
1756 cs_real_t *array)
1757 {
1758 assert(pty != NULL && array != NULL);
1759
1760 const cs_cdo_quantities_t *quant = cs_cdo_quant;
1761
1762 if (pty->type & CS_PROPERTY_BY_PRODUCT) {
1763
1764 assert(pty->related_properties != NULL);
1765 const cs_property_t *a = pty->related_properties[0];
1766 const cs_property_t *b = pty->related_properties[1];
1767
1768 if (pty->type & CS_PROPERTY_ISO) {
1769
1770 cs_real_t *val_a = NULL;
1771 BFT_MALLOC(val_a, quant->n_cells, cs_real_t);
1772 memset(val_a, 0, quant->n_cells*sizeof(cs_real_t));
1773
1774 /* 1. Evaluates the property A */
1775 for (int i = 0; i < a->n_definitions; i++) {
1776
1777 cs_xdef_t *def = a->defs[i];
1778 const cs_zone_t *z = cs_volume_zone_by_id(def->z_id);
1779
1780 a->get_eval_at_cell[i](z->n_elts,
1781 z->elt_ids,
1782 false, /* without dense output */
1783 cs_glob_mesh,
1784 cs_cdo_connect,
1785 quant,
1786 t_eval,
1787 def->context,
1788 val_a);
1789
1790 } /* Loop on definitions */
1791
1792 /* 2. Evaluates the property B and operates the product */
1793 for (int i = 0; i < b->n_definitions; i++) {
1794
1795 cs_xdef_t *def = b->defs[i];
1796 const cs_zone_t *z = cs_volume_zone_by_id(def->z_id);
1797
1798 b->get_eval_at_cell[i](z->n_elts,
1799 z->elt_ids,
1800 false, /* without dense output */
1801 cs_glob_mesh,
1802 cs_cdo_connect,
1803 quant,
1804 t_eval,
1805 def->context,
1806 array);
1807
1808 for (cs_lnum_t j = 0; j < z->n_elts; j++)
1809 array[z->elt_ids[j]] *= val_a[z->elt_ids[j]];
1810
1811 } /* Loop on definitions */
1812
1813 BFT_FREE(val_a);
1814
1815 }
1816 else {
1817
1818 if (a->type & CS_PROPERTY_ISO) {
1819
1820 cs_real_t *val_a = NULL;
1821 BFT_MALLOC(val_a, quant->n_cells, cs_real_t);
1822 memset(val_a, 0, quant->n_cells*sizeof(cs_real_t));
1823
1824 /* 1. Evaluates the property A */
1825 for (int i = 0; i < a->n_definitions; i++) {
1826
1827 cs_xdef_t *def = a->defs[i];
1828 const cs_zone_t *z = cs_volume_zone_by_id(def->z_id);
1829
1830 a->get_eval_at_cell[i](z->n_elts,
1831 z->elt_ids,
1832 false, /* without dense output */
1833 cs_glob_mesh,
1834 cs_cdo_connect,
1835 quant,
1836 t_eval,
1837 def->context,
1838 val_a);
1839
1840 } /* Loop on definitions */
1841
1842 int b_dim = cs_property_get_dim(b);
1843
1844 /* 2. Evaluates the property B and operates the product */
1845 for (int i = 0; i < b->n_definitions; i++) {
1846
1847 cs_xdef_t *def = b->defs[i];
1848 const cs_zone_t *z = cs_volume_zone_by_id(def->z_id);
1849
1850 b->get_eval_at_cell[i](z->n_elts,
1851 z->elt_ids,
1852 false, /* without dense output */
1853 cs_glob_mesh,
1854 cs_cdo_connect,
1855 quant,
1856 t_eval,
1857 def->context,
1858 array);
1859
1860 for (cs_lnum_t j = 0; j < z->n_elts; j++) {
1861 const cs_real_t acoef = val_a[z->elt_ids[j]];
1862 cs_real_t *_a = array + b_dim*z->elt_ids[j];
1863 for (int k = 0; k < b_dim; k++)
1864 _a[k] *= acoef;
1865 }
1866
1867 } /* Loop on definitions */
1868
1869 BFT_FREE(val_a);
1870
1871 }
1872 else if (b->type & CS_PROPERTY_ISO) {
1873
1874 cs_real_t *val_b = NULL;
1875 BFT_MALLOC(val_b, quant->n_cells, cs_real_t);
1876 memset(val_b, 0, quant->n_cells*sizeof(cs_real_t));
1877
1878 /* 1. Evaluates the property B */
1879 for (int i = 0; i < b->n_definitions; i++) {
1880
1881 cs_xdef_t *def = b->defs[i];
1882 const cs_zone_t *z = cs_volume_zone_by_id(def->z_id);
1883
1884 b->get_eval_at_cell[i](z->n_elts,
1885 z->elt_ids,
1886 false, /* without dense output */
1887 cs_glob_mesh,
1888 cs_cdo_connect,
1889 quant,
1890 t_eval,
1891 def->context,
1892 val_b);
1893
1894 } /* Loop on definitions */
1895
1896 int a_dim = cs_property_get_dim(a);
1897
1898 /* 2. Evaluates the property A and operates the product */
1899 for (int i = 0; i < a->n_definitions; i++) {
1900
1901 cs_xdef_t *def = a->defs[i];
1902 const cs_zone_t *z = cs_volume_zone_by_id(def->z_id);
1903
1904 a->get_eval_at_cell[i](z->n_elts,
1905 z->elt_ids,
1906 false, /* without dense output */
1907 cs_glob_mesh,
1908 cs_cdo_connect,
1909 quant,
1910 t_eval,
1911 def->context,
1912 array);
1913
1914 for (cs_lnum_t j = 0; j < z->n_elts; j++) {
1915 const cs_real_t bcoef = val_b[z->elt_ids[j]];
1916 cs_real_t *_a = array + a_dim*z->elt_ids[j];
1917 for (int k = 0; k < a_dim; k++)
1918 _a[k] *= bcoef;
1919 }
1920
1921 } /* Loop on definitions */
1922
1923 BFT_FREE(val_b);
1924
1925 }
1926 else
1927 bft_error(__FILE__, __LINE__, 0,
1928 " %s: Case not handled yet.\n", __func__);
1929
1930 } /* Either a or b is an isotropic property */
1931
1932 }
1933 else { /* Simple case: One has to evaluate the property */
1934
1935 if ((pty->type & CS_PROPERTY_ISO) && cs_property_is_constant(pty)) {
1936
1937 # pragma omp parallel for if (cs_cdo_connect->n_cells > CS_THR_MIN)
1938 for (cs_lnum_t i = 0; i < cs_cdo_connect->n_cells; i++)
1939 array[i] = pty->ref_value;
1940
1941 }
1942 else {
1943
1944 for (int i = 0; i < pty->n_definitions; i++) {
1945
1946 cs_xdef_t *def = pty->defs[i];
1947 const cs_zone_t *z = cs_volume_zone_by_id(def->z_id);
1948
1949 pty->get_eval_at_cell[i](z->n_elts,
1950 z->elt_ids,
1951 false, /* without dense output */
1952 cs_glob_mesh,
1953 cs_cdo_connect,
1954 quant,
1955 t_eval,
1956 def->context,
1957 array);
1958
1959 } /* Loop on definitions */
1960
1961 } /* Not isotropic and not constant */
1962
1963 } /* Not defined as the product of two existing properties */
1964 }
1965
1966 /*----------------------------------------------------------------------------*/
1967 /*!
1968 * \brief Compute the value of the tensor attached a property at the cell
1969 * center
1970 *
1971 * \param[in] c_id id of the current cell
1972 * \param[in] t_eval physical time at which one evaluates the term
1973 * \param[in] pty pointer to a cs_property_t structure
1974 * \param[in] do_inversion true or false
1975 * \param[in, out] tensor 3x3 matrix
1976 */
1977 /*----------------------------------------------------------------------------*/
1978
1979 void
cs_property_get_cell_tensor(cs_lnum_t c_id,cs_real_t t_eval,const cs_property_t * pty,bool do_inversion,cs_real_t tensor[3][3])1980 cs_property_get_cell_tensor(cs_lnum_t c_id,
1981 cs_real_t t_eval,
1982 const cs_property_t *pty,
1983 bool do_inversion,
1984 cs_real_t tensor[3][3])
1985 {
1986 if (pty == NULL)
1987 return;
1988
1989 /* Initialize extra-diag. values of the tensor */
1990 tensor[0][1] = tensor[1][0] = tensor[2][0] = 0;
1991 tensor[0][2] = tensor[1][2] = tensor[2][1] = 0;
1992
1993 if (pty->type & CS_PROPERTY_BY_PRODUCT)
1994 _get_cell_tensor_by_property_product(c_id, t_eval, pty, tensor);
1995 else
1996 _get_cell_tensor(c_id, t_eval, pty, tensor);
1997
1998 if (do_inversion)
1999 _invert_tensor(tensor, pty->type);
2000
2001 #if defined(DEBUG) && !defined(NDEBUG) && CS_PROPERTY_DBG > 1
2002 cs_log_printf(CS_LOG_DEFAULT, "\n pty: %s, cell_id: %d\n", pty->name, c_id);
2003 _print_tensor(tensor);
2004 #endif
2005 }
2006
2007 /*----------------------------------------------------------------------------*/
2008 /*!
2009 * \brief Compute the value of a property at the cell center
2010 *
2011 * \param[in] c_id id of the current cell
2012 * \param[in] t_eval physical time at which one evaluates the term
2013 * \param[in] pty pointer to a cs_property_t structure
2014 *
2015 * \return the value of the property for the given cell
2016 */
2017 /*----------------------------------------------------------------------------*/
2018
2019 cs_real_t
cs_property_get_cell_value(cs_lnum_t c_id,cs_real_t t_eval,const cs_property_t * pty)2020 cs_property_get_cell_value(cs_lnum_t c_id,
2021 cs_real_t t_eval,
2022 const cs_property_t *pty)
2023 {
2024 cs_real_t result = 0;
2025
2026 if (pty == NULL)
2027 return result;
2028
2029 if ((pty->type & CS_PROPERTY_ISO) == 0)
2030 bft_error(__FILE__, __LINE__, 0,
2031 " %s: Invalid type of property for this function.\n"
2032 " Property %s has to be isotropic.", __func__, pty->name);
2033
2034 if (pty->type & CS_PROPERTY_BY_PRODUCT) {
2035
2036 assert(pty->related_properties != NULL);
2037 const cs_real_t result_a = _get_cell_value(c_id, t_eval,
2038 pty->related_properties[0]);
2039 const cs_real_t result_b = _get_cell_value(c_id, t_eval,
2040 pty->related_properties[1]);
2041
2042 return result_a * result_b;
2043
2044 }
2045 else {
2046
2047 if (cs_property_is_constant(pty))
2048 return pty->ref_value;
2049 else
2050 return _get_cell_value(c_id, t_eval, pty);
2051
2052 }
2053 }
2054
2055 /*----------------------------------------------------------------------------*/
2056 /*!
2057 * \brief Compute the value of the tensor attached to a property at the cell
2058 * center
2059 * Version using a cs_cell_mesh_t structure
2060 *
2061 * \param[in] cm pointer to a cs_cell_mesh_t structure
2062 * \param[in] pty pointer to a cs_property_t structure
2063 * \param[in] t_eval physical time at which one evaluates the term
2064 * \param[in] do_inversion true or false
2065 * \param[in, out] tensor 3x3 matrix
2066 */
2067 /*----------------------------------------------------------------------------*/
2068
2069 void
cs_property_tensor_in_cell(const cs_cell_mesh_t * cm,const cs_property_t * pty,cs_real_t t_eval,bool do_inversion,cs_real_t tensor[3][3])2070 cs_property_tensor_in_cell(const cs_cell_mesh_t *cm,
2071 const cs_property_t *pty,
2072 cs_real_t t_eval,
2073 bool do_inversion,
2074 cs_real_t tensor[3][3])
2075 {
2076 if (pty == NULL)
2077 return;
2078
2079 /* Initialize extra-diag. values of the tensor */
2080 tensor[0][1] = tensor[1][0] = tensor[2][0] = 0;
2081 tensor[0][2] = tensor[1][2] = tensor[2][1] = 0;
2082
2083 if (pty->type & CS_PROPERTY_BY_PRODUCT)
2084 _tensor_in_cell_by_property_product(cm, pty, t_eval, tensor);
2085 else
2086 _tensor_in_cell(cm, pty, t_eval, tensor);
2087
2088 if (do_inversion)
2089 _invert_tensor(tensor, pty->type);
2090
2091 #if defined(DEBUG) && !defined(NDEBUG) && CS_PROPERTY_DBG > 1
2092 cs_log_printf(CS_LOG_DEFAULT, "\n pty: %s, cell_id: %d\n",
2093 pty->name, cm->c_id);
2094 _print_tensor(tensor);
2095 #endif
2096 }
2097
2098 /*----------------------------------------------------------------------------*/
2099 /*!
2100 * \brief Compute the value of a property at the cell center
2101 * Version using a cs_cell_mesh_t structure
2102 *
2103 * \param[in] cm pointer to a cs_cell_mesh_t structure
2104 * \param[in] pty pointer to a cs_property_t structure
2105 * \param[in] t_eval physical time at which one evaluates the term
2106 *
2107 * \return the value of the property for the given cell
2108 */
2109 /*----------------------------------------------------------------------------*/
2110
2111 cs_real_t
cs_property_value_in_cell(const cs_cell_mesh_t * cm,const cs_property_t * pty,cs_real_t t_eval)2112 cs_property_value_in_cell(const cs_cell_mesh_t *cm,
2113 const cs_property_t *pty,
2114 cs_real_t t_eval)
2115 {
2116 cs_real_t result = 0;
2117
2118 if (pty == NULL)
2119 return result;
2120
2121 if ((pty->type & CS_PROPERTY_ISO) == 0)
2122 bft_error(__FILE__, __LINE__, 0,
2123 " Invalid type of property for this function.\n"
2124 " Property %s has to be isotropic.", pty->name);
2125
2126 if (pty->type & CS_PROPERTY_BY_PRODUCT) {
2127 assert(pty->related_properties != NULL);
2128 const cs_real_t result_a = _value_in_cell(cm, pty->related_properties[0],
2129 t_eval);
2130 const cs_real_t result_b = _value_in_cell(cm, pty->related_properties[1],
2131 t_eval);
2132 return result_a * result_b;
2133 }
2134 else {
2135
2136 if (cs_property_is_constant(pty))
2137 return pty->ref_value;
2138 else
2139 return _value_in_cell(cm, pty, t_eval);
2140
2141 }
2142 }
2143
2144 /*----------------------------------------------------------------------------*/
2145 /*!
2146 * \brief Compute the Fourier number in each cell
2147 *
2148 * \param[in] pty pointer to the diffusive property struct.
2149 * \param[in] t_eval physical time at which one evaluates the term
2150 * \param[in] dt value of the current time step
2151 * \param[in, out] fourier pointer to an array storing Fourier numbers
2152 */
2153 /*----------------------------------------------------------------------------*/
2154
2155 void
cs_property_get_fourier(const cs_property_t * pty,cs_real_t t_eval,double dt,cs_real_t fourier[])2156 cs_property_get_fourier(const cs_property_t *pty,
2157 cs_real_t t_eval,
2158 double dt,
2159 cs_real_t fourier[])
2160 {
2161 assert(fourier != NULL); /* Sanity check */
2162 assert(dt > 0.);
2163
2164 const bool pty_uniform = cs_property_is_uniform(pty);
2165 const cs_cdo_quantities_t *cdoq = cs_cdo_quant;
2166
2167 if (pty->type & CS_PROPERTY_ISO) {
2168
2169 cs_real_t ptyval = 0.;
2170 if (pty_uniform)
2171 ptyval = cs_property_get_cell_value(0, t_eval, pty);
2172
2173 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++) {
2174
2175 const cs_real_t hc = cbrt(cdoq->cell_vol[c_id]);
2176 if (!pty_uniform)
2177 ptyval = cs_property_get_cell_value(c_id, t_eval, pty);
2178
2179 fourier[c_id] = dt * ptyval / (hc*hc);
2180
2181 }
2182
2183 }
2184 else { /* Property is orthotropic or anisotropic */
2185
2186 cs_real_t eig_max, eig_ratio;
2187 cs_real_t ptymat[3][3];
2188
2189 /* Get the value of the material property at the first cell center */
2190 if (pty_uniform) {
2191 cs_property_get_cell_tensor(0, t_eval, pty, false, ptymat);
2192 cs_math_33_eigen((const cs_real_t (*)[3])ptymat, &eig_ratio, &eig_max);
2193 }
2194
2195 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++) {
2196
2197 const cs_real_t hc = cbrt(cdoq->cell_vol[c_id]);
2198
2199 /* Get the value of the material property at the cell center */
2200 if (!pty_uniform) {
2201 cs_property_get_cell_tensor(c_id, t_eval, pty, false, ptymat);
2202 cs_math_33_eigen((const cs_real_t (*)[3])ptymat, &eig_ratio, &eig_max);
2203 }
2204
2205 fourier[c_id] = dt * eig_max / (hc*hc);
2206
2207 }
2208
2209 } /* Type of property */
2210
2211 }
2212
2213 /*----------------------------------------------------------------------------*/
2214 /*!
2215 * \brief Print a summary of the settings for all defined cs_property_t
2216 * structures
2217 */
2218 /*----------------------------------------------------------------------------*/
2219
2220 void
cs_property_log_setup(void)2221 cs_property_log_setup(void)
2222 {
2223 if (_n_properties == 0)
2224 return;
2225 assert(_properties != NULL);
2226
2227 cs_log_printf(CS_LOG_SETUP, "\nSummary of the definition of properties\n");
2228 cs_log_printf(CS_LOG_SETUP, "%s", cs_sep_h1);
2229
2230 char prefix[256];
2231
2232 for (int i = 0; i < _n_properties; i++) {
2233
2234 bool is_uniform = false, is_steady = false;
2235 const cs_property_t *pty = _properties[i];
2236
2237 if (pty == NULL)
2238 continue;
2239 assert(strlen(pty->name) < 200); /* Check that prefix is large enough */
2240
2241 if (pty->state_flag & CS_FLAG_STATE_UNIFORM) is_uniform = true;
2242 if (pty->state_flag & CS_FLAG_STATE_STEADY) is_steady = true;
2243
2244 cs_log_printf(CS_LOG_SETUP, "\n * %s | Uniform %s Steady %s\n", pty->name,
2245 cs_base_strtf(is_uniform), cs_base_strtf(is_steady));
2246 cs_log_printf(CS_LOG_SETUP, " * %s | Reference value % -8.4e\n",
2247 pty->name, pty->ref_value);
2248
2249 if (pty->type & CS_PROPERTY_ISO)
2250 cs_log_printf(CS_LOG_SETUP, " * %s | Type: isotropic", pty->name);
2251 else if (pty->type & CS_PROPERTY_ORTHO)
2252 cs_log_printf(CS_LOG_SETUP, " * %s | Type: orthotropic", pty->name);
2253 else if (pty->type & CS_PROPERTY_ANISO_SYM)
2254 cs_log_printf(CS_LOG_SETUP,
2255 " * %s | Type: anisotropic with a symmetric storage",
2256 pty->name);
2257 else if (pty->type & CS_PROPERTY_ANISO)
2258 cs_log_printf(CS_LOG_SETUP, " * %s | Type: anisotropic", pty->name);
2259 else
2260 bft_error(__FILE__, __LINE__, 0, _("%s: Invalid type of property."),
2261 __func__);
2262
2263 if (pty->type & CS_PROPERTY_BY_PRODUCT)
2264 cs_log_printf(CS_LOG_SETUP, " | by product\n");
2265 else
2266 cs_log_printf(CS_LOG_SETUP, "\n");
2267
2268 cs_log_printf(CS_LOG_SETUP, " * %s | Number of definitions: %d\n\n",
2269 pty->name, pty->n_definitions);
2270
2271 for (int j = 0; j < pty->n_definitions; j++) {
2272 sprintf(prefix, " Definition %3d", j);
2273 cs_xdef_log(prefix, pty->defs[j]);
2274 }
2275
2276 } /* Loop on properties */
2277
2278 }
2279
2280 /*----------------------------------------------------------------------------*/
2281
2282 END_C_DECLS
2283