1 /*============================================================================
2 * Functions and structures to deal with source term computations
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 #include "cs_defs.h"
26
27 /*----------------------------------------------------------------------------
28 * Standard C library headers
29 *----------------------------------------------------------------------------*/
30
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include <string.h>
34 #include <assert.h>
35 #include <float.h>
36
37 /*----------------------------------------------------------------------------
38 * Local headers
39 *----------------------------------------------------------------------------*/
40
41 #include <bft_mem.h>
42
43 #include "cs_basis_func.h"
44 #include "cs_evaluate.h"
45 #include "cs_hho_builder.h"
46 #include "cs_hodge.h"
47 #include "cs_log.h"
48 #include "cs_math.h"
49 #include "cs_scheme_geometry.h"
50 #include "cs_volume_zone.h"
51
52 /*----------------------------------------------------------------------------
53 * Header for the current file
54 *----------------------------------------------------------------------------*/
55
56 #include "cs_source_term.h"
57
58 /*----------------------------------------------------------------------------*/
59
60 BEGIN_C_DECLS
61
62 /*=============================================================================
63 * Local Macro definitions and structure definitions
64 *============================================================================*/
65
66 #define CS_SOURCE_TERM_DBG 0
67
68 /*============================================================================
69 * Private variables
70 *============================================================================*/
71
72 static const char _err_empty_st[] =
73 " Stop setting an empty cs_xdef_t structure.\n"
74 " Please check your settings.\n";
75
76 /* Pointer to shared structures (owned by a cs_domain_t structure) */
77 static const cs_cdo_quantities_t *cs_cdo_quant;
78 static const cs_cdo_connect_t *cs_cdo_connect;
79
80 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
81
82 /*============================================================================
83 * Private function prototypes
84 *============================================================================*/
85
86 /*----------------------------------------------------------------------------*/
87 /*!
88 * \brief Allocate and initialize a name (copy or generic name)
89 *
90 * \param[in] name input name
91 * \param[in] base_name generic name by default if input name is NULL
92 * \param[in] id id related to this name
93 *
94 * \return a pointer to a new allocated string
95 */
96 /*----------------------------------------------------------------------------*/
97
98 static inline char *
_get_name(const char * name,const char * base_name,int id)99 _get_name(const char *name,
100 const char *base_name,
101 int id)
102 {
103 char *n = NULL;
104
105 if (name == NULL) { /* Define a name by default */
106 assert(id < 100);
107 int len = strlen(base_name) + 4; // 1 + 3 = "_00\n"
108 BFT_MALLOC(n, len, char);
109 sprintf(n, "%s_%2d", base_name, id);
110 }
111 else { /* Copy name */
112 size_t len = strlen(name);
113 BFT_MALLOC(n, len + 1, char);
114 strncpy(n, name, len);
115 n[len] = '\0';
116 }
117
118 return n;
119 }
120
121 /*----------------------------------------------------------------------------*/
122 /*!
123 * \brief Update the mask associated to each cell from the mask related to
124 * the given source term structure
125 *
126 * \param[in] st pointer to a cs_xdef_t structure
127 * \param[in] st_id id related to this source term
128 * \param[in, out] cell_mask mask related to each cell to be updated
129 */
130 /*----------------------------------------------------------------------------*/
131
132 static void
_set_mask(const cs_xdef_t * st,int st_id,cs_mask_t * cell_mask)133 _set_mask(const cs_xdef_t *st,
134 int st_id,
135 cs_mask_t *cell_mask)
136 {
137 if (st == NULL)
138 bft_error(__FILE__, __LINE__, 0, _(_err_empty_st));
139
140 /* value of the mask for the source term */
141 const cs_mask_t mask = (1 << st_id);
142
143 if (st->meta & CS_FLAG_FULL_LOC) /* All cells are selected */
144 # pragma omp parallel for if (cs_cdo_quant->n_cells > CS_THR_MIN)
145 for (cs_lnum_t i = 0; i < cs_cdo_quant->n_cells; i++) cell_mask[i] |= mask;
146
147 else {
148
149 /* Retrieve information from the volume zone structure */
150 const cs_zone_t *z = cs_volume_zone_by_id(st->z_id);
151 for (cs_lnum_t i = 0; i < z->n_elts; i++)
152 cell_mask[z->elt_ids[i]] |= mask;
153
154 }
155
156 }
157
158 /*----------------------------------------------------------------------------*/
159 /*!
160 * \brief Compute the reduction onto the cell polynomial space of a function
161 * defined by a constant value
162 *
163 * \param[in] const_val constant value
164 * \param[in] cbf pointer to a structure for face basis functions
165 * \param[in] xv1 first vertex
166 * \param[in] xv2 second vertex
167 * \param[in] xv3 third vertex
168 * \param[in] xv4 fourth vertex
169 * \param[in] vol volume of the tetrahedron
170 * \param[in, out] cb pointer to a cs_cell_builder_structure_t
171 * \param[in, out] array array storing values to compute
172 */
173 /*----------------------------------------------------------------------------*/
174
175 static void
_hho_add_tetra_by_val(cs_real_t const_val,const cs_basis_func_t * cbf,const cs_real_3_t xv1,const cs_real_3_t xv2,const cs_real_3_t xv3,const cs_real_3_t xv4,const double vol,cs_cell_builder_t * cb,cs_real_t array[])176 _hho_add_tetra_by_val(cs_real_t const_val,
177 const cs_basis_func_t *cbf,
178 const cs_real_3_t xv1,
179 const cs_real_3_t xv2,
180 const cs_real_3_t xv3,
181 const cs_real_3_t xv4,
182 const double vol,
183 cs_cell_builder_t *cb,
184 cs_real_t array[])
185 {
186 cs_real_3_t *gpts = cb->vectors;
187 cs_real_t *gw = cb->values;
188 cs_real_t *phi_eval = cb->values + 15;
189
190 /* Compute Gauss points and related weights */
191 cs_quadrature_tet_15pts(xv1, xv2, xv3, xv4, vol, gpts, gw);
192
193 for (short int gp = 0; gp < 15; gp++) {
194
195 cbf->eval_all_at_point(cbf, gpts[gp], phi_eval);
196
197 const cs_real_t w = gw[gp] * const_val;
198 for (short int i = 0; i < cbf->size; i++)
199 array[i] += w * phi_eval[i];
200
201 } /* End of loop on Gauss points */
202 }
203
204 /*----------------------------------------------------------------------------*/
205 /*!
206 * \brief Compute the reduction onto the cell polynomial space of a function
207 * defined by an analytical expression depending on the location and
208 * the current time
209 *
210 * \param[in] ac pointer to an analytical definition
211 * \param[in] cbf pointer to a structure for face basis functions
212 * \param[in] xv1 first vertex
213 * \param[in] xv2 second vertex
214 * \param[in] xv3 third vertex
215 * \param[in] xv4 fourth vertex
216 * \param[in] vol volume of the tetrahedron
217 * \param[in] time_eval physical time at which one evaluates the term
218 * \param[in, out] cb pointer to a cs_cell_builder_structure_t
219 * \param[in, out] array array storing values to compute
220 */
221 /*----------------------------------------------------------------------------*/
222
223 static void
_hho_add_tetra_by_ana(const cs_xdef_analytic_context_t * ac,const cs_basis_func_t * cbf,const cs_real_3_t xv1,const cs_real_3_t xv2,const cs_real_3_t xv3,const cs_real_3_t xv4,const double vol,cs_real_t time_eval,cs_cell_builder_t * cb,cs_real_t array[])224 _hho_add_tetra_by_ana(const cs_xdef_analytic_context_t *ac,
225 const cs_basis_func_t *cbf,
226 const cs_real_3_t xv1,
227 const cs_real_3_t xv2,
228 const cs_real_3_t xv3,
229 const cs_real_3_t xv4,
230 const double vol,
231 cs_real_t time_eval,
232 cs_cell_builder_t *cb,
233 cs_real_t array[])
234 {
235 cs_real_3_t *gpts = cb->vectors;
236 cs_real_t *gw = cb->values;
237 cs_real_t *ana_eval = cb->values + 15;
238 cs_real_t *phi_eval = cb->values + 30;
239
240 /* Compute Gauss points and related weights */
241 cs_quadrature_tet_15pts(xv1, xv2, xv3, xv4, vol, gpts, gw);
242
243 /* Evaluate the analytical function at the Gauss points */
244 ac->func(time_eval, 15, NULL, (const cs_real_t *)gpts, true,
245 ac->input, ana_eval);
246
247 for (short int gp = 0; gp < 15; gp++) {
248
249 cbf->eval_all_at_point(cbf, gpts[gp], phi_eval);
250
251 const cs_real_t w = gw[gp] * ana_eval[gp];
252 for (short int i = 0; i < cbf->size; i++)
253 array[i] += w * phi_eval[i];
254
255 } /* End of loop on Gauss points */
256 }
257
258
259 /*----------------------------------------------------------------------------*/
260 /*!
261 * \brief Compute the reduction onto the cell polynomial space of a function
262 * defined by an analytical expression depending on the location and
263 * the current time. This is the vector case.
264 *
265 * \param[in] ac pointer to an analytical definition
266 * \param[in] cbf pointer to a structure for cell basis functions
267 * \param[in] xv1 first vertex
268 * \param[in] xv2 second vertex
269 * \param[in] xv3 third vertex
270 * \param[in] xv4 third vertex
271 * \param[in] vol volume of the tetrahedron
272 * \param[in] time_eval physical time at which one evaluates the term
273 * \param[in, out] cb pointer to a cs_cell_builder_structure_t
274 * \param[in, out] array array storing values to compute
275 */
276 /*----------------------------------------------------------------------------*/
277
278 static void
_hho_add_tetra_by_ana_vd(const cs_xdef_analytic_context_t * ac,const cs_basis_func_t * cbf,const cs_real_3_t xv1,const cs_real_3_t xv2,const cs_real_3_t xv3,const cs_real_3_t xv4,const double vol,cs_real_t time_eval,cs_cell_builder_t * cb,cs_real_t array[])279 _hho_add_tetra_by_ana_vd(const cs_xdef_analytic_context_t *ac,
280 const cs_basis_func_t *cbf,
281 const cs_real_3_t xv1,
282 const cs_real_3_t xv2,
283 const cs_real_3_t xv3,
284 const cs_real_3_t xv4,
285 const double vol,
286 cs_real_t time_eval,
287 cs_cell_builder_t *cb,
288 cs_real_t array[])
289 {
290 cs_real_3_t *gpts = cb->vectors;
291
292 /* cb->values is big enough to store:
293 gw+ana_eval+phi_eval= 15+3*15+cell_basis
294 = 64 < 78 + 12 for k=1
295 = 70 < 465 + 30 for k=2
296 */
297 cs_real_t *gw = cb->values;
298 cs_real_t *ana_eval = cb->values + 15;
299 cs_real_t *phi_eval = cb->values + 15 + 3*15;
300
301 /* Compute Gauss points and related weights */
302 cs_quadrature_tet_15pts(xv1, xv2, xv3, xv4, vol, gpts, gw);
303
304 /* Evaluate the analytical function at the Gauss points */
305 ac->func(time_eval, 15, NULL, (const cs_real_t *)gpts, true, ac->input,
306 ana_eval);
307
308 for (short int gp = 0; gp < 15; gp++) {
309
310 cbf->eval_all_at_point(cbf, gpts[gp], phi_eval);
311
312 for (short int i = 0; i < cbf->size; i++) {
313
314 const double gcoef = gw[gp] * phi_eval[i];
315
316 /* x-component */
317 array[i ] += gcoef * ana_eval[3*gp];
318 /* y-component */
319 array[i + cbf->size] += gcoef * ana_eval[3*gp+1];
320 /* z-component */
321 array[i + 2*cbf->size] += gcoef * ana_eval[3*gp+2];
322
323 }
324
325 } /* End of loop on Gauss points */
326
327 }
328
329 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
330
331 /*============================================================================
332 * Public function prototypes
333 *============================================================================*/
334
335 /*----------------------------------------------------------------------------*/
336 /*!
337 * \brief Set shared pointers to main domain members
338 *
339 * \param[in] quant additional mesh quantities struct.
340 * \param[in] connect pointer to a cs_cdo_connect_t struct.
341 */
342 /*----------------------------------------------------------------------------*/
343
344 void
cs_source_term_set_shared_pointers(const cs_cdo_quantities_t * quant,const cs_cdo_connect_t * connect)345 cs_source_term_set_shared_pointers(const cs_cdo_quantities_t *quant,
346 const cs_cdo_connect_t *connect)
347 {
348 /* Assign static const pointers */
349 cs_cdo_quant = quant;
350 cs_cdo_connect = connect;
351 }
352
353 /*----------------------------------------------------------------------------*/
354 /*!
355 * \brief Set the default flag related to a source term according to the
356 * numerical scheme chosen for discretizing an equation
357 *
358 * \param[in] scheme numerical scheme used for the discretization
359 * \param[out] state_flag flag describing the status of the source term
360 * \param[out] meta_flag additional flags associated to a source term
361 */
362 /*----------------------------------------------------------------------------*/
363
364 void
cs_source_term_set_default_flag(cs_param_space_scheme_t scheme,cs_flag_t * state_flag,cs_flag_t * meta_flag)365 cs_source_term_set_default_flag(cs_param_space_scheme_t scheme,
366 cs_flag_t *state_flag,
367 cs_flag_t *meta_flag)
368 {
369 switch (scheme) {
370
371 case CS_SPACE_SCHEME_CDOVB:
372 *state_flag = CS_FLAG_STATE_DENSITY;
373 *meta_flag = cs_flag_dual_cell; /* Predefined mask */
374 break;
375
376 case CS_SPACE_SCHEME_CDOEB:
377 *state_flag = CS_FLAG_STATE_FLUX;
378 *meta_flag = cs_flag_dual_face; /* Predefined mask */
379 break;
380
381 case CS_SPACE_SCHEME_CDOFB:
382 *state_flag = CS_FLAG_STATE_DENSITY;
383 *meta_flag = cs_flag_primal_cell; /* Predefined mask */
384 break;
385
386 case CS_SPACE_SCHEME_CDOVCB:
387 case CS_SPACE_SCHEME_HHO_P0:
388 case CS_SPACE_SCHEME_HHO_P1:
389 case CS_SPACE_SCHEME_HHO_P2:
390 *state_flag = CS_FLAG_STATE_DENSITY;
391 *meta_flag = CS_FLAG_PRIMAL;
392 break;
393
394 default:
395 bft_error(__FILE__, __LINE__, 0,
396 " %s: Invalid numerical scheme to set a source term.",
397 __func__);
398
399 }
400 }
401
402 /*----------------------------------------------------------------------------*/
403 /*!
404 * \brief Set advanced parameters which are defined by default in a
405 * source term structure.
406 *
407 * \param[in, out] st pointer to a cs_xdef_t structure
408 * \param[in] flag CS_FLAG_DUAL or CS_FLAG_PRIMAL
409 */
410 /*----------------------------------------------------------------------------*/
411
412 void
cs_source_term_set_reduction(cs_xdef_t * st,cs_flag_t flag)413 cs_source_term_set_reduction(cs_xdef_t *st,
414 cs_flag_t flag)
415 {
416 if (st == NULL)
417 bft_error(__FILE__, __LINE__, 0, _(_err_empty_st));
418
419 if (st->meta & flag)
420 return; /* Nothing to do */
421
422 cs_flag_t save_meta = st->meta;
423
424 st->meta = 0;
425 /* Set unchanged parts of the existing flag */
426 if (save_meta & CS_FLAG_SCALAR) st->meta |= CS_FLAG_SCALAR;
427 if (save_meta & CS_FLAG_VECTOR) st->meta |= CS_FLAG_VECTOR;
428 if (save_meta & CS_FLAG_TENSOR) st->meta |= CS_FLAG_TENSOR;
429 if (save_meta & CS_FLAG_BORDER) st->meta |= CS_FLAG_BORDER;
430 if (save_meta & CS_FLAG_BY_CELL) st->meta |= CS_FLAG_BY_CELL;
431 if (save_meta & CS_FLAG_FULL_LOC) st->meta |= CS_FLAG_FULL_LOC;
432
433 if (flag & CS_FLAG_DUAL) {
434 assert(save_meta & CS_FLAG_PRIMAL);
435 if (save_meta & CS_FLAG_VERTEX)
436 st->meta |= CS_FLAG_DUAL | CS_FLAG_CELL;
437 else
438 bft_error(__FILE__, __LINE__, 0,
439 " %s: Stop modifying the source term flag.\n"
440 " This case is not handled.", __func__);
441 }
442 else if (flag & CS_FLAG_PRIMAL) {
443 assert(save_meta & CS_FLAG_DUAL);
444 if (save_meta & CS_FLAG_CELL)
445 st->meta |= CS_FLAG_PRIMAL | CS_FLAG_VERTEX;
446 else
447 bft_error(__FILE__, __LINE__, 0,
448 " Stop modifying the source term flag.\n"
449 " This case is not handled.");
450 }
451 else
452 bft_error(__FILE__, __LINE__, 0,
453 " Stop modifying the source term flag.\n"
454 " This case is not handled.");
455 }
456
457 /*----------------------------------------------------------------------------*/
458 /*!
459 * \brief Get metadata related to the given source term structure
460 *
461 * \param[in, out] st pointer to a cs_xdef_t structure
462 *
463 * \return the value of the flag related to this source term
464 */
465 /*----------------------------------------------------------------------------*/
466
467 cs_flag_t
cs_source_term_get_flag(const cs_xdef_t * st)468 cs_source_term_get_flag(const cs_xdef_t *st)
469 {
470 if (st == NULL)
471 bft_error(__FILE__, __LINE__, 0, _(_err_empty_st));
472
473 return st->meta;
474 }
475
476 /*----------------------------------------------------------------------------*/
477 /*!
478 * \brief Initialize data to build the source terms
479 *
480 * \param[in] space_scheme scheme used to discretize in space
481 * \param[in] n_source_terms number of source terms
482 * \param[in] source_terms pointer to the definitions of source terms
483 * \param[in, out] compute_source array of function pointers
484 * \param[in, out] sys_flag metadata about the algebraic system
485 * \param[in, out] source_mask pointer to an array storing in a compact way
486 * which source term is defined in a given cell
487 *
488 * \return a flag which indicates what to build in a cell mesh structure
489 */
490 /*----------------------------------------------------------------------------*/
491
492 cs_eflag_t
cs_source_term_init(cs_param_space_scheme_t space_scheme,const int n_source_terms,cs_xdef_t * const * source_terms,cs_source_term_cellwise_t * compute_source[],cs_flag_t * sys_flag,cs_mask_t * source_mask[])493 cs_source_term_init(cs_param_space_scheme_t space_scheme,
494 const int n_source_terms,
495 cs_xdef_t *const *source_terms,
496 cs_source_term_cellwise_t *compute_source[],
497 cs_flag_t *sys_flag,
498 cs_mask_t *source_mask[])
499 {
500 if (n_source_terms > CS_N_MAX_SOURCE_TERMS)
501 bft_error(__FILE__, __LINE__, 0,
502 " Limitation to %d source terms has been reached!",
503 CS_N_MAX_SOURCE_TERMS);
504
505 cs_eflag_t msh_flag = 0;
506 *source_mask = NULL;
507 for (short int i = 0; i < CS_N_MAX_SOURCE_TERMS; i++)
508 compute_source[i] = NULL;
509
510 if (n_source_terms == 0)
511 return msh_flag;
512
513 bool need_mask = false;
514
515 for (int st_id = 0; st_id < n_source_terms; st_id++) {
516
517 const cs_xdef_t *st_def = source_terms[st_id];
518
519 if (st_def->meta & CS_FLAG_PRIMAL) {
520 if (space_scheme == CS_SPACE_SCHEME_CDOVB ||
521 space_scheme == CS_SPACE_SCHEME_CDOVCB) {
522 msh_flag |= CS_FLAG_COMP_PVQ | CS_FLAG_COMP_DEQ | CS_FLAG_COMP_PFQ |
523 CS_FLAG_COMP_EV | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_HFQ;
524 *sys_flag |= CS_FLAG_SYS_MASS_MATRIX | CS_FLAG_SYS_SOURCES_HLOC;
525 }
526 else if (space_scheme == CS_SPACE_SCHEME_CDOEB) {
527 msh_flag |= CS_FLAG_COMP_DFQ | CS_FLAG_COMP_EV;
528 *sys_flag |= CS_FLAG_SYS_MASS_MATRIX | CS_FLAG_SYS_SOURCES_HLOC;
529 }
530 }
531
532 /* Not defined on the whole mesh */
533 if ((st_def->meta & CS_FLAG_FULL_LOC) == 0)
534 need_mask = true;
535
536 switch (space_scheme) {
537
538 case CS_SPACE_SCHEME_CDOVB:
539
540 msh_flag |= CS_FLAG_COMP_PV;
541 if (st_def->meta & CS_FLAG_DUAL) {
542
543 switch (st_def->type) {
544
545 case CS_XDEF_BY_VALUE:
546 msh_flag |= CS_FLAG_COMP_PVQ;
547 if ((*sys_flag) & CS_FLAG_SYS_VECTOR)
548 compute_source[st_id] = cs_source_term_dcvd_by_value;
549 else
550 compute_source[st_id] = cs_source_term_dcsd_by_value;
551 break;
552
553 case CS_XDEF_BY_ANALYTIC_FUNCTION:
554
555 switch (st_def->qtype) {
556
557 case CS_QUADRATURE_BARY:
558 msh_flag |=
559 CS_FLAG_COMP_PVQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_PFC |
560 CS_FLAG_COMP_FE | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV |
561 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PEC | CS_FLAG_COMP_DEQ;
562 compute_source[st_id] = cs_source_term_dcsd_bary_by_analytic;
563 break;
564
565 case CS_QUADRATURE_BARY_SUBDIV:
566 msh_flag |=
567 CS_FLAG_COMP_EV | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_PFC |
568 CS_FLAG_COMP_FE | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_DEQ;
569 compute_source[st_id] = cs_source_term_dcsd_q1o1_by_analytic;
570 break;
571
572 case CS_QUADRATURE_HIGHER:
573 msh_flag |=
574 CS_FLAG_COMP_PFQ | CS_FLAG_COMP_PFC | CS_FLAG_COMP_FE |
575 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV | CS_FLAG_COMP_PVQ |
576 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DEQ;
577 compute_source[st_id] = cs_source_term_dcsd_q10o2_by_analytic;
578 break;
579
580 case CS_QUADRATURE_HIGHEST:
581 msh_flag |=
582 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
583 CS_FLAG_COMP_EV | CS_FLAG_COMP_PFC | CS_FLAG_COMP_FEQ |
584 CS_FLAG_COMP_DEQ;
585 compute_source[st_id] = cs_source_term_dcsd_q5o3_by_analytic;
586 break;
587
588 default:
589 bft_error(__FILE__, __LINE__, 0,
590 " Invalid type of quadrature for computing a source term"
591 " with CDOVB schemes");
592 } /* quad_type */
593 break;
594
595 case CS_XDEF_BY_ARRAY:
596 msh_flag |= CS_FLAG_COMP_PVQ;
597 compute_source[st_id] = cs_source_term_dcsd_by_array;
598 break;
599
600 case CS_XDEF_BY_DOF_FUNCTION:
601 msh_flag |= CS_FLAG_COMP_PVQ;
602 compute_source[st_id] = cs_source_term_dcsd_by_dof_func;
603 break;
604
605 default:
606 bft_error(__FILE__, __LINE__, 0,
607 "%s: Invalid type of definition for a source term in CDOVB",
608 __func__);
609 break;
610 } /* switch def_type */
611
612 }
613 else {
614 assert(st_def->meta & CS_FLAG_PRIMAL);
615
616 switch (st_def->type) {
617
618 case CS_XDEF_BY_VALUE:
619 msh_flag |= CS_FLAG_COMP_PV;
620 compute_source[st_id] = cs_source_term_pvsp_by_value;
621 break;
622
623 case CS_XDEF_BY_ANALYTIC_FUNCTION:
624 msh_flag |= CS_FLAG_COMP_PV;
625 compute_source[st_id] = cs_source_term_pvsp_by_analytic;
626 break;
627
628 default:
629 bft_error(__FILE__, __LINE__, 0,
630 "%s: Invalid type of definition for a source term in CDOVB",
631 __func__);
632 break;
633
634 } /* switch def_type */
635
636 } /* flag PRIMAL or DUAL */
637 break; /* CDOVB */
638
639 case CS_SPACE_SCHEME_CDOVCB:
640 if (st_def->meta & CS_FLAG_DUAL) {
641
642 bft_error(__FILE__, __LINE__, 0,
643 "%s: Invalid type of definition for a source term in CDOVB",
644 __func__);
645
646 /* TODO:
647 case CS_XDEF_BY_VALUE:
648 cs_source_term_vcsd_by_value; --> case CS_QUADRATURE_BARY:
649
650 case CS_XDEF_BY_ANALYTIC_FUNCTION:
651 cs_source_term_vcsd_q1o1_by_analytic; --> case CS_QUADRATURE_BARY:
652 cs_source_term_vcsd_q10o2_by_analytic; --> case CS_QUADRATURE_HIGHER:
653 cs_source_term_vcsd_q5o3_by_analytic; --> case CS_QUADRATURE_HIGHEST:
654 */
655
656 }
657 else {
658 assert(st_def->meta & CS_FLAG_PRIMAL);
659
660 switch (st_def->type) {
661
662 case CS_XDEF_BY_VALUE:
663 msh_flag |= CS_FLAG_COMP_PV;
664 compute_source[st_id] = cs_source_term_vcsp_by_value;
665 break;
666
667 case CS_XDEF_BY_ANALYTIC_FUNCTION:
668 msh_flag |= CS_FLAG_COMP_PV;
669 compute_source[st_id] = cs_source_term_vcsp_by_analytic;
670 break;
671
672 default:
673 bft_error(__FILE__, __LINE__, 0,
674 " Invalid type of definition for a source term in CDOVB");
675 break;
676
677 } /* switch def_type */
678
679 }
680 break; /* CDOVCB */
681
682 case CS_SPACE_SCHEME_CDOEB:
683 switch (st_def->type) {
684
685 case CS_XDEF_BY_VALUE:
686 compute_source[st_id] = cs_source_term_dfsf_by_value;
687 break;
688
689 default:
690 bft_error(__FILE__, __LINE__, 0,
691 "%s: Invalid type of definition for a source term in CDOEB",
692 __func__);
693 break;
694 }
695 break; /* CDOEB */
696
697 case CS_SPACE_SCHEME_CDOFB:
698 case CS_SPACE_SCHEME_HHO_P0:
699 switch (st_def->type) {
700
701 case CS_XDEF_BY_VALUE:
702 if ((*sys_flag) & CS_FLAG_SYS_VECTOR)
703 compute_source[st_id] = cs_source_term_pcvd_by_value;
704 else
705 compute_source[st_id] = cs_source_term_pcsd_by_value;
706 break;
707
708 case CS_XDEF_BY_DOF_FUNCTION:
709 if ((*sys_flag) & CS_FLAG_SYS_VECTOR)
710 compute_source[st_id] = cs_source_term_pcvd_by_dof_func;
711 else
712 compute_source[st_id] = cs_source_term_pcsd_by_dof_func;
713 break;
714
715 case CS_XDEF_BY_ANALYTIC_FUNCTION:
716 msh_flag |= CS_FLAG_COMP_PV;
717 if ((*sys_flag) & CS_FLAG_SYS_VECTOR) {
718
719 /* Switch only to allow more precise quadratures */
720 if (st_def->qtype == CS_QUADRATURE_BARY)
721 compute_source[st_id] = cs_source_term_pcvd_bary_by_analytic;
722
723 else {
724
725 /* TODO: Are all these flags really necessary? Check in the */
726 /* integration */
727 msh_flag |= CS_FLAG_COMP_EV |CS_FLAG_COMP_EF | CS_FLAG_COMP_PFQ |
728 CS_FLAG_COMP_FEQ |CS_FLAG_COMP_HFQ| CS_FLAG_COMP_PEQ |
729 CS_FLAG_COMP_FE;
730 compute_source[st_id] = cs_source_term_pcvd_by_analytic;
731
732 }
733 }
734 else { /* Scalar-valued case */
735
736 /* Switch only to allow more precise quadratures */
737 if (st_def->qtype == CS_QUADRATURE_BARY)
738 compute_source[st_id] = cs_source_term_pcsd_bary_by_analytic;
739
740 else {
741
742 /* TODO: Are all these flags really necessary? Check in the */
743 /* integration */
744 msh_flag |= CS_FLAG_COMP_EV |CS_FLAG_COMP_EF | CS_FLAG_COMP_PFQ |
745 CS_FLAG_COMP_FEQ |CS_FLAG_COMP_HFQ| CS_FLAG_COMP_PEQ |
746 CS_FLAG_COMP_FE;
747 compute_source[st_id] = cs_source_term_pcsd_by_analytic;
748
749 }
750
751 }
752 break;
753
754 case CS_XDEF_BY_ARRAY:
755 if ((*sys_flag) & CS_FLAG_SYS_VECTOR)
756 compute_source[st_id] = cs_source_term_pcvd_by_array;
757 else
758 compute_source[st_id] = cs_source_term_pcsd_by_array;
759 break;
760
761 default:
762 bft_error(__FILE__, __LINE__, 0,
763 "%s: Invalid type of definition for a source term in CDOFB",
764 __func__);
765 break;
766
767 } /* switch def_type */
768 break;
769
770 case CS_SPACE_SCHEME_HHO_P1:
771 case CS_SPACE_SCHEME_HHO_P2:
772 switch (st_def->type) {
773
774 case CS_XDEF_BY_VALUE:
775 if ((*sys_flag) & CS_FLAG_SYS_VECTOR)
776 bft_error(__FILE__, __LINE__, 0,
777 " %s: Invalid type of definition for a source term in HHO",
778 __func__);
779 else
780 compute_source[st_id] = cs_source_term_hhosd_by_value;
781 break;
782
783 case CS_XDEF_BY_ANALYTIC_FUNCTION:
784 if ((*sys_flag) & CS_FLAG_SYS_VECTOR)
785 compute_source[st_id] = cs_source_term_hhovd_by_analytic;
786 else
787 compute_source[st_id] = cs_source_term_hhosd_by_analytic;
788 break;
789
790 default:
791 bft_error(__FILE__, __LINE__, 0,
792 " %s: Invalid type of definition for a source term in HHO",
793 __func__);
794 break;
795
796 } /* switch def_type */
797 break;
798
799 default:
800 bft_error(__FILE__, __LINE__, 0,
801 "%s: Invalid space scheme for setting the source term.",
802 __func__);
803 break;
804
805 } /* Switch on space scheme */
806
807 } /* Loop on source terms */
808
809 if (need_mask) {
810
811 const cs_lnum_t n_cells = cs_cdo_quant->n_cells;
812
813 /* Initialize mask buffer */
814 cs_mask_t *mask = NULL;
815 BFT_MALLOC(mask, n_cells, cs_mask_t);
816 # pragma omp parallel for if (n_cells > CS_THR_MIN)
817 for (int i = 0; i < n_cells; i++) mask[i] = 0;
818
819 for (int st_id = 0; st_id < n_source_terms; st_id++)
820 _set_mask(source_terms[st_id], st_id, mask);
821
822 *source_mask = mask;
823
824 } /* Build a tag related to the source terms defined in each cell */
825
826 return msh_flag;
827 }
828
829 /*----------------------------------------------------------------------------*/
830 /*!
831 * \brief Compute the local contributions of source terms in a cell
832 *
833 * \param[in] n_source_terms number of source terms
834 * \param[in] source_terms pointer to the definitions of source terms
835 * \param[in] cm pointer to a cs_cell_mesh_t structure
836 * \param[in] source_mask array storing in a compact way which source
837 * term is defined in a given cell
838 * \param[in] compute_source array of function pointers
839 * \param[in] time_eval physical time at which one evaluates the term
840 * \param[in, out] input pointer to an element cast on-the-fly
841 * \param[in, out] cb pointer to a cs_cell_builder_t structure
842 * \param[in, out] result array storing the result of the evaluation
843 */
844 /*----------------------------------------------------------------------------*/
845
846 void
cs_source_term_compute_cellwise(const int n_source_terms,cs_xdef_t * const * source_terms,const cs_cell_mesh_t * cm,const cs_mask_t * source_mask,cs_source_term_cellwise_t * compute_source[],cs_real_t time_eval,void * input,cs_cell_builder_t * cb,cs_real_t * result)847 cs_source_term_compute_cellwise(const int n_source_terms,
848 cs_xdef_t *const *source_terms,
849 const cs_cell_mesh_t *cm,
850 const cs_mask_t *source_mask,
851 cs_source_term_cellwise_t *compute_source[],
852 cs_real_t time_eval,
853 void *input,
854 cs_cell_builder_t *cb,
855 cs_real_t *result)
856 {
857 if (n_source_terms < 1)
858 return;
859
860 if (source_mask == NULL) { /* Source terms are defined on the whole mesh */
861
862 for (short int st_id = 0; st_id < n_source_terms; st_id++) {
863
864 cs_source_term_cellwise_t *compute = compute_source[st_id];
865
866 /* result is updated inside */
867 compute(source_terms[st_id], cm, time_eval, cb, input, result);
868
869 } /* Loop on source terms */
870
871 }
872 else { /* Some source terms are only defined on a selection of cells */
873
874 for (short int st_id = 0; st_id < n_source_terms; st_id++) {
875
876 const cs_mask_t st_mask = (1 << st_id);
877 if (source_mask[cm->c_id] & st_mask) {
878
879 cs_source_term_cellwise_t *compute = compute_source[st_id];
880
881 /* result is updated inside */
882 compute(source_terms[st_id], cm, time_eval, cb, input, result);
883
884 } /* Compute the source term on this cell */
885
886 } /* Loop on source terms */
887
888 } /* Source terms are defined on the whole domain or not ? */
889
890 }
891
892 /*----------------------------------------------------------------------------*/
893 /*!
894 * \brief Compute the contribution for a cell related to a source term and
895 * add it to the given array of values.
896 * Case of a scalar potential defined at primal vertices by a constant
897 * value.
898 * A discrete Hodge operator has to be computed before this call and
899 * given as input parameter
900 *
901 * \param[in] source pointer to a cs_xdef_t structure
902 * \param[in] cm pointer to a cs_cell_mesh_t structure
903 * \param[in] time_eval physical time at which one evaluates the term
904 * \param[in, out] cb pointer to a cs_cell_builder_t structure
905 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
906 * \param[in, out] values pointer to the computed values
907 */
908 /*----------------------------------------------------------------------------*/
909
910 void
cs_source_term_pvsp_by_value(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)911 cs_source_term_pvsp_by_value(const cs_xdef_t *source,
912 const cs_cell_mesh_t *cm,
913 cs_real_t time_eval,
914 cs_cell_builder_t *cb,
915 void *input,
916 double *values)
917 {
918 CS_UNUSED(time_eval);
919
920 if (source == NULL)
921 return;
922
923 cs_hodge_t *mass_hodge = (cs_hodge_t *)input;
924
925 /* Sanity checks */
926 assert(values != NULL && cm != NULL && cb != NULL);
927 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PV));
928 assert(mass_hodge != NULL);
929 assert(mass_hodge->matrix != NULL);
930
931 const cs_real_t *s_values = (const cs_real_t *)source->context;
932 const cs_real_t pot_value = s_values[0];
933
934 /* Retrieve the values of the potential at each cell vertices */
935 double *eval = cb->values;
936 for (short int v = 0; v < cm->n_vc; v++)
937 eval[v] = pot_value;
938
939 /* Multiply these values by a cellwise Hodge operator previously computed */
940 double *hdg_eval = cb->values + cm->n_vc;
941 cs_sdm_square_matvec(mass_hodge->matrix, eval, hdg_eval);
942
943 for (short int v = 0; v < cm->n_vc; v++)
944 values[v] += hdg_eval[v];
945 }
946
947 /*----------------------------------------------------------------------------*/
948 /*!
949 * \brief Compute the contribution for a cell related to a source term and
950 * add it to the given array of values.
951 * Case of a scalar potential defined at primal vertices by an
952 * analytical function.
953 * A discrete Hodge operator has to be computed before this call and
954 * given as input parameter
955 *
956 * \param[in] source pointer to a cs_xdef_t structure
957 * \param[in] cm pointer to a cs_cell_mesh_t structure
958 * \param[in] time_eval physical time at which one evaluates the term
959 * \param[in, out] cb pointer to a cs_cell_builder_t structure
960 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
961 * \param[in, out] values pointer to the computed values
962 */
963 /*----------------------------------------------------------------------------*/
964
965 void
cs_source_term_pvsp_by_analytic(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)966 cs_source_term_pvsp_by_analytic(const cs_xdef_t *source,
967 const cs_cell_mesh_t *cm,
968 cs_real_t time_eval,
969 cs_cell_builder_t *cb,
970 void *input,
971 double *values)
972 {
973 if (source == NULL)
974 return;
975
976 cs_hodge_t *mass_hodge = (cs_hodge_t *)input;
977
978 /* Sanity checks */
979 assert(values != NULL && cm != NULL && cb != NULL);
980 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PV));
981 assert(mass_hodge != NULL);
982 assert(mass_hodge->matrix != NULL);
983
984 cs_xdef_analytic_context_t *ac =
985 (cs_xdef_analytic_context_t *)source->context;
986
987 /* Retrieve the values of the potential at each cell vertices */
988 double *eval = cb->values;
989 ac->func(time_eval, cm->n_vc, NULL, cm->xv,
990 true, /* compacted output ? */
991 ac->input,
992 eval);
993
994 /* Multiply these values by a cellwise Hodge operator previously computed */
995 double *hdg_eval = cb->values + cm->n_vc;
996 cs_sdm_square_matvec(mass_hodge->matrix, eval, hdg_eval);
997
998 for (short int v = 0; v < cm->n_vc; v++)
999 values[v] += hdg_eval[v];
1000 }
1001
1002 /*----------------------------------------------------------------------------*/
1003 /*!
1004 * \brief Compute the contribution for a cell related to a source term and
1005 * add it to the given array of values.
1006 * Case of a scalar density defined at dual cells by a value.
1007 *
1008 * \param[in] source pointer to a cs_xdef_t structure
1009 * \param[in] cm pointer to a cs_cell_mesh_t structure
1010 * \param[in] time_eval physical time at which one evaluates the term
1011 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1012 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1013 * \param[in, out] values pointer to the computed values
1014 */
1015 /*----------------------------------------------------------------------------*/
1016
1017 void
cs_source_term_dcsd_by_value(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1018 cs_source_term_dcsd_by_value(const cs_xdef_t *source,
1019 const cs_cell_mesh_t *cm,
1020 cs_real_t time_eval,
1021 cs_cell_builder_t *cb,
1022 void *input,
1023 double *values)
1024 {
1025 CS_UNUSED(cb);
1026 CS_UNUSED(input);
1027 CS_UNUSED(time_eval);
1028
1029 if (source == NULL)
1030 return;
1031
1032 /* Sanity checks */
1033 assert(values != NULL && cm != NULL);
1034 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PVQ));
1035
1036 const cs_real_t *s_values = (const cs_real_t *)source->context;
1037 const cs_real_t density_value = s_values[0];
1038
1039 for (int v = 0; v < cm->n_vc; v++)
1040 values[v] += density_value * cm->wvc[v] * cm->vol_c;
1041 }
1042
1043 /*----------------------------------------------------------------------------*/
1044 /*!
1045 * \brief Compute the contribution for a cell related to a source term and
1046 * add it to the given array of values.
1047 * Case of a vector-valued density defined at dual cells by a value.
1048 *
1049 * \param[in] source pointer to a cs_xdef_t structure
1050 * \param[in] cm pointer to a cs_cell_mesh_t structure
1051 * \param[in] time_eval physical time at which one evaluates the term
1052 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1053 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1054 * \param[in, out] values pointer to the computed values
1055 */
1056 /*----------------------------------------------------------------------------*/
1057
1058 void
cs_source_term_dcvd_by_value(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1059 cs_source_term_dcvd_by_value(const cs_xdef_t *source,
1060 const cs_cell_mesh_t *cm,
1061 cs_real_t time_eval,
1062 cs_cell_builder_t *cb,
1063 void *input,
1064 double *values)
1065 {
1066 CS_UNUSED(cb);
1067 CS_UNUSED(input);
1068 CS_UNUSED(time_eval);
1069
1070 if (source == NULL)
1071 return;
1072
1073 /* Sanity checks */
1074 assert(values != NULL && cm != NULL);
1075 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PVQ));
1076
1077 const cs_real_t *st_vect = (const cs_real_t *)source->context;
1078
1079 for (int v = 0; v < cm->n_vc; v++)
1080 for (int k = 0; k < 3; k++)
1081 values[3*v+k] += st_vect[k] * cm->wvc[v] * cm->vol_c;
1082 }
1083
1084 /*----------------------------------------------------------------------------*/
1085 /*!
1086 * \brief Compute the contribution for a cell related to a source term and
1087 * add it to the given array of values.
1088 * Case of a scalar density defined at dual cells by an array.
1089 *
1090 * \param[in] source pointer to a cs_xdef_t structure
1091 * \param[in] cm pointer to a cs_cell_mesh_t structure
1092 * \param[in] time_eval physical time at which one evaluates the term
1093 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1094 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1095 * \param[in, out] values pointer to the computed values
1096 */
1097 /*----------------------------------------------------------------------------*/
1098
1099 void
cs_source_term_dcsd_by_array(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1100 cs_source_term_dcsd_by_array(const cs_xdef_t *source,
1101 const cs_cell_mesh_t *cm,
1102 cs_real_t time_eval,
1103 cs_cell_builder_t *cb,
1104 void *input,
1105 double *values)
1106 {
1107 CS_UNUSED(cb);
1108 CS_UNUSED(input);
1109 CS_UNUSED(time_eval);
1110
1111 if (source == NULL)
1112 return;
1113
1114 /* Sanity checks */
1115 assert(values != NULL && cm != NULL);
1116 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PVQ));
1117
1118 const cs_xdef_array_context_t *ac =
1119 (const cs_xdef_array_context_t *)source->context;
1120
1121 assert(cs_flag_test(ac->loc, cs_flag_primal_vtx));
1122 for (int v = 0; v < cm->n_vc; v++)
1123 values[v] += ac->values[cm->v_ids[v]] * cm->wvc[v] * cm->vol_c;
1124 }
1125
1126 /*----------------------------------------------------------------------------*/
1127 /*!
1128 * \brief Compute the contribution for a cell related to a source term and
1129 * add it to the given array of values.
1130 * Case of a scalar density defined at dual cells by a DoF function.
1131 *
1132 * \param[in] source pointer to a cs_xdef_t structure
1133 * \param[in] cm pointer to a cs_cell_mesh_t structure
1134 * \param[in] time_eval physical time at which one evaluates the term
1135 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1136 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1137 * \param[in, out] values pointer to the computed values
1138 */
1139 /*----------------------------------------------------------------------------*/
1140
1141 void
cs_source_term_dcsd_by_dof_func(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1142 cs_source_term_dcsd_by_dof_func(const cs_xdef_t *source,
1143 const cs_cell_mesh_t *cm,
1144 cs_real_t time_eval,
1145 cs_cell_builder_t *cb,
1146 void *input,
1147 double *values)
1148 {
1149 CS_UNUSED(cb);
1150 CS_UNUSED(time_eval);
1151 CS_UNUSED(input);
1152
1153 if (source == NULL)
1154 return;
1155
1156 /* Sanity checks */
1157 assert(values != NULL && cm != NULL);
1158 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PVQ));
1159
1160 const cs_xdef_dof_context_t *dc = (cs_xdef_dof_context_t *)source->context;
1161
1162 /* Up to now this should be the only location allowed */
1163 assert(cs_flag_test(dc->loc, cs_flag_primal_cell));
1164
1165 /* Call the DoF function to evaluate the function at xc */
1166 double cell_eval;
1167 dc->func(1, &(cm->c_id), true, /* compacted output ? */
1168 dc->input,
1169 &cell_eval);
1170
1171 cell_eval *= cm->vol_c;
1172 for (int v = 0; v < cm->n_vc; v++)
1173 values[v] += cell_eval * cm->wvc[v];
1174 }
1175
1176 /*----------------------------------------------------------------------------*/
1177 /*!
1178 * \brief Compute the contribution for a cell related to a source term and
1179 * add it to the given array of values.
1180 * Case of a scalar density defined at dual cells by an analytical
1181 * function.
1182 * Use the barycentric approximation as quadrature to evaluate the
1183 * integral. Exact for linear function.
1184 *
1185 * \param[in] source pointer to a cs_xdef_t structure
1186 * \param[in] cm pointer to a cs_cell_mesh_t structure
1187 * \param[in] time_eval physical time at which one evaluates the term
1188 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1189 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1190 * \param[in, out] values pointer to the computed values
1191 */
1192 /*----------------------------------------------------------------------------*/
1193
1194 void
cs_source_term_dcsd_bary_by_analytic(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1195 cs_source_term_dcsd_bary_by_analytic(const cs_xdef_t *source,
1196 const cs_cell_mesh_t *cm,
1197 cs_real_t time_eval,
1198 cs_cell_builder_t *cb,
1199 void *input,
1200 double *values)
1201 {
1202 CS_UNUSED(input);
1203
1204 if (source == NULL)
1205 return;
1206
1207 /* Sanity checks */
1208 assert(values != NULL && cm != NULL);
1209 assert(cs_eflag_test(cm->flag,
1210 CS_FLAG_COMP_PVQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_PFC |
1211 CS_FLAG_COMP_FE | CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV |
1212 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PEC));
1213
1214 cs_xdef_analytic_context_t *ac =
1215 (cs_xdef_analytic_context_t *)source->context;
1216
1217 /* Compute the barycenter of each portion of dual cells */
1218 double *vol_vc = cb->values;
1219 for (short int v = 0; v < cm->n_vc; v++)
1220 vol_vc[v] = cm->vol_c * cm->wvc[v];
1221
1222 /* cell and vertex contribution */
1223 cs_real_3_t *xgv = cb->vectors;
1224 for (short int v = 0; v < cm->n_vc; v++) {
1225 xgv[v][0] = 0.25 * vol_vc[v] * (cm->xc[0] + cm->xv[0]);
1226 xgv[v][1] = 0.25 * vol_vc[v] * (cm->xc[1] + cm->xv[1]);
1227 xgv[v][2] = 0.25 * vol_vc[v] * (cm->xc[2] + cm->xv[2]);
1228 }
1229
1230 /* edge contribution */
1231 for (short int e = 0; e < cm->n_ec; e++) {
1232
1233 cs_real_t *xgv1 = xgv[cm->e2v_ids[2*e]];
1234 cs_real_t *xgv2 = xgv[cm->e2v_ids[2*e+1]];
1235
1236 const cs_real_t *xe = cm->edge[e].center;
1237 const double e_coef = 0.125 * cm->pvol_e[e]; /* 0.25* (0.5*|pvol_ec|) */
1238 for (int k = 0; k < 3; k++) {
1239 xgv1[k] += e_coef * xe[k];
1240 xgv2[k] += e_coef * xe[k];
1241 }
1242
1243 } /* Loop on cell edges */
1244
1245 /* face contribution */
1246 cs_real_t *wvf = cb->values + cm->n_vc;
1247 for (short int f = 0; f < cm->n_fc; f++) {
1248
1249 cs_compute_wvf(f, cm, wvf);
1250
1251 const double *xf = cm->face[f].center;
1252 const double f_coef = 0.25 * cm->pvol_f[f];
1253
1254 for (short int v = 0; v < cm->n_vc; v++) {
1255 if (wvf[v] > 0) {
1256 const double vf_coef = f_coef * wvf[v];
1257 xgv[v][0] += vf_coef * xf[0];
1258 xgv[v][1] += vf_coef * xf[1];
1259 xgv[v][2] += vf_coef * xf[2];
1260 }
1261
1262 }
1263
1264 } /* Loop on faces */
1265
1266 for (short int v = 0; v < cm->n_vc; v++) {
1267 const double invvol = 1/vol_vc[v];
1268 for (int k = 0; k < 3; k++) xgv[v][k] *= invvol;
1269 }
1270
1271 /* Call the analytic function to evaluate the function at xgv */
1272 double *eval_xgv = vol_vc + cm->n_vc;
1273 ac->func(time_eval, cm->n_vc, NULL, (const cs_real_t *)xgv,
1274 true, /* compacted output ? */
1275 ac->input,
1276 eval_xgv);
1277
1278 for (short int v = 0; v < cm->n_vc; v++)
1279 values[v] = vol_vc[v] * eval_xgv[v];
1280 }
1281
1282 /*----------------------------------------------------------------------------*/
1283 /*!
1284 * \brief Compute the contribution for a cell related to a source term and
1285 * add it to the given array of values.
1286 * Case of a scalar density defined at dual cells by an analytical
1287 * function.
1288 * Use a the barycentric approximation as quadrature to evaluate the
1289 * integral. Exact for linear function.
1290 *
1291 * \param[in] source pointer to a cs_xdef_t structure
1292 * \param[in] cm pointer to a cs_cell_mesh_t structure
1293 * \param[in] time_eval physical time at which one evaluates the term
1294 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1295 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1296 * \param[in, out] values pointer to the computed values
1297 */
1298 /*----------------------------------------------------------------------------*/
1299
1300 void
cs_source_term_dcsd_q1o1_by_analytic(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1301 cs_source_term_dcsd_q1o1_by_analytic(const cs_xdef_t *source,
1302 const cs_cell_mesh_t *cm,
1303 cs_real_t time_eval,
1304 cs_cell_builder_t *cb,
1305 void *input,
1306 double *values)
1307 {
1308 CS_UNUSED(cb);
1309 CS_UNUSED(input);
1310
1311 if (source == NULL)
1312 return;
1313
1314 /* Sanity checks */
1315 assert(values != NULL && cm != NULL);
1316 assert(cs_eflag_test(cm->flag,
1317 CS_FLAG_COMP_PFQ | CS_FLAG_COMP_PFC | CS_FLAG_COMP_FE |
1318 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
1319
1320 cs_xdef_analytic_context_t *ac =
1321 (cs_xdef_analytic_context_t *)source->context;
1322
1323 for (short int f = 0; f < cm->n_fc; f++) {
1324
1325 cs_real_3_t xg[2], xfc;
1326 cs_real_t eval_xg[2];
1327
1328 const double *xf = cm->face[f].center;
1329 const double hf_coef = 0.5 * cm->pvol_f[f]/cm->face[f].meas;
1330
1331 for (int k = 0; k < 3; k++) xfc[k] = 0.25*(xf[k] + cm->xc[k]);
1332
1333 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1334
1335 const short int e = cm->f2e_ids[i];
1336 const short int v1 = cm->e2v_ids[2*e];
1337 const short int v2 = cm->e2v_ids[2*e+1];
1338 const double *xv1 = cm->xv + 3*v1, *xv2 = cm->xv + 3*v2;
1339
1340 /* xg = 0.25(xv1 + xe + xf + xc) where xe = 0.5*(xv1 + xv2) */
1341 for (int k = 0; k < 3; k++)
1342 xg[0][k] = xfc[k] + 0.375*xv1[k] + 0.125*xv2[k];
1343
1344 /* xg = 0.25(xv2 + xe + xf + xc) where xe = 0.5*(xv1 + xv2) */
1345 for (int k = 0; k < 3; k++)
1346 xg[1][k] = xfc[k] + 0.375*xv2[k] + 0.125*xv1[k];
1347
1348 ac->func(time_eval, 2, NULL, (const cs_real_t *)xg,
1349 true, /* compacted output ? */
1350 ac->input,
1351 eval_xg);
1352
1353 const double half_pef_vol = cm->tef[i]*hf_coef;
1354 values[v1] += half_pef_vol * eval_xg[0];
1355 values[v2] += half_pef_vol * eval_xg[1];
1356
1357 } /* Loop on face edges */
1358
1359 } /* Loop on cell faces */
1360
1361 }
1362
1363 /*----------------------------------------------------------------------------*/
1364 /*!
1365 * \brief Compute the contribution for a cell related to a source term and
1366 * add it to the given array of values.
1367 * Case of a scalar density defined at dual cells by an analytical
1368 * function.
1369 * Use a ten-point quadrature rule to evaluate the integral.
1370 * Exact for quadratic function.
1371 *
1372 * \param[in] source pointer to a cs_xdef_t structure
1373 * \param[in] cm pointer to a cs_cell_mesh_t structure
1374 * \param[in] time_eval physical time at which one evaluates the term
1375 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1376 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1377 * \param[in, out] values pointer to the computed values
1378 */
1379 /*----------------------------------------------------------------------------*/
1380
1381 void
cs_source_term_dcsd_q10o2_by_analytic(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1382 cs_source_term_dcsd_q10o2_by_analytic(const cs_xdef_t *source,
1383 const cs_cell_mesh_t *cm,
1384 cs_real_t time_eval,
1385 cs_cell_builder_t *cb,
1386 void *input,
1387 double *values)
1388 {
1389 CS_UNUSED(input);
1390
1391 if (source == NULL)
1392 return;
1393
1394 /* Sanity checks */
1395 assert(values != NULL && cm != NULL);
1396 assert(cb != NULL);
1397 assert(cs_eflag_test(cm->flag,
1398 CS_FLAG_COMP_PFQ | CS_FLAG_COMP_PFC | CS_FLAG_COMP_FE |
1399 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV | CS_FLAG_COMP_PVQ |
1400 CS_FLAG_COMP_PEQ));
1401
1402 cs_xdef_analytic_context_t *ac =
1403 (cs_xdef_analytic_context_t *)source->context;
1404
1405 /* Temporary buffers */
1406 double *contrib = cb->values; /* size n_vc */
1407
1408 /* 1) Compute the contributions seen by the whole portion of dual cell */
1409
1410 /* Cell evaluation */
1411 double eval_c;
1412 ac->func(time_eval, 1, NULL, cm->xc,
1413 true, /* compacted output ? */
1414 ac->input,
1415 &eval_c);
1416
1417 /* Contributions related to vertices */
1418 double *eval_v = cb->values + cm->n_vc; /* size n_vc */
1419 ac->func(time_eval, cm->n_vc, NULL, cm->xv,
1420 true, /* compacted output ? */
1421 ac->input,
1422 eval_v);
1423
1424 cs_real_3_t *xvc = cb->vectors;
1425 for (short int v = 0; v < cm->n_vc; v++) {
1426 const double *xv = cm->xv + 3*v;
1427 for (int k = 0; k < 3; k++) xvc[v][k] = 0.5*(cm->xc[k] + xv[k]);
1428 }
1429
1430 double *eval_vc = cb->values + 2*cm->n_vc; /* size n_vc */
1431 ac->func(time_eval, cm->n_vc, NULL, (const cs_real_t *)xvc,
1432 true, /* compacted output ? */
1433 ac->input,
1434 eval_vc);
1435
1436 for (short int v = 0; v < cm->n_vc; v++) {
1437
1438 /* Set the initial values
1439 -1/20 on extremity points and 1/5 on midpoints */
1440 contrib[v] = cm->wvc[v]*cm->vol_c
1441 * (-0.05*(eval_c + eval_v[v]) + 0.2*eval_vc[v]);
1442
1443 } /* Loop on vertices */
1444
1445 /* 2) Compute the contribution related to edge
1446 The portion of dual cell seen by each vertex is 1/2 |pec| */
1447 cs_real_3_t *x_e = cb->vectors;
1448 cs_real_3_t *xec = cb->vectors + cm->n_ec; /* size = n_ec (overwrite xvc) */
1449
1450 for (short int e = 0; e < cm->n_ec; e++) {
1451 for (int k = 0; k < 3; k++) {
1452 x_e[e][k] = cm->edge[e].center[k];
1453 xec[e][k] = 0.5*(cm->xc[k] + x_e[e][k]);
1454 }
1455 }
1456
1457 /* Evaluate the analytic function at xe and xec */
1458 double *eval_e = cb->values + cm->n_vc; /* size=n_ec (overwrite eval_v) */
1459 double *eval_ec = eval_e + cm->n_ec; /* size=n_ec (overwrite eval_vc) */
1460 ac->func(time_eval, 2*cm->n_ec, NULL, (const cs_real_t *)cb->vectors,
1461 true, /* compacted output ? */
1462 ac->input,
1463 eval_e);
1464
1465 /* xev (size = 2*n_ec) */
1466 cs_real_3_t *xve = cb->vectors; /* size=2*n_ec (overwrite xe and xec) */
1467 for (short int e = 0; e < cm->n_ec; e++) {
1468
1469 const cs_real_t *xe = cm->edge[e].center;
1470 const short int v1 = cm->e2v_ids[2*e];
1471 const double *xv1 = cm->xv + 3*v1;
1472 const short int v2 = cm->e2v_ids[2*e+1];
1473 const double *xv2 = cm->xv + 3*v2;
1474
1475 for (int k = 0; k < 3; k++) {
1476 xve[2*e ][k] = 0.5*(xv1[k] + xe[k]);
1477 xve[2*e+1][k] = 0.5*(xv2[k] + xe[k]);
1478 }
1479
1480 } /* Loop on edges */
1481
1482 double *eval_ve = eval_ec + cm->n_ec; /* size = 2*n_ec */
1483 ac->func(time_eval, 2*cm->n_ec, NULL, (const cs_real_t *)cb->vectors,
1484 true, /* compacted output ? */
1485 ac->input,
1486 eval_ve);
1487
1488 /* 3) Main loop on faces */
1489 double *pvf_vol = eval_ve + 2*cm->n_ec; /* size n_vc */
1490
1491 for (short int f = 0; f < cm->n_fc; f++) {
1492
1493 const double *xf = cm->face[f].center;
1494 const double hf_coef = 0.5* cm->pvol_f[f]/cm->face[f].meas;
1495
1496 /* Reset volume of the face related to a vertex */
1497 for (short int v = 0; v < cm->n_vc; v++) pvf_vol[v] = 0;
1498
1499 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1500
1501 const short int e = cm->f2e_ids[i];
1502 const short int v1 = cm->e2v_ids[2*e];
1503 const short int v2 = cm->e2v_ids[2*e+1];
1504 const double half_pef_vol = hf_coef * cm->tef[i];
1505
1506 pvf_vol[v1] += half_pef_vol;
1507 pvf_vol[v2] += half_pef_vol;
1508
1509 cs_real_3_t xef;
1510 cs_real_t eval_ef;
1511 for (int k = 0; k < 3; k++) xef[k] = 0.5*(cm->edge[e].center[k] + xf[k]);
1512 ac->func(time_eval, 1, NULL, xef,
1513 true, /* compacted output ? */
1514 ac->input,
1515 &eval_ef);
1516
1517 /* 1/5 (EF + EC) -1/20 * (E) */
1518 const double common_ef_contrib =
1519 0.2*(eval_ef + eval_ec[e]) -0.05*eval_e[e];
1520
1521 contrib[v1] += half_pef_vol*(common_ef_contrib + 0.2*eval_ve[2*e]);
1522 contrib[v2] += half_pef_vol*(common_ef_contrib + 0.2*eval_ve[2*e+1]);
1523
1524 } /* Loop on face edges */
1525
1526 /* Contributions related to this face */
1527 cs_real_3_t *xvfc = cb->vectors; /* size=2+n_vc (overwrite xev) */
1528 for (int k = 0; k < 3; k++) {
1529 xvfc[0][k] = xf[k]; /* xf */
1530 xvfc[1][k] = 0.5*(xf[k] + cm->xc[k]); /* xfc */
1531 }
1532
1533 short int n_vf = 0;
1534 for (short int v = 0; v < cm->n_vc; v++) {
1535 if (pvf_vol[v] > 0) {
1536 cb->ids[n_vf] = v;
1537 for (int k = 0; k < 3; k++)
1538 xvfc[2+n_vf][k] = 0.5*(xf[k] + cm->xv[3*v+k]);
1539 n_vf++;
1540 }
1541 }
1542
1543 double *eval_vfc = pvf_vol + cm->n_vc; /* size=n_vf + 2 */
1544 ac->func(time_eval, 2+n_vf, NULL, (const cs_real_t *)xvfc,
1545 true, /* compacted output ? */
1546 ac->input,
1547 eval_vfc);
1548
1549 for (short int i = 0; i < n_vf; i++) {
1550 short int v = cb->ids[i];
1551 const double val_vfc = -0.05*eval_vfc[0] + 0.2*eval_vfc[1];
1552 contrib[v] += pvf_vol[v] * (val_vfc + 0.2*eval_vfc[2+i]);
1553 }
1554
1555 } /* Loop on cell faces */
1556
1557 /* Add the computed contributions to the return values */
1558 for (short int v = 0; v < cm->n_vc; v++)
1559 values[v] += contrib[v];
1560
1561 }
1562
1563 /*----------------------------------------------------------------------------*/
1564 /*!
1565 * \brief Compute the contribution for a cell related to a source term and
1566 * add it to the given array of values.
1567 * Case of a scalar density defined at dual cells by an analytical
1568 * function.
1569 * Use a five-point quadrature rule to evaluate the integral.
1570 * Exact for cubic function.
1571 * This function may be expensive since many evaluations are needed.
1572 * Please use it with care.
1573 *
1574 * \param[in] source pointer to a cs_xdef_t structure
1575 * \param[in] cm pointer to a cs_cell_mesh_t structure
1576 * \param[in] time_eval physical time at which one evaluates the term
1577 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1578 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1579 * \param[in, out] values pointer to the computed values
1580 */
1581 /*----------------------------------------------------------------------------*/
1582
1583 void
cs_source_term_dcsd_q5o3_by_analytic(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1584 cs_source_term_dcsd_q5o3_by_analytic(const cs_xdef_t *source,
1585 const cs_cell_mesh_t *cm,
1586 cs_real_t time_eval,
1587 cs_cell_builder_t *cb,
1588 void *input,
1589 double *values)
1590 {
1591 CS_UNUSED(input);
1592
1593 if (source == NULL)
1594 return;
1595
1596 double sum, weights[5], results[5];
1597 cs_real_3_t gauss_pts[5];
1598
1599 /* Sanity checks */
1600 assert(values != NULL && cm != NULL);
1601 assert(cb != NULL);
1602 assert(cs_eflag_test(cm->flag,
1603 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
1604 CS_FLAG_COMP_EV | CS_FLAG_COMP_PFC | CS_FLAG_COMP_FEQ));
1605
1606 cs_xdef_analytic_context_t *ac =
1607 (cs_xdef_analytic_context_t *)source->context;
1608
1609 /* Temporary buffers */
1610 double *contrib = cb->values;
1611 memset(contrib, 0, cm->n_vc*sizeof(double));
1612
1613 /* Main loop on faces */
1614 for (short int f = 0; f < cm->n_fc; f++) {
1615
1616 const double *xf = cm->face[f].center;
1617 const double hf_coef = 0.5* cm->pvol_f[f]/cm->face[f].meas;
1618
1619 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1620
1621 const short int e = cm->f2e_ids[i];
1622 const short int v1 = cm->e2v_ids[2*e];
1623 const short int v2 = cm->e2v_ids[2*e+1];
1624 const double half_pef_vol = hf_coef * cm->tef[i];
1625
1626 /* Compute Gauss points and its weights */
1627 cs_quadrature_tet_5pts(cm->xv + 3*v1, cm->edge[e].center, xf, cm->xc,
1628 half_pef_vol,
1629 gauss_pts, weights);
1630
1631 ac->func(time_eval, 5, NULL, (const cs_real_t *)gauss_pts,
1632 true, /* compacted output ? */
1633 ac->input,
1634 results);
1635
1636 sum = 0.;
1637 for (int p = 0; p < 5; p++) sum += results[p] * weights[p];
1638 contrib[v1] += sum;
1639
1640 /* Compute Gauss points and its weights */
1641 cs_quadrature_tet_5pts(cm->xv + 3*v2, cm->edge[e].center, xf, cm->xc,
1642 half_pef_vol,
1643 gauss_pts, weights);
1644
1645 ac->func(time_eval, 5, NULL, (const cs_real_t *)gauss_pts,
1646 true, /* compacted output ? */
1647 ac->input,
1648 results);
1649
1650 sum = 0.;
1651 for (int p = 0; p < 5; p++) sum += results[p] * weights[p];
1652 contrib[v2] += sum;
1653
1654 } /* Loop on face edges */
1655
1656 } /* Loop on cell faces */
1657
1658 /* Add the computed contributions to the return values */
1659 for (short int v = 0; v < cm->n_vc; v++)
1660 values[v] += contrib[v];
1661 }
1662
1663 /*----------------------------------------------------------------------------*/
1664 /*!
1665 * \brief Compute the contribution for a cell related to a source term and
1666 * add it to the given array of values.
1667 * Case of a scalar potential defined at primal vertices and cells
1668 * by a constant value.
1669 * A discrete Hodge operator has to be computed before this call and
1670 * given as input parameter
1671 *
1672 * \param[in] source pointer to a cs_xdef_t structure
1673 * \param[in] cm pointer to a cs_cell_mesh_t structure
1674 * \param[in] time_eval physical time at which one evaluates the term
1675 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1676 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1677 * \param[in, out] values pointer to the computed values
1678 */
1679 /*----------------------------------------------------------------------------*/
1680
1681 void
cs_source_term_vcsp_by_value(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1682 cs_source_term_vcsp_by_value(const cs_xdef_t *source,
1683 const cs_cell_mesh_t *cm,
1684 cs_real_t time_eval,
1685 cs_cell_builder_t *cb,
1686 void *input,
1687 double *values)
1688 {
1689 CS_UNUSED(time_eval);
1690
1691 if (source == NULL)
1692 return;
1693
1694 cs_hodge_t *mass_hodge = (cs_hodge_t *)input;
1695
1696 /* Sanity checks */
1697 assert(values != NULL && cm != NULL && cb != NULL);
1698 assert(mass_hodge != NULL);
1699 assert(mass_hodge->matrix != NULL);
1700
1701 const cs_real_t *s_values = (const cs_real_t *)source->context;
1702 const cs_real_t pot_value = s_values[0];
1703
1704 /* Retrieve the values of the potential at each cell vertices */
1705 double *eval = cb->values;
1706 for (short int v = 0; v < cm->n_vc; v++)
1707 eval[v] = pot_value;
1708 eval[cm->n_vc] = pot_value;
1709
1710 /* Multiply these values by a cellwise Hodge operator previously computed */
1711 double *hdg_eval = cb->values + cm->n_vc + 1;
1712 cs_sdm_square_matvec(mass_hodge->matrix, eval, hdg_eval);
1713
1714 for (short int v = 0; v < cm->n_vc + 1; v++)
1715 values[v] += hdg_eval[v];
1716 }
1717
1718 /*----------------------------------------------------------------------------*/
1719 /*!
1720 * \brief Compute the contribution for a cell related to a source term and
1721 * add it to the given array of values.
1722 * Case of a scalar potential defined at primal vertices and cells by
1723 * an analytical function.
1724 * A discrete Hodge operator has to be computed before this call and
1725 * given as input parameter
1726 *
1727 * \param[in] source pointer to a cs_xdef_t structure
1728 * \param[in] cm pointer to a cs_cell_mesh_t structure
1729 * \param[in] time_eval physical time at which one evaluates the term
1730 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1731 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1732 * \param[in, out] values pointer to the computed values
1733 */
1734 /*----------------------------------------------------------------------------*/
1735
1736 void
cs_source_term_vcsp_by_analytic(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1737 cs_source_term_vcsp_by_analytic(const cs_xdef_t *source,
1738 const cs_cell_mesh_t *cm,
1739 cs_real_t time_eval,
1740 cs_cell_builder_t *cb,
1741 void *input,
1742 double *values)
1743 {
1744 if (source == NULL)
1745 return;
1746
1747 cs_hodge_t *mass_hodge = (cs_hodge_t *)input;
1748
1749 /* Sanity checks */
1750 assert(values != NULL && cm != NULL && cb != NULL);
1751 assert(mass_hodge != NULL);
1752 assert(mass_hodge->matrix != NULL);
1753
1754 cs_xdef_analytic_context_t *ac =
1755 (cs_xdef_analytic_context_t *)source->context;
1756
1757 /* Retrieve the values of the potential at each cell vertices */
1758 double *eval = cb->values;
1759
1760 ac->func(time_eval, cm->n_vc, NULL, cm->xv,
1761 true, /* compacted output ? */
1762 ac->input,
1763 eval);
1764
1765 /* Retrieve the value at the cell center */
1766 ac->func(time_eval, 1, NULL, cm->xc,
1767 true, /* compacted output ? */
1768 ac->input,
1769 eval + cm->n_vc);
1770
1771 /* Multiply these values by a cellwise Hodge operator previously computed */
1772 double *hdg_eval = cb->values + cm->n_vc + 1;
1773 cs_sdm_square_matvec(mass_hodge->matrix, eval, hdg_eval);
1774
1775 for (short int v = 0; v < cm->n_vc + 1; v++)
1776 values[v] += hdg_eval[v];
1777 }
1778
1779 /*----------------------------------------------------------------------------*/
1780 /*!
1781 * \brief Compute the contribution for a cell related to a source term and
1782 * add it to the given array of values.
1783 * Case of a scalar density (sd) defined on primal cells by a value.
1784 * Case of face-based schemes
1785 *
1786 * \param[in] source pointer to a cs_xdef_t structure
1787 * \param[in] cm pointer to a cs_cell_mesh_t structure
1788 * \param[in] time_eval physical time at which one evaluates the term
1789 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1790 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1791 * \param[in, out] values pointer to the computed value
1792 */
1793 /*----------------------------------------------------------------------------*/
1794
1795 void
cs_source_term_pcsd_by_value(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1796 cs_source_term_pcsd_by_value(const cs_xdef_t *source,
1797 const cs_cell_mesh_t *cm,
1798 cs_real_t time_eval,
1799 cs_cell_builder_t *cb,
1800 void *input,
1801 double *values)
1802 {
1803 CS_UNUSED(cb);
1804 CS_UNUSED(input);
1805 CS_UNUSED(time_eval);
1806
1807 if (source == NULL)
1808 return;
1809
1810 /* Sanity checks */
1811 assert(values != NULL && cm != NULL);
1812
1813 const cs_real_t *s_values = (const cs_real_t *)source->context;
1814
1815 values[cm->n_fc] += s_values[0] * cm->vol_c;
1816 }
1817
1818 /*----------------------------------------------------------------------------*/
1819 /*!
1820 * \brief Compute the contribution for a cell related to a source term and
1821 * add it to the given array of values.
1822 * Case of a density defined on primal cells by a DoF function.
1823 * Case of scalar-valued face-based schemes
1824 *
1825 * \param[in] source pointer to a cs_xdef_t structure
1826 * \param[in] cm pointer to a cs_cell_mesh_t structure
1827 * \param[in] time_eval physical time at which one evaluates the term
1828 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1829 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1830 * \param[in, out] values pointer to the computed value
1831 */
1832 /*----------------------------------------------------------------------------*/
1833
1834 void
cs_source_term_pcsd_by_dof_func(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1835 cs_source_term_pcsd_by_dof_func(const cs_xdef_t *source,
1836 const cs_cell_mesh_t *cm,
1837 cs_real_t time_eval,
1838 cs_cell_builder_t *cb,
1839 void *input,
1840 double *values)
1841 {
1842 CS_UNUSED(cb);
1843 CS_UNUSED(time_eval);
1844 CS_UNUSED(input);
1845
1846 if (source == NULL)
1847 return;
1848
1849 cs_xdef_dof_context_t *context = (cs_xdef_dof_context_t *)source->context;
1850
1851 /* Sanity checks */
1852 assert(values != NULL && cm != NULL);
1853
1854 /* Up to now this should be the only location allowed */
1855 assert(cs_flag_test(context->loc, cs_flag_primal_cell));
1856
1857 /* Call the DoF function to evaluate the function at xc */
1858 double cell_eval;
1859 context->func(1, &(cm->c_id), true, /* compacted output ? */
1860 context->input,
1861 &cell_eval);
1862
1863 values[cm->n_fc] += cell_eval * cm->vol_c;
1864 }
1865
1866 /*----------------------------------------------------------------------------*/
1867 /*!
1868 * \brief Compute the contribution for a cell related to a source term and
1869 * add it to the given array of values.
1870 * Case of a density defined on primal cells by a DoF function.
1871 * Case of vector-valued face-based schemes
1872 *
1873 * \param[in] source pointer to a cs_xdef_t structure
1874 * \param[in] cm pointer to a cs_cell_mesh_t structure
1875 * \param[in] time_eval physical time at which one evaluates the term
1876 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1877 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1878 * \param[in, out] values pointer to the computed value
1879 */
1880 /*----------------------------------------------------------------------------*/
1881
1882 void
cs_source_term_pcvd_by_dof_func(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1883 cs_source_term_pcvd_by_dof_func(const cs_xdef_t *source,
1884 const cs_cell_mesh_t *cm,
1885 cs_real_t time_eval,
1886 cs_cell_builder_t *cb,
1887 void *input,
1888 double *values)
1889 {
1890 CS_UNUSED(cb);
1891 CS_UNUSED(time_eval);
1892 CS_UNUSED(input);
1893
1894 if (source == NULL)
1895 return;
1896
1897 cs_xdef_dof_context_t *context = (cs_xdef_dof_context_t *)source->context;
1898
1899 /* Sanity checks */
1900 assert(values != NULL && cm != NULL);
1901
1902 /* Up to now this should be the only location allowed */
1903 assert(cs_flag_test(context->loc, cs_flag_primal_cell));
1904
1905 /* Call the DoF function to evaluate the function at xc */
1906 cs_real_t cell_eval[3];
1907 context->func(1, &(cm->c_id), true, /* compacted output ? */
1908 context->input,
1909 cell_eval);
1910
1911 for (int k = 0; k < 3; k++)
1912 values[3*cm->n_fc+k] += cell_eval[k] * cm->vol_c;
1913 }
1914
1915 /*----------------------------------------------------------------------------*/
1916 /*!
1917 * \brief Compute the contribution for a cell related to a source term and
1918 * add it the given array of values.
1919 * Case of a vector-valued density (vd) defined on primal cells
1920 * by a value.
1921 * Case of face-based schemes
1922 *
1923 * \param[in] source pointer to a cs_xdef_t structure
1924 * \param[in] cm pointer to a cs_cell_mesh_t structure
1925 * \param[in] time_eval physical time at which one evaluates the term
1926 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1927 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
1928 * \param[in, out] values pointer to the computed value
1929 */
1930 /*----------------------------------------------------------------------------*/
1931
1932 void
cs_source_term_pcvd_by_value(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1933 cs_source_term_pcvd_by_value(const cs_xdef_t *source,
1934 const cs_cell_mesh_t *cm,
1935 cs_real_t time_eval,
1936 cs_cell_builder_t *cb,
1937 void *input,
1938 double *values)
1939 {
1940 CS_UNUSED(cb);
1941 CS_UNUSED(input);
1942 CS_UNUSED(time_eval);
1943
1944 if (source == NULL)
1945 return;
1946
1947 /* Sanity checks */
1948 assert(values != NULL && cm != NULL);
1949
1950 const cs_real_t *v_input = (const cs_real_t *)source->context;
1951 for (int k = 0; k < source->dim; k++)
1952 values[source->dim*cm->n_fc + k] = v_input[k] * cm->vol_c;
1953 }
1954
1955 /*----------------------------------------------------------------------------*/
1956 /*!
1957 * \brief Compute the contribution for a cell related to a source term and
1958 * add it to the given array of values.
1959 * Case of a scalar density defined at primal cells by an analytical
1960 * function.
1961 * Use the barycentric approximation as quadrature to evaluate the
1962 * integral. Exact for linear function.
1963 * Case of face-based schemes
1964 *
1965 * \param[in] source pointer to a cs_xdef_t structure
1966 * \param[in] cm pointer to a cs_cell_mesh_t structure
1967 * \param[in] time_eval physical time at which one evaluates the term
1968 * \param[in, out] cb pointer to a cs_cell_builder_t structure
1969 * \param[in] input NULL or pointer to a structure cast on-the-fly
1970 * \param[in, out] values pointer to the computed values
1971 */
1972 /*----------------------------------------------------------------------------*/
1973
1974 void
cs_source_term_pcsd_bary_by_analytic(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)1975 cs_source_term_pcsd_bary_by_analytic(const cs_xdef_t *source,
1976 const cs_cell_mesh_t *cm,
1977 cs_real_t time_eval,
1978 cs_cell_builder_t *cb,
1979 void *input,
1980 double *values)
1981 {
1982 CS_UNUSED(cb);
1983 CS_UNUSED(input);
1984
1985 if (source == NULL)
1986 return;
1987
1988 /* Sanity checks */
1989 assert(values != NULL && cm != NULL);
1990
1991 cs_xdef_analytic_context_t *ac =
1992 (cs_xdef_analytic_context_t *)source->context;
1993
1994 /* Call the analytic function to evaluate the function at xc */
1995 double eval_xc;
1996 ac->func(time_eval, 1, NULL, (const cs_real_t *)cm->xc,
1997 true, /* compacted output ? */
1998 ac->input,
1999 &eval_xc);
2000
2001 values[cm->n_fc] += cm->vol_c * eval_xc;
2002 }
2003
2004 /*----------------------------------------------------------------------------*/
2005 /*!
2006 * \brief Compute the contribution of a source term for a cell and add it to
2007 * the given array of values.
2008 * Case of a scalar density (sd) defined on primal cells by an analytic
2009 * function.
2010 * Case of face-based schemes
2011 *
2012 * \param[in] source pointer to a cs_xdef_t structure
2013 * \param[in] cm pointer to a cs_cell_mesh_t structure
2014 * \param[in] time_eval physical time at which one evaluates the term
2015 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2016 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
2017 * \param[in, out] values pointer to the computed value
2018 */
2019 /*----------------------------------------------------------------------------*/
2020
2021 void
cs_source_term_pcsd_by_analytic(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)2022 cs_source_term_pcsd_by_analytic(const cs_xdef_t *source,
2023 const cs_cell_mesh_t *cm,
2024 cs_real_t time_eval,
2025 cs_cell_builder_t *cb,
2026 void *input,
2027 double *values)
2028 {
2029 CS_UNUSED(input);
2030 if (source == NULL)
2031 return;
2032
2033 /* Sanity checks */
2034 assert(values != NULL && cm != NULL);
2035 assert(cs_eflag_test(cm->flag,
2036 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
2037 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
2038
2039 assert(source->dim == 1);
2040
2041 if (source->qtype == CS_QUADRATURE_BARY)
2042 cs_source_term_pcsd_bary_by_analytic(source, cm, time_eval, cb, input,
2043 values);
2044
2045 else {
2046
2047 const cs_real_t *xv = cm->xv;
2048
2049 double cell_values = 0.0;
2050
2051 cs_quadrature_tetra_integral_t *qfunc =
2052 cs_quadrature_get_tetra_integral(1, source->qtype);
2053
2054 cs_xdef_analytic_context_t *ac =
2055 (cs_xdef_analytic_context_t *)source->context;
2056
2057 /* Switch according to the cell type: optimised version for tetra */
2058 switch (cm->type) {
2059
2060 case FVM_CELL_TETRA:
2061 {
2062 assert(cm->n_fc == 4 && cm->n_vc == 4);
2063 qfunc(time_eval, xv, xv+3, xv+6, xv+9, cm->vol_c, ac->func, ac->input,
2064 &cell_values);
2065 }
2066 break;
2067
2068 case FVM_CELL_PYRAM:
2069 case FVM_CELL_PRISM:
2070 case FVM_CELL_HEXA:
2071 case FVM_CELL_POLY:
2072 {
2073 for (short int f = 0; f < cm->n_fc; f++) {
2074
2075 const cs_quant_t pfq = cm->face[f];
2076 const double hf_coef = cs_math_1ov3 * cm->hfc[f];
2077 const int start = cm->f2e_idx[f];
2078 const int end = cm->f2e_idx[f+1];
2079 const short int n_vf = end - start; /* #vertices (=#edges) */
2080 const short int *f2e_ids = cm->f2e_ids + start;
2081
2082 assert(n_vf > 2);
2083 switch(n_vf){
2084
2085 case 3: /* triangle (optimized version, no subdivision) */
2086 {
2087 short int v0, v1, v2;
2088 cs_cell_mesh_get_next_3_vertices(f2e_ids, cm->e2v_ids,
2089 &v0, &v1, &v2);
2090
2091 qfunc(time_eval, cm->xc, xv+3*v0, xv+3*v1, xv+3*v2,
2092 hf_coef*pfq.meas,
2093 ac->func, ac->input,
2094 &cell_values);
2095 }
2096 break;
2097
2098 default:
2099 {
2100 const double *tef = cm->tef + start;
2101
2102 for (short int e = 0; e < n_vf; e++) { /* Loop on face edges */
2103
2104 /* Edge-related variables */
2105 const short int e0 = f2e_ids[e];
2106 const double *xv0 = xv + 3*cm->e2v_ids[2*e0];
2107 const double *xv1 = xv + 3*cm->e2v_ids[2*e0+1];
2108
2109 qfunc(time_eval, cm->xc, pfq.center, xv0, xv1,
2110 hf_coef*tef[e], ac->func, ac->input, &cell_values);
2111
2112 }
2113 }
2114 break;
2115
2116 } /* End of switch */
2117
2118 } /* End of loop on faces */
2119
2120 }
2121 break;
2122
2123 default:
2124 bft_error(__FILE__, __LINE__, 0, "%s: Unknown cell-type.\n", __func__);
2125 break;
2126
2127 } /* End of switch on the cell-type */
2128
2129 values[cm->n_fc] += cell_values;
2130
2131 } /* If not a barycentric quadrature */
2132 }
2133
2134 /*----------------------------------------------------------------------------*/
2135 /*!
2136 * \brief Compute the contribution for a cell related to a source term and
2137 * add it to the given array of values.
2138 * Case of a scalar density defined at primal cells by an array.
2139 *
2140 * \param[in] source pointer to a cs_xdef_t structure
2141 * \param[in] cm pointer to a cs_cell_mesh_t structure
2142 * \param[in] time_eval physical time at which one evaluates the term
2143 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2144 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
2145 * \param[in, out] values pointer to the computed values
2146 */
2147 /*----------------------------------------------------------------------------*/
2148
2149 void
cs_source_term_pcsd_by_array(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)2150 cs_source_term_pcsd_by_array(const cs_xdef_t *source,
2151 const cs_cell_mesh_t *cm,
2152 cs_real_t time_eval,
2153 cs_cell_builder_t *cb,
2154 void *input,
2155 double *values)
2156 {
2157 CS_UNUSED(cb);
2158 CS_UNUSED(input);
2159 CS_UNUSED(time_eval);
2160
2161 if (source == NULL)
2162 return;
2163
2164 /* Sanity checks */
2165 assert(values != NULL && cm != NULL);
2166
2167 const cs_xdef_array_context_t *ai =
2168 (cs_xdef_array_context_t *)source->context;
2169
2170 assert(cs_flag_test(ai->loc, cs_flag_primal_cell));
2171 values[cm->n_fc] += ai->values[cm->c_id];
2172 }
2173
2174 /*----------------------------------------------------------------------------*/
2175 /*!
2176 * \brief Compute the contribution for a cell related to a source term and
2177 * add it to the given array of values.
2178 * Case of a vector-valued density defined at primal cells by an
2179 * analytical function.
2180 * Use the barycentric approximation as quadrature to evaluate the
2181 * integral. Exact for linear function.
2182 * Case of face-based schemes
2183 *
2184 * \param[in] source pointer to a cs_xdef_t structure
2185 * \param[in] cm pointer to a cs_cell_mesh_t structure
2186 * \param[in] time_eval physical time at which one evaluates the term
2187 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2188 * \param[in] input NULL or pointer to a structure cast on-the-fly
2189 * \param[in, out] values pointer to the computed values
2190 */
2191 /*----------------------------------------------------------------------------*/
2192
2193 void
cs_source_term_pcvd_bary_by_analytic(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)2194 cs_source_term_pcvd_bary_by_analytic(const cs_xdef_t *source,
2195 const cs_cell_mesh_t *cm,
2196 cs_real_t time_eval,
2197 cs_cell_builder_t *cb,
2198 void *input,
2199 double *values)
2200 {
2201 CS_UNUSED(cb);
2202 CS_UNUSED(input);
2203
2204 if (source == NULL)
2205 return;
2206
2207 /* Sanity checks */
2208 assert(values != NULL && cm != NULL);
2209 assert(source->dim == 3);
2210
2211 cs_xdef_analytic_context_t *ac =
2212 (cs_xdef_analytic_context_t *)source->context;
2213
2214 /* Call the analytic function to evaluate the function at xc */
2215 double eval_xc[3];
2216 ac->func(time_eval, 1, NULL, (const cs_real_t *)cm->xc,
2217 true, /* compacted output ? */
2218 ac->input,
2219 eval_xc);
2220
2221 cs_real_t *c_val = values + 3*cm->n_fc;
2222 for (int k = 0; k < source->dim; k++)
2223 c_val[k] += cm->vol_c * eval_xc[k];
2224 }
2225
2226 /*----------------------------------------------------------------------------*/
2227 /*!
2228 * \brief Compute the contribution of a source term for a cell and add it to
2229 * the given array of values.
2230 * Case of a vector density (vd) defined on primal cells by an analytic
2231 * function.
2232 * Case of face-based schemes
2233 *
2234 * \param[in] source pointer to a cs_xdef_t structure
2235 * \param[in] cm pointer to a cs_cell_mesh_t structure
2236 * \param[in] time_eval physical time at which one evaluates the term
2237 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2238 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
2239 * \param[in, out] values pointer to the computed value
2240 */
2241 /*----------------------------------------------------------------------------*/
2242
2243 void
cs_source_term_pcvd_by_analytic(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)2244 cs_source_term_pcvd_by_analytic(const cs_xdef_t *source,
2245 const cs_cell_mesh_t *cm,
2246 cs_real_t time_eval,
2247 cs_cell_builder_t *cb,
2248 void *input,
2249 double *values)
2250 {
2251 CS_UNUSED(input);
2252
2253 if (source == NULL)
2254 return;
2255
2256 /* Sanity checks */
2257 assert(values != NULL && cm != NULL);
2258 assert(cs_eflag_test(cm->flag,
2259 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
2260 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
2261 assert(source->dim == 3);
2262
2263 if (source->qtype == CS_QUADRATURE_BARY)
2264 cs_source_term_pcvd_bary_by_analytic(source, cm, time_eval, cb, input,
2265 values);
2266
2267 else {
2268
2269 const cs_real_t *xv = cm->xv;
2270
2271 cs_real_3_t cell_values = {0.0, 0.0, 0.0};
2272
2273 cs_quadrature_tetra_integral_t *qfunc =
2274 cs_quadrature_get_tetra_integral(3, source->qtype);
2275
2276 cs_xdef_analytic_context_t *ac =
2277 (cs_xdef_analytic_context_t *)source->context;
2278
2279 /* Switch according to the cell type: optimized version for tetra */
2280 switch (cm->type) {
2281
2282 case FVM_CELL_TETRA:
2283 {
2284 assert(cm->n_fc == 4 && cm->n_vc == 4);
2285 qfunc(time_eval, xv, xv+3, xv+6, xv+9, cm->vol_c,
2286 ac->func, ac->input,
2287 cell_values);
2288 }
2289 break;
2290
2291 case FVM_CELL_PYRAM:
2292 case FVM_CELL_PRISM:
2293 case FVM_CELL_HEXA:
2294 case FVM_CELL_POLY:
2295 {
2296 for (short int f = 0; f < cm->n_fc; f++) {
2297
2298 const cs_quant_t pfq = cm->face[f];
2299 const double hf_coef = cs_math_1ov3 * cm->hfc[f];
2300 const int start = cm->f2e_idx[f];
2301 const int end = cm->f2e_idx[f+1];
2302 const short int n_vf = end - start; // #vertices (=#edges)
2303 const short int *f2e_ids = cm->f2e_ids + start;
2304
2305 assert(n_vf > 2);
2306 switch(n_vf){
2307
2308 case 3: /* triangle (optimized version, no subdivision) */
2309 {
2310 short int v0, v1, v2;
2311 cs_cell_mesh_get_next_3_vertices(f2e_ids, cm->e2v_ids,
2312 &v0, &v1, &v2);
2313
2314 qfunc(time_eval, cm->xc, xv+3*v0, xv+3*v1, xv+3*v2,
2315 hf_coef*pfq.meas,
2316 ac->func, ac->input,
2317 cell_values);
2318 }
2319 break;
2320
2321 default:
2322 {
2323 const double *tef = cm->tef + start;
2324
2325 for (short int e = 0; e < n_vf; e++) { /* Loop on face edges */
2326
2327 /* Edge-related variables */
2328 const short int e0 = f2e_ids[e];
2329 const double *xv0 = xv + 3*cm->e2v_ids[2*e0];
2330 const double *xv1 = xv + 3*cm->e2v_ids[2*e0+1];
2331
2332 qfunc(time_eval, cm->xc, pfq.center, xv0, xv1,
2333 hf_coef*tef[e],
2334 ac->func, ac->input,
2335 cell_values);
2336 }
2337 }
2338 break;
2339
2340 } /* End of switch */
2341
2342 } /* End of loop on faces */
2343
2344 }
2345 break;
2346
2347 default:
2348 bft_error(__FILE__, __LINE__, 0, "%s: Unknown cell-type.\n", __func__);
2349 break;
2350
2351 } /* End of switch on the cell-type */
2352
2353 cs_real_t *c_val = values + 3*cm->n_fc;
2354 c_val[0] += cell_values[0];
2355 c_val[1] += cell_values[1];
2356 c_val[2] += cell_values[2];
2357
2358 } /* If not a barycentric quadrature */
2359 }
2360
2361 /*----------------------------------------------------------------------------*/
2362 /*!
2363 * \brief Compute the contribution of a source term for a cell and add it to
2364 * the given array of values.
2365 * Case of a vector density (vd) defined on primal cells by an array
2366 *
2367 * \param[in] source pointer to a cs_xdef_t structure
2368 * \param[in] cm pointer to a cs_cell_mesh_t structure
2369 * \param[in] time_eval physical time at which one evaluates the term
2370 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2371 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
2372 * \param[in, out] values pointer to the computed value
2373 */
2374 /*----------------------------------------------------------------------------*/
2375
2376 void
cs_source_term_pcvd_by_array(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)2377 cs_source_term_pcvd_by_array(const cs_xdef_t *source,
2378 const cs_cell_mesh_t *cm,
2379 cs_real_t time_eval,
2380 cs_cell_builder_t *cb,
2381 void *input,
2382 double *values)
2383 {
2384 CS_UNUSED(cb);
2385 CS_UNUSED(input);
2386 CS_UNUSED(time_eval);
2387
2388 if (source == NULL)
2389 return;
2390
2391 /* Sanity checks */
2392 assert(values != NULL && cm != NULL);
2393 assert(source->dim == 3);
2394
2395 const cs_xdef_array_context_t *ai =
2396 (cs_xdef_array_context_t *)source->context;
2397 const double *arr = ai->values + 3*cm->c_id;
2398
2399 assert(cs_flag_test(ai->loc, cs_flag_primal_cell));
2400
2401 double *val_c = values + 3*cm->n_fc;
2402 val_c[0] += arr[0];
2403 val_c[1] += arr[1];
2404 val_c[2] += arr[2];
2405 }
2406
2407 /*----------------------------------------------------------------------------*/
2408 /*!
2409 * \brief Compute the contribution of a source term for a cell and add it to
2410 * the given array of values.
2411 * Case of a scalar density (sd) defined on primal cells by a value.
2412 * Case of HHO schemes
2413 *
2414 * \param[in] source pointer to a cs_xdef_t structure
2415 * \param[in] cm pointer to a cs_cell_mesh_t structure
2416 * \param[in] time_eval physical time at which one evaluates the term
2417 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2418 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
2419 * \param[in, out] values pointer to the computed value
2420 */
2421 /*----------------------------------------------------------------------------*/
2422
2423 void
cs_source_term_hhosd_by_value(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)2424 cs_source_term_hhosd_by_value(const cs_xdef_t *source,
2425 const cs_cell_mesh_t *cm,
2426 cs_real_t time_eval,
2427 cs_cell_builder_t *cb,
2428 void *input,
2429 double *values)
2430 {
2431 CS_UNUSED(cb);
2432 CS_UNUSED(time_eval);
2433
2434 if (source == NULL)
2435 return;
2436
2437 /* Sanity checks */
2438 assert(values != NULL && cm != NULL && input != NULL);
2439
2440 const cs_real_t *const_val = (const cs_real_t *)source->context;
2441
2442 cs_hho_builder_t *hhob = (cs_hho_builder_t *)input;
2443 cs_real_t *cell_values = values + cm->n_fc * hhob->face_basis[0]->size;
2444
2445 const cs_basis_func_t *cbf = hhob->cell_basis;
2446
2447 switch (cbf->poly_order) {
2448
2449 case 0:
2450 case 1:
2451 cbf->eval_all_at_point(cbf, cm->xc, cell_values);
2452 for (int i = 0; i < cbf->size; i++)
2453 cell_values[i] *= cm->vol_c * const_val[0];
2454 break;
2455
2456 default:
2457 /* Reset cell values */
2458 memset(cell_values, 0, sizeof(cs_real_t)*cbf->size);
2459
2460 /* Switch according to the cell type: optimised version for tetra */
2461 switch (cm->type) {
2462
2463 case FVM_CELL_TETRA:
2464 assert(cm->n_fc == 4 && cm->n_vc == 4);
2465 _hho_add_tetra_by_val(const_val[0], cbf,
2466 cm->xv, cm->xv + 3, cm->xv + 6, cm->xv + 9,
2467 cm->vol_c, cb, cell_values);
2468 break;
2469
2470 case FVM_CELL_PYRAM:
2471 case FVM_CELL_PRISM:
2472 case FVM_CELL_HEXA:
2473 case FVM_CELL_POLY:
2474 {
2475 for (short int f = 0; f < cm->n_fc; f++) {
2476
2477 const cs_quant_t pfq = cm->face[f];
2478 const double hf_coef = cs_math_1ov3 * cm->hfc[f];
2479 const int start = cm->f2e_idx[f];
2480 const int end = cm->f2e_idx[f+1];
2481 const short int n_vf = end - start; /* #vertices (=#edges) */
2482 const short int *f2e_ids = cm->f2e_ids + start;
2483
2484 assert(n_vf > 2);
2485 switch(n_vf){
2486
2487 case CS_TRIANGLE_CASE: /* Optimized version, no subdivision */
2488 {
2489 short int v0, v1, v2;
2490 cs_cell_mesh_get_next_3_vertices(f2e_ids, cm->e2v_ids,
2491 &v0, &v1, &v2);
2492
2493 _hho_add_tetra_by_val(const_val[0], cbf,
2494 cm->xv+3*v0, cm->xv+3*v1, cm->xv+3*v2,
2495 cm->xc,
2496 hf_coef * pfq.meas, cb, cell_values);
2497 }
2498 break;
2499
2500 default:
2501 {
2502 const double *tef = cm->tef + start;
2503
2504 for (short int e = 0; e < n_vf; e++) { /* Loop on face edges */
2505
2506 /* Edge-related variables */
2507 const short int e0 = f2e_ids[e];
2508 const double *xv0 = cm->xv + 3*cm->e2v_ids[2*e0];
2509 const double *xv1 = cm->xv + 3*cm->e2v_ids[2*e0+1];
2510
2511 _hho_add_tetra_by_val(const_val[0], cbf,
2512 xv0, xv1, pfq.center, cm->xc,
2513 hf_coef*tef[e], cb, cell_values);
2514 }
2515 }
2516 break;
2517
2518 } /* End of switch */
2519
2520 } /* End of loop on faces */
2521
2522 }
2523 break;
2524
2525 default:
2526 bft_error(__FILE__, __LINE__, 0, _(" Unknown cell-type.\n"));
2527 break;
2528
2529 } /* End of switch on the cell-type */
2530 break;
2531
2532 } /* Switch on polynomial order */
2533 }
2534
2535 /*----------------------------------------------------------------------------*/
2536 /*!
2537 * \brief Compute the contribution of a source term for a cell and add it to
2538 * the given array of values.
2539 * Case of a scalar density (sd) defined on primal cells by an analytic
2540 * function.
2541 * Case of HHO schemes
2542 *
2543 * \param[in] source pointer to a cs_xdef_t structure
2544 * \param[in] cm pointer to a cs_cell_mesh_t structure
2545 * \param[in] time_eval physical time at which one evaluates the term
2546 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2547 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
2548 * \param[in, out] values pointer to the computed value
2549 */
2550 /*----------------------------------------------------------------------------*/
2551
2552 void
cs_source_term_hhosd_by_analytic(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)2553 cs_source_term_hhosd_by_analytic(const cs_xdef_t *source,
2554 const cs_cell_mesh_t *cm,
2555 cs_real_t time_eval,
2556 cs_cell_builder_t *cb,
2557 void *input,
2558 double *values)
2559 {
2560 CS_UNUSED(cb);
2561 CS_UNUSED(time_eval);
2562
2563 if (source == NULL)
2564 return;
2565
2566 /* Sanity checks */
2567 assert(values != NULL && cm != NULL && input != NULL);
2568 assert(cs_eflag_test(cm->flag,
2569 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
2570 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
2571
2572 cs_hho_builder_t *hhob = (cs_hho_builder_t *)input;
2573 cs_xdef_analytic_context_t *ac =
2574 (cs_xdef_analytic_context_t *)source->context;
2575 cs_real_t *cell_values = values + cm->n_fc * hhob->face_basis[0]->size;
2576
2577 const cs_basis_func_t *cbf = hhob->cell_basis;
2578
2579 /* Reset cell values */
2580 memset(cell_values, 0, sizeof(cs_real_t)*cbf->size);
2581
2582 /* Switch according to the cell type: optimised version for tetra */
2583 switch (cm->type) {
2584
2585 case FVM_CELL_TETRA:
2586 {
2587 assert(cm->n_fc == 4 && cm->n_vc == 4);
2588 _hho_add_tetra_by_ana(ac, cbf,
2589 cm->xv, cm->xv+3, cm->xv+6, cm->xv+9,
2590 cm->vol_c, time_eval,
2591 cb, cell_values);
2592 }
2593 break;
2594
2595 case FVM_CELL_PYRAM:
2596 case FVM_CELL_PRISM:
2597 case FVM_CELL_HEXA:
2598 case FVM_CELL_POLY:
2599 {
2600 for (short int f = 0; f < cm->n_fc; f++) {
2601
2602 const cs_quant_t pfq = cm->face[f];
2603 const double hf_coef = cs_math_1ov3 * cm->hfc[f];
2604 const int start = cm->f2e_idx[f];
2605 const int end = cm->f2e_idx[f+1];
2606 const short int n_vf = end - start; /* #vertices (=#edges) */
2607 const short int *f2e_ids = cm->f2e_ids + start;
2608
2609 assert(n_vf > 2);
2610 switch(n_vf){
2611
2612 case CS_TRIANGLE_CASE: /* triangle (optimized version, no subdivision) */
2613 {
2614 short int v0, v1, v2;
2615 cs_cell_mesh_get_next_3_vertices(f2e_ids, cm->e2v_ids, &v0, &v1, &v2);
2616
2617 _hho_add_tetra_by_ana(ac, cbf,
2618 cm->xv+3*v0, cm->xv+3*v1, cm->xv+3*v2, cm->xc,
2619 hf_coef * pfq.meas, time_eval,
2620 cb, cell_values);
2621 }
2622 break;
2623
2624 default:
2625 {
2626 const double *tef = cm->tef + start;
2627
2628 for (short int e = 0; e < n_vf; e++) { /* Loop on face edges */
2629
2630 /* Edge-related variables */
2631 const short int e0 = f2e_ids[e];
2632 const double *xv0 = cm->xv + 3*cm->e2v_ids[2*e0];
2633 const double *xv1 = cm->xv + 3*cm->e2v_ids[2*e0+1];
2634
2635 _hho_add_tetra_by_ana(ac, cbf,
2636 xv0, xv1, pfq.center, cm->xc,
2637 hf_coef*tef[e], time_eval,
2638 cb, cell_values);
2639 }
2640 }
2641 break;
2642
2643 } /* End of switch */
2644
2645 } /* End of loop on faces */
2646
2647 }
2648 break;
2649
2650 default:
2651 bft_error(__FILE__, __LINE__, 0, _(" Unknown cell-type.\n"));
2652 break;
2653
2654 } /* End of switch on the cell-type */
2655 }
2656
2657 /*----------------------------------------------------------------------------*/
2658 /*!
2659 * \brief Compute the contribution of a source term for a cell and add it to
2660 * the given array of values.
2661 * Case of a vector field (vd) defined on primal cells by an analytic
2662 * function.
2663 * Case of HHO schemes
2664 *
2665 * \param[in] source pointer to a cs_xdef_t structure
2666 * \param[in] cm pointer to a cs_cell_mesh_t structure
2667 * \param[in] time_eval physical time at which one evaluates the term
2668 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2669 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
2670 * \param[in, out] values pointer to the computed value
2671 */
2672 /*----------------------------------------------------------------------------*/
2673
2674 void
cs_source_term_hhovd_by_analytic(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)2675 cs_source_term_hhovd_by_analytic(const cs_xdef_t *source,
2676 const cs_cell_mesh_t *cm,
2677 cs_real_t time_eval,
2678 cs_cell_builder_t *cb,
2679 void *input,
2680 double *values)
2681 {
2682 CS_UNUSED(cb);
2683 if (source == NULL)
2684 return;
2685
2686 /* Sanity checks */
2687 assert(values != NULL && cm != NULL && input != NULL);
2688 assert(cs_eflag_test(cm->flag,
2689 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
2690 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
2691
2692 cs_hho_builder_t *hhob = (cs_hho_builder_t *)input;
2693 cs_xdef_analytic_context_t *ac =
2694 (cs_xdef_analytic_context_t *)source->context;
2695 cs_real_t *cell_values = values + 3*cm->n_fc * hhob->face_basis[0]->size;
2696
2697 const cs_basis_func_t *cbf = hhob->cell_basis;
2698
2699 /* Reset cell values */
2700 memset(cell_values, 0, 3*sizeof(cs_real_t)*cbf->size);
2701
2702 /* Switch according to the cell type: optimised version for tetra */
2703 switch (cm->type) {
2704
2705 case FVM_CELL_TETRA:
2706 {
2707 assert(cm->n_fc == 4 && cm->n_vc == 4);
2708 _hho_add_tetra_by_ana_vd(ac, cbf,
2709 cm->xv, cm->xv+3, cm->xv+6, cm->xv+9,
2710 cm->vol_c, time_eval,
2711 cb, cell_values);
2712 }
2713 break;
2714
2715 case FVM_CELL_PYRAM:
2716 case FVM_CELL_PRISM:
2717 case FVM_CELL_HEXA:
2718 case FVM_CELL_POLY:
2719 {
2720 for (short int f = 0; f < cm->n_fc; f++) {
2721
2722 const cs_quant_t pfq = cm->face[f];
2723 const double hf_coef = cs_math_1ov3 * cm->hfc[f];
2724 const int start = cm->f2e_idx[f];
2725 const int end = cm->f2e_idx[f+1];
2726 const short int n_vf = end - start; /* #vertices (=#edges) */
2727 const short int *f2e_ids = cm->f2e_ids + start;
2728
2729 assert(n_vf > 2);
2730 switch(n_vf){
2731
2732 case 3: /* triangle (optimized version, no subdivision) */
2733 {
2734 short int v0, v1, v2;
2735 cs_cell_mesh_get_next_3_vertices(f2e_ids, cm->e2v_ids, &v0, &v1, &v2);
2736
2737 _hho_add_tetra_by_ana_vd(ac, cbf,
2738 cm->xv+3*v0, cm->xv+3*v1, cm->xv+3*v2,
2739 cm->xc, hf_coef * pfq.meas, time_eval,
2740 cb, cell_values);
2741 }
2742 break;
2743
2744 default:
2745 {
2746 const double *tef = cm->tef + start;
2747
2748 for (short int e = 0; e < n_vf; e++) { /* Loop on face edges */
2749
2750 /* Edge-related variables */
2751 const short int e0 = f2e_ids[e];
2752 const double *xv0 = cm->xv + 3*cm->e2v_ids[2*e0];
2753 const double *xv1 = cm->xv + 3*cm->e2v_ids[2*e0+1];
2754
2755 _hho_add_tetra_by_ana_vd(ac, cbf,
2756 xv0, xv1, pfq.center, cm->xc,
2757 hf_coef*tef[e], time_eval,
2758 cb, cell_values);
2759 }
2760 }
2761 break;
2762
2763 } /* End of switch */
2764
2765 } /* End of loop on faces */
2766
2767 }
2768 break;
2769
2770 default:
2771 bft_error(__FILE__, __LINE__, 0, _(" Unknown cell-type.\n"));
2772 break;
2773
2774 } /* End of switch on the cell-type */
2775 }
2776
2777 /*----------------------------------------------------------------------------*/
2778 /*!
2779 * \brief Compute the contribution for a cell related to a source term and
2780 * add it to the given array of values.
2781 * Case of a scalar flux defined at dual faces by a constant value.
2782 *
2783 * \param[in] source pointer to a cs_xdef_t structure
2784 * \param[in] cm pointer to a cs_cell_mesh_t structure
2785 * \param[in] time_eval physical time at which one evaluates the term
2786 * \param[in, out] cb pointer to a cs_cell_builder_t structure
2787 * \param[in, out] input pointer to an element cast on-the-fly (or NULL)
2788 * \param[in, out] values pointer to the computed values
2789 */
2790 /*----------------------------------------------------------------------------*/
2791
2792 void
cs_source_term_dfsf_by_value(const cs_xdef_t * source,const cs_cell_mesh_t * cm,cs_real_t time_eval,cs_cell_builder_t * cb,void * input,double * values)2793 cs_source_term_dfsf_by_value(const cs_xdef_t *source,
2794 const cs_cell_mesh_t *cm,
2795 cs_real_t time_eval,
2796 cs_cell_builder_t *cb,
2797 void *input,
2798 double *values)
2799 {
2800 if (source == NULL)
2801 return;
2802
2803 CS_UNUSED(input);
2804 CS_UNUSED(cb);
2805 CS_UNUSED(time_eval);
2806
2807 /* Sanity checks */
2808 assert(values != NULL && cm != NULL);
2809 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_DFQ));
2810
2811 const cs_real_t *vector = (const cs_real_t *)source->context;
2812
2813 /* Retrieve the values of the normal flux for each dual face */
2814 for (short int e = 0; e < cm->n_ec; e++)
2815 values[e] = cm->dface[e].meas *
2816 cs_math_3_dot_product(vector, cm->dface[e].unitv);
2817 }
2818
2819 /*----------------------------------------------------------------------------*/
2820
2821 END_C_DECLS
2822