1 /*============================================================================
2 * Manage the (generic) cellwise evaluation of extended definitions
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <ctype.h>
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <string.h>
38 #include <math.h>
39
40 /*----------------------------------------------------------------------------
41 * Local headers
42 *----------------------------------------------------------------------------*/
43
44 #include <bft_mem.h>
45
46 #include "cs_defs.h"
47 #include "cs_field.h"
48 #include "cs_mesh_location.h"
49 #include "cs_reco.h"
50
51 /*----------------------------------------------------------------------------
52 * Header for the current file
53 *----------------------------------------------------------------------------*/
54
55 #include "cs_xdef_cw_eval.h"
56
57 /*----------------------------------------------------------------------------*/
58
59 BEGIN_C_DECLS
60
61 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
62
63 /*=============================================================================
64 * Local macro definition (unset at the end of file)
65 *============================================================================*/
66
67 /* Redefined the name of functions from cs_math to get shorter names */
68 #define _dp3 cs_math_3_dot_product
69
70 /*=============================================================================
71 * Local variables
72 *============================================================================*/
73
74 static const char _err_empty_array[] =
75 " %s: Array storing the evaluation should be allocated before the call"
76 " to this function.";
77
78 /*============================================================================
79 * Private function prototypes
80 *============================================================================*/
81
82 /*! \endcond DOXYGEN_SHOULD_SKIP_THIS */
83
84 /*============================================================================
85 * Public function prototypes
86 *============================================================================*/
87
88 /*----------------------------------------------------------------------------*/
89 /*!
90 * \brief Integrate an analytic function over a face
91 *
92 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
93 * \param[in] t_eval time at which the function is evaluated
94 * \param[in] f face id in the local cell numbering
95 * \param[in] ana analytic function to integrate
96 * \param[in] input pointer to an input structure
97 * \param[in] qfunc quadrature function to use
98 * \param[in, out] eval result of the evaluation
99 */
100 /*----------------------------------------------------------------------------*/
101
102 void
cs_xdef_cw_eval_f_int_by_analytic(const cs_cell_mesh_t * cm,double t_eval,short int f,cs_analytic_func_t * ana,void * input,cs_quadrature_tria_integral_t * qfunc,cs_real_t * eval)103 cs_xdef_cw_eval_f_int_by_analytic(const cs_cell_mesh_t *cm,
104 double t_eval,
105 short int f,
106 cs_analytic_func_t *ana,
107 void *input,
108 cs_quadrature_tria_integral_t *qfunc,
109 cs_real_t *eval)
110 {
111 const cs_quant_t pfq = cm->face[f];
112 const int start = cm->f2e_idx[f];
113 const int end = cm->f2e_idx[f+1];
114 const short int n_vf = end - start; /* #vertices (=#edges) */
115 const short int *f2e_ids = cm->f2e_ids + start;
116
117 switch (n_vf) {
118 case CS_TRIANGLE_CASE:
119 {
120 short int v0, v1, v2;
121 cs_cell_mesh_get_next_3_vertices(f2e_ids, cm->e2v_ids, &v0, &v1, &v2);
122
123 qfunc(t_eval, cm->xv+3*v0, cm->xv+3*v1, cm->xv+3*v2, pfq.meas,
124 ana, input, eval);
125 }
126 break;
127
128 default:
129 {
130 const double *tef = cm->tef + start;
131 for (short int e = 0; e < n_vf; e++) { /* Loop on face edges */
132
133 /* Edge-related variables */
134 const short int e0 = f2e_ids[e];
135 const double *xv0 = cm->xv + 3*cm->e2v_ids[2*e0];
136 const double *xv1 = cm->xv + 3*cm->e2v_ids[2*e0+1];
137
138 qfunc(t_eval, xv0, xv1, pfq.center, tef[e], ana, input, eval);
139 }
140 }
141
142 } /* Switch */
143 }
144
145 /*----------------------------------------------------------------------------*/
146 /*!
147 * \brief Integrate an analytic function over a cell
148 *
149 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
150 * \param[in] t_eval time at which the function is evaluated
151 * \param[in] ana analytic function to integrate
152 * \param[in] input pointer to an input structure
153 * \param[in] qfunc quadrature function to use
154 * \param[in, out] eval result of the evaluation
155 */
156 /*----------------------------------------------------------------------------*/
157
158 void
cs_xdef_cw_eval_c_int_by_analytic(const cs_cell_mesh_t * cm,double t_eval,cs_analytic_func_t * ana,void * input,cs_quadrature_tetra_integral_t * qfunc,cs_real_t * eval)159 cs_xdef_cw_eval_c_int_by_analytic(const cs_cell_mesh_t *cm,
160 double t_eval,
161 cs_analytic_func_t *ana,
162 void *input,
163 cs_quadrature_tetra_integral_t *qfunc,
164 cs_real_t *eval)
165 {
166 switch (cm->type) {
167
168 case FVM_CELL_TETRA:
169 {
170 assert(cm->n_fc == 4 && cm->n_vc == 4);
171 qfunc(t_eval, cm->xv, cm->xv+3, cm->xv+6, cm->xv+9, cm->vol_c,
172 ana, input, eval);
173 }
174 break;
175
176 case FVM_CELL_PYRAM:
177 case FVM_CELL_PRISM:
178 case FVM_CELL_HEXA:
179 case FVM_CELL_POLY:
180 {
181 for (short int f = 0; f < cm->n_fc; ++f) {
182
183 const cs_quant_t pfq = cm->face[f];
184 const double hf_coef = cs_math_1ov3 * cm->hfc[f];
185 const int start = cm->f2e_idx[f];
186 const int end = cm->f2e_idx[f+1];
187 const short int n_vf = end - start; /* #vertices (=#edges) */
188 const short int *f2e_ids = cm->f2e_ids + start;
189
190 assert(n_vf > 2);
191 switch(n_vf){
192
193 case CS_TRIANGLE_CASE: /* Optimized version, no subdivision */
194 {
195 short int v0, v1, v2;
196 cs_cell_mesh_get_next_3_vertices(f2e_ids, cm->e2v_ids,
197 &v0, &v1, &v2);
198
199 const double *xv0 = cm->xv + 3*v0;
200 const double *xv1 = cm->xv + 3*v1;
201 const double *xv2 = cm->xv + 3*v2;
202
203 qfunc(t_eval, xv0, xv1, xv2, cm->xc, hf_coef * pfq.meas,
204 ana, input, eval);
205 }
206 break;
207
208 default:
209 {
210 const double *tef = cm->tef + start;
211
212 for (short int e = 0; e < n_vf; e++) { /* Loop on face edges */
213
214 /* Edge-related variables */
215 const short int e0 = f2e_ids[e];
216 const double *xv0 = cm->xv + 3*cm->e2v_ids[2*e0];
217 const double *xv1 = cm->xv + 3*cm->e2v_ids[2*e0+1];
218
219 qfunc(t_eval, xv0, xv1, pfq.center, cm->xc, hf_coef * tef[e],
220 ana, input, eval);
221 }
222 }
223 break;
224
225 } /* End of switch */
226 } /* End of loop on faces */
227
228 }
229 break;
230
231 default:
232 bft_error(__FILE__, __LINE__, 0, _(" Unknown cell-type.\n"));
233 break;
234
235 } /* End of switch on the cell-type */
236 }
237
238 /*----------------------------------------------------------------------------*/
239 /*!
240 * \brief Routine to integrate an analytic function over a cell and its faces
241 *
242 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
243 * \param[in] t_eval physical time at which one evaluates the term
244 * \param[in] ana analytic function to integrate
245 * \param[in] input pointer to an input structure
246 * \param[in] dim dimension of the function
247 * \param[in] q_tet quadrature function to use on tetrahedra
248 * \param[in] q_tri quadrature function to use on triangles
249 * \param[out] c_int result of the evaluation on the cell
250 * \param[out] f_int result of the evaluation on the faces
251 */
252 /*----------------------------------------------------------------------------*/
253
254 void
cs_xdef_cw_eval_fc_int_by_analytic(const cs_cell_mesh_t * cm,cs_real_t t_eval,cs_analytic_func_t * ana,void * input,const short int dim,cs_quadrature_tetra_integral_t * q_tet,cs_quadrature_tria_integral_t * q_tri,cs_real_t * c_int,cs_real_t * f_int)255 cs_xdef_cw_eval_fc_int_by_analytic(const cs_cell_mesh_t *cm,
256 cs_real_t t_eval,
257 cs_analytic_func_t *ana,
258 void *input,
259 const short int dim,
260 cs_quadrature_tetra_integral_t *q_tet,
261 cs_quadrature_tria_integral_t *q_tri,
262 cs_real_t *c_int,
263 cs_real_t *f_int)
264 {
265 short int v0, v1, v2;
266
267 const short int nf = cm->n_fc;
268
269 /* Switch on the cell-type: optimised version for tetra */
270 switch (cm->type) {
271
272 case FVM_CELL_TETRA:
273 {
274 assert(cm->n_fc == 4 && cm->n_vc == 4);
275 q_tet(t_eval, cm->xv, cm->xv+3, cm->xv+6, cm->xv+9, cm->vol_c,
276 ana, input, c_int);
277
278 for (short int f = 0; f < nf; ++f) {
279
280 const cs_quant_t pfq = cm->face[f];
281 const short int *f2e_ids = cm->f2e_ids + cm->f2e_idx[f];
282
283 cs_cell_mesh_get_next_3_vertices(f2e_ids, cm->e2v_ids, &v0, &v1, &v2);
284 q_tri(t_eval, cm->xv+3*v0, cm->xv+3*v1, cm->xv+3*v2, pfq.meas,
285 ana, input, f_int + dim*f);
286 }
287 }
288 break;
289
290 case FVM_CELL_PYRAM:
291 case FVM_CELL_PRISM:
292 case FVM_CELL_HEXA:
293 case FVM_CELL_POLY:
294 {
295 for (short int f = 0; f < nf; ++f) {
296
297 const cs_quant_t pfq = cm->face[f];
298 const double hf_coef = cs_math_1ov3 * cm->hfc[f];
299 const int start = cm->f2e_idx[f];
300 const int end = cm->f2e_idx[f+1];
301 const short int n_vf = end - start; /* #vertices (=#edges) */
302 const short int *f2e_ids = cm->f2e_ids + start;
303
304 assert(n_vf > 2);
305 switch(n_vf){
306
307 case CS_TRIANGLE_CASE: /* triangle (optimized version, no subdivision) */
308 {
309 cs_cell_mesh_get_next_3_vertices(f2e_ids, cm->e2v_ids, &v0, &v1, &v2);
310
311 const double *xv0 = cm->xv + 3*v0;
312 const double *xv1 = cm->xv + 3*v1;
313 const double *xv2 = cm->xv + 3*v2;
314
315 q_tet(t_eval, xv0, xv1, xv2, cm->xc, hf_coef * pfq.meas,
316 ana, input, c_int);
317 q_tri(t_eval, cm->xv+3*v0, cm->xv+3*v1, cm->xv+3*v2, pfq.meas,
318 ana, input, f_int + dim*f);
319 }
320 break;
321
322 default:
323 {
324 const double *tef = cm->tef + start;
325
326 for (short int e = 0; e < n_vf; e++) { /* Loop on face edges */
327
328 /* Edge-related variables */
329 const short int e0 = f2e_ids[e];
330 const double *xv0 = cm->xv + 3*cm->e2v_ids[2*e0];
331 const double *xv1 = cm->xv + 3*cm->e2v_ids[2*e0+1];
332
333 q_tet(t_eval, xv0, xv1, pfq.center, cm->xc, hf_coef*tef[e],
334 ana, input, c_int);
335 q_tri(t_eval, xv0, xv1, pfq.center, tef[e],
336 ana, input, f_int + dim*f);
337
338 }
339 }
340 break;
341
342 } /* End of switch */
343
344 } /* End of loop on faces */
345
346 }
347 break;
348
349 default:
350 bft_error(__FILE__, __LINE__, 0, _(" Unknown cell-type.\n"));
351 break;
352
353 } /* End of switch on the cell-type */
354 }
355
356 /*----------------------------------------------------------------------------*/
357 /*!
358 * \brief Function pointer for evaluating the average on a face of a scalar
359 * function defined through a descriptor (\ref cs_xdef_t structure) by
360 * a cellwise process (usage of a \ref cs_cell_mesh_t structure)
361 *
362 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
363 * \param[in] f local face id
364 * \param[in] time_eval physical time at which one evaluates the term
365 * \param[in] context pointer to a context structure
366 * \param[in] qtype level of quadrature to use
367 * \param[in, out] eval result of the evaluation
368 */
369 /*----------------------------------------------------------------------------*/
370
371 void
cs_xdef_cw_eval_scalar_face_avg_by_analytic(const cs_cell_mesh_t * cm,short int f,cs_real_t time_eval,void * context,cs_quadrature_type_t qtype,cs_real_t * eval)372 cs_xdef_cw_eval_scalar_face_avg_by_analytic(const cs_cell_mesh_t *cm,
373 short int f,
374 cs_real_t time_eval,
375 void *context,
376 cs_quadrature_type_t qtype,
377 cs_real_t *eval)
378 {
379 if (eval == NULL)
380 bft_error(__FILE__, __LINE__, 0, _err_empty_array, __func__);
381
382 assert(context != NULL);
383 assert(cs_eflag_test(cm->flag,
384 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
385 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
386
387 cs_quadrature_tria_integral_t
388 *qfunc = cs_quadrature_get_tria_integral(1, qtype);
389 cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
390
391 cs_xdef_cw_eval_f_int_by_analytic(cm, time_eval, f,
392 ac->func, ac->input, qfunc, eval);
393
394 /* Average */
395 eval[0] /= cm->face[f].meas;
396 }
397
398 /*----------------------------------------------------------------------------*/
399 /*!
400 * \brief Function pointer for evaluating the average on a face of a vector
401 * function defined through a descriptor (\ref cs_xdef_t structure) by
402 * a cellwise process (usage of a \ref cs_cell_mesh_t structure)
403 *
404 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
405 * \param[in] f local face id
406 * \param[in] t_eval physical time at which one evaluates the term
407 * \param[in] context pointer to a context structure
408 * \param[in] qtype level of quadrature to use
409 * \param[in, out] eval result of the evaluation
410 */
411 /*----------------------------------------------------------------------------*/
412
413 void
cs_xdef_cw_eval_vector_face_avg_by_analytic(const cs_cell_mesh_t * cm,short int f,cs_real_t t_eval,void * context,cs_quadrature_type_t qtype,cs_real_t * eval)414 cs_xdef_cw_eval_vector_face_avg_by_analytic(const cs_cell_mesh_t *cm,
415 short int f,
416 cs_real_t t_eval,
417 void *context,
418 cs_quadrature_type_t qtype,
419 cs_real_t *eval)
420 {
421 if (eval == NULL)
422 bft_error(__FILE__, __LINE__, 0, _err_empty_array, __func__);
423
424 assert(context != NULL);
425 assert(cs_eflag_test(cm->flag,
426 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
427 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
428
429 cs_quadrature_tria_integral_t
430 *qfunc = cs_quadrature_get_tria_integral(3, qtype);
431 cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
432
433 cs_xdef_cw_eval_f_int_by_analytic(cm, t_eval, f,
434 ac->func, ac->input, qfunc, eval);
435
436 /* Average */
437 const double _oversurf = 1./cm->face[f].meas;
438 for (short int xyz = 0; xyz < 3; xyz++)
439 eval[xyz] *= _oversurf;
440 }
441
442 /*----------------------------------------------------------------------------*/
443 /*!
444 * \brief Function pointer for evaluating the average on a face of a tensor
445 * function defined through a descriptor (\ref cs_xdef_t structure) by
446 * a cellwise process (usage of a \ref cs_cell_mesh_t structure)
447 *
448 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
449 * \param[in] f local face id
450 * \param[in] t_eval physical time at which one evaluates the term
451 * \param[in] context pointer to a context structure
452 * \param[in] qtype level of quadrature to use
453 * \param[in, out] eval result of the evaluation
454 */
455 /*----------------------------------------------------------------------------*/
456
457 void
cs_xdef_cw_eval_tensor_face_avg_by_analytic(const cs_cell_mesh_t * cm,short int f,cs_real_t t_eval,void * context,cs_quadrature_type_t qtype,cs_real_t * eval)458 cs_xdef_cw_eval_tensor_face_avg_by_analytic(const cs_cell_mesh_t *cm,
459 short int f,
460 cs_real_t t_eval,
461 void *context,
462 cs_quadrature_type_t qtype,
463 cs_real_t *eval)
464 {
465 if (eval == NULL)
466 bft_error(__FILE__, __LINE__, 0, _err_empty_array, __func__);
467
468 assert(context != NULL);
469 assert(cs_eflag_test(cm->flag,
470 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
471 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
472
473 cs_quadrature_tria_integral_t
474 *qfunc = cs_quadrature_get_tria_integral(9, qtype);;
475 cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
476
477 cs_xdef_cw_eval_f_int_by_analytic(cm, t_eval, f,
478 ac->func, ac->input, qfunc, eval);
479
480 /* Average */
481 const double _oversurf = 1./cm->face[f].meas;
482 for (short int i = 0; i < 9; i++)
483 eval[i] *= _oversurf;
484 }
485
486 /*----------------------------------------------------------------------------*/
487 /*!
488 * \brief Function pointer for evaluating a quantity defined through a
489 * descriptor (\ref cs_xdef_t structure) by a cellwise process (usage
490 * of a \ref cs_cell_mesh_t structure).
491 * This evaluation hinges on the computation of integrals
492 *
493 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
494 * \param[in] t_eval physical time at which one evaluates the term
495 * \param[in] qtype quadrature type
496 * \param[in] context pointer to a context structure
497 * \param[in, out] eval result of the evaluation
498 */
499 /*----------------------------------------------------------------------------*/
500
501 void
cs_xdef_cw_eval_scalar_avg_by_analytic(const cs_cell_mesh_t * cm,cs_real_t t_eval,void * context,cs_quadrature_type_t qtype,cs_real_t * eval)502 cs_xdef_cw_eval_scalar_avg_by_analytic(const cs_cell_mesh_t *cm,
503 cs_real_t t_eval,
504 void *context,
505 cs_quadrature_type_t qtype,
506 cs_real_t *eval)
507 {
508 if (eval == NULL)
509 bft_error(__FILE__, __LINE__, 0, _err_empty_array, __func__);
510
511 assert(context != NULL);
512 assert(cs_eflag_test(cm->flag,
513 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
514 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
515
516 cs_quadrature_tetra_integral_t
517 *qfunc = cs_quadrature_get_tetra_integral(1, qtype);
518 cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
519
520 cs_xdef_cw_eval_c_int_by_analytic(cm, t_eval,
521 ac->func, ac->input, qfunc,
522 eval);
523
524 /* Average */
525 eval[0] /= cm->vol_c;
526 }
527
528 /*----------------------------------------------------------------------------*/
529 /*!
530 * \brief Function pointer for evaluating a quantity defined through a
531 * descriptor (\ref cs_xdef_t structure) by a cellwise process (usage
532 * of a \ref cs_cell_mesh_t structure).
533 * This evaluation hinges on the computation of integrals
534 * Vector-valued case.
535 *
536 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
537 * \param[in] t_eval physical time at which one evaluates the term
538 * \param[in] qtype quadrature type
539 * \param[in] context pointer to a context structure
540 * \param[in, out] eval result of the evaluation
541 */
542 /*----------------------------------------------------------------------------*/
543
544 void
cs_xdef_cw_eval_vector_avg_by_analytic(const cs_cell_mesh_t * cm,cs_real_t t_eval,void * context,cs_quadrature_type_t qtype,cs_real_t * eval)545 cs_xdef_cw_eval_vector_avg_by_analytic(const cs_cell_mesh_t *cm,
546 cs_real_t t_eval,
547 void *context,
548 cs_quadrature_type_t qtype,
549 cs_real_t *eval)
550 {
551 if (eval == NULL)
552 bft_error(__FILE__, __LINE__, 0, _err_empty_array, __func__);
553
554 assert(context != NULL);
555 assert(cs_eflag_test(cm->flag,
556 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
557 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
558
559 cs_quadrature_tetra_integral_t
560 *qfunc = cs_quadrature_get_tetra_integral(3, qtype);
561 cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
562
563 cs_xdef_cw_eval_c_int_by_analytic(cm, t_eval,
564 ac->func, ac->input,
565 qfunc, eval);
566
567 /* Average */
568 const double _overvol = 1./cm->vol_c;
569 for (short int xyz = 0; xyz < 3; xyz++)
570 eval[xyz] *= _overvol;
571 }
572
573 /*----------------------------------------------------------------------------*/
574 /*!
575 * \brief Function pointer for evaluating a quantity defined through a
576 * descriptor (\ref cs_xdef_t structure) by a cellwise process (usage
577 * of a \ref cs_cell_mesh_t structure).
578 * This evaluation hinges on the computation of integrals
579 * Tensor-valued case.
580 *
581 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
582 * \param[in] t_eval physical time at which one evaluates the term
583 * \param[in] qtype quadrature type
584 * \param[in] context pointer to a context structure
585 * \param[in, out] eval result of the evaluation
586 */
587 /*----------------------------------------------------------------------------*/
588
589 void
cs_xdef_cw_eval_tensor_avg_by_analytic(const cs_cell_mesh_t * cm,cs_real_t t_eval,void * context,cs_quadrature_type_t qtype,cs_real_t * eval)590 cs_xdef_cw_eval_tensor_avg_by_analytic(const cs_cell_mesh_t *cm,
591 cs_real_t t_eval,
592 void *context,
593 cs_quadrature_type_t qtype,
594 cs_real_t *eval)
595 {
596 if (eval == NULL)
597 bft_error(__FILE__, __LINE__, 0, _err_empty_array, __func__);
598
599 assert(context != NULL);
600 assert(cs_eflag_test(cm->flag,
601 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
602 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
603
604 cs_quadrature_tetra_integral_t
605 *qfunc = cs_quadrature_get_tetra_integral(9, qtype);
606 cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
607
608 cs_xdef_cw_eval_c_int_by_analytic(cm, t_eval,
609 ac->func, ac->input, qfunc,
610 eval);
611
612 /* Average */
613 const double _overvol = 1./cm->vol_c;
614 for (short int xyz = 0; xyz < 9; xyz++)
615 eval[xyz] *= _overvol;
616 }
617
618 /*----------------------------------------------------------------------------*/
619 /*!
620 * \brief Evaluate a quantity by a cellwise process using a definition by
621 * time function
622 *
623 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
624 * \param[in] time_eval physical time at which one evaluates the term
625 * \param[in] context pointer to a context structure
626 * \param[out] eval result of the evaluation
627 */
628 /*----------------------------------------------------------------------------*/
629
630 void
cs_xdef_cw_eval_by_time_func(const cs_cell_mesh_t * cm,cs_real_t time_eval,void * context,cs_real_t * eval)631 cs_xdef_cw_eval_by_time_func(const cs_cell_mesh_t *cm,
632 cs_real_t time_eval,
633 void *context,
634 cs_real_t *eval)
635 {
636 CS_UNUSED(cm);
637 cs_xdef_time_func_context_t *tfc = (cs_xdef_time_func_context_t *)context;
638
639 /* Evaluate the quantity */
640 tfc->func(time_eval, tfc->input, eval);
641 }
642
643 /*----------------------------------------------------------------------------*/
644 /*!
645 * \brief Evaluate a quantity defined using an analytic function by a
646 * cellwise process (usage of a \ref cs_cell_mesh_t structure)
647 *
648 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
649 * \param[in] time_eval physical time at which one evaluates the term
650 * \param[in] context pointer to a context structure
651 * \param[in, out] eval result of the evaluation at cell center
652 */
653 /*----------------------------------------------------------------------------*/
654
655 void
cs_xdef_cw_eval_by_analytic(const cs_cell_mesh_t * cm,cs_real_t time_eval,void * context,cs_real_t * eval)656 cs_xdef_cw_eval_by_analytic(const cs_cell_mesh_t *cm,
657 cs_real_t time_eval,
658 void *context,
659 cs_real_t *eval)
660 {
661 cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
662
663 /* Evaluate the function for this time at the cell center */
664 ac->func(time_eval,
665 1, NULL, cm->xc, true, /* compacted output ? */
666 ac->input,
667 eval);
668 }
669
670 /*----------------------------------------------------------------------------*/
671 /*!
672 * \brief Evaluate a quantity at cells defined by an array.
673 * Array is assumed to be interlaced.
674 * Variation using a \ref cs_cell_mesh_t structure
675 *
676 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
677 * \param[in] time_eval physical time at which one evaluates the term
678 * \param[in] context pointer to a context structure
679 * \param[in, out] eval result of the evaluation at cell center
680 */
681 /*----------------------------------------------------------------------------*/
682
683 void
cs_xdef_cw_eval_by_array(const cs_cell_mesh_t * cm,cs_real_t time_eval,void * context,cs_real_t * eval)684 cs_xdef_cw_eval_by_array(const cs_cell_mesh_t *cm,
685 cs_real_t time_eval,
686 void *context,
687 cs_real_t *eval)
688 {
689 CS_UNUSED(time_eval);
690
691 cs_xdef_array_context_t *ac = (cs_xdef_array_context_t *)context;
692
693 const int stride = ac->stride;
694
695 /* array is assumed to be interlaced */
696 if (cs_flag_test(ac->loc, cs_flag_primal_cell)) {
697
698 for (int k = 0; k < stride; k++)
699 eval[k] = ac->values[stride*cm->c_id + k];
700
701 }
702 else if (cs_flag_test(ac->loc, cs_flag_primal_vtx)) {
703
704 /* Sanity checks */
705 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PVQ));
706
707 /* Reconstruct (or interpolate) value at the current cell center */
708 for (short int v = 0; v < cm->n_vc; v++) {
709 for (int k = 0; k < stride; k++)
710 eval[k] += cm->wvc[v] * ac->values[stride*cm->v_ids[v] + k];
711 }
712
713 }
714 else if (cs_flag_test(ac->loc, cs_flag_dual_face_byc)) {
715
716 assert(ac->index != NULL);
717
718 /* Reconstruct (or interpolate) value at the current cell center */
719 cs_reco_dfbyc_in_cell(cm,
720 ac->values + ac->index[cm->c_id],
721 eval);
722
723 }
724 else
725 bft_error(__FILE__, __LINE__, 0,
726 " %s: Invalid support for the input array", __func__);
727
728 }
729
730 /*----------------------------------------------------------------------------*/
731 /*!
732 * \brief Evaluate a quantity inside a cell defined using a field
733 * Variation using a \ref cs_cell_mesh_t structure
734 *
735 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
736 * \param[in] time_eval physical time at which one evaluates the term
737 * \param[in] context pointer to a context structure
738 * \param[in, out] eval value of the property at the cell center
739 */
740 /*----------------------------------------------------------------------------*/
741
742 void
cs_xdef_cw_eval_by_field(const cs_cell_mesh_t * cm,cs_real_t time_eval,void * context,cs_real_t * eval)743 cs_xdef_cw_eval_by_field(const cs_cell_mesh_t *cm,
744 cs_real_t time_eval,
745 void *context,
746 cs_real_t *eval)
747 {
748 CS_UNUSED(time_eval);
749
750 cs_field_t *field = (cs_field_t *)context;
751 assert(field != NULL);
752 cs_real_t *values = field->val;
753 const int c_ml_id = cs_mesh_location_get_id_by_name(N_("cells"));
754 const int v_ml_id = cs_mesh_location_get_id_by_name(N_("vertices"));
755
756 if (field->location_id == c_ml_id) {
757
758 for (int k = 0; k < field->dim; k++)
759 eval[k] = values[field->dim*cm->c_id + k];
760
761 }
762 else if (field->location_id == v_ml_id) {
763
764 /* Sanity checks */
765 assert(field->dim == 1);
766 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PVQ));
767
768 /* Reconstruct (or interpolate) value at the current cell center */
769 for (short int v = 0; v < cm->n_vc; v++)
770 eval[0] += cm->wvc[v] * values[cm->v_ids[v]];
771
772
773 }
774 else
775 bft_error(__FILE__, __LINE__, 0,
776 " %s: Invalid support for the input array", __func__);
777
778 }
779
780 /*----------------------------------------------------------------------------*/
781 /*!
782 * \brief Function pointer for evaluating a quantity defined by analytic
783 * function at a precise location (x, y, z) inside a cell
784 * Use of a \ref cs_cell_mesh_t structure.
785 *
786 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
787 * \param[in] n_points number of points where to compute the evaluation
788 * \param[in] xyz where to compute the evaluation
789 * \param[in] time_eval physical time at which one evaluates the term
790 * \param[in] context pointer to a context structure
791 * \param[in, out] eval result of the evaluation
792 */
793 /*----------------------------------------------------------------------------*/
794
795 void
cs_xdef_cw_eval_at_xyz_by_analytic(const cs_cell_mesh_t * cm,cs_lnum_t n_points,const cs_real_t * xyz,cs_real_t time_eval,void * context,cs_real_t * eval)796 cs_xdef_cw_eval_at_xyz_by_analytic(const cs_cell_mesh_t *cm,
797 cs_lnum_t n_points,
798 const cs_real_t *xyz,
799 cs_real_t time_eval,
800 void *context,
801 cs_real_t *eval)
802 {
803 CS_UNUSED(cm);
804
805 cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
806
807 /* Evaluate the function for this time at the given coordinates */
808 ac->func(time_eval,
809 n_points, NULL, xyz, true, /* compacted output ? */
810 ac->input,
811 eval);
812 }
813
814 /*----------------------------------------------------------------------------*/
815 /*!
816 * \brief Function pointer for evaluating a quantity defined by analytic
817 * function at a precise location inside a cell
818 * Use of a \ref cs_cell_mesh_t structure.
819 * Vector-valued case.
820 *
821 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
822 * \param[in] n_points number of points where to compute the evaluation
823 * \param[in] xyz where to compute the evaluation
824 * \param[in] time_eval physical time at which one evaluates the term
825 * \param[in] context pointer to a context structure
826 * \param[in, out] eval result of the evaluation
827 */
828 /*----------------------------------------------------------------------------*/
829
830 void
cs_xdef_cw_eval_vector_at_xyz_by_array(const cs_cell_mesh_t * cm,cs_lnum_t n_points,const cs_real_t * xyz,cs_real_t time_eval,void * context,cs_real_t * eval)831 cs_xdef_cw_eval_vector_at_xyz_by_array(const cs_cell_mesh_t *cm,
832 cs_lnum_t n_points,
833 const cs_real_t *xyz,
834 cs_real_t time_eval,
835 void *context,
836 cs_real_t *eval)
837 {
838 CS_UNUSED(xyz);
839 CS_UNUSED(time_eval);
840
841 cs_xdef_array_context_t *ac = (cs_xdef_array_context_t *)context;
842
843 const int stride = ac->stride;
844
845 /* array is assumed to be interlaced */
846 if (cs_flag_test(ac->loc, cs_flag_primal_cell)) {
847
848 assert(stride == 3);
849 cs_real_3_t cell_vector;
850 for (int k = 0; k < stride; k++)
851 cell_vector[k] = ac->values[stride*cm->c_id + k];
852 for (int i = 0; i < n_points; i++) {
853 eval[3*i ] = cell_vector[0];
854 eval[3*i + 1] = cell_vector[1];
855 eval[2*i + 2] = cell_vector[2];
856 }
857
858 }
859 else if (cs_flag_test(ac->loc, cs_flag_primal_vtx)) {
860
861 /* Sanity checks */
862 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PVQ));
863 assert(stride == 3);
864
865 /* Reconstruct (or interpolate) value at the current cell center */
866 for (int k = 0; k < stride; k++) {
867 for (short int v = 0; v < cm->n_vc; v++)
868 eval[k] += cm->wvc[v] * ac->values[stride*cm->v_ids[v] + k];
869 }
870
871 }
872 else if (cs_flag_test(ac->loc, cs_flag_dual_face_byc)) {
873
874 assert(ac->index != NULL);
875
876 /* Reconstruct (or interpolate) value at the current cell center */
877 cs_real_3_t cell_vector;
878 cs_reco_dfbyc_in_cell(cm,
879 ac->values + ac->index[cm->c_id],
880 cell_vector);
881
882 for (int i = 0; i < n_points; i++) {
883 eval[3*i ] = cell_vector[0];
884 eval[3*i + 1] = cell_vector[1];
885 eval[3*i + 2] = cell_vector[2];
886 }
887
888 }
889 else
890 bft_error(__FILE__, __LINE__, 0,
891 " %s: Invalid support for the input array", __func__);
892
893 }
894
895 /*----------------------------------------------------------------------------*/
896 /*!
897 * \brief Function pointer for evaluating a quantity defined by a field
898 * at a precise location inside a cell
899 * Use of a \ref cs_cell_mesh_t structure.
900 * Vector-valued case.
901 *
902 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
903 * \param[in] n_points number of points where to compute the evaluation
904 * \param[in] xyz where to compute the evaluation
905 * \param[in] time_eval physical time at which one evaluates the term
906 * \param[in] context pointer to a context structure
907 * \param[in, out] eval result of the evaluation
908 */
909 /*----------------------------------------------------------------------------*/
910
911 void
cs_xdef_cw_eval_vector_at_xyz_by_field(const cs_cell_mesh_t * cm,cs_lnum_t n_points,const cs_real_t * xyz,cs_real_t time_eval,void * context,cs_real_t * eval)912 cs_xdef_cw_eval_vector_at_xyz_by_field(const cs_cell_mesh_t *cm,
913 cs_lnum_t n_points,
914 const cs_real_t *xyz,
915 cs_real_t time_eval,
916 void *context,
917 cs_real_t *eval)
918 {
919 CS_UNUSED(xyz);
920 CS_UNUSED(time_eval);
921
922 cs_field_t *field = (cs_field_t *)context;
923 const cs_real_t *values = field->val;
924
925 assert(field != NULL);
926 assert(field->dim == 3);
927
928 const int c_ml_id = cs_mesh_location_get_id_by_name(N_("cells"));
929 const int v_ml_id = cs_mesh_location_get_id_by_name(N_("vertices"));
930
931 /* array is assumed to be interlaced */
932 if (field->location_id == c_ml_id) {
933
934 cs_real_3_t cell_vector;
935 for (int k = 0; k < 3; k++)
936 cell_vector[k] = values[3*cm->c_id + k];
937 for (int i = 0; i < n_points; i++) { /* No interpolation */
938 eval[3*i ] = cell_vector[0];
939 eval[3*i + 1] = cell_vector[1];
940 eval[3*i + 2] = cell_vector[2];
941 }
942
943 }
944 else if (field->location_id == v_ml_id) {
945
946 /* Sanity check */
947 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_PVQ));
948
949 /* Reconstruct (or interpolate) value at the current cell center */
950 for (int k = 0; k < 3; k++)
951 for (short int v = 0; v < cm->n_vc; v++)
952 eval[k] += cm->wvc[v] * values[3*cm->v_ids[v] + k];
953
954 }
955 else
956 bft_error(__FILE__, __LINE__, 0,
957 " %s: Invalid support for the input array", __func__);
958
959 }
960
961 /*----------------------------------------------------------------------------*/
962 /*!
963 * \brief Function pointer for evaluating the normal flux of a quantity
964 * defined by values. The normal flux is then added to each portion of
965 * face related to a vertex.
966 * Use of a \ref cs_cell_mesh_t structure.
967 *
968 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
969 * \param[in] f local face id
970 * \param[in] time_eval physical time at which one evaluates the term
971 * \param[in] context pointer to a context structure
972 * \param[in, out] eval result of the evaluation (updated inside)
973 */
974 /*----------------------------------------------------------------------------*/
975
976 void
cs_xdef_cw_eval_flux_at_vtx_by_val(const cs_cell_mesh_t * cm,short int f,cs_real_t time_eval,void * context,cs_real_t * eval)977 cs_xdef_cw_eval_flux_at_vtx_by_val(const cs_cell_mesh_t *cm,
978 short int f,
979 cs_real_t time_eval,
980 void *context,
981 cs_real_t *eval)
982 {
983 CS_UNUSED(time_eval);
984 assert(cs_eflag_test(cm->flag, CS_FLAG_COMP_EV | CS_FLAG_COMP_FE));
985
986 const cs_real_t *flux = (cs_real_t *)context;
987
988 if (cs_eflag_test(cm->flag, CS_FLAG_COMP_FEQ)) {
989
990 /* Loop on face edges */
991 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
992
993 const double _flx = 0.5 * cm->tef[i] * _dp3(flux, cm->face[f].unitv);
994 const short int ee = 2*cm->f2e_ids[i];
995
996 eval[cm->e2v_ids[ee ]] += _flx;
997 eval[cm->e2v_ids[ee+1]] += _flx;
998
999 }
1000
1001 }
1002 else {
1003
1004 /* Loop on face edges */
1005 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1006
1007 const short int e = cm->f2e_ids[i];
1008 const double tef = cs_compute_area_from_quant(cm->edge[e],
1009 cm->face[f].center);
1010 const double _flx = 0.5 * tef * _dp3(flux, cm->face[f].unitv);
1011
1012 eval[cm->e2v_ids[2*e ]] += _flx;
1013 eval[cm->e2v_ids[2*e+1]] += _flx;
1014
1015 }
1016
1017 }
1018
1019 }
1020
1021 /*----------------------------------------------------------------------------*/
1022 /*!
1023 * \brief Function pointer for evaluating the normal flux of a quantity
1024 * defined by analytic function. The normal flux is then added to each
1025 * portion of face related to a vertex.
1026 * Use of a \ref cs_cell_mesh_t structure.
1027 *
1028 * \param[in] cm pointer to a cs_cell_mesh_t structure
1029 * \param[in] f local face id
1030 * \param[in] time_eval physical time at which one evaluates the term
1031 * \param[in] context pointer to a context structure
1032 * \param[in] qtype level of quadrature to use
1033 * \param[in, out] eval result of the evaluation (updated inside)
1034 */
1035 /*----------------------------------------------------------------------------*/
1036
1037 void
cs_xdef_cw_eval_flux_at_vtx_by_analytic(const cs_cell_mesh_t * cm,short int f,cs_real_t time_eval,void * context,cs_quadrature_type_t qtype,cs_real_t * eval)1038 cs_xdef_cw_eval_flux_at_vtx_by_analytic(const cs_cell_mesh_t *cm,
1039 short int f,
1040 cs_real_t time_eval,
1041 void *context,
1042 cs_quadrature_type_t qtype,
1043 cs_real_t *eval)
1044 {
1045 assert(cs_eflag_test(cm->flag,
1046 CS_FLAG_COMP_PFQ | CS_FLAG_COMP_EV | CS_FLAG_COMP_FE));
1047
1048 const cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
1049 const cs_quant_t fq = cm->face[f];
1050
1051 switch (qtype) {
1052
1053 case CS_QUADRATURE_NONE:
1054 case CS_QUADRATURE_BARY:
1055 {
1056 cs_real_3_t flux_xf = {0, 0, 0};
1057
1058 /* Evaluate the function for this time at the given coordinates */
1059 ac->func(time_eval, 1, NULL, fq.center, true, /* compacted output ? */
1060 ac->input,
1061 flux_xf);
1062
1063 /* Plug into the evaluation by value now */
1064 cs_xdef_cw_eval_flux_at_vtx_by_val(cm, f, time_eval, flux_xf, eval);
1065 }
1066 break;
1067
1068 case CS_QUADRATURE_BARY_SUBDIV:
1069 {
1070 assert(cs_flag_test(cm->flag, CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ));
1071 cs_real_3_t _val[2], _xyz[2];
1072
1073 if (cs_flag_test(cm->flag, CS_FLAG_COMP_FEQ)) {
1074
1075 /* Loop on face edges */
1076 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1077
1078 const short int e = cm->f2e_ids[i];
1079 const short int v1 = cm->e2v_ids[2*e];
1080 const short int v2 = cm->e2v_ids[2*e+1];
1081
1082 for (int k = 0; k < 3; k++) {
1083 const double xef = cm->edge[e].center[k] + fq.center[k];
1084 _xyz[0][k] = cs_math_1ov3 * (xef + cm->xv[3*v1+k]);
1085 _xyz[1][k] = cs_math_1ov3 * (xef + cm->xv[3*v2+k]);
1086 }
1087
1088 /* Evaluate the function for this time at the given coordinates */
1089 ac->func(time_eval, 2, NULL,
1090 (const cs_real_t *)_xyz, true, /* compacted output ? */
1091 ac->input,
1092 (cs_real_t *)_val);
1093
1094 eval[v1] += 0.5*cm->tef[i] * _dp3(_val[0], fq.unitv);
1095 eval[v2] += 0.5*cm->tef[i] * _dp3(_val[1], fq.unitv);
1096
1097 }
1098
1099 }
1100 else {
1101
1102 /* Loop on face edges */
1103 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1104
1105 const short int e = cm->f2e_ids[i];
1106 const short int v1 = cm->e2v_ids[2*e];
1107 const short int v2 = cm->e2v_ids[2*e+1];
1108
1109 for (int k = 0; k < 3; k++) {
1110 const double xef = cm->edge[e].center[k] + fq.center[k];
1111 _xyz[0][k] = cs_math_1ov3 * (xef + cm->xv[3*v1+k]);
1112 _xyz[1][k] = cs_math_1ov3 * (xef + cm->xv[3*v2+k]);
1113 }
1114
1115 /* Evaluate the function for this time at the given coordinates */
1116 ac->func(time_eval, 2, NULL,
1117 (const cs_real_t *)_xyz, true, /* compacted output ? */
1118 ac->input,
1119 (cs_real_t *)_val);
1120
1121 const double tef = cs_compute_area_from_quant(cm->edge[e], fq.center);
1122
1123 eval[v1] += 0.5 * tef * _dp3(_val[0], fq.unitv);
1124 eval[v2] += 0.5 * tef * _dp3(_val[1], fq.unitv);
1125
1126 }
1127
1128 }
1129
1130 }
1131 break; /* BARY_SUBDIV */
1132
1133 case CS_QUADRATURE_HIGHER:
1134 {
1135 assert(cs_flag_test(cm->flag, CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ));
1136
1137 /* Two triangles s_{vef} related to a vertex and four values by triangle
1138 * --> 2*4 = 8 Gauss points
1139 * The flux returns by the analytic function is a vector. So the size
1140 * of _val is 24 = 8*3
1141 */
1142 cs_real_t _val[24], w[8];
1143 cs_real_3_t gpts[8];
1144
1145 if (cs_flag_test(cm->flag, CS_FLAG_COMP_FEQ)) { /* tef is pre-computed */
1146
1147 /* Loop on face edges */
1148 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1149
1150 const short int e = cm->f2e_ids[i];
1151 const short int v1 = cm->e2v_ids[2*e];
1152 const short int v2 = cm->e2v_ids[2*e+1];
1153 const cs_real_t svef = 0.5 * cm->tef[i];
1154
1155 /* Two triangles composing the portion of face related to a vertex
1156 Evaluate the field at the four quadrature points */
1157 cs_quadrature_tria_4pts(cm->edge[e].center, fq.center, cm->xv + 3*v1,
1158 svef,
1159 gpts, w);
1160
1161 cs_quadrature_tria_4pts(cm->edge[e].center, fq.center, cm->xv + 3*v2,
1162 svef,
1163 gpts + 4, w + 4);
1164
1165 /* Evaluate the function for this time at the given coordinates */
1166 ac->func(time_eval, 8, NULL,
1167 (const cs_real_t *)gpts, true, /* compacted output ? */
1168 ac->input,
1169 _val);
1170
1171 cs_real_t add0 = 0, add1 = 0;
1172 for (int p = 0; p < 4; p++) {
1173 add0 += w[ p] * _dp3(_val + 3*p, fq.unitv);
1174 add1 += w[4 + p] * _dp3(_val + 3*(4+p), fq.unitv);
1175 }
1176
1177 eval[v1] += add0;
1178 eval[v2] += add1;
1179
1180 }
1181
1182 }
1183 else {
1184
1185 /* Loop on face edges */
1186 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1187
1188 const short int e = cm->f2e_ids[i];
1189 const short int v1 = cm->e2v_ids[2*e];
1190 const short int v2 = cm->e2v_ids[2*e+1];
1191 const double svef = 0.5 * cs_compute_area_from_quant(cm->edge[e],
1192 fq.center);
1193
1194 /* Two triangles composing the portion of face related to a vertex
1195 Evaluate the field at the four quadrature points */
1196 cs_quadrature_tria_4pts(cm->edge[e].center, fq.center, cm->xv + 3*v1,
1197 svef,
1198 gpts, w);
1199
1200 cs_quadrature_tria_4pts(cm->edge[e].center, fq.center, cm->xv + 3*v2,
1201 svef,
1202 gpts + 4, w + 4);
1203
1204 /* Evaluate the function for this time at the given coordinates */
1205 ac->func(time_eval, 8, NULL,
1206 (const cs_real_t *)gpts, true, /* compacted output ? */
1207 ac->input,
1208 _val);
1209
1210 cs_real_t add0 = 0, add1 = 0;
1211 for (int p = 0; p < 4; p++) {
1212 add0 += w[ p] * _dp3(_val + 3*p, fq.unitv);
1213 add1 += w[4 + p] * _dp3(_val + 3*(4+p), fq.unitv);
1214 }
1215
1216 eval[v1] += add0;
1217 eval[v2] += add1;
1218
1219 }
1220
1221 } /* Is tef already computed ? */
1222
1223 }
1224 break;
1225
1226 case CS_QUADRATURE_HIGHEST:
1227 {
1228 assert(cs_flag_test(cm->flag, CS_FLAG_COMP_PV | CS_FLAG_COMP_PEQ));
1229
1230 /* Two triangles s_{vef} related to a vertex and seven values by triangle
1231 * --> 2*7 = 14 Gauss points
1232 * The flux returns by the analytic function is a vector. So the size
1233 * of _val is 42 = 14*3
1234 */
1235 cs_real_t _val[42], w[14];
1236 cs_real_3_t gpts[14];
1237
1238 if (cs_flag_test(cm->flag, CS_FLAG_COMP_FEQ)) { /* tef is pre-computed */
1239
1240 /* Loop on face edges */
1241 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1242
1243 const short int e = cm->f2e_ids[i];
1244 const short int v1 = cm->e2v_ids[2*e];
1245 const short int v2 = cm->e2v_ids[2*e+1];
1246 const cs_real_t svef = 0.5 * cm->tef[i];
1247
1248 /* Two triangles composing the portion of face related to a vertex
1249 Evaluate the field at the seven quadrature points */
1250 cs_quadrature_tria_7pts(cm->edge[e].center, fq.center, cm->xv + 3*v1,
1251 svef,
1252 gpts, w);
1253
1254 cs_quadrature_tria_7pts(cm->edge[e].center, fq.center, cm->xv + 3*v2,
1255 svef,
1256 gpts + 7, w + 7);
1257
1258 /* Evaluate the function for this time at the given coordinates */
1259 ac->func(time_eval, 14, NULL,
1260 (const cs_real_t *)gpts, true, /* compacted output ? */
1261 ac->input,
1262 _val);
1263
1264 cs_real_t add0 = 0, add1 = 0;
1265 for (int p = 0; p < 7; p++) {
1266 add0 += w[ p] * _dp3(_val + 3*p, fq.unitv);
1267 add1 += w[7 + p] * _dp3(_val + 3*(7+p), fq.unitv);
1268 }
1269
1270 eval[v1] += add0;
1271 eval[v2] += add1;
1272
1273 }
1274
1275 }
1276 else {
1277
1278 /* Loop on face edges */
1279 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1280
1281 const short int e = cm->f2e_ids[i];
1282 const short int v1 = cm->e2v_ids[2*e];
1283 const short int v2 = cm->e2v_ids[2*e+1];
1284 const double svef = 0.5 * cs_compute_area_from_quant(cm->edge[e],
1285 fq.center);
1286
1287 /* Two triangles composing the portion of face related to a vertex
1288 Evaluate the field at the seven quadrature points */
1289 cs_quadrature_tria_7pts(cm->edge[e].center, fq.center, cm->xv + 3*v1,
1290 svef,
1291 gpts, w);
1292
1293 cs_quadrature_tria_7pts(cm->edge[e].center, fq.center, cm->xv + 3*v2,
1294 svef,
1295 gpts + 7, w + 7);
1296
1297 /* Evaluate the function for this time at the given coordinates */
1298 ac->func(time_eval, 14, NULL,
1299 (const cs_real_t *)gpts, true, /* compacted output ? */
1300 ac->input,
1301 _val);
1302
1303 cs_real_t add0 = 0, add1 = 0;
1304 for (int p = 0; p < 7; p++) {
1305 add0 += w[ p] * _dp3(_val + 3*p, fq.unitv);
1306 add1 += w[7 + p] * _dp3(_val + 3*(7+p), fq.unitv);
1307 }
1308
1309 eval[v1] += add0;
1310 eval[v2] += add1;
1311
1312 }
1313
1314 } /* Is tef already computed ? */
1315
1316 }
1317 break;
1318
1319 default:
1320 bft_error(__FILE__, __LINE__, 0,
1321 " %s: Invalid type of quadrature.", __func__);
1322 break;
1323
1324 } /* switch type of quadrature */
1325
1326 }
1327
1328 /*----------------------------------------------------------------------------*/
1329 /*!
1330 * \brief Function pointer for evaluating the normal flux of a quantity
1331 * defined by analytic function.
1332 * Use of a \ref cs_cell_mesh_t structure.
1333 *
1334 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1335 * \param[in] f local face id
1336 * \param[in] time_eval physical time at which one evaluates the term
1337 * \param[in] context pointer to a context structure
1338 * \param[in] qtype level of quadrature to use
1339 * \param[in, out] eval result of the evaluation (set inside)
1340 */
1341 /*----------------------------------------------------------------------------*/
1342
1343 void
cs_xdef_cw_eval_flux_by_analytic(const cs_cell_mesh_t * cm,short int f,cs_real_t time_eval,void * context,cs_quadrature_type_t qtype,cs_real_t * eval)1344 cs_xdef_cw_eval_flux_by_analytic(const cs_cell_mesh_t *cm,
1345 short int f,
1346 cs_real_t time_eval,
1347 void *context,
1348 cs_quadrature_type_t qtype,
1349 cs_real_t *eval)
1350 {
1351 cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
1352
1353 switch (qtype) {
1354
1355 case CS_QUADRATURE_NONE:
1356 case CS_QUADRATURE_BARY:
1357 {
1358 cs_real_3_t flux_xf = {0, 0, 0};
1359
1360 /* Evaluate the function for this time at the given coordinates */
1361 ac->func(time_eval, 1, NULL, cm->face[f].center,
1362 true, /* compacted output ? */
1363 ac->input,
1364 flux_xf);
1365
1366 /* Plug into the evaluation by value now */
1367 cs_xdef_cw_eval_flux_by_val(cm, f, time_eval, flux_xf, eval);
1368 }
1369 break;
1370
1371 case CS_QUADRATURE_BARY_SUBDIV:
1372 {
1373 assert(cs_flag_test(cm->flag,
1374 CS_FLAG_COMP_EV |CS_FLAG_COMP_FE |CS_FLAG_COMP_FEQ));
1375
1376 const cs_quant_t fq = cm->face[f];
1377
1378 cs_real_3_t _val, _xyz;
1379
1380 eval[f] = 0.; /* Reset value */
1381
1382 /* Loop on face edges */
1383 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1384
1385 const short int e = cm->f2e_ids[i];
1386 const short int v1 = cm->e2v_ids[2*e];
1387 const short int v2 = cm->e2v_ids[2*e+1];
1388
1389 for (int k = 0; k < 3; k++)
1390 _xyz[k] = cs_math_1ov3 *
1391 (fq.center[k] + cm->xv[3*v1+k] + cm->xv[3*v2+k]);
1392
1393 /* Evaluate the function for this time at the given coordinates */
1394 ac->func(time_eval, 1, NULL,
1395 (const cs_real_t *)_xyz, true, /* compacted output ? */
1396 ac->input,
1397 (cs_real_t *)_val);
1398
1399 eval[f] += cm->tef[i] * _dp3(_val, fq.unitv);
1400
1401 }
1402
1403 }
1404 break; /* BARY_SUBDIV */
1405
1406 case CS_QUADRATURE_HIGHER:
1407 {
1408 assert(cs_flag_test(cm->flag,
1409 CS_FLAG_COMP_EV |CS_FLAG_COMP_FE | CS_FLAG_COMP_FEQ));
1410
1411 /* Four values by triangle --> 4 Gauss points
1412 * The flux returns by the analytic function is a vector. So the
1413 * size of _val is 12=4*3 */
1414 cs_real_t w[4], _val[12];
1415 cs_real_3_t gpts[4];
1416
1417 const cs_quant_t fq = cm->face[f];
1418
1419 eval[f] = 0.; /* Reset value */
1420
1421 /* Loop on face edges */
1422 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1423
1424 const short int e = cm->f2e_ids[i];
1425 const short int v1 = cm->e2v_ids[2*e];
1426 const short int v2 = cm->e2v_ids[2*e+1];
1427
1428 /* Evaluate the field at the three quadrature points */
1429 cs_quadrature_tria_4pts(fq.center, cm->xv + 3*v1, cm->xv + 3*v2,
1430 cm->tef[i],
1431 gpts, w);
1432
1433 /* Evaluate the function for this time at the given coordinates */
1434 ac->func(time_eval, 4, NULL,
1435 (const cs_real_t *)gpts, true, /* compacted output ? */
1436 ac->input,
1437 _val);
1438
1439 cs_real_t add = 0;
1440 for (int p = 0; p < 4; p++)
1441 add += w[p] * _dp3(_val+3*p, fq.unitv);
1442
1443 eval[f] += add;
1444
1445 } /* Loop on face edges */
1446
1447 }
1448 break;
1449
1450 case CS_QUADRATURE_HIGHEST:
1451 {
1452 assert(cs_flag_test(cm->flag,
1453 CS_FLAG_COMP_EV |CS_FLAG_COMP_FE | CS_FLAG_COMP_FEQ));
1454
1455 /* Seven values by triangle --> 7 Gauss points
1456 * The flux returns by the analytic function is a vector. So the
1457 * size of _val is 21=7*3
1458 */
1459 cs_real_t w[7], _val[21];
1460 cs_real_3_t gpts[7];
1461
1462 const cs_quant_t fq = cm->face[f];
1463
1464 eval[f] = 0.; /* Reset value */
1465
1466 /* Loop on face edges */
1467 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1468
1469 const short int e = cm->f2e_ids[i];
1470 const short int v1 = cm->e2v_ids[2*e];
1471 const short int v2 = cm->e2v_ids[2*e+1];
1472
1473 /* Evaluate the field at the three quadrature points */
1474 cs_quadrature_tria_7pts(fq.center, cm->xv + 3*v1, cm->xv + 3*v2,
1475 cm->tef[i],
1476 gpts, w);
1477
1478 /* Evaluate the function for this time at the given coordinates */
1479 ac->func(time_eval, 7, NULL,
1480 (const cs_real_t *)gpts, true, /* compacted output ? */
1481 ac->input,
1482 _val);
1483
1484 cs_real_t add = 0;
1485 for (int p = 0; p < 7; p++)
1486 add += w[p] * _dp3(_val+3*p, fq.unitv);
1487
1488 eval[f] += add;
1489
1490 }
1491
1492 }
1493 break;
1494
1495 default:
1496 bft_error(__FILE__, __LINE__, 0,
1497 "%s: Invalid type of quadrature.", __func__);
1498 break;
1499
1500 } /* Switch type of quadrature */
1501
1502 }
1503
1504 /*----------------------------------------------------------------------------*/
1505 /*!
1506 * \brief Function pointer for evaluating the normal flux of a quantity
1507 * defined by analytic function.
1508 * Use of a \ref cs_cell_mesh_t structure.
1509 * Case of tensor-valued quantities.
1510 *
1511 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1512 * \param[in] f local face id
1513 * \param[in] time_eval physical time at which one evaluates the term
1514 * \param[in] context pointer to a context structure
1515 * \param[in] qtype level of quadrature to use
1516 * \param[in, out] eval result of the evaluation (set inside)
1517 */
1518 /*----------------------------------------------------------------------------*/
1519
1520 void
cs_xdef_cw_eval_tensor_flux_by_analytic(const cs_cell_mesh_t * cm,short int f,cs_real_t time_eval,void * context,cs_quadrature_type_t qtype,cs_real_t * eval)1521 cs_xdef_cw_eval_tensor_flux_by_analytic(const cs_cell_mesh_t *cm,
1522 short int f,
1523 cs_real_t time_eval,
1524 void *context,
1525 cs_quadrature_type_t qtype,
1526 cs_real_t *eval)
1527 {
1528 cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
1529
1530 switch (qtype) {
1531
1532 case CS_QUADRATURE_NONE:
1533 case CS_QUADRATURE_BARY:
1534 {
1535 cs_real_33_t flux_xc = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
1536
1537 /* Evaluate the function for this time at the given coordinates */
1538 ac->func(time_eval, 1, NULL, cm->xc, true, /* compacted output ? */
1539 ac->input,
1540 (cs_real_t *)flux_xc);
1541
1542 /* Plug into the evaluation by value now */
1543 cs_xdef_cw_eval_tensor_flux_by_val(cm, f, time_eval, flux_xc, eval);
1544 }
1545 break;
1546
1547 case CS_QUADRATURE_BARY_SUBDIV:
1548 {
1549 assert(cs_flag_test(cm->flag,
1550 CS_FLAG_COMP_EV| CS_FLAG_COMP_FE| CS_FLAG_COMP_FEQ));
1551
1552 const cs_quant_t fq = cm->face[f];
1553
1554 cs_real_3_t _xyz, _val;
1555 cs_real_33_t _eval;
1556
1557 for (int k = 0; k < 3; k++)
1558 eval[3*f+k] = 0.; /* Reset value */
1559
1560 /* Loop on face edges */
1561 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1562
1563 const short int _2e = cm->f2e_ids[i];
1564 const short int v1 = cm->e2v_ids[_2e];
1565 const short int v2 = cm->e2v_ids[_2e+1];
1566
1567 for (int k = 0; k < 3; k++)
1568 _xyz[k] = cs_math_1ov3 *
1569 (fq.center[k] + cm->xv[3*v1+k] + cm->xv[3*v2+k]);
1570
1571 /* Evaluate the function for this time at the given coordinates */
1572 ac->func(time_eval, 1, NULL,
1573 (const cs_real_t *)_xyz, true, /* compacted output ? */
1574 ac->input,
1575 (cs_real_t *)_eval);
1576
1577 cs_math_33_3_product((const cs_real_t (*)[3])_eval, fq.unitv, _val);
1578 for (int k = 0; k < 3; k++)
1579 eval[3*f+k] += cm->tef[i] * _val[k];
1580
1581 }
1582
1583 }
1584 break; /* BARY_SUBDIV */
1585
1586 case CS_QUADRATURE_HIGHER:
1587 {
1588 assert(cs_flag_test(cm->flag,
1589 CS_FLAG_COMP_EV| CS_FLAG_COMP_FE | CS_FLAG_COMP_FEQ));
1590
1591 /* Four values by triangle --> 4 Gauss points
1592 * The flux returns by the analytic function is a 3x3 tensor. */
1593 cs_real_t w[4];
1594 cs_real_3_t gpts[4], _val;
1595 cs_real_33_t _eval[4];
1596
1597 const cs_quant_t fq = cm->face[f];
1598
1599 for (int k = 0; k < 3; k++)
1600 eval[3*f+k] = 0.; /* Reset value */
1601
1602 /* Loop on face edges */
1603 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1604
1605 const short int e = cm->f2e_ids[i];
1606 const short int v1 = cm->e2v_ids[2*e];
1607 const short int v2 = cm->e2v_ids[2*e+1];
1608
1609 /* Evaluate the field at the three quadrature points */
1610 cs_quadrature_tria_4pts(fq.center, cm->xv + 3*v1, cm->xv + 3*v2,
1611 cm->tef[i],
1612 gpts, w);
1613
1614 /* Evaluate the function for this time at the given coordinates */
1615 ac->func(time_eval, 4, NULL,
1616 (const cs_real_t *)gpts, true, /* compacted output ? */
1617 ac->input,
1618 (cs_real_t *)_eval);
1619
1620 for (int p = 0; p < 4; p++) {
1621 cs_math_33_3_product((const cs_real_t (*)[3])_eval[p], fq.unitv,
1622 _val);
1623 for (int k = 0; k < 3; k++)
1624 eval[3*f+k] += w[p] * _val[k];
1625 }
1626
1627 }
1628
1629 }
1630 break;
1631
1632 case CS_QUADRATURE_HIGHEST:
1633 {
1634 assert(cs_flag_test(cm->flag,
1635 CS_FLAG_COMP_EV| CS_FLAG_COMP_FE | CS_FLAG_COMP_FEQ));
1636
1637 /* Seven values by triangle --> 7 Gauss points
1638 * The flux returns by the analytic function is a 3x3 tensor. */
1639 cs_real_t w[7];
1640 cs_real_3_t gpts[7], _val;
1641 cs_real_33_t _eval[7];
1642
1643 const cs_quant_t fq = cm->face[f];
1644
1645 for (int k = 0; k < 3; k++)
1646 eval[3*f+k] = 0.; /* Reset value */
1647
1648 /* Loop on face edges */
1649 for (int i = cm->f2e_idx[f]; i < cm->f2e_idx[f+1]; i++) {
1650
1651 const short int e = cm->f2e_ids[i];
1652 const short int v1 = cm->e2v_ids[2*e];
1653 const short int v2 = cm->e2v_ids[2*e+1];
1654
1655 /* Evaluate the field at the three quadrature points */
1656 cs_quadrature_tria_7pts(fq.center, cm->xv + 3*v1, cm->xv + 3*v2,
1657 cm->tef[i],
1658 gpts, w);
1659
1660 /* Evaluate the function for this time at the given coordinates */
1661 ac->func(time_eval, 7, NULL,
1662 (const cs_real_t *)gpts, true, /* compacted output ? */
1663 ac->input,
1664 (cs_real_t *)_eval);
1665
1666 for (int p = 0; p < 7; p++) {
1667 cs_math_33_3_product((const cs_real_t (*)[3])_eval[p], fq.unitv,
1668 _val);
1669 for (int k = 0; k < 3; k++)
1670 eval[3*f+k] += w[p] * _val[k];
1671 }
1672
1673 }
1674
1675 }
1676 break;
1677
1678 default:
1679 bft_error(__FILE__, __LINE__, 0,
1680 " %s: Invalid type of quadrature.", __func__);
1681 break;
1682
1683 } /* Switch type of quadrature */
1684
1685 }
1686
1687 /*----------------------------------------------------------------------------*/
1688 /*!
1689 * \brief Function pointer for evaluating the reduction by averages of a
1690 * analytic function by a cellwise process (usage of a
1691 * \ref cs_cell_mesh_t structure).
1692 * This evaluation hinges on the computation of integrals (faces first,
1693 * then cell)
1694 * Scalar-valued case.
1695 *
1696 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1697 * \param[in] t_eval physical time at which one evaluates the term
1698 * \param[in] qtype quadrature type
1699 * \param[in] context pointer to a context structure
1700 * \param[in, out] eval result of the evaluation
1701 */
1702 /*----------------------------------------------------------------------------*/
1703
1704 void
cs_xdef_cw_eval_scal_avg_reduction_by_analytic(const cs_cell_mesh_t * cm,cs_real_t t_eval,void * context,cs_quadrature_type_t qtype,cs_real_t * eval)1705 cs_xdef_cw_eval_scal_avg_reduction_by_analytic(const cs_cell_mesh_t *cm,
1706 cs_real_t t_eval,
1707 void *context,
1708 cs_quadrature_type_t qtype,
1709 cs_real_t *eval)
1710 {
1711 if (eval == NULL)
1712 bft_error(__FILE__, __LINE__, 0, _err_empty_array, __func__);
1713
1714 assert(context != NULL);
1715 assert(cs_flag_test(cm->flag,
1716 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
1717 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
1718
1719 const int dim = 1;
1720 const short int nf = cm->n_fc;
1721
1722 cs_quadrature_tetra_integral_t
1723 *q_tet = cs_quadrature_get_tetra_integral(dim, qtype);
1724 cs_quadrature_tria_integral_t
1725 *q_tri = cs_quadrature_get_tria_integral(dim, qtype);
1726 cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
1727 cs_real_t *c_eval = eval + nf;
1728
1729 cs_xdef_cw_eval_fc_int_by_analytic(cm, t_eval,
1730 ac->func, ac->input,
1731 dim,
1732 q_tet, q_tri,
1733 c_eval, eval);
1734
1735 /* Compute the averages */
1736 for (short int f = 0; f < nf; f++)
1737 eval[f] /= cm->face[f].meas;
1738 eval[nf] /= cm->vol_c;
1739 }
1740
1741 /*----------------------------------------------------------------------------*/
1742 /*!
1743 * \brief Function pointer for evaluating the reduction by averages of a
1744 * analytic function by a cellwise process (usage of a
1745 * \ref cs_cell_mesh_t structure).
1746 * This evaluation hinges on the computation of integrals (faces first,
1747 * then cell)
1748 * Vector-valued case.
1749 *
1750 * \param[in] cm pointer to a \ref cs_cell_mesh_t structure
1751 * \param[in] t_eval physical time at which one evaluates the term
1752 * \param[in] qtype quadrature type
1753 * \param[in] context pointer to a context structure
1754 * \param[in, out] eval result of the evaluation
1755 */
1756 /*----------------------------------------------------------------------------*/
1757
1758 void
cs_xdef_cw_eval_vect_avg_reduction_by_analytic(const cs_cell_mesh_t * cm,cs_real_t t_eval,void * context,cs_quadrature_type_t qtype,cs_real_t * eval)1759 cs_xdef_cw_eval_vect_avg_reduction_by_analytic(const cs_cell_mesh_t *cm,
1760 cs_real_t t_eval,
1761 void *context,
1762 cs_quadrature_type_t qtype,
1763 cs_real_t *eval)
1764 {
1765 if (eval == NULL)
1766 bft_error(__FILE__, __LINE__, 0, _err_empty_array, __func__);
1767
1768 assert(context != NULL);
1769 assert(cs_flag_test(cm->flag,
1770 CS_FLAG_COMP_PEQ | CS_FLAG_COMP_PFQ | CS_FLAG_COMP_FE |
1771 CS_FLAG_COMP_FEQ | CS_FLAG_COMP_EV));
1772
1773 const int dim = 3;
1774 const short int nf = cm->n_fc;
1775
1776 cs_quadrature_tetra_integral_t
1777 *q_tet = cs_quadrature_get_tetra_integral(3, qtype);
1778 cs_quadrature_tria_integral_t
1779 *q_tri = cs_quadrature_get_tria_integral(3, qtype);
1780 cs_xdef_analytic_context_t *ac = (cs_xdef_analytic_context_t *)context;
1781 cs_real_t *c_eval = eval + dim*nf;
1782
1783 cs_xdef_cw_eval_fc_int_by_analytic(cm, t_eval,
1784 ac->func, ac->input,
1785 dim,
1786 q_tet, q_tri,
1787 c_eval, eval);
1788
1789 /* Compute the averages */
1790 for (short int f = 0; f < nf; f++) {
1791 const cs_real_t _os = 1. / cm->face[f].meas;
1792 cs_real_t *f_eval = eval + dim*f;
1793 f_eval[0] *= _os, f_eval[1] *= _os, f_eval[2] *= _os;
1794 }
1795
1796 const cs_real_t _ov = 1. / cm->vol_c;
1797 c_eval[0] *= _ov, c_eval[1] *= _ov, c_eval[2] *= _ov;
1798 }
1799
1800 /*----------------------------------------------------------------------------*/
1801
1802 #undef _dp3
1803
1804 END_C_DECLS
1805