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