1 /*============================================================================
2 * Functions to handle the reconstruction of fields
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <math.h>
35
36 /*----------------------------------------------------------------------------
37 * Local headers
38 *----------------------------------------------------------------------------*/
39
40 #include <bft_mem.h>
41
42 #include "cs_math.h"
43 #include "cs_scheme_geometry.h"
44
45 /*----------------------------------------------------------------------------
46 * Header for the current file
47 *----------------------------------------------------------------------------*/
48
49 #include "cs_reco.h"
50
51 /*----------------------------------------------------------------------------*/
52
53 BEGIN_C_DECLS
54
55 /*=============================================================================
56 * Local macro and structure definitions
57 *============================================================================*/
58
59 /* Redefined the name of functions from cs_math to get shorter names */
60 #define _dp3 cs_math_3_dot_product
61
62 /*============================================================================
63 * Public function prototypes
64 *============================================================================*/
65
66 /*----------------------------------------------------------------------------*/
67 /*!
68 * \brief Reconstruct at cell centers and face centers a vertex-based field
69 * Linear interpolation. If p_crec and/or p_frec are not allocated, this
70 * done in this subroutine.
71 *
72 * \param[in] connect pointer to the connectivity struct.
73 * \param[in] quant pointer to the additional quantities struct.
74 * \param[in] dof pointer to the field of vtx-based DoFs
75 * \param[in, out] p_crec reconstructed values at cell centers
76 * \param[in, out] p_frec reconstructed values at face centers
77 */
78 /*----------------------------------------------------------------------------*/
79
80 void
cs_reco_conf_vtx_dofs(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const double * dof,double * p_crec[],double * p_frec[])81 cs_reco_conf_vtx_dofs(const cs_cdo_connect_t *connect,
82 const cs_cdo_quantities_t *quant,
83 const double *dof,
84 double *p_crec[],
85 double *p_frec[])
86 {
87 double *crec = *p_crec, *frec = *p_frec;
88
89 const cs_adjacency_t *c2v = connect->c2v;
90 const double *dcv = quant->dcell_vol;
91 const cs_adjacency_t *f2e = connect->f2e;
92 const cs_adjacency_t *e2v = connect->e2v;
93
94 if (dof == NULL)
95 return;
96
97 /* Allocate reconstruction arrays if necessary */
98 if (crec == NULL)
99 BFT_MALLOC(crec, quant->n_cells, double);
100 if (frec == NULL)
101 BFT_MALLOC(frec, quant->n_faces, double);
102
103 /* Reconstruction at cell centers */
104 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
105
106 crec[c_id] = 0;
107 for (cs_lnum_t j = c2v->idx[c_id]; j < c2v->idx[c_id+1]; j++)
108 crec[c_id] += dcv[j]*dof[c2v->ids[j]];
109 crec[c_id] /= quant->cell_vol[c_id];
110
111 }
112
113 /* Reconstruction at face centers */
114 for (cs_lnum_t f_id = 0; f_id < quant->n_faces; f_id++) {
115
116 const cs_real_t *xf = cs_quant_get_face_center(f_id, quant);
117
118 frec[f_id] = 0;
119 double f_surf = 0.;
120 for (cs_lnum_t j = f2e->idx[f_id]; j < f2e->idx[f_id+1]; j++) {
121
122 const cs_lnum_t e_id = f2e->ids[j];
123 const cs_lnum_t v1_id = e2v->ids[2*e_id];
124 const cs_lnum_t v2_id = e2v->ids[2*e_id+1];
125 const cs_real_t *xv1 = quant->vtx_coord + 3*v1_id;
126 const cs_real_t *xv2 = quant->vtx_coord + 3*v2_id;
127
128 cs_real_3_t xe;
129 for (int k = 0; k < 3; k++)
130 xe[k] = 0.5 * (xv1[k] + xv2[k]);
131
132 double lef, lve;
133 cs_real_3_t uef, uve, cp;
134 cs_math_3_length_unitv(xe, xf, &lef, uef);
135 cs_math_3_length_unitv(xv1, xv2, &lve, uve);
136 cs_math_3_cross_product(uve, uef, cp);
137
138 const double tef = 0.5 * lve * lef * cs_math_3_norm(cp);
139
140 f_surf += tef;
141 frec[f_id] += 0.5 * tef * (dof[v1_id] + dof[v2_id]);
142
143 } /* End of loop on face edges */
144
145 frec[f_id] /= f_surf;
146
147 } /* End of loop on faces */
148
149 /* Return pointers */
150 *p_crec = crec;
151 *p_frec = frec;
152 }
153
154 /*----------------------------------------------------------------------------*/
155 /*!
156 * \brief Reconstruct the value at all cell centers from an array of values
157 * defined on primal vertices.
158 *
159 * \param[in] c2v cell -> vertices connectivity
160 * \param[in] quant pointer to the additional quantities struct.
161 * \param[in] array pointer to the array of values
162 * \param[in, out] val_xc values of the reconstruction at the cell center
163 */
164 /*----------------------------------------------------------------------------*/
165
166 void
cs_reco_pv_at_cell_centers(const cs_adjacency_t * c2v,const cs_cdo_quantities_t * quant,const double * array,cs_real_t * val_xc)167 cs_reco_pv_at_cell_centers(const cs_adjacency_t *c2v,
168 const cs_cdo_quantities_t *quant,
169 const double *array,
170 cs_real_t *val_xc)
171 {
172 if (array == NULL)
173 return;
174
175 /* Sanity checks */
176 assert(c2v != NULL && quant != NULL);
177
178 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
179 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
180
181 const double invvol = 1/quant->cell_vol[c_id];
182 const cs_real_t *dcvol = quant->dcell_vol;
183
184 cs_real_t reco_val = 0;
185 for (cs_lnum_t jv = c2v->idx[c_id]; jv < c2v->idx[c_id+1]; jv++) {
186
187 const cs_lnum_t v_id = c2v->ids[jv];
188
189 reco_val += dcvol[jv] * array[v_id];
190
191 } /* Loop on cell vertices; */
192
193 val_xc[c_id] = invvol * reco_val;
194
195 } /* Loop on cells */
196
197 }
198
199 /*----------------------------------------------------------------------------*/
200 /*!
201 * \brief Reconstruct the value at all cell centers from an array of values
202 * defined on primal vertices.
203 * Case of vector-valued fields.
204 *
205 * \param[in] c2v cell -> vertices connectivity
206 * \param[in] quant pointer to the additional quantities struct.
207 * \param[in] array pointer to the array of values
208 * \param[in, out] val_xc values of the reconstruction at the cell center
209 */
210 /*----------------------------------------------------------------------------*/
211
212 void
cs_reco_vect_pv_at_cell_centers(const cs_adjacency_t * c2v,const cs_cdo_quantities_t * quant,const double * array,cs_real_t * val_xc)213 cs_reco_vect_pv_at_cell_centers(const cs_adjacency_t *c2v,
214 const cs_cdo_quantities_t *quant,
215 const double *array,
216 cs_real_t *val_xc)
217 {
218 if (array == NULL)
219 return;
220
221 /* Sanity checks */
222 assert(c2v != NULL && quant != NULL);
223
224 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
225 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
226
227 const double invvol = 1/quant->cell_vol[c_id];
228 const cs_real_t *dcvol = quant->dcell_vol;
229
230 cs_real_t reco_val[3] = {0, 0, 0};
231 for (cs_lnum_t jv = c2v->idx[c_id]; jv < c2v->idx[c_id+1]; jv++) {
232
233 const cs_real_t *_array = array + 3*c2v->ids[jv];
234 const cs_real_t vc_vol = dcvol[jv];
235
236 reco_val[0] += vc_vol * _array[0];
237 reco_val[1] += vc_vol * _array[1];
238 reco_val[2] += vc_vol * _array[2];
239
240 } /* Loop on cell vertices */
241
242 for (int k = 0; k < 3; k++)
243 val_xc[3*c_id+k] = invvol * reco_val[k];
244
245 } /* Loop on cells */
246
247 }
248
249 /*----------------------------------------------------------------------------*/
250 /*!
251 * \brief Reconstruct the vector-valued quantity inside each cell from the face
252 * DoFs (interior and boundary). Scalar-valued face DoFs are related to
253 * the normal flux across faces.
254 *
255 * \param[in] c2f cell -> faces connectivity
256 * \param[in] quant pointer to the additional quantities struct.
257 * \param[in] i_face_vals array of DoF values for interior faces
258 * \param[in] b_face_vals array of DoF values for border faces
259 * \param[out] cell_reco vector-valued reconstruction inside cells
260 */
261 /*----------------------------------------------------------------------------*/
262
263 void
cs_reco_cell_vectors_by_ib_face_dofs(const cs_adjacency_t * c2f,const cs_cdo_quantities_t * cdoq,const cs_real_t i_face_vals[],const cs_real_t b_face_vals[],cs_real_t * cell_reco)264 cs_reco_cell_vectors_by_ib_face_dofs(const cs_adjacency_t *c2f,
265 const cs_cdo_quantities_t *cdoq,
266 const cs_real_t i_face_vals[],
267 const cs_real_t b_face_vals[],
268 cs_real_t *cell_reco)
269 {
270 assert(c2f != NULL && i_face_vals != NULL && b_face_vals != NULL);
271
272 /* Initialization */
273
274 memset(cell_reco, 0, 3*cdoq->n_cells*sizeof(cs_real_t));
275
276 # pragma omp parallel for if (cdoq->n_cells > CS_THR_MIN)
277 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++) {
278
279 cs_real_t *cval = cell_reco + 3*c_id;
280
281 /* Loop on cell faces */
282
283 for (cs_lnum_t j = c2f->idx[c_id]; j < c2f->idx[c_id+1]; j++) {
284
285 const cs_lnum_t f_id = c2f->ids[j];
286 const cs_lnum_t bf_id = f_id - cdoq->n_i_faces;
287 const cs_real_t *dedge_vect = cdoq->dedge_vector + 3*j;
288
289 if (bf_id > -1) { /* Border face */
290
291 for (int k = 0; k < 3; k++)
292 cval[k] += b_face_vals[bf_id] * dedge_vect[k];
293
294 }
295 else { /* Interior face */
296
297 for (int k = 0; k < 3; k++)
298 cval[k] += i_face_vals[f_id] * dedge_vect[k];
299
300 }
301
302 } /* Loop on cell faces */
303
304 const cs_real_t invvol = 1./cdoq->cell_vol[c_id];
305 for (int k =0; k < 3; k++) cval[k] *= invvol;
306
307 } /* Loop on cells */
308 }
309
310 /*----------------------------------------------------------------------------*/
311 /*!
312 * \brief Reconstruct the vector-valued quantity inside a cell from the face
313 * DoFs (interior and boundary). Scalar-valued face DoFs are related to
314 * the normal flux across faces.
315 *
316 * \param[in] c_id id of the cell to handle
317 * \param[in] c2f cell -> faces connectivity
318 * \param[in] quant pointer to the additional quantities struct.
319 * \param[in] face_dofs array of DoF values at faces
320 * \param[in] local_input true means that face_dofs is of size n_cell_faces
321 * \param[out] cell_reco vector-valued reconstruction inside cells. This
322 * quantity should have been allocated before calling
323 * this function
324 */
325 /*----------------------------------------------------------------------------*/
326
327 void
cs_reco_cell_vector_by_face_dofs(cs_lnum_t c_id,const cs_adjacency_t * c2f,const cs_cdo_quantities_t * cdoq,const cs_real_t face_dofs[],bool local_input,cs_real_t * cell_reco)328 cs_reco_cell_vector_by_face_dofs(cs_lnum_t c_id,
329 const cs_adjacency_t *c2f,
330 const cs_cdo_quantities_t *cdoq,
331 const cs_real_t face_dofs[],
332 bool local_input,
333 cs_real_t *cell_reco)
334 {
335 assert(c2f != NULL && cdoq != NULL && face_dofs != NULL);
336
337 /* Initialization */
338
339 cell_reco[0] = cell_reco[1] = cell_reco[2] = 0;
340
341 /* Loop on cell faces */
342
343 if (local_input) {
344
345 const cs_lnum_t s = c2f->idx[c_id], e = c2f->idx[c_id+1];
346 const cs_real_t *_dedge_vector = cdoq->dedge_vector + 3*s;
347
348 for (cs_lnum_t j = 0; j < e-s; j++) {
349
350 const cs_real_t *dedge_vect = _dedge_vector + 3*j;
351 for (int k = 0; k < 3; k++)
352 cell_reco[k] += face_dofs[j] * dedge_vect[k];
353
354 } /* Loop on cell faces */
355
356 }
357 else { /* face_dofs is accessed with f_id */
358
359 for (cs_lnum_t j = c2f->idx[c_id]; j < c2f->idx[c_id+1]; j++) {
360
361 const cs_lnum_t f_id = c2f->ids[j];
362 const cs_real_t *dedge_vect = cdoq->dedge_vector + 3*j;
363
364 for (int k = 0; k < 3; k++)
365 cell_reco[k] += face_dofs[f_id] * dedge_vect[k];
366
367 } /* Loop on cell faces */
368
369 }
370
371 const cs_real_t invvol = 1./cdoq->cell_vol[c_id];
372 for (int k =0; k < 3; k++) cell_reco[k] *= invvol;
373 }
374
375 /*----------------------------------------------------------------------------*/
376 /*!
377 * \brief Reconstruct the vector-valued quantity inside each cell from the face
378 * DoFs (interior and boundary). Scalar-valued face DoFs are related to
379 * the normal flux across faces.
380 *
381 * \param[in] c2f cell -> faces connectivity
382 * \param[in] quant pointer to the additional quantities struct.
383 * \param[in] face_dofs array of DoF values at faces
384 * \param[out] cell_reco vector-valued reconstruction inside cells. This
385 * quantity should have been allocated before calling
386 * this function
387 */
388 /*----------------------------------------------------------------------------*/
389
390 void
cs_reco_cell_vectors_by_face_dofs(const cs_adjacency_t * c2f,const cs_cdo_quantities_t * cdoq,const cs_real_t face_dofs[],cs_real_t * cell_reco)391 cs_reco_cell_vectors_by_face_dofs(const cs_adjacency_t *c2f,
392 const cs_cdo_quantities_t *cdoq,
393 const cs_real_t face_dofs[],
394 cs_real_t *cell_reco)
395 {
396 assert(c2f != NULL && cdoq != NULL && face_dofs != NULL);
397
398 /* Initialization */
399
400 memset(cell_reco, 0, 3*cdoq->n_cells*sizeof(cs_real_t));
401
402 # pragma omp parallel for if (cdoq->n_cells > CS_THR_MIN)
403 for (cs_lnum_t c_id = 0; c_id < cdoq->n_cells; c_id++) {
404
405 cs_real_t *cval = cell_reco + 3*c_id;
406
407 for (cs_lnum_t j = c2f->idx[c_id]; j < c2f->idx[c_id+1]; j++) {
408
409 const cs_lnum_t f_id = c2f->ids[j];
410 const cs_real_t *dedge_vect = cdoq->dedge_vector + 3*j;
411
412 for (int k = 0; k < 3; k++)
413 cval[k] += face_dofs[f_id] * dedge_vect[k];
414
415 } /* Loop on cell faces */
416
417 const cs_real_t invvol = 1./cdoq->cell_vol[c_id];
418 for (int k =0; k < 3; k++) cval[k] *= invvol;
419
420 } /* Loop on cells */
421 }
422
423 /*----------------------------------------------------------------------------*/
424 /*!
425 * \brief Reconstruct the value at the cell center from an array of values
426 * defined on primal vertices.
427 *
428 * \param[in] c_id cell id
429 * \param[in] c2v cell -> vertices connectivity
430 * \param[in] quant pointer to the additional quantities struct.
431 * \param[in] array pointer to the array of values
432 * \param[in, out] val_xc value of the reconstruction at the cell center
433 */
434 /*----------------------------------------------------------------------------*/
435
436 void
cs_reco_pv_at_cell_center(cs_lnum_t c_id,const cs_adjacency_t * c2v,const cs_cdo_quantities_t * quant,const double * array,cs_real_t * val_xc)437 cs_reco_pv_at_cell_center(cs_lnum_t c_id,
438 const cs_adjacency_t *c2v,
439 const cs_cdo_quantities_t *quant,
440 const double *array,
441 cs_real_t *val_xc)
442 {
443 cs_real_t reco_val = 0;
444
445 if (array == NULL) {
446 *val_xc = reco_val;
447 return;
448 }
449
450 assert(c2v != NULL && quant != NULL && c_id > -1);
451
452 const double invvol = 1/quant->cell_vol[c_id];
453 const cs_real_t *dcvol = quant->dcell_vol;
454
455 for (cs_lnum_t jv = c2v->idx[c_id]; jv < c2v->idx[c_id+1]; jv++) {
456
457 const cs_lnum_t v_id = c2v->ids[jv];
458
459 reco_val += dcvol[jv] * array[v_id];
460
461 } /* Loop on cell vertices; */
462
463 *val_xc = invvol * reco_val;
464 }
465
466 /*----------------------------------------------------------------------------*/
467 /*!
468 * \brief Reconstruct a vector-valued array at vertices from a vector-valued
469 * array at cells.
470 *
471 * \param[in] c2v cell -> vertices connectivity
472 * \param[in] quant pointer to the additional quantities struct.
473 * \param[in] val pointer to the array of values
474 * \param[in, out] reco_val values of the reconstruction at vertices
475 */
476 /*----------------------------------------------------------------------------*/
477
478 void
cs_reco_vect_pv_from_pc(const cs_adjacency_t * c2v,const cs_cdo_quantities_t * quant,const double * val,cs_real_t * reco_val)479 cs_reco_vect_pv_from_pc(const cs_adjacency_t *c2v,
480 const cs_cdo_quantities_t *quant,
481 const double *val,
482 cs_real_t *reco_val)
483 {
484 if (val == NULL || reco_val == NULL)
485 return;
486
487 memset(reco_val, 0, 3*quant->n_vertices*sizeof(cs_real_t));
488
489 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
490 for (cs_lnum_t j = c2v->idx[c_id]; j < c2v->idx[c_id+1]; j++) {
491
492 const cs_real_t vc_vol = quant->dcell_vol[j];
493 cs_real_t *rval = reco_val + 3*c2v->ids[j];
494
495 rval[0] += vc_vol * val[3*c_id];
496 rval[1] += vc_vol * val[3*c_id + 1];
497 rval[2] += vc_vol * val[3*c_id + 2];
498
499 }
500 } /* Loop on cells */
501
502 cs_real_t *dual_vol = NULL;
503 BFT_MALLOC(dual_vol, quant->n_vertices, cs_real_t);
504 cs_cdo_quantities_compute_dual_volumes(quant, c2v, dual_vol);
505
506 # pragma omp parallel for if (quant->n_vertices > CS_THR_MIN)
507 for (cs_lnum_t v_id = 0; v_id < quant->n_vertices; v_id++) {
508 const cs_real_t invvol = 1./dual_vol[v_id];
509 for (int k = 0; k < 3; k++)
510 reco_val[3*v_id+k] *= invvol;
511 }
512
513 BFT_FREE(dual_vol);
514 }
515
516 /*----------------------------------------------------------------------------*/
517 /*!
518 * \brief Reconstruct the value at the face center from an array of values
519 * defined on primal vertices.
520 *
521 * \param[in] f_id face id (interior and border faces)
522 * \param[in] connect pointer to a cs_cdo_connect_t structure
523 * \param[in] quant pointer to the additional quantities struct.
524 * \param[in] pdi pointer to the array of values
525 * \param[in, out] pdi_f value of the reconstruction at the face center
526 */
527 /*----------------------------------------------------------------------------*/
528
529 void
cs_reco_pf_from_pv(cs_lnum_t f_id,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const double * pdi,cs_real_t * pdi_f)530 cs_reco_pf_from_pv(cs_lnum_t f_id,
531 const cs_cdo_connect_t *connect,
532 const cs_cdo_quantities_t *quant,
533 const double *pdi,
534 cs_real_t *pdi_f)
535 {
536 *pdi_f = 0.;
537
538 if (pdi == NULL)
539 return;
540
541 const cs_real_t *xf = cs_quant_get_face_center(f_id, quant);
542 const cs_real_t *xyz = quant->vtx_coord;
543 const cs_adjacency_t *e2v = connect->e2v;
544 const cs_adjacency_t *f2e = connect->f2e;
545
546 double f_surf = 0.;
547 for (cs_lnum_t i = f2e->idx[f_id]; i < f2e->idx[f_id+1]; i++) {
548
549 const cs_lnum_t e_id = f2e->ids[i];
550 const cs_lnum_t shift_e = 2*e_id;
551 const cs_lnum_t v1_id = e2v->ids[shift_e];
552 const cs_lnum_t v2_id = e2v->ids[shift_e+1];
553 const double pdi_e = 0.5*(pdi[v1_id] + pdi[v2_id]);
554 const cs_real_t *xv1 = xyz + 3*v1_id;
555 const cs_real_t *xv2 = xyz + 3*v2_id;
556 const cs_real_t tef = cs_math_surftri(xv1, xv2, xf);
557
558 f_surf += tef;
559 *pdi_f += pdi_e * tef;
560
561 } /* Loop on face edges */
562
563 *pdi_f /= f_surf;
564 }
565
566 /*----------------------------------------------------------------------------*/
567 /*!
568 * \brief Reconstruct a constant vector at the cell center from an array of
569 * values defined on dual faces lying inside each cell.
570 * This array is scanned thanks to the c2e connectivity.
571 *
572 * \param[in] c_id cell id
573 * \param[in] c2e cell -> edges connectivity
574 * \param[in] quant pointer to the additional quantities struct.
575 * \param[in] array pointer to the array of values
576 * \param[in, out] val_xc value of the reconstruction at the cell center
577 */
578 /*----------------------------------------------------------------------------*/
579
580 void
cs_reco_dfbyc_at_cell_center(cs_lnum_t c_id,const cs_adjacency_t * c2e,const cs_cdo_quantities_t * quant,const double * array,cs_real_3_t val_xc)581 cs_reco_dfbyc_at_cell_center(cs_lnum_t c_id,
582 const cs_adjacency_t *c2e,
583 const cs_cdo_quantities_t *quant,
584 const double *array,
585 cs_real_3_t val_xc)
586 {
587 /* Initialization */
588
589 val_xc[0] = val_xc[1] = val_xc[2] = 0.;
590
591 if (array == NULL)
592 return;
593
594 for (cs_lnum_t j = c2e->idx[c_id]; j < c2e->idx[c_id+1]; j++) {
595
596 const cs_real_t *e_vect = quant->edge_vector + 3*c2e->ids[j];
597 for (int k = 0; k < 3; k++)
598 val_xc[k] += array[j] * e_vect[k];
599
600 } /* Loop on cell edges */
601
602 const double invvol = 1/quant->cell_vol[c_id];
603 for (int k = 0; k < 3; k++)
604 val_xc[k] *= invvol;
605 }
606
607 /*----------------------------------------------------------------------------*/
608 /*!
609 * \brief Reconstruct a constant vector inside the cell c.
610 * array is scanned thanks to the c2e connectivity. Pointer is already
611 * located at the beginning of the cell sequence.
612 * Reconstruction used is based on DGA (stabilization = 1/d where d is
613 * the space dimension)
614 *
615 * \param[in] cm pointer to a cs_cell_mesh_t structure
616 * \param[in] array local pointer to the array of values
617 * \param[in, out] val_c value of the reconstructed vector in the cell
618 */
619 /*----------------------------------------------------------------------------*/
620
621 void
cs_reco_dfbyc_in_cell(const cs_cell_mesh_t * cm,const cs_real_t * array,cs_real_3_t val_c)622 cs_reco_dfbyc_in_cell(const cs_cell_mesh_t *cm,
623 const cs_real_t *array,
624 cs_real_3_t val_c)
625 {
626 /* Initialization */
627
628 val_c[0] = val_c[1] = val_c[2] = 0.;
629
630 if (array == NULL)
631 return;
632
633 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PEQ));
634
635 const double invvol = 1/cm->vol_c;
636
637 for (short int e = 0; e < cm->n_ec; e++) {
638
639 const cs_quant_t peq = cm->edge[e];
640 const cs_real_t edge_contrib = array[e]*peq.meas;
641
642 for (int k = 0; k < 3; k++)
643 val_c[k] += edge_contrib * peq.unitv[k];
644
645 } /* Loop on cell edges */
646
647 for (int k = 0; k < 3; k++)
648 val_c[k] *= invvol;
649 }
650
651 /*----------------------------------------------------------------------------*/
652 /*!
653 * \brief Reconstruct a constant vector inside pec which is a volume
654 * surrounding the edge e inside the cell c.
655 * array is scanned thanks to the c2e connectivity. Pointer is already
656 * located at the beginning of the cell sequence.
657 * Reconstruction used is based on DGA (stabilization = 1/d where d is
658 * the space dimension)
659 *
660 * \param[in] cm pointer to a cs_cell_mesh_t structure
661 * \param[in] e local edge id
662 * \param[in] array local pointer to the array of values
663 * \param[in, out] val_pec value of the reconstruction in pec
664 */
665 /*----------------------------------------------------------------------------*/
666
667 void
cs_reco_dfbyc_in_pec(const cs_cell_mesh_t * cm,short int e,const cs_real_t * array,cs_real_3_t val_pec)668 cs_reco_dfbyc_in_pec(const cs_cell_mesh_t *cm,
669 short int e,
670 const cs_real_t *array,
671 cs_real_3_t val_pec)
672 {
673 /* Initialize values */
674
675 val_pec[0] = val_pec[1] = val_pec[2] = 0.;
676
677 if (array == NULL)
678 return;
679
680 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PEQ | CS_FLAG_COMP_DFQ));
681
682 cs_real_3_t val_c = {0., 0., 0.};
683
684 /* Compute val_c */
685
686 for (short int _e = 0; _e < cm->n_ec; _e++) {
687
688 const cs_quant_t _peq = cm->edge[_e];
689
690 for (int k = 0; k < 3; k++)
691 val_c[k] += array[_e] * _peq.meas * _peq.unitv[k];
692
693 } /* Loop on cell edges */
694
695 const double invvol = 1/cm->vol_c;
696
697 /* Compute the consistency part related to this cell */
698
699 for (int k = 0; k < 3; k++)
700 val_c[k] *= invvol;
701
702 /* Compute the reconstruction inside pec */
703
704 const cs_quant_t peq = cm->edge[e];
705 const cs_nvec3_t dfq = cm->dface[e];
706 const double ecoef = (array[e] - dfq.meas * _dp3(dfq.unitv, val_c))
707 / (dfq.meas * _dp3(dfq.unitv, peq.unitv));
708
709 for (int k = 0; k < 3; k++)
710 val_pec[k] = val_c[k] + ecoef * peq.unitv[k];
711 }
712
713 /*----------------------------------------------------------------------------*/
714 /*!
715 * \brief Reconstruct at the cell center a field of edge-based DoFs
716 *
717 * \param[in] c_id cell id
718 * \param[in] c2e cell -> edges connectivity
719 * \param[in] quant pointer to the additional quantities struct.
720 * \param[in] dof pointer to the field of edge-based DoFs
721 * \param[in, out] reco value of the reconstructed field at cell center
722 */
723 /*----------------------------------------------------------------------------*/
724
725 void
cs_reco_ccen_edge_dof(cs_lnum_t c_id,const cs_adjacency_t * c2e,const cs_cdo_quantities_t * quant,const double * dof,double reco[])726 cs_reco_ccen_edge_dof(cs_lnum_t c_id,
727 const cs_adjacency_t *c2e,
728 const cs_cdo_quantities_t *quant,
729 const double *dof,
730 double reco[])
731 {
732 if (dof == NULL)
733 return;
734
735 /* Initialize value */
736
737 reco[0] = reco[1] = reco[2] = 0.0;
738
739 const cs_lnum_t *c2e_idx = c2e->idx + c_id;
740 const cs_lnum_t *c2e_ids = c2e->ids + c2e_idx[0];
741 const cs_real_t *dface = quant->dface_normal + 3*c2e_idx[0];
742 for (cs_lnum_t i = 0; i < c2e_idx[1] - c2e_idx[0]; i++) {
743
744 /* Dual face quantities */
745
746 const double val = dof[c2e_ids[i]]; /* Edge value */
747 for (int k = 0; k < 3; k++)
748 reco[k] += val * dface[3*i+k];
749
750 } /* End of loop on cell edges */
751
752 /* Divide by the cell volume */
753
754 const double invvol = 1/quant->cell_vol[c_id];
755 for (int k = 0; k < 3; k++)
756 reco[k] *= invvol;
757 }
758
759 /*----------------------------------------------------------------------------*/
760 /*!
761 * \brief Reconstruct at each cell center a field of edge-based DoFs
762 *
763 * \param[in] connect pointer to the connectivity struct.
764 * \param[in] quant pointer to the additional quantities struct.
765 * \param[in] dof pointer to the field of edge-based DoFs
766 * \param[in, out] p_ccrec pointer to the reconstructed values
767 */
768 /*----------------------------------------------------------------------------*/
769
770 void
cs_reco_ccen_edge_dofs(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const double * dof,double * p_ccrec[])771 cs_reco_ccen_edge_dofs(const cs_cdo_connect_t *connect,
772 const cs_cdo_quantities_t *quant,
773 const double *dof,
774 double *p_ccrec[])
775 {
776 assert(connect->c2e != NULL); /* Sanity check */
777
778 double *ccrec = *p_ccrec;
779
780 if (dof == NULL)
781 return;
782
783 /* Allocate reconstructed vector field at each cell barycenter */
784
785 if (ccrec == NULL)
786 BFT_MALLOC(ccrec, 3*quant->n_cells, double);
787
788 # pragma omp parallel for if (quant->n_cells > CS_THR_MIN)
789 for (int c_id = 0; c_id < quant->n_cells; c_id++)
790 cs_reco_ccen_edge_dof(c_id,
791 connect->c2e,
792 quant,
793 dof,
794 ccrec + 3*c_id);
795
796 /* Return pointer */
797
798 *p_ccrec = ccrec;
799 }
800
801 /*----------------------------------------------------------------------------*/
802 /*!
803 * \brief Reconstruct a cell-wise constant curl from the knowledge of the
804 * circulation at primal edges
805 *
806 * \param[in] connect pointer to a cs_cdo_connect_t structure
807 * \param[in] quant pointer to the additional quantities struct.
808 * \param[in] circ pointer to the array of circulations at edges
809 * \param[in, out] p_curl pointer to value of the reconstructed curl inside
810 * cells (allocated if set to NULL)
811 */
812 /*----------------------------------------------------------------------------*/
813
814 void
cs_reco_cell_curl_by_edge_dofs(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_real_t * circ,cs_real_t ** p_curl)815 cs_reco_cell_curl_by_edge_dofs(const cs_cdo_connect_t *connect,
816 const cs_cdo_quantities_t *quant,
817 const cs_real_t *circ,
818 cs_real_t **p_curl)
819 {
820 if (circ == NULL)
821 return;
822
823 assert(connect != NULL && quant != NULL);
824
825 cs_real_t *curl_vectors = *p_curl;
826 if (curl_vectors == NULL)
827 BFT_MALLOC(curl_vectors, 3*quant->n_cells, cs_real_t);
828
829 cs_real_t *face_curl = NULL;
830 cs_cdo_connect_discrete_curl(connect, circ, &face_curl);
831
832 cs_reco_cell_vectors_by_face_dofs(connect->c2f, quant, face_curl,
833 curl_vectors);
834
835 BFT_FREE(face_curl);
836
837 /* Returns pointer */
838
839 *p_curl = curl_vectors;
840 }
841
842 /*----------------------------------------------------------------------------*/
843 /*!
844 * \brief Reconstruct the mean-value of the gradient field with DoFs arising
845 * from a face-based scheme (values at face center and cell center)
846 * The reconstruction only deals with the consistent part so that there
847 * is no distinction between Fb schemes
848 *
849 * \param[in] c_id cell id
850 * \param[in] connect pointer to a cs_cdo_connect_t structure
851 * \param[in] quant pointer to the additional quantities struct.
852 * \param[in] p_c pointer to the array of values in cells
853 * \param[in] p_f pointer to the array of values on faces
854 * \param[in, out] grd_c value of the reconstructed gradient at cell center
855 */
856 /*----------------------------------------------------------------------------*/
857
858 void
cs_reco_grad_cell_from_fb_dofs(cs_lnum_t c_id,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_real_t * p_c,const cs_real_t * p_f,cs_real_t grd_c[])859 cs_reco_grad_cell_from_fb_dofs(cs_lnum_t c_id,
860 const cs_cdo_connect_t *connect,
861 const cs_cdo_quantities_t *quant,
862 const cs_real_t *p_c,
863 const cs_real_t *p_f,
864 cs_real_t grd_c[])
865 {
866 /* Initialize the value to return */
867
868 grd_c[0] = grd_c[1] = grd_c[2] = 0.;
869
870 if (p_c == NULL || p_f == NULL)
871 return;
872 assert(connect != NULL && quant != NULL); /* sanity checks */
873
874 const cs_real_t _p_c = p_c[c_id];
875 const cs_lnum_t s = connect->c2f->idx[c_id], e = connect->c2f->idx[c_id+1];
876 const cs_lnum_t *c2f_ids = connect->c2f->ids + s;
877 const short int *c2f_sgn = connect->c2f->sgn + s;
878
879 for (cs_lnum_t i = 0; i < e-s; i++) {
880
881 const cs_lnum_t f_id = c2f_ids[i];
882 const cs_real_t *fq = cs_quant_get_face_vector_area(f_id, quant);
883 for (int k = 0; k < 3; k++)
884 grd_c[k] += c2f_sgn[i] * (p_f[f_id] - _p_c) * fq[k];
885
886 } /* Loop on cell faces */
887
888 const cs_real_t invvol = 1./quant->cell_vol[c_id];
889 for (int k = 0; k < 3; k++)
890 grd_c[k] *= invvol;
891 }
892
893 /*----------------------------------------------------------------------------*/
894 /*!
895 * \brief Reconstruct the mean-value of the tensor gradient field with DoFs
896 * arising from a face-based scheme (vector-valued at face center and
897 * cell center) The reconstruction only deals with the consistent part
898 * so that there is no distinction between Fb schemes
899 *
900 * \param[in] c_id cell id
901 * \param[in] connect pointer to a cs_cdo_connect_t structure
902 * \param[in] quant pointer to the additional quantities struct.
903 * \param[in] u_c pointer to the array of values in cells
904 * \param[in] u_f pointer to the array of values on faces
905 * \param[in, out] grd_c value of the reconstructed gradient at cell center
906 */
907 /*----------------------------------------------------------------------------*/
908
909 void
cs_reco_grad_33_cell_from_fb_dofs(cs_lnum_t c_id,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_real_t * u_c,const cs_real_t * u_f,cs_real_t grd_c[])910 cs_reco_grad_33_cell_from_fb_dofs(cs_lnum_t c_id,
911 const cs_cdo_connect_t *connect,
912 const cs_cdo_quantities_t *quant,
913 const cs_real_t *u_c,
914 const cs_real_t *u_f,
915 cs_real_t grd_c[])
916 {
917 /* Initialize the value to return */
918
919 grd_c[0] = grd_c[1] = grd_c[2] = 0.;
920 grd_c[3] = grd_c[4] = grd_c[5] = 0.;
921 grd_c[6] = grd_c[7] = grd_c[8] = 0.;
922
923 if (u_c == NULL || u_f == NULL)
924 return;
925 assert(connect != NULL && quant != NULL); /* sanity checks */
926
927 const cs_real_t *_u_c = u_c + 3*c_id;
928 const cs_lnum_t s = connect->c2f->idx[c_id], e = connect->c2f->idx[c_id+1];
929 const cs_lnum_t *c2f_ids = connect->c2f->ids + s;
930 const short int *c2f_sgn = connect->c2f->sgn + s;
931
932 for (cs_lnum_t i = 0; i < e-s; i++) {
933
934 const cs_lnum_t f_id = c2f_ids[i];
935 const cs_real_t *fq = cs_quant_get_face_vector_area(f_id, quant);
936 const cs_real_t *_u_f = u_f + 3*f_id;
937
938 for (int ki = 0; ki < 3; ki++) {
939 const cs_real_t gki = c2f_sgn[i] * (_u_f[ki] - _u_c[ki]);
940 for (int kj = 0; kj < 3; kj++)
941 grd_c[3*ki+kj] += gki * fq[kj];
942 }
943
944 } /* Loop on cell faces */
945
946 const cs_real_t invvol = 1./quant->cell_vol[c_id];
947 for (int k = 0; k < 9; k++)
948 grd_c[k] *= invvol;
949 }
950
951 /*----------------------------------------------------------------------------*/
952 /*!
953 * \brief Reconstruct the value at the cell center of the gradient of a field
954 * defined on primal vertices.
955 *
956 * \param[in] c_id cell id
957 * \param[in] connect pointer to a cs_cdo_connect_t structure
958 * \param[in] quant pointer to the additional quantities struct.
959 * \param[in] pdi pointer to the array of values
960 * \param[in, out] val_xc value of the reconstructed gradient at cell center
961 */
962 /*----------------------------------------------------------------------------*/
963
964 void
cs_reco_grad_cell_from_pv(cs_lnum_t c_id,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_real_t * pdi,cs_real_t val_xc[])965 cs_reco_grad_cell_from_pv(cs_lnum_t c_id,
966 const cs_cdo_connect_t *connect,
967 const cs_cdo_quantities_t *quant,
968 const cs_real_t *pdi,
969 cs_real_t val_xc[])
970 {
971 val_xc[0] = val_xc[1] = val_xc[2] = 0.;
972
973 if (pdi == NULL)
974 return;
975
976 const cs_adjacency_t *e2v = connect->e2v;
977 const cs_adjacency_t *c2e = connect->c2e;
978 const cs_lnum_t *c2e_idx = c2e->idx + c_id;
979 const cs_lnum_t *c2e_ids = c2e->ids + c2e_idx[0];
980 const cs_real_t *dface = quant->dface_normal + 3*c2e_idx[0];
981
982 for (cs_lnum_t i = 0; i < c2e_idx[1] - c2e_idx[0]; i++) {
983
984 const cs_lnum_t shift_e = 2*c2e_ids[i];
985 const short int sgn_v1 = e2v->sgn[shift_e];
986 const cs_real_t pv1 = pdi[e2v->ids[shift_e]];
987 const cs_real_t pv2 = pdi[e2v->ids[shift_e+1]];
988 const cs_real_t gdi_e = sgn_v1*(pv1 - pv2);
989
990 for (int k = 0; k < 3; k++)
991 val_xc[k] += gdi_e * dface[3*i+k];
992
993 } /* Loop on cell edges */
994
995 /* Divide by cell volume */
996
997 const double invvol = 1/quant->cell_vol[c_id];
998 for (int k = 0; k < 3; k++)
999 val_xc[k] *= invvol;
1000 }
1001
1002 /*----------------------------------------------------------------------------*/
1003 /*!
1004 * \brief Reconstruct the vector-valued quantity inside each cell from the
1005 * given flux array. This array is stored in the same order as cm->f_ids
1006 * Scalar-valued face DoFs are related to the normal flux across primal
1007 * faces.
1008 *
1009 * \param[in] cm pointer to a cs_cell_mesh_t structure
1010 * \param[in] fluxes array of normal fluxes on primal faces
1011 * \param[out] cell_reco vector-valued reconstruction inside the cell
1012 */
1013 /*----------------------------------------------------------------------------*/
1014
1015 void
cs_reco_cw_cell_vect_from_flux(const cs_cell_mesh_t * cm,const cs_real_t * fluxes,cs_real_t * cell_reco)1016 cs_reco_cw_cell_vect_from_flux(const cs_cell_mesh_t *cm,
1017 const cs_real_t *fluxes,
1018 cs_real_t *cell_reco)
1019 {
1020 if (fluxes == NULL)
1021 return;
1022
1023 /* Sanity checks */
1024 assert(cell_reco != NULL && cm != NULL);
1025 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ));
1026
1027 /* Initialization */
1028 cell_reco[0] = cell_reco[1] = cell_reco[2] = 0.;
1029
1030 /* Loop on cell faces */
1031 for (short int f = 0; f < cm->n_fc; f++) {
1032
1033 /* Related dual edge quantity */
1034 const cs_nvec3_t deq = cm->dedge[f];
1035 const cs_real_t coef = fluxes[f] * deq.meas;
1036
1037 cell_reco[0] += coef*deq.unitv[0];
1038 cell_reco[1] += coef*deq.unitv[1];
1039 cell_reco[2] += coef*deq.unitv[2];
1040
1041 } /* Loop on cell faces */
1042
1043 const cs_real_t invvol = 1./cm->vol_c;
1044 for (int k =0; k < 3; k++) cell_reco[k] *= invvol;
1045 }
1046
1047 /*----------------------------------------------------------------------------*/
1048 /*!
1049 * \brief Reconstruct the vector-valued quantity inside each cell from the face
1050 * DoFs (interior and boundary). Scalar-valued face DoFs are related to
1051 * the normal flux across faces.
1052 *
1053 * \param[in] cm pointer to a cs_cell_mesh_t structure
1054 * \param[in] i_face_vals array of DoF values for interior faces
1055 * \param[in] b_face_vals array of DoF values for border faces
1056 * \param[out] cell_reco vector-valued reconstruction inside the cell
1057 */
1058 /*----------------------------------------------------------------------------*/
1059
1060 void
cs_reco_cw_cell_vect_from_face_dofs(const cs_cell_mesh_t * cm,const cs_real_t i_face_vals[],const cs_real_t b_face_vals[],cs_real_t * cell_reco)1061 cs_reco_cw_cell_vect_from_face_dofs(const cs_cell_mesh_t *cm,
1062 const cs_real_t i_face_vals[],
1063 const cs_real_t b_face_vals[],
1064 cs_real_t *cell_reco)
1065 {
1066 /* Sanity checks */
1067 assert(cm != NULL && i_face_vals != NULL && b_face_vals != NULL);
1068 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ));
1069
1070 /* Initialization */
1071 cell_reco[0] = cell_reco[1] = cell_reco[2] = 0.;
1072
1073 /* Loop on cell faces */
1074 for (short int f = 0; f < cm->n_fc; f++) {
1075
1076 /* Related dual edge quantity */
1077 const cs_lnum_t f_id = cm->f_ids[f];
1078 const cs_nvec3_t deq = cm->dedge[f];
1079 const cs_real_t coef = (cm->f_ids[f] < cm->bface_shift) ?
1080 i_face_vals[f_id] * deq.meas :
1081 b_face_vals[f_id - cm->bface_shift] * deq.meas;
1082
1083 cell_reco[0] += coef*deq.unitv[0];
1084 cell_reco[1] += coef*deq.unitv[1];
1085 cell_reco[2] += coef*deq.unitv[2];
1086
1087 } /* Loop on cell faces */
1088
1089 const cs_real_t invvol = 1./cm->vol_c;
1090 for (int k =0; k < 3; k++) cell_reco[k] *= invvol;
1091 }
1092
1093 /*----------------------------------------------------------------------------*/
1094 /*!
1095 * \brief Reconstruct the value of a scalar potential at a point inside a cell
1096 * The scalar potential has DoFs located at primal vertices
1097 *
1098 * \param[in] cm pointer to a cs_cell_mesh_t structure
1099 * \param[in] pdi array of DoF values at vertices
1100 * \param[out] cell_gradient gradient inside the cell
1101 */
1102 /*----------------------------------------------------------------------------*/
1103
1104 void
cs_reco_cw_cell_grad_from_scalar_pv(const cs_cell_mesh_t * cm,const cs_real_t pdi[],cs_real_t * cell_gradient)1105 cs_reco_cw_cell_grad_from_scalar_pv(const cs_cell_mesh_t *cm,
1106 const cs_real_t pdi[],
1107 cs_real_t *cell_gradient)
1108 {
1109 /* Sanity checks */
1110 assert(cm != NULL && pdi != NULL);
1111 assert(cs_eflag_test(cm->flag,
1112 CS_FLAG_COMP_PVQ | CS_FLAG_COMP_EV | CS_FLAG_COMP_DFQ));
1113
1114 /* Reconstruct a constant gradient inside the current cell */
1115 cell_gradient[0] = cell_gradient[1] = cell_gradient[2] = 0;
1116 for (short int e = 0; e < cm->n_ec; e++) {
1117
1118 const cs_lnum_t v0 = cm->v_ids[cm->e2v_ids[2*e]];
1119 const cs_lnum_t v1 = cm->v_ids[cm->e2v_ids[2*e+1]];
1120 const cs_real_t coeff = cm->e2v_sgn[e]*(pdi[v0]-pdi[v1])*cm->dface[e].meas;
1121 for (int k = 0; k < 3; k++)
1122 cell_gradient[k] += coeff * cm->dface[e].unitv[k];
1123
1124 } /* Loop on cell edges */
1125
1126 const cs_real_t invcell = 1./cm->vol_c;
1127 for (int k = 0; k < 3; k++) cell_gradient[k] *= invcell;
1128 }
1129
1130 /*----------------------------------------------------------------------------*/
1131 /*!
1132 * \brief Reconstruct the value of a scalar potential at a point inside a cell
1133 * The scalar potential has DoFs located at primal vertices
1134 *
1135 * \param[in] cm pointer to a cs_cell_mesh_t structure
1136 * \param[in] pdi array of DoF values at vertices
1137 * \param[in] length_xcxp lenght of the segment [x_c, x_p]
1138 * \param[in] unitv_xcxp unitary vector pointed from x_c to x_p
1139 * \param[in, out] wbuf pointer to a temporary buffer
1140 *
1141 * \return the value of the reconstructed potential at the evaluation point
1142 */
1143 /*----------------------------------------------------------------------------*/
1144
1145 cs_real_t
cs_reco_cw_scalar_pv_inside_cell(const cs_cell_mesh_t * cm,const cs_real_t pdi[],const cs_real_t length_xcxp,const cs_real_t unitv_xcxp[],cs_real_t wbuf[])1146 cs_reco_cw_scalar_pv_inside_cell(const cs_cell_mesh_t *cm,
1147 const cs_real_t pdi[],
1148 const cs_real_t length_xcxp,
1149 const cs_real_t unitv_xcxp[],
1150 cs_real_t wbuf[])
1151 {
1152 /* Sanity checks */
1153 assert(cm != NULL && pdi != NULL && wbuf != NULL);
1154 assert(cs_eflag_test(cm->flag,
1155 CS_FLAG_COMP_PVQ | CS_FLAG_COMP_EV | CS_FLAG_COMP_DFQ));
1156
1157 cs_real_t *_pv = wbuf; /* Local value of the potential field */
1158
1159 /* Reconstruct the value at the cell center */
1160 cs_real_t pc = 0.;
1161 for (short int v = 0; v < cm->n_vc; v++) {
1162 _pv[v] = pdi[cm->v_ids[v]];
1163 pc += cm->wvc[v] * _pv[v];
1164 }
1165
1166 /* Reconstruct a constant gradient inside the current cell */
1167 cs_real_3_t gcell = {0., 0., 0.};
1168 for (short int e = 0; e < cm->n_ec; e++) {
1169 const cs_real_t ge =
1170 cm->e2v_sgn[e]*( _pv[cm->e2v_ids[2*e]] - _pv[cm->e2v_ids[2*e+1]]);
1171 const cs_real_t coef_e = ge * cm->dface[e].meas;
1172 for (int k = 0; k < 3; k++)
1173 gcell[k] += coef_e * cm->dface[e].unitv[k];
1174 }
1175
1176 const cs_real_t invcell = 1./cm->vol_c;
1177 for (int k = 0; k < 3; k++) gcell[k] *= invcell;
1178
1179 /* Evaluation at the given point */
1180 cs_real_t p_rec = pc;
1181 p_rec += length_xcxp * cs_math_3_dot_product(gcell, unitv_xcxp);
1182
1183 return p_rec;
1184 }
1185
1186 /*----------------------------------------------------------------------------*/
1187 /*!
1188 * \brief Compute the weighted (by volume) gradient inside a given primal
1189 * cell for the related vertices.
1190 * Use the WBS algo. for approximating the gradient.
1191 * The computation takes into account a subdivision into tetrahedra of
1192 * the current cell based on p_{ef,c}
1193 *
1194 * \param[in] cm pointer to a cs_cell_mesh_t structure
1195 * \param[in] pot values of the potential fields at vertices + cell
1196 * \param[in, out] cb auxiliary structure for computing the flux
1197 * \param[in, out] vgrd gradient at vertices inside this cell
1198 */
1199 /*----------------------------------------------------------------------------*/
1200
1201 void
cs_reco_cw_vgrd_wbs_from_pvc(const cs_cell_mesh_t * cm,const cs_real_t * pot,cs_cell_builder_t * cb,cs_real_t * vgrd)1202 cs_reco_cw_vgrd_wbs_from_pvc(const cs_cell_mesh_t *cm,
1203 const cs_real_t *pot,
1204 cs_cell_builder_t *cb,
1205 cs_real_t *vgrd)
1206 {
1207 /* Sanity checks */
1208 assert(cs_eflag_test(cm->flag,
1209 CS_FLAG_COMP_PV | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ |
1210 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV | CS_FLAG_COMP_HFQ));
1211
1212 cs_real_3_t grd_c, grd_v1, grd_v2;
1213
1214 /* Temporary buffers */
1215 cs_real_3_t *u_vc = cb->vectors;
1216 double *l_vc = cb->values;
1217
1218 const double *p_v = pot;
1219 const double p_c = pot[cm->n_vc];
1220
1221 /* Reset local fluxes */
1222 for (int i = 0; i < 3*cm->n_vc; i++) vgrd[i] = 0.;
1223
1224 /* Store segments xv --> xc for this cell */
1225 for (short int v = 0; v < cm->n_vc; v++)
1226 cs_math_3_length_unitv(cm->xc, cm->xv + 3*v, l_vc + v, u_vc[v]);
1227
1228 /* Loop on cell faces */
1229 for (short int f = 0; f < cm->n_fc; f++) {
1230
1231 const cs_quant_t pfq = cm->face[f];
1232 const cs_nvec3_t deq = cm->dedge[f];
1233
1234 /* Compute geometrical quantities for the current face:
1235 - the weighting of each triangle defined by a base e and an apex f
1236 - volume of each sub-tetrahedron pef_c
1237 - the gradient of the Lagrange function related xc in p_{f,c} */
1238 const cs_real_t ohf = -cm->f_sgn[f]/cm->hfc[f];
1239 for (int k = 0; k < 3; k++) grd_c[k] = ohf * pfq.unitv[k];
1240
1241 /* Compute the reconstructed value of the potential at p_f */
1242 double p_f = 0.;
1243 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1244
1245 const short int *_v_ids = cm->e2v_ids + 2*cm->f2e_ids[i];
1246
1247 p_f += cm->tef[i]*( p_v[_v_ids[0]] + p_v[_v_ids[1]] );
1248
1249 }
1250 p_f *= 0.5/pfq.meas;
1251
1252 const double dp_cf = p_c - p_f;
1253
1254 /* Loop on face edges to scan p_{ef,c} subvolumes */
1255 const cs_real_t hf_coef = cs_math_1ov3 * cm->hfc[f];
1256 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1257
1258 const short int *_v_ids = cm->e2v_ids + 2*cm->f2e_ids[i];
1259 const short int v1 = _v_ids[0], v2 = _v_ids[1];
1260
1261 cs_compute_grd_ve(v1, v2, deq, (const cs_real_t (*)[3])u_vc, l_vc,
1262 grd_v1, grd_v2);
1263
1264 /* Gradient of the Lagrange function related to a face.
1265 grd_f = -(grd_c + grd_v1 + grd_v2)
1266 This formula is a consequence of the Partition of the Unity.
1267 This yields the following formula for grd(Lv^conf)|_p_{ef,c} */
1268 const cs_real_t sefv_vol = 0.5 * hf_coef * cm->tef[i]; /* 0.5 pef_vol */
1269 cs_real_t _grd;
1270 for (int k = 0; k < 3; k++) {
1271 _grd = sefv_vol * ( dp_cf * grd_c[k] +
1272 (p_v[v1] - p_f) * grd_v1[k] +
1273 (p_v[v2] - p_f) * grd_v2[k]);
1274 vgrd[3*v1 + k] += _grd;
1275 vgrd[3*v2 + k] += _grd;
1276 }
1277
1278 } /* Loop on face edges */
1279
1280 } /* Loop on cell faces */
1281
1282 }
1283
1284 /*----------------------------------------------------------------------------*/
1285 /*!
1286 * \brief Compute the mean value of a gradient inside a given primal cell.
1287 * Use the WBS algo. for approximating the gradient.
1288 * The computation takes into account a subdivision into tetrahedra of
1289 * the current cell based on p_{ef,c}
1290 *
1291 * \param[in] cm pointer to a cs_cell_mesh_t structure
1292 * \param[in] pot values of the potential fields at vertices + cell
1293 * \param[in, out] cb auxiliary structure for computing the flux
1294 * \param[in, out] cgrd mean-value of the cell gradient
1295 */
1296 /*----------------------------------------------------------------------------*/
1297
1298 void
cs_reco_cw_cgrd_wbs_from_pvc(const cs_cell_mesh_t * cm,const cs_real_t * pot,cs_cell_builder_t * cb,cs_real_t * cgrd)1299 cs_reco_cw_cgrd_wbs_from_pvc(const cs_cell_mesh_t *cm,
1300 const cs_real_t *pot,
1301 cs_cell_builder_t *cb,
1302 cs_real_t *cgrd)
1303 {
1304 /* Sanity checks */
1305 assert(cs_eflag_test(cm->flag,
1306 CS_FLAG_COMP_PV | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_DEQ |
1307 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV | CS_FLAG_COMP_HFQ));
1308
1309 cs_real_3_t grd_c, grd_v1, grd_v2;
1310
1311 /* Temporary buffers */
1312 cs_real_3_t *u_vc = cb->vectors;
1313 double *l_vc = cb->values;
1314
1315 cgrd[0] = cgrd[1] = cgrd[2] = 0.;
1316
1317 const double *p_v = pot;
1318 const double p_c = pot[cm->n_vc];
1319
1320 /* Store segments xv --> xc for this cell */
1321 for (short int v = 0; v < cm->n_vc; v++)
1322 cs_math_3_length_unitv(cm->xc, cm->xv + 3*v, l_vc + v, u_vc[v]);
1323
1324 /* Loop on cell faces */
1325 for (short int f = 0; f < cm->n_fc; f++) {
1326
1327 const cs_quant_t pfq = cm->face[f];
1328 const cs_nvec3_t deq = cm->dedge[f];
1329
1330 /* Compute geometrical quantities for the current face:
1331 - the weighting of each triangle defined by a base e and an apex f
1332 - volume of each sub-tetrahedron pef_c
1333 - the gradient of the Lagrange function related xc in p_{f,c} */
1334 const cs_real_t ohf = -cm->f_sgn[f]/cm->hfc[f];
1335 for (int k = 0; k < 3; k++) grd_c[k] = ohf * pfq.unitv[k];
1336
1337 /* Compute the reconstructed value of the potential at p_f */
1338 double p_f = 0.;
1339 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1340
1341 const short int ee = 2*cm->f2e_ids[i];
1342
1343 p_f += cm->tef[i]*( p_v[cm->e2v_ids[ee]] /* p_v1 */
1344 + p_v[cm->e2v_ids[ee+1]] ); /* p_v2 */
1345 }
1346 p_f *= 0.5/pfq.meas;
1347
1348 const double dp_cf = p_c - p_f;
1349
1350 /* Loop on face edges to scan p_{ef,c} subvolumes */
1351 const cs_real_t hf_coef = cs_math_1ov3 * cm->hfc[f];
1352 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1353
1354 const short int ee = 2*cm->f2e_ids[i];
1355 const short int v1 = cm->e2v_ids[ee];
1356 const short int v2 = cm->e2v_ids[ee+1];
1357
1358 cs_compute_grd_ve(v1, v2, deq, (const cs_real_t (*)[3])u_vc, l_vc,
1359 grd_v1, grd_v2);
1360
1361 /* Gradient of the Lagrange function related to a face.
1362 grd_f = -(grd_c + grd_v1 + grd_v2)
1363 This formula is a consequence of the Partition of the Unity.
1364 This yields the following formula for grd(Lv^conf)|_p_{ef,c} */
1365 const cs_real_t pefv_vol = hf_coef * cm->tef[i];
1366 for (int k = 0; k < 3; k++)
1367 cgrd[k] += pefv_vol * ( dp_cf * grd_c[k] +
1368 (p_v[v1] - p_f) * grd_v1[k] +
1369 (p_v[v2] - p_f) * grd_v2[k]);
1370
1371 } /* Loop on face edges */
1372
1373 } /* Loop on cell faces */
1374
1375 const double invvol = 1/cm->vol_c;
1376 for (int k = 0; k < 3; k++) cgrd[k] *= invvol;
1377 }
1378
1379 /*----------------------------------------------------------------------------*/
1380
1381 #undef _dp3
1382
1383 END_C_DECLS
1384