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