1 /*============================================================================
2 * Functions to handle the evaluation of boundary conditions when building the
3 * algebraic system in CDO/HHO schemes
4 *============================================================================*/
5
6 /*
7 This file is part of Code_Saturne, a general-purpose CFD tool.
8
9 Copyright (C) 1998-2021 EDF S.A.
10
11 This program is free software; you can redistribute it and/or modify it under
12 the terms of the GNU General Public License as published by the Free Software
13 Foundation; either version 2 of the License, or (at your option) any later
14 version.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
19 details.
20
21 You should have received a copy of the GNU General Public License along with
22 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
23 Street, Fifth Floor, Boston, MA 02110-1301, USA.
24 */
25
26 /*----------------------------------------------------------------------------*/
27
28 #include "cs_defs.h"
29
30 /*----------------------------------------------------------------------------
31 * Standard C library headers
32 *----------------------------------------------------------------------------*/
33
34 #include <assert.h>
35 #include <string.h>
36
37 /*----------------------------------------------------------------------------
38 * Local headers
39 *----------------------------------------------------------------------------*/
40
41 #include <bft_mem.h>
42
43 #include "cs_boundary_zone.h"
44 #include "cs_equation_common.h"
45 #include "cs_evaluate.h"
46 #include "cs_xdef.h"
47
48 #if defined(DEBUG) && !defined(NDEBUG)
49 #include "cs_dbg.h"
50 #endif
51
52 /*----------------------------------------------------------------------------*/
53
54 BEGIN_C_DECLS
55
56 /*----------------------------------------------------------------------------
57 * Header for the current file
58 *----------------------------------------------------------------------------*/
59
60 #include "cs_equation_bc.h"
61
62 /*----------------------------------------------------------------------------*/
63
64 BEGIN_C_DECLS
65
66 /*============================================================================
67 * Type definitions and macros
68 *============================================================================*/
69
70 #define CS_EQUATION_BC_DBG 0
71
72 /*============================================================================
73 * Local private variables
74 *============================================================================*/
75
76 /*============================================================================
77 * Private function prototypes
78 *============================================================================*/
79
80 /*----------------------------------------------------------------------------*/
81 /*!
82 * \brief Set the Dirichlet BC values to the face vertices.
83 * Case of vertex-based schemes
84 *
85 * \param[in] dim number of values to assign to each vertex
86 * \param[in] n_vf number of vertices in a face
87 * \param[in] lst list of vertex numbering
88 * \param[in] eval result of the evaluation to set
89 * \param[in] is_constant same value for all vertices ?
90 * \param[in, out] vvals vertex values to update
91 * \param[in, out] counter counter to update
92 */
93 /*----------------------------------------------------------------------------*/
94
95 static inline void
_assign_vb_dirichlet_values(int dim,int n_vf,const cs_lnum_t * lst,const cs_real_t * eval,bool is_constant,cs_real_t * vvals,int counter[])96 _assign_vb_dirichlet_values(int dim,
97 int n_vf,
98 const cs_lnum_t *lst,
99 const cs_real_t *eval,
100 bool is_constant,
101 cs_real_t *vvals,
102 int counter[])
103 {
104 if (dim == 1) {
105
106 for (short int v = 0; v < n_vf; v++) {
107 const cs_lnum_t v_id = lst[v];
108 const short int _v = is_constant ? 0 : v;
109 counter[v_id] += 1;
110 vvals[v_id] += eval[_v];
111 }
112
113 }
114 else { /* Not a scalar-valued quantity */
115
116 for (short int v = 0; v < n_vf; v++) {
117 const cs_lnum_t v_id = lst[v];
118 const short int _v = is_constant ? 0 : v;
119 counter[v_id] += 1;
120 for (int k = 0; k < dim; k++)
121 vvals[dim*v_id + k] += eval[dim*_v + k];
122 }
123
124 }
125 }
126
127 /*----------------------------------------------------------------------------*/
128 /*!
129 * \brief Set members of the cs_cell_sys_t structure related to the boundary
130 * conditions. Only the generic part is done here. The remaining part
131 * is performed specifically for each scheme
132 *
133 * \param[in] face_bc pointer to a cs_cdo_bc_face_t structure
134 * \param[in] cm pointer to a cs_cell_mesh_t structure
135 * \param[in, out] csys pointer to a cs_cell_system_t structure
136 */
137 /*----------------------------------------------------------------------------*/
138
139 static inline void
_init_cell_sys_bc(const cs_cdo_bc_face_t * face_bc,const cs_cell_mesh_t * cm,cs_cell_sys_t * csys)140 _init_cell_sys_bc(const cs_cdo_bc_face_t *face_bc,
141 const cs_cell_mesh_t *cm,
142 cs_cell_sys_t *csys)
143 {
144 for (short int f = 0; f < cm->n_fc; f++) {
145
146 const cs_lnum_t bf_id = cm->f_ids[f] - cm->bface_shift;
147
148 csys->bf_ids[f] = bf_id;
149
150 if (bf_id > -1) { /* This is a boundary face */
151
152 csys->bf_flag[f] = face_bc->flag[bf_id];
153 csys->_f_ids[csys->n_bc_faces] = f;
154 csys->n_bc_faces++;
155
156 }
157
158 } /* Loop on cell faces */
159 }
160
161 /*----------------------------------------------------------------------------*/
162 /*!
163 * \brief Synchronize the boundary definitions related to the enforcement of
164 * a circulation along a boundary edge
165 *
166 * \param[in] connect pointer to a cs_cdo_connect_t structure
167 * \param[in] n_defs number of definitions
168 * \param[in] defs number of times the values has been updated
169 * \param[in, out] def2e_idx index array to define
170 * \param[in, out] def2e_ids array of ids to define
171 */
172 /*----------------------------------------------------------------------------*/
173
174 static void
_sync_circulation_def_at_edges(const cs_cdo_connect_t * connect,int n_defs,cs_xdef_t ** defs,cs_lnum_t def2e_idx[],cs_lnum_t def2e_ids[])175 _sync_circulation_def_at_edges(const cs_cdo_connect_t *connect,
176 int n_defs,
177 cs_xdef_t **defs,
178 cs_lnum_t def2e_idx[],
179 cs_lnum_t def2e_ids[])
180 {
181 if (n_defs == 0)
182 return;
183
184 const cs_lnum_t n_edges = connect->n_edges;
185 const cs_adjacency_t *f2e = connect->f2e;
186
187 int *e2def_ids = NULL;
188 BFT_MALLOC(e2def_ids, n_edges, int);
189 # pragma omp parallel for if (n_edges > CS_THR_MIN)
190 for (cs_lnum_t e = 0; e < n_edges; e++)
191 e2def_ids[e] = -1; /* default: not associated to a definition */
192
193 const cs_lnum_t face_shift = connect->n_faces[CS_INT_FACES];
194
195 for (int def_id = 0; def_id < n_defs; def_id++) {
196
197 /* Get and then set the definition of the initial condition */
198
199 const cs_xdef_t *def = defs[def_id];
200 assert(def->support == CS_XDEF_SUPPORT_BOUNDARY);
201
202 if ((def->meta & CS_CDO_BC_TANGENTIAL_DIRICHLET) ||
203 (def->meta & CS_CDO_BC_DIRICHLET)) {
204
205 const cs_zone_t *z = cs_boundary_zone_by_id(def->z_id);
206
207 for (cs_lnum_t i = 0; i < z->n_elts; i++) { /* Loop on selected faces */
208 const cs_lnum_t f_id = face_shift + z->elt_ids[i];
209 for (cs_lnum_t j = f2e->idx[f_id]; j < f2e->idx[f_id+1]; j++)
210 e2def_ids[f2e->ids[j]] = def_id;
211 }
212
213 } /* Enforcement of the tangential component */
214
215 } /* Loop on definitions */
216
217 if (connect->interfaces[CS_CDO_CONNECT_EDGE_SCAL] != NULL) {
218
219 /* Last definition is used if there is a conflict between several
220 definitions */
221
222 cs_interface_set_max(connect->interfaces[CS_CDO_CONNECT_EDGE_SCAL],
223 n_edges,
224 1, /* stride */
225 false, /* interlace (not useful here) */
226 CS_INT_TYPE, /* int */
227 e2def_ids);
228
229 }
230
231 /* 0. Initialization */
232
233 cs_lnum_t *count = NULL;
234 BFT_MALLOC(count, n_defs, cs_lnum_t);
235 memset(count, 0, n_defs*sizeof(cs_lnum_t));
236 memset(def2e_idx, 0, (n_defs+1)*sizeof(cs_lnum_t));
237
238 /* 1. Count the number of edges related to each definition */
239
240 for (cs_lnum_t e = 0; e < n_edges; e++)
241 if (e2def_ids[e] > -1)
242 def2e_idx[e2def_ids[e]+1] += 1;
243
244 /* 2. Build the index */
245
246 for (int def_id = 0; def_id < n_defs; def_id++)
247 def2e_idx[def_id+1] += def2e_idx[def_id];
248
249 /* 3. Build the list */
250
251 for (cs_lnum_t e = 0; e < n_edges; e++) {
252 const int def_id = e2def_ids[e];
253 if (def_id > -1) {
254 def2e_ids[def2e_idx[def_id] + count[def_id]] = e;
255 count[def_id] += 1;
256 }
257 }
258
259 /* Free memory */
260
261 BFT_FREE(e2def_ids);
262 BFT_FREE(count);
263 }
264
265 /*============================================================================
266 * Public function prototypes
267 *============================================================================*/
268
269 /*----------------------------------------------------------------------------*/
270 /*!
271 * \brief Set the values for the normal boundary flux stemming from the
272 * Neumann boundary conditions (zero is left where a Dirichlet is
273 * set. This can be updated later on)
274 *
275 * \param[in] t_eval time at which one performs the evaluation
276 * \param[in] cdoq pointer to a cs_cdo_quantities_t structure
277 * \param[in] eqp pointer to a cs_equation_param_t structure
278 * \param[in, out] values pointer to the array of values to set
279 */
280 /*----------------------------------------------------------------------------*/
281
282 void
cs_equation_init_boundary_flux_from_bc(cs_real_t t_eval,const cs_cdo_quantities_t * cdoq,const cs_equation_param_t * eqp,cs_real_t * values)283 cs_equation_init_boundary_flux_from_bc(cs_real_t t_eval,
284 const cs_cdo_quantities_t *cdoq,
285 const cs_equation_param_t *eqp,
286 cs_real_t *values)
287 {
288 /* We assume a homogeneous Neumann boundary condition as a default */
289
290 memset(values, 0, sizeof(cs_real_t)*cdoq->n_b_faces);
291
292 for (int def_id = 0; def_id < eqp->n_bc_defs; def_id++) {
293
294 const cs_xdef_t *def = eqp->bc_defs[def_id];
295 const cs_zone_t *z = cs_boundary_zone_by_id(def->z_id);
296 assert(def->support == CS_XDEF_SUPPORT_BOUNDARY);
297
298 if (cs_flag_test(def->meta, CS_CDO_BC_NEUMANN)) {
299
300 switch (def->type) {
301
302 case CS_XDEF_BY_VALUE:
303 {
304 const cs_real_t *constant_val = (cs_real_t *)def->context;
305
306 switch (eqp->dim) {
307
308 case 1: /* scalar-valued equation */
309 # pragma omp parallel for if (z->n_elts > CS_THR_MIN)
310 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
311 const cs_lnum_t elt_id =
312 (z->elt_ids != NULL) ? z->elt_ids[i] : i;
313 values[elt_id] = constant_val[0];
314 }
315 break;
316
317 default:
318 # pragma omp parallel for if (z->n_elts > CS_THR_MIN)
319 for (cs_lnum_t i = 0; i < z->n_elts; i++) {
320 const cs_lnum_t elt_id =
321 (z->elt_ids != NULL) ? z->elt_ids[i] : i;
322 for (int k = 0; k < eqp->dim; k++)
323 values[eqp->dim*elt_id + k] = constant_val[k];
324 }
325 break;
326
327 } /* switch on dimension */
328
329 }
330 break; /* definition by value */
331
332 case CS_XDEF_BY_ANALYTIC_FUNCTION:
333 {
334 cs_xdef_analytic_context_t *ac =
335 (cs_xdef_analytic_context_t *)def->context;
336
337 ac->func(t_eval,
338 z->n_elts, z->elt_ids, cdoq->b_face_center,
339 false, /* dense output ? */
340 ac->input,
341 values);
342 }
343 break;
344
345 default:
346 bft_error(__FILE__, __LINE__, 0, " %s: Invalid case.", __func__);
347 }
348
349 } /* Neumann boundary conditions */
350
351 } /* Loop on boundary conditions applied to the Richards equation */
352
353 }
354
355 /*----------------------------------------------------------------------------*/
356 /*!
357 * \brief Set the BC into a cellwise view of the current system.
358 * Case of vertex-based schemes
359 *
360 * \param[in] cm pointer to a cellwise view of the mesh
361 * \param[in] eqp pointer to a cs_equation_param_t structure
362 * \param[in] face_bc pointer to a cs_cdo_bc_face_t structure
363 * \param[in] vtx_bc_flag BC flags associated to vertices
364 * \param[in] dir_values Dirichlet values associated to each vertex
365 * \param[in] t_eval time at which one performs the evaluation
366 * \param[in, out] csys pointer to a cellwise view of the system
367 * \param[in, out] cb pointer to a cellwise builder
368 */
369 /*----------------------------------------------------------------------------*/
370
371 void
cs_equation_vb_set_cell_bc(const cs_cell_mesh_t * cm,const cs_equation_param_t * eqp,const cs_cdo_bc_face_t * face_bc,const cs_flag_t vtx_bc_flag[],const cs_real_t dir_values[],cs_real_t t_eval,cs_cell_sys_t * csys,cs_cell_builder_t * cb)372 cs_equation_vb_set_cell_bc(const cs_cell_mesh_t *cm,
373 const cs_equation_param_t *eqp,
374 const cs_cdo_bc_face_t *face_bc,
375 const cs_flag_t vtx_bc_flag[],
376 const cs_real_t dir_values[],
377 cs_real_t t_eval,
378 cs_cell_sys_t *csys,
379 cs_cell_builder_t *cb)
380 {
381 CS_UNUSED(cb);
382
383 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_EV | CS_FLAG_COMP_FE));
384
385 /* Initialize the common part */
386
387 _init_cell_sys_bc(face_bc, cm, csys);
388
389 const int d = eqp->dim;
390
391 /* First pass: Set the DoF flag and the Dirichlet values */
392
393 for (short int v = 0; v < cm->n_vc; v++) {
394
395 const cs_flag_t bc_flag = vtx_bc_flag[cm->v_ids[v]];
396 for (int k = 0; k < d; k++)
397 csys->dof_flag[d*v+k] = bc_flag;
398
399 if (cs_cdo_bc_is_dirichlet(bc_flag)) {
400 csys->has_dirichlet = true;
401 if (bc_flag & CS_CDO_BC_HMG_DIRICHLET)
402 continue; /* Nothing else to do for this vertex */
403 else {
404 assert(bc_flag & CS_CDO_BC_DIRICHLET);
405 for (int k = 0; k < d; k++)
406 csys->dir_values[d*v+k] = dir_values[d*cm->v_ids[v]+k];
407 }
408 }
409
410 } /* Loop on cell vertices */
411
412 /* Second pass: BC related to faces. Neumann or Robin or sliding condition
413 (for vector-valued equations). Identify which face is a boundary face */
414
415 for (short int f = 0; f < cm->n_fc; f++) {
416 if (csys->bf_ids[f] > -1) { /* This a boundary face */
417
418 switch(csys->bf_flag[f]) {
419
420 case CS_CDO_BC_NEUMANN:
421 csys->has_nhmg_neumann = true;
422 cs_equation_compute_neumann_sv(t_eval,
423 face_bc->def_ids[csys->bf_ids[f]],
424 f,
425 eqp,
426 cm,
427 csys->neu_values);
428 break;
429
430 case CS_CDO_BC_ROBIN:
431 csys->has_robin = true;
432
433 /* The values which define the Robin BC are stored for each boundary
434 face. This is different from Neumann and Dirichlet where the values
435 are defined at each vertices. Be aware of that when computing */
436
437 cs_equation_compute_robin(t_eval,
438 face_bc->def_ids[csys->bf_ids[f]],
439 f,
440 eqp,
441 cm,
442 csys->rob_values);
443 break;
444
445 case CS_CDO_BC_SLIDING:
446 csys->has_sliding = true;
447 break;
448
449 default: /* Nothing to do for */
450 /* case CS_CDO_BC_HMG_DIRICHLET: */
451 /* case CS_CDO_BC_DIRICHLET: */
452 /* case CS_CDO_BC_HMG_NEUMANN: */
453 break;
454
455 } /* End of switch */
456
457 } /* Boundary face */
458 } /* Loop on cell faces */
459 }
460
461 /*----------------------------------------------------------------------------*/
462 /*!
463 * \brief Set the BC into a cellwise view of the current system.
464 * Case of edge-based schemes
465 *
466 * \param[in] cm pointer to a cellwise view of the mesh
467 * \param[in] eqp pointer to a cs_equation_param_t structure
468 * \param[in] face_bc pointer to a cs_cdo_bc_face_t structure
469 * \param[in] dir_values Dirichlet values associated to each vertex
470 * \param[in, out] csys pointer to a cellwise view of the system
471 * \param[in, out] cb pointer to a cellwise builder
472 */
473 /*----------------------------------------------------------------------------*/
474
475 void
cs_equation_eb_set_cell_bc(const cs_cell_mesh_t * cm,const cs_equation_param_t * eqp,const cs_cdo_bc_face_t * face_bc,const cs_real_t dir_values[],cs_cell_sys_t * csys,cs_cell_builder_t * cb)476 cs_equation_eb_set_cell_bc(const cs_cell_mesh_t *cm,
477 const cs_equation_param_t *eqp,
478 const cs_cdo_bc_face_t *face_bc,
479 const cs_real_t dir_values[],
480 cs_cell_sys_t *csys,
481 cs_cell_builder_t *cb)
482 {
483 CS_UNUSED(cb);
484 CS_UNUSED(eqp);
485
486 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_FE));
487
488 /* Initialize the common part */
489
490 _init_cell_sys_bc(face_bc, cm, csys);
491
492 /* From BC related to faces to edges. */
493
494 for (short int f = 0; f < cm->n_fc; f++) {
495 if (csys->bf_ids[f] > -1) { /* This a boundary face */
496
497 switch(csys->bf_flag[f]) {
498
499 case CS_CDO_BC_DIRICHLET:
500 case CS_CDO_BC_TANGENTIAL_DIRICHLET:
501 csys->has_dirichlet = true;
502 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
503 const short int e = cm->f2e_ids[i];
504 csys->dof_flag[e] |= CS_CDO_BC_DIRICHLET;
505 csys->dir_values[e] = dir_values[cm->e_ids[e]];
506 }
507 break;
508
509 case CS_CDO_BC_HMG_DIRICHLET:
510 csys->has_dirichlet = true;
511 for (int i = cm->f2e_idx[f]; f < cm->f2e_idx[f+1]; f++) {
512 const short int e = cm->f2e_ids[i];
513 csys->dof_flag[e] |= CS_CDO_BC_HMG_DIRICHLET;
514 csys->dir_values[e] = 0.;
515 }
516 break;
517
518 default: /* Nothing to do for */
519 /* case CS_CDO_BC_HMG_DIRICHLET: */
520 /* case CS_CDO_BC_HMG_NEUMANN: */
521 break;
522
523 } /* End of switch */
524
525 } /* Boundary face */
526 } /* Loop on cell faces */
527 }
528
529 /*----------------------------------------------------------------------------*/
530 /*!
531 * \brief Set the BC into a cellwise view of the current system.
532 * Case of Face-based schemes
533 *
534 * \param[in] cm pointer to a cellwise view of the mesh
535 * \param[in] eqp pointer to a cs_equation_param_t structure
536 * \param[in] face_bc pointer to a cs_cdo_bc_face_t structure
537 * \param[in] dir_values Dirichlet values associated to each vertex
538 * \param[in, out] csys pointer to a cellwise view of the system
539 * \param[in, out] cb pointer to a cellwise builder
540 */
541 /*----------------------------------------------------------------------------*/
542
543 void
cs_equation_fb_set_cell_bc(const cs_cell_mesh_t * cm,const cs_equation_param_t * eqp,const cs_cdo_bc_face_t * face_bc,const cs_real_t dir_values[],cs_cell_sys_t * csys,cs_cell_builder_t * cb)544 cs_equation_fb_set_cell_bc(const cs_cell_mesh_t *cm,
545 const cs_equation_param_t *eqp,
546 const cs_cdo_bc_face_t *face_bc,
547 const cs_real_t dir_values[],
548 cs_cell_sys_t *csys,
549 cs_cell_builder_t *cb)
550 {
551 /* Initialize the common part */
552
553 _init_cell_sys_bc(face_bc, cm, csys);
554
555 const int d = eqp->dim;
556
557 /* Identify which face is a boundary face */
558
559 for (short int f = 0; f < cm->n_fc; f++) {
560 if (csys->bf_ids[f] > -1) { /* This a boundary face */
561
562 switch(csys->bf_flag[f]) {
563
564 case CS_CDO_BC_HMG_DIRICHLET:
565 csys->has_dirichlet = true;
566 for (int k = 0; k < d; k++)
567 csys->dof_flag[d*f + k] |= CS_CDO_BC_HMG_DIRICHLET;
568 break;
569
570 case CS_CDO_BC_DIRICHLET:
571 csys->has_dirichlet = true;
572 for (int k = 0; k < d; k++) {
573 csys->dof_flag[d*f + k] |= CS_CDO_BC_DIRICHLET;
574 csys->dir_values[d*f + k] = dir_values[d*csys->bf_ids[f] + k];
575 }
576 break;
577
578 case CS_CDO_BC_NEUMANN:
579 csys->has_nhmg_neumann = true;
580 for (int k = 0; k < d; k++)
581 csys->dof_flag[d*f + k] |= CS_CDO_BC_NEUMANN;
582
583 cs_equation_compute_neumann_fb(cb->t_bc_eval,
584 face_bc->def_ids[csys->bf_ids[f]],
585 f,
586 eqp,
587 cm,
588 csys->neu_values);
589 break;
590
591 case CS_CDO_BC_ROBIN:
592 csys->has_robin = true;
593 for (int k = 0; k < d; k++)
594 csys->dof_flag[d*f + k] |= CS_CDO_BC_ROBIN;
595
596 cs_equation_compute_robin(cb->t_bc_eval,
597 face_bc->def_ids[csys->bf_ids[f]],
598 f,
599 eqp,
600 cm,
601 csys->rob_values);
602 break;
603
604 case CS_CDO_BC_SLIDING:
605 csys->has_sliding = true;
606 break;
607
608 default:
609 /* Nothing to do for instance in case of homogeneous Neumann */
610 break;
611
612 } /* End of switch */
613
614 } /* Boundary face */
615 } /* Loop on cell faces */
616 }
617
618 /*----------------------------------------------------------------------------*/
619 /*!
620 * \brief Compute the values of the Dirichlet BCs when DoFs are attached to
621 * vertices
622 *
623 * \param[in] t_eval time at which one performs the evaluation
624 * \param[in] mesh pointer to a cs_mesh_t structure
625 * \param[in] quant pointer to a cs_cdo_quantities_t structure
626 * \param[in] connect pointer to a cs_cdo_connect_t struct.
627 * \param[in] eqp pointer to a cs_equation_param_t
628 * \param[in] face_bc pointer to a cs_cdo_bc_face_t structure
629 * \param[in, out] cb pointer to a cs_cell_builder_t structure
630 * \param[in, out] bcflag pointer to an array storing type of BC
631 * \param[in, out] values pointer to the array of values to set
632 */
633 /*----------------------------------------------------------------------------*/
634
635 void
cs_equation_compute_dirichlet_vb(cs_real_t t_eval,const cs_mesh_t * mesh,const cs_cdo_quantities_t * quant,const cs_cdo_connect_t * connect,const cs_equation_param_t * eqp,const cs_cdo_bc_face_t * face_bc,cs_cell_builder_t * cb,cs_flag_t * bcflag,cs_real_t * values)636 cs_equation_compute_dirichlet_vb(cs_real_t t_eval,
637 const cs_mesh_t *mesh,
638 const cs_cdo_quantities_t *quant,
639 const cs_cdo_connect_t *connect,
640 const cs_equation_param_t *eqp,
641 const cs_cdo_bc_face_t *face_bc,
642 cs_cell_builder_t *cb,
643 cs_flag_t *bcflag,
644 cs_real_t *values)
645 {
646 assert(face_bc != NULL && bcflag != NULL && values != NULL);
647
648 const cs_lnum_t *bf2v_idx = mesh->b_face_vtx_idx;
649 const cs_lnum_t *bf2v_lst = mesh->b_face_vtx_lst;
650
651 /* Initialization */
652
653 cs_real_t *bcvals = cs_equation_get_tmpbuf();
654 memset(bcvals, 0, eqp->dim*quant->n_vertices*sizeof(cs_real_t));
655
656 /* Number of faces with a Dir. related to a vertex */
657
658 int *counter = NULL;
659 BFT_MALLOC(counter, quant->n_vertices, int);
660 memset(counter, 0, quant->n_vertices*sizeof(int));
661
662 if (face_bc->is_steady == false) /* Update bcflag if needed */
663 cs_equation_set_vertex_bc_flag(connect, face_bc, bcflag);
664
665 /* Define array storing the Dirichlet values */
666
667 for (cs_lnum_t i = 0; i < face_bc->n_nhmg_dir_faces; i++) {
668
669 const cs_lnum_t bf_id = face_bc->nhmg_dir_ids[i];
670 const short int def_id = face_bc->def_ids[bf_id];
671 const cs_xdef_t *def = eqp->bc_defs[def_id];
672 const cs_lnum_t *idx = bf2v_idx + bf_id;
673 const cs_lnum_t *lst = bf2v_lst + idx[0];
674 const int n_vf = idx[1] - idx[0];
675
676 switch(def->type) {
677
678 case CS_XDEF_BY_VALUE:
679 _assign_vb_dirichlet_values(eqp->dim, n_vf, lst,
680 (const cs_real_t *)def->context,
681 true, /* is constant for all vertices ? */
682 bcvals,
683 counter);
684 break;
685
686 case CS_XDEF_BY_ARRAY:
687 {
688 cs_real_t *eval = cb->values;
689
690 /* Evaluate the boundary condition at each boundary vertex */
691
692 cs_xdef_eval_at_vertices_by_array(n_vf,
693 lst,
694 true, /* dense output */
695 mesh,
696 connect,
697 quant,
698 t_eval,
699 def->context,
700 eval);
701
702 _assign_vb_dirichlet_values(eqp->dim, n_vf, lst,
703 eval,
704 false, /* is constant for all vertices ? */
705 bcvals,
706 counter);
707 }
708 break; /* By array */
709
710 case CS_XDEF_BY_ANALYTIC_FUNCTION:
711 {
712 cs_real_t *eval = cb->values;
713
714 /* Evaluate the boundary condition at each boundary vertex */
715
716 cs_xdef_eval_at_vertices_by_analytic(n_vf,
717 lst,
718 true, /* dense output */
719 mesh,
720 connect,
721 quant,
722 t_eval,
723 def->context,
724 eval);
725
726 _assign_vb_dirichlet_values(eqp->dim, n_vf, lst,
727 eval,
728 false, /* is constant for all vertices ? */
729 bcvals,
730 counter);
731 }
732 break;
733
734 default:
735 bft_error(__FILE__, __LINE__, 0,
736 _(" %s: Invalid type of definition.\n"
737 " Stop computing the Dirichlet value.\n"), __func__);
738
739 } /* End of switch: def_type */
740
741 } /* Loop on faces with a non-homogeneous Dirichlet BC */
742
743 cs_equation_sync_vertex_mean_values(connect, eqp->dim, counter, bcvals);
744
745 /* Homogeneous Dirichlet are always enforced (even in case of multiple BCs).
746 If multi-valued Dirichlet BCs are set, a weighted sum is used to set the
747 Dirichlet value at each corresponding vertex */
748
749 if (eqp->dim == 1) {
750
751 # pragma omp parallel if (quant->n_vertices > CS_THR_MIN)
752 {
753 # pragma omp for
754 for (cs_lnum_t v_id = 0; v_id < quant->n_vertices; v_id++) {
755 if (bcflag[v_id] & CS_CDO_BC_HMG_DIRICHLET)
756 bcvals[v_id] = 0.;
757 }
758
759 /* BC value overwrites the initial value */
760
761 # pragma omp for
762 for (cs_lnum_t v = 0; v < quant->n_vertices; v++)
763 if (cs_cdo_bc_is_dirichlet(bcflag[v]))
764 values[v] = bcvals[v];
765 }
766
767 }
768 else { /* eqp->dim > 1 */
769
770 # pragma omp parallel if (quant->n_vertices > CS_THR_MIN)
771 {
772 # pragma omp for
773 for (cs_lnum_t v_id = 0; v_id < quant->n_vertices; v_id++) {
774 if (bcflag[v_id] & CS_CDO_BC_HMG_DIRICHLET)
775 memset(bcvals + eqp->dim*v_id, 0, eqp->dim*sizeof(cs_real_t));
776 }
777
778 /* BC value overwrites the initial value */
779
780 # pragma omp for
781 for (cs_lnum_t v = 0; v < quant->n_vertices; v++) {
782 if (cs_cdo_bc_is_dirichlet(bcflag[v])) {
783 for (int k = 0; k < 3; k++)
784 values[3*v+k] = bcvals[3*v+k];
785 }
786 }
787
788 }
789
790 } /* eqp->dim ? */
791
792 /* Free temporary buffers */
793
794 BFT_FREE(counter);
795
796 #if defined(DEBUG) && !defined(NDEBUG) && CS_EQUATION_BC_DBG > 1
797 cs_dbg_darray_to_listing("DIRICHLET_VALUES",
798 eqp->dim*quant->n_vertices, bcvals, 6*eqp->dim);
799 #endif
800 }
801
802 /*----------------------------------------------------------------------------*/
803 /*!
804 * \brief Compute the values of the Dirichlet BCs when DoFs are attached to
805 * CDO face-based schemes
806 *
807 * \param[in] mesh pointer to a cs_mesh_t structure
808 * \param[in] quant pointer to a cs_cdo_quantities_t structure
809 * \param[in] connect pointer to a cs_cdo_connect_t struct.
810 * \param[in] eqp pointer to a cs_equation_param_t
811 * \param[in] face_bc pointer to a cs_cdo_bc_face_t structure
812 * \param[in] t_eval time at which one evaluates the boundary cond.
813 * \param[in, out] cb pointer to a cs_cell_builder_t structure
814 * \param[in, out] values pointer to the array of values to set
815 */
816 /*----------------------------------------------------------------------------*/
817
818 void
cs_equation_compute_dirichlet_fb(const cs_mesh_t * mesh,const cs_cdo_quantities_t * quant,const cs_cdo_connect_t * connect,const cs_equation_param_t * eqp,const cs_cdo_bc_face_t * face_bc,cs_real_t t_eval,cs_cell_builder_t * cb,cs_real_t * values)819 cs_equation_compute_dirichlet_fb(const cs_mesh_t *mesh,
820 const cs_cdo_quantities_t *quant,
821 const cs_cdo_connect_t *connect,
822 const cs_equation_param_t *eqp,
823 const cs_cdo_bc_face_t *face_bc,
824 cs_real_t t_eval,
825 cs_cell_builder_t *cb,
826 cs_real_t *values)
827 {
828 assert(face_bc != NULL);
829
830 CS_UNUSED(cb);
831
832 /* Define the array storing the Dirichlet values */
833
834 for (int def_id = 0; def_id < eqp->n_bc_defs; def_id++) {
835
836 const cs_xdef_t *def = eqp->bc_defs[def_id];
837
838 if (def->meta & CS_CDO_BC_DIRICHLET) {
839
840 assert(eqp->dim == def->dim);
841
842 const cs_zone_t *bz = cs_boundary_zone_by_id(def->z_id);
843 const cs_lnum_t *elt_ids = bz->elt_ids;
844
845 switch(def->type) {
846
847 case CS_XDEF_BY_VALUE:
848 {
849 const cs_real_t *constant_val = (cs_real_t *)def->context;
850
851 if (def->dim == 1) {
852
853 # pragma omp parallel for if (bz->n_elts > CS_THR_MIN)
854 for (cs_lnum_t i = 0; i < bz->n_elts; i++) {
855 const cs_lnum_t elt_id = (elt_ids == NULL) ? i : elt_ids[i];
856 values[elt_id] = constant_val[0];
857 }
858
859 }
860 else {
861
862 # pragma omp parallel for if (bz->n_elts > CS_THR_MIN)
863 for (cs_lnum_t i = 0; i < bz->n_elts; i++) {
864 const cs_lnum_t elt_id = (elt_ids == NULL) ? i : elt_ids[i];
865 for (int k = 0; k < def->dim; k++)
866 values[def->dim*elt_id+k] = constant_val[k];
867 }
868
869 }
870 }
871 break;
872
873 case CS_XDEF_BY_ARRAY:
874 {
875 cs_xdef_array_context_t *ac =
876 (cs_xdef_array_context_t *)def->context;
877
878 assert(ac->stride == eqp->dim);
879 assert(cs_flag_test(ac->loc, cs_flag_primal_face) ||
880 cs_flag_test(ac->loc, cs_flag_boundary_face));
881
882 if (eqp->n_bc_defs == 1) { /* Only one definition */
883
884 assert(bz->n_elts == quant->n_b_faces);
885 memcpy(values, ac->values, sizeof(cs_real_t)*bz->n_elts*eqp->dim);
886
887 }
888 else { /* Only a selection has to be considered */
889
890 assert(elt_ids != NULL);
891 # pragma omp parallel for if (bz->n_elts > CS_THR_MIN)
892 for (cs_lnum_t i = 0; i < bz->n_elts; i++) {
893 const cs_lnum_t shift = def->dim*elt_ids[i];
894 for (int k = 0; k < def->dim; k++)
895 values[shift+k] = ac->values[shift+k];
896 }
897
898 }
899 }
900 break;
901
902 case CS_XDEF_BY_ANALYTIC_FUNCTION:
903 /* Evaluate the boundary condition at each boundary face */
904 switch(eqp->dof_reduction) {
905
906 case CS_PARAM_REDUCTION_DERHAM:
907 cs_xdef_eval_at_b_faces_by_analytic(bz->n_elts,
908 bz->elt_ids,
909 false, /* dense output */
910 mesh,
911 connect,
912 quant,
913 t_eval,
914 def->context,
915 values);
916 break;
917
918 case CS_PARAM_REDUCTION_AVERAGE:
919 cs_xdef_eval_avg_at_b_faces_by_analytic(bz->n_elts,
920 bz->elt_ids,
921 false, /* dense output */
922 mesh,
923 connect,
924 quant,
925 t_eval,
926 def->context,
927 def->qtype,
928 def->dim,
929 values);
930 break;
931
932 default:
933 bft_error(__FILE__, __LINE__, 0,
934 _(" %s: Invalid type of reduction.\n"
935 " Stop computing the Dirichlet value.\n"), __func__);
936
937 } /* switch on reduction */
938 break;
939
940 case CS_XDEF_BY_DOF_FUNCTION:
941 cs_xdef_eval_at_b_faces_by_dof_func(bz->n_elts,
942 bz->elt_ids,
943 false, /* dense output */
944 mesh,
945 connect,
946 quant,
947 t_eval,
948 def->context,
949 values);
950 break;
951
952 default:
953 bft_error(__FILE__, __LINE__, 0,
954 _(" %s: Invalid type of definition.\n"
955 " Stop computing the Dirichlet value.\n"), __func__);
956
957 } /* switch on def_type */
958
959 } /* Definition based on Dirichlet BC */
960 } /* Loop on definitions */
961
962 /* Set the values to zero for all faces attached to a homogeneous
963 Diriclet BC */
964
965 # pragma omp parallel for if (quant->n_b_faces > CS_THR_MIN)
966 for (cs_lnum_t bf_id = 0; bf_id < quant->n_b_faces; bf_id++)
967 if (face_bc->flag[bf_id] & CS_CDO_BC_HMG_DIRICHLET)
968 for (int k = 0; k < eqp->dim; k++)
969 values[eqp->dim*bf_id + k] = 0;
970
971 #if defined(DEBUG) && !defined(NDEBUG) && CS_EQUATION_BC_DBG > 1
972 cs_dbg_darray_to_listing("DIRICHLET_VALUES",
973 eqp->dim*quant->n_b_faces, values, 9);
974 #endif
975 }
976
977 /*----------------------------------------------------------------------------*/
978 /*!
979 * \brief Define an array of flags for each vertex collecting the flags
980 * of associated boundary faces
981 *
982 * \param[in] connect pointer to a \ref cs_cdo_connect_t struct.
983 * \param[in] face_bc pointer to a structure collecting boundary
984 * conditions applied to faces
985 * \param[in, out] vflag BC flag on vertices to define
986 */
987 /*----------------------------------------------------------------------------*/
988
989 void
cs_equation_set_vertex_bc_flag(const cs_cdo_connect_t * connect,const cs_cdo_bc_face_t * face_bc,cs_flag_t * vflag)990 cs_equation_set_vertex_bc_flag(const cs_cdo_connect_t *connect,
991 const cs_cdo_bc_face_t *face_bc,
992 cs_flag_t *vflag)
993 {
994 if (vflag == NULL)
995 return;
996
997 assert(connect->bf2v != NULL);
998
999 const cs_adjacency_t *bf2v = connect->bf2v;
1000 const cs_lnum_t n_vertices = connect->n_vertices;
1001 const cs_lnum_t n_b_faces = connect->n_faces[CS_BND_FACES];
1002
1003 /* Initialization */
1004
1005 memset(vflag, 0, n_vertices*sizeof(cs_flag_t));
1006
1007 for (cs_lnum_t bf_id = 0; bf_id < n_b_faces; bf_id++) {
1008
1009 const cs_flag_t bc_flag = face_bc->flag[bf_id];
1010 for (cs_lnum_t j = bf2v->idx[bf_id]; j < bf2v->idx[bf_id+1]; j++) {
1011 vflag[bf2v->ids[j]] |= bc_flag;
1012 }
1013
1014 } /* Loop on border faces */
1015
1016 #if defined(DEBUG) && !defined(NDEBUG)
1017 for (cs_lnum_t bf_id = 0; bf_id < n_b_faces; bf_id++) {
1018 for (cs_lnum_t j = bf2v->idx[bf_id]; j < bf2v->idx[bf_id+1]; j++) {
1019 const cs_lnum_t v_id = bf2v->ids[j];
1020 if (vflag[v_id] == 0)
1021 bft_error(__FILE__, __LINE__, 0,
1022 "%s: Border vertices %ld without any boundary conditions.",
1023 __func__, (long)v_id);
1024 }
1025 } /* Loop on border faces */
1026 #endif
1027
1028 if (connect->interfaces[CS_CDO_CONNECT_VTX_SCAL] != NULL)
1029 cs_interface_set_inclusive_or(connect->interfaces[CS_CDO_CONNECT_VTX_SCAL],
1030 n_vertices,
1031 1, /* stride */
1032 false, /* interlace */
1033 CS_FLAG_TYPE, /* unsigned short int */
1034 vflag);
1035 }
1036
1037 /*----------------------------------------------------------------------------*/
1038 /*!
1039 * \brief Define an array of flags for each edge collecting the flags
1040 * of associated boundary faces
1041 *
1042 * \param[in] connect pointer to a \ref cs_cdo_connect_t struct.
1043 * \param[in] face_bc pointer to a structure collecting boundary
1044 * conditions applied to faces
1045 * \param[in, out] edge_flag BC flag on edges to define
1046 */
1047 /*----------------------------------------------------------------------------*/
1048
1049 void
cs_equation_set_edge_bc_flag(const cs_cdo_connect_t * connect,const cs_cdo_bc_face_t * face_bc,cs_flag_t * edge_flag)1050 cs_equation_set_edge_bc_flag(const cs_cdo_connect_t *connect,
1051 const cs_cdo_bc_face_t *face_bc,
1052 cs_flag_t *edge_flag)
1053 {
1054 if (edge_flag == NULL)
1055 return;
1056
1057 const cs_adjacency_t *f2e = connect->f2e;
1058 const cs_lnum_t n_edges = connect->n_edges;
1059 const cs_lnum_t n_i_faces = connect->n_faces[CS_INT_FACES];
1060 const cs_lnum_t n_faces = connect->n_faces[CS_ALL_FACES];
1061
1062 /* Initialization */
1063
1064 memset(edge_flag, 0, n_edges*sizeof(cs_flag_t));
1065
1066 for (cs_lnum_t bf_id = 0, f_id = n_i_faces; f_id < n_faces;
1067 f_id++, bf_id++) {
1068
1069 const cs_flag_t bc_flag = face_bc->flag[bf_id];
1070 for (cs_lnum_t j = f2e->idx[f_id]; j < f2e->idx[f_id+1]; j++) {
1071 edge_flag[f2e->ids[j]] |= bc_flag;
1072 }
1073
1074 } /* Loop on border faces */
1075
1076 #if defined(DEBUG) && !defined(NDEBUG)
1077 for (cs_lnum_t bf_id = n_i_faces; bf_id < n_faces; bf_id++) {
1078 for (cs_lnum_t j = f2e->idx[bf_id]; j < f2e->idx[bf_id+1]; j++) {
1079 const cs_lnum_t e_id = f2e->ids[j];
1080 if (edge_flag[e_id] == 0)
1081 bft_error(__FILE__, __LINE__, 0,
1082 "%s: Border edge %ld without any boundary conditions.",
1083 __func__, (long)e_id);
1084 }
1085 } /* Loop on border faces */
1086 #endif
1087
1088 if (connect->interfaces[CS_CDO_CONNECT_EDGE_SCAL] != NULL)
1089 cs_interface_set_inclusive_or(connect->interfaces[CS_CDO_CONNECT_EDGE_SCAL],
1090 n_edges,
1091 1, /* stride */
1092 false, /* interlace */
1093 CS_FLAG_TYPE, /* unsigned short int */
1094 edge_flag);
1095 }
1096
1097 /*----------------------------------------------------------------------------*/
1098 /*!
1099 * \brief Compute the values of the Neumann BCs when DoFs are scalar-valued
1100 * and attached to vertices.
1101 *
1102 * \param[in] t_eval time at which one performs the evaluation
1103 * \param[in] def_id id of the definition for setting the Neumann BC
1104 * \param[in] f local face number in the cs_cell_mesh_t
1105 * \param[in] eqp pointer to a cs_equation_param_t
1106 * \param[in] cm pointer to a cs_cell_mesh_t structure
1107 * \param[in, out] neu_values array storing the Neumann values
1108 */
1109 /*----------------------------------------------------------------------------*/
1110
1111 void
cs_equation_compute_neumann_sv(cs_real_t t_eval,short int def_id,short int f,const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,double * neu_values)1112 cs_equation_compute_neumann_sv(cs_real_t t_eval,
1113 short int def_id,
1114 short int f,
1115 const cs_equation_param_t *eqp,
1116 const cs_cell_mesh_t *cm,
1117 double *neu_values)
1118 {
1119 assert(neu_values != NULL && cm != NULL && eqp != NULL);
1120 assert(def_id > -1);
1121 assert(cs_eflag_test(cm->flag,
1122 CS_FLAG_COMP_EV | CS_FLAG_COMP_FE | CS_FLAG_COMP_FV));
1123
1124 const cs_xdef_t *def = eqp->bc_defs[def_id];
1125
1126 assert(def->dim == 3); /* flux is a vector in the scalar case */
1127 assert(def->meta & CS_CDO_BC_NEUMANN); /* Neuman BC */
1128
1129 /* Evaluate the boundary condition at each boundary vertex */
1130
1131 switch(def->type) {
1132
1133 case CS_XDEF_BY_VALUE:
1134 cs_xdef_cw_eval_flux_at_vtx_by_val(cm, f, t_eval, def->context, neu_values);
1135 break;
1136
1137 case CS_XDEF_BY_ANALYTIC_FUNCTION:
1138 cs_xdef_cw_eval_flux_at_vtx_by_analytic(cm,
1139 f,
1140 t_eval,
1141 def->context,
1142 def->qtype,
1143 neu_values);
1144 break;
1145
1146 case CS_XDEF_BY_ARRAY:
1147 {
1148 cs_xdef_array_context_t *ac
1149 = (cs_xdef_array_context_t *)def->context;
1150
1151 assert(eqp->n_bc_defs == 1); /* Only one definition allowed */
1152 assert(ac->stride == 3);
1153
1154 cs_lnum_t bf_id = cm->f_ids[f] - cm->bface_shift;
1155 assert(bf_id > -1);
1156
1157 if (cs_flag_test(ac->loc, cs_flag_primal_face) ||
1158 cs_flag_test(ac->loc, cs_flag_boundary_face))
1159 cs_xdef_cw_eval_flux_at_vtx_by_val(cm, f, t_eval,
1160 ac->values + 3*bf_id,
1161 neu_values);
1162
1163 else if (cs_flag_test(ac->loc, cs_flag_dual_closure_byf)) {
1164
1165 assert(ac->index != NULL);
1166
1167 /* Retrieve the bf2v->idx stored in the cs_cdo_connect_t structure */
1168
1169 cs_lnum_t shift = ac->index[bf_id];
1170 for (short int iv = cm->f2v_idx[f]; iv < cm->f2v_idx[f+1]; iv++)
1171 neu_values[cm->f2v_ids[iv]] = ac->values[shift++];
1172
1173 }
1174 else
1175 bft_error(__FILE__, __LINE__, 0,
1176 " %s: Invalid array location.", __func__);
1177
1178 }
1179 break;
1180
1181 default:
1182 bft_error(__FILE__, __LINE__, 0,
1183 _(" Invalid type of definition.\n"
1184 " Stop computing the Neumann value.\n"));
1185
1186 } /* switch def_type */
1187 }
1188
1189 /*----------------------------------------------------------------------------*/
1190 /*!
1191 * \brief Compute the values of the Neumann BCs when DoFs are attached to
1192 * faces.
1193 *
1194 * \param[in] t_eval time at which one performs the evaluation
1195 * \param[in] def_id id of the definition for setting the Neumann BC
1196 * \param[in] f local face number in the cs_cell_mesh_t
1197 * \param[in] eqp pointer to a cs_equation_param_t
1198 * \param[in] cm pointer to a cs_cell_mesh_t structure
1199 * \param[in, out] neu_values array storing Neumann values to use
1200 */
1201 /*----------------------------------------------------------------------------*/
1202
1203 void
cs_equation_compute_neumann_fb(cs_real_t t_eval,short int def_id,short int f,const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,double * neu_values)1204 cs_equation_compute_neumann_fb(cs_real_t t_eval,
1205 short int def_id,
1206 short int f,
1207 const cs_equation_param_t *eqp,
1208 const cs_cell_mesh_t *cm,
1209 double *neu_values)
1210 {
1211 assert(neu_values != NULL && cm != NULL && eqp != NULL);
1212 assert(def_id > -1);
1213 assert(eqp->dim == 1 || eqp->dim == 3);
1214
1215 const cs_xdef_t *def = eqp->bc_defs[def_id];
1216
1217 /* Flux is a vector in the scalar-valued case and a tensor in the
1218 vector-valued case */
1219
1220 assert(def->meta & CS_CDO_BC_NEUMANN); /* Neuman BC */
1221
1222 /* Evaluate the boundary condition at each boundary face */
1223
1224 switch(def->type) {
1225
1226 case CS_XDEF_BY_VALUE:
1227 if (eqp->dim == 1)
1228 cs_xdef_cw_eval_flux_by_val(cm, f, t_eval, def->context, neu_values);
1229 else if (eqp->dim == 3)
1230 cs_xdef_cw_eval_tensor_flux_by_val(cm, f, t_eval,
1231 def->context,
1232 neu_values);
1233 break;
1234
1235 case CS_XDEF_BY_ANALYTIC_FUNCTION:
1236 if (eqp->dim == 1)
1237 cs_xdef_cw_eval_flux_by_analytic(cm,
1238 f,
1239 t_eval,
1240 def->context,
1241 def->qtype,
1242 neu_values);
1243 else if (eqp->dim == 3)
1244 cs_xdef_cw_eval_tensor_flux_by_analytic(cm,
1245 f,
1246 t_eval,
1247 def->context,
1248 def->qtype,
1249 neu_values);
1250 break;
1251
1252 case CS_XDEF_BY_ARRAY:
1253 {
1254 cs_xdef_array_context_t *ac
1255 = (cs_xdef_array_context_t *)def->context;
1256
1257 assert(eqp->n_bc_defs == 1); /* Only one definition allowed */
1258 assert(ac->stride == 3);
1259 assert(cs_flag_test(ac->loc, cs_flag_primal_face) ||
1260 cs_flag_test(ac->loc, cs_flag_boundary_face));
1261
1262 cs_lnum_t bf_id = cm->f_ids[f] - cm->bface_shift;
1263 assert(bf_id > -1);
1264
1265 cs_real_t *face_val = ac->values + 3*bf_id;
1266
1267 cs_xdef_cw_eval_flux_by_val(cm, f, t_eval, face_val, neu_values);
1268 }
1269 break;
1270
1271 default:
1272 bft_error(__FILE__, __LINE__, 0,
1273 _(" Invalid type of definition.\n"
1274 " Stop computing the Neumann value.\n"));
1275
1276 } /* switch def_type */
1277 }
1278
1279 /*----------------------------------------------------------------------------*/
1280 /*!
1281 * \brief Compute the values of the Robin BCs
1282 *
1283 * \param[in] t_eval time at which one performs the evaluation
1284 * \param[in] def_id id of the definition for setting the Neumann BC
1285 * \param[in] f local face number in the cs_cell_mesh_t
1286 * \param[in] eqp pointer to a cs_equation_param_t
1287 * \param[in] cm pointer to a cs_cell_mesh_t structure
1288 * \param[in, out] rob_values array storing Robin values to use
1289 */
1290 /*----------------------------------------------------------------------------*/
1291
1292 void
cs_equation_compute_robin(cs_real_t t_eval,short int def_id,short int f,const cs_equation_param_t * eqp,const cs_cell_mesh_t * cm,double * rob_values)1293 cs_equation_compute_robin(cs_real_t t_eval,
1294 short int def_id,
1295 short int f,
1296 const cs_equation_param_t *eqp,
1297 const cs_cell_mesh_t *cm,
1298 double *rob_values)
1299 {
1300 assert(rob_values != NULL && cm != NULL && eqp != NULL);
1301 assert(def_id > -1);
1302 assert(eqp->dim == 1);
1303
1304 const cs_xdef_t *def = eqp->bc_defs[def_id];
1305
1306 /* Flux is a vector in the scalar-valued case and a tensor in the
1307 vector-valued case */
1308 assert(def->meta & CS_CDO_BC_ROBIN); /* Robin BC */
1309
1310 /* Evaluate the boundary condition at each boundary face */
1311 switch(def->type) {
1312
1313 case CS_XDEF_BY_VALUE:
1314 {
1315 const cs_real_t *parameters = (cs_real_t *)def->context;
1316
1317 rob_values[3*f ] = parameters[0];
1318 rob_values[3*f+1] = parameters[1];
1319 rob_values[3*f+2] = parameters[2];
1320 }
1321 break;
1322
1323 case CS_XDEF_BY_ANALYTIC_FUNCTION:
1324 {
1325 cs_real_3_t parameters = {0, 0, 0};
1326 cs_xdef_analytic_context_t *ac =
1327 (cs_xdef_analytic_context_t *)def->context;
1328
1329 ac->func(t_eval, 1, NULL,
1330 cm->face[f].center, true, /* dense output ? */
1331 ac->input,
1332 (cs_real_t *)parameters);
1333
1334 rob_values[3*f ] = parameters[0];
1335 rob_values[3*f+1] = parameters[1];
1336 rob_values[3*f+2] = parameters[2];
1337 }
1338 break;
1339
1340 case CS_XDEF_BY_ARRAY:
1341 {
1342 cs_xdef_array_context_t *c = (cs_xdef_array_context_t *)def->context;
1343
1344 assert(eqp->n_bc_defs == 1); /* Only one definition allowed */
1345 assert(c->stride == 3);
1346 assert(cs_flag_test(c->loc, cs_flag_primal_face) ||
1347 cs_flag_test(c->loc, cs_flag_boundary_face));
1348
1349 cs_lnum_t bf_id = cm->f_ids[f] - cm->bface_shift;
1350 assert(bf_id > -1);
1351 cs_real_t *parameters = c->values + 3*bf_id;
1352
1353 rob_values[3*f ] = parameters[0];
1354 rob_values[3*f+1] = parameters[1];
1355 rob_values[3*f+2] = parameters[2];
1356 }
1357 break;
1358
1359 default:
1360 bft_error(__FILE__, __LINE__, 0,
1361 _(" Invalid type of definition.\n"
1362 " Stop computing the Robin value.\n"));
1363
1364 } /* switch def_type */
1365
1366 }
1367
1368 /*----------------------------------------------------------------------------*/
1369 /*!
1370 * \brief Compute the values of the tangential component lying on the domain
1371 * boundary. Kind of BCs used when DoFs are attached to CDO (primal)
1372 * edge-based schemes. One sets the values of the circulation.
1373 *
1374 * \param[in] t_eval time at which one evaluates the boundary cond.
1375 * \param[in] mesh pointer to a cs_mesh_t structure
1376 * \param[in] quant pointer to a cs_cdo_quantities_t structure
1377 * \param[in] connect pointer to a cs_cdo_connect_t struct.
1378 * \param[in] eqp pointer to a cs_equation_param_t
1379 * \param[in, out] values pointer to the array of values to set
1380 */
1381 /*----------------------------------------------------------------------------*/
1382
1383 void
cs_equation_compute_circulation_eb(cs_real_t t_eval,const cs_mesh_t * mesh,const cs_cdo_quantities_t * quant,const cs_cdo_connect_t * connect,const cs_equation_param_t * eqp,cs_real_t * values)1384 cs_equation_compute_circulation_eb(cs_real_t t_eval,
1385 const cs_mesh_t *mesh,
1386 const cs_cdo_quantities_t *quant,
1387 const cs_cdo_connect_t *connect,
1388 const cs_equation_param_t *eqp,
1389 cs_real_t *values)
1390 {
1391 CS_UNUSED(mesh);
1392 CS_UNUSED(quant);
1393 assert(values != NULL);
1394
1395 /* Synchronization of the definition of the circulation if needed */
1396 cs_lnum_t *def2e_ids = (cs_lnum_t *)cs_equation_get_tmpbuf();
1397 cs_lnum_t *def2e_idx = NULL;
1398 BFT_MALLOC(def2e_idx, eqp->n_bc_defs + 1, cs_lnum_t);
1399
1400 _sync_circulation_def_at_edges(connect,
1401 eqp->n_bc_defs,
1402 eqp->bc_defs,
1403 def2e_idx,
1404 def2e_ids);
1405
1406 /* Define the array storing the circulation values */
1407 for (int def_id = 0; def_id < eqp->n_bc_defs; def_id++) {
1408
1409 const cs_xdef_t *def = eqp->bc_defs[def_id];
1410
1411 if ((def->meta & CS_CDO_BC_TANGENTIAL_DIRICHLET) ||
1412 (def->meta & CS_CDO_BC_DIRICHLET)) {
1413
1414 const cs_lnum_t n_elts = def2e_idx[def_id+1] - def2e_idx[def_id];
1415 const cs_lnum_t *elt_ids = def2e_ids + def2e_idx[def_id];
1416
1417 switch(def->type) {
1418
1419 case CS_XDEF_BY_VALUE:
1420 cs_evaluate_circulation_along_edges_by_value(def,
1421 n_elts,
1422 elt_ids,
1423 values);
1424 break;
1425
1426 case CS_XDEF_BY_ARRAY:
1427 cs_evaluate_circulation_along_edges_by_array(def,
1428 n_elts,
1429 elt_ids,
1430 values);
1431 break;
1432
1433 case CS_XDEF_BY_ANALYTIC_FUNCTION:
1434 cs_evaluate_circulation_along_edges_by_analytic(def,
1435 t_eval,
1436 n_elts,
1437 elt_ids,
1438 values);
1439 break;
1440
1441 default:
1442 bft_error(__FILE__, __LINE__, 0,
1443 _(" %s: Invalid type of definition.\n"
1444 " Stop computing the circulation.\n"), __func__);
1445
1446 } /* Switch on def_type */
1447
1448 } /* Definition related to a circulation */
1449
1450 } /* Loop on definitions */
1451
1452 BFT_FREE(def2e_idx);
1453
1454 #if defined(DEBUG) && !defined(NDEBUG) && CS_EQUATION_BC_DBG > 1
1455 cs_dbg_darray_to_listing("CIRCULATION_VALUES", quant->n_edges, values, 9);
1456 #endif
1457 }
1458
1459 /*----------------------------------------------------------------------------*/
1460
1461 END_C_DECLS
1462