8.3
general documentation
cs_quadrature.h
Go to the documentation of this file.
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-2024 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
42BEGIN_C_DECLS
43
44/*============================================================================
45 * Macro definitions
46 *============================================================================*/
47
48/*============================================================================
49 * Type definitions
50 *============================================================================*/
51
84typedef enum {
85
86 CS_QUADRATURE_NONE,
87 CS_QUADRATURE_BARY, /* Value at the barycenter * meas */
88 CS_QUADRATURE_BARY_SUBDIV, /* Value at the barycenter * meas on a sub-mesh */
89 CS_QUADRATURE_HIGHER, /* Unique weight but several Gauss points */
90 CS_QUADRATURE_HIGHEST, /* Specific weight for each Gauss points */
91 CS_QUADRATURE_N_TYPES
92
93} cs_quadrature_type_t;
94
95/*----------------------------------------------------------------------------*/
106/*----------------------------------------------------------------------------*/
107
108typedef void
109(cs_quadrature_edge_t) (const cs_real_3_t v1,
110 const cs_real_3_t v2,
111 double len,
112 cs_real_3_t gpts[],
113 double *weights);
114
115/*----------------------------------------------------------------------------*/
127/*----------------------------------------------------------------------------*/
128
129typedef void
130(cs_quadrature_tria_t) (const cs_real_3_t v1,
131 const cs_real_3_t v2,
132 const cs_real_3_t v3,
133 double area,
134 cs_real_3_t gpts[],
135 double *weights);
136
137/*----------------------------------------------------------------------------*/
149/*----------------------------------------------------------------------------*/
150
151typedef void
152(cs_quadrature_tet_t) (const cs_real_3_t v1,
153 const cs_real_3_t v2,
154 const cs_real_3_t v3,
155 const cs_real_3_t v4,
156 double vol,
157 cs_real_3_t gpts[],
158 double weights[]);
159
160/*----------------------------------------------------------------------------*/
173/*----------------------------------------------------------------------------*/
174
175typedef void
176(cs_quadrature_edge_integral_t)(double tcur,
177 const cs_real_3_t v1,
178 const cs_real_3_t v2,
179 double len,
180 cs_analytic_func_t *ana,
181 void *input,
182 double results[]);
183
184/*----------------------------------------------------------------------------*/
198/*----------------------------------------------------------------------------*/
199
200typedef void
201(cs_quadrature_tria_integral_t)(double tcur,
202 const cs_real_3_t v1,
203 const cs_real_3_t v2,
204 const cs_real_3_t v3,
205 double area,
206 cs_analytic_func_t *ana,
207 void *input,
208 double results[]);
209
210/*----------------------------------------------------------------------------*/
225/*----------------------------------------------------------------------------*/
226
227typedef void
228(cs_quadrature_tetra_integral_t)(double tcur,
229 const cs_real_3_t v1,
230 const cs_real_3_t v2,
231 const cs_real_3_t v3,
232 const cs_real_3_t v4,
233 double vol,
234 cs_analytic_func_t *ana,
235 void *input,
236 double results[]);
237
238/*----------------------------------------------------------------------------*/
252/*----------------------------------------------------------------------------*/
253
254typedef void(cs_quadrature_hexa_integral_t)(double tcur,
255 const cs_real_3_t vb,
256 const double hx,
257 const double hy,
258 const double hz,
259 cs_analytic_func_t *ana,
260 void *input,
261 double results[]);
262
263/*============================================================================
264 * Public function prototypes
265 *============================================================================*/
266
267/*----------------------------------------------------------------------------*/
271/*----------------------------------------------------------------------------*/
272
273void
274cs_quadrature_setup(void);
275
276/*----------------------------------------------------------------------------*/
284/*----------------------------------------------------------------------------*/
285
286const char *
287cs_quadrature_get_type_name(const cs_quadrature_type_t type);
288
289/*----------------------------------------------------------------------------*/
300/*----------------------------------------------------------------------------*/
301
302static inline void
303cs_quadrature_edge_1pt(const cs_real_3_t v1,
304 const cs_real_3_t v2,
305 double len,
306 cs_real_3_t gpts[],
307 double *w)
308{
309 gpts[0][0] = 0.5*(v1[0] + v2[0]);
310 gpts[0][1] = 0.5*(v1[1] + v2[1]);
311 gpts[0][2] = 0.5*(v1[2] + v2[2]);
312 w[0] = len;
313}
314
315/*----------------------------------------------------------------------------*/
326/*----------------------------------------------------------------------------*/
327
328void
329cs_quadrature_edge_2pts(const cs_real_3_t v1,
330 const cs_real_3_t v2,
331 double len,
332 cs_real_3_t gpts[],
333 double *w);
334
335/*----------------------------------------------------------------------------*/
346/*----------------------------------------------------------------------------*/
347
348void
349cs_quadrature_edge_3pts(const cs_real_3_t v1,
350 const cs_real_3_t v2,
351 double len,
352 cs_real_3_t gpts[],
353 double w[]);
354
355/*----------------------------------------------------------------------------*/
367/*----------------------------------------------------------------------------*/
368
369static inline void
370cs_quadrature_hex_1pt(const cs_real_3_t vb,
371 const double hx,
372 const double hy,
373 const double hz,
374 cs_real_3_t gpts[],
375 double *w)
376{
377 gpts[0][0] = vb[0];
378 gpts[0][1] = vb[1];
379 gpts[0][2] = vb[2];
380 w[0] = hx * hy * hz;
381}
382
383/*----------------------------------------------------------------------------*/
395/*----------------------------------------------------------------------------*/
396
397void cs_quadrature_hex_8pts(const cs_real_3_t vb,
398 const double hx,
399 const double hy,
400 const double hz,
401 cs_real_3_t gpts[],
402 double *w);
403
404/*----------------------------------------------------------------------------*/
416/*----------------------------------------------------------------------------*/
417
418void cs_quadrature_hex_27pts(const cs_real_3_t vb,
419 const double hx,
420 const double hy,
421 const double hz,
422 cs_real_3_t gpts[],
423 double w[]);
424
425/*----------------------------------------------------------------------------*/
437/*----------------------------------------------------------------------------*/
438
439static inline void
440cs_quadrature_tria_1pt(const cs_real_3_t v1,
441 const cs_real_3_t v2,
442 const cs_real_3_t v3,
443 double area,
444 cs_real_3_t gpts[],
445 double *w)
446{
447#ifdef __cplusplus
448 constexpr cs_real_t c_1ov3 = 1./3.;
449 gpts[0][0] = c_1ov3 * (v1[0] + v2[0] + v3[0]);
450 gpts[0][1] = c_1ov3 * (v1[1] + v2[1] + v3[1]);
451 gpts[0][2] = c_1ov3 * (v1[2] + v2[2] + v3[2]);
452 w[0] = area;
453#else
454 gpts[0][0] = cs_math_1ov3 * (v1[0] + v2[0] + v3[0]);
455 gpts[0][1] = cs_math_1ov3 * (v1[1] + v2[1] + v3[1]);
456 gpts[0][2] = cs_math_1ov3 * (v1[2] + v2[2] + v3[2]);
457 w[0] = area;
458#endif
459}
460
461/*----------------------------------------------------------------------------*/
473/*----------------------------------------------------------------------------*/
474
475void
476cs_quadrature_tria_3pts(const cs_real_3_t v1,
477 const cs_real_3_t v2,
478 const cs_real_3_t v3,
479 double area,
480 cs_real_3_t gpts[],
481 double *w);
482
483/*----------------------------------------------------------------------------*/
495/*----------------------------------------------------------------------------*/
496
497void
498cs_quadrature_tria_4pts(const cs_real_3_t v1,
499 const cs_real_3_t v2,
500 const cs_real_3_t v3,
501 double area,
502 cs_real_3_t gpts[],
503 double w[]);
504
505/*----------------------------------------------------------------------------*/
517/*----------------------------------------------------------------------------*/
518
519void
520cs_quadrature_tria_7pts(const cs_real_3_t v1,
521 const cs_real_3_t v2,
522 const cs_real_3_t v3,
523 double area,
524 cs_real_3_t gpts[],
525 double w[]);
526
527/*----------------------------------------------------------------------------*/
540/*----------------------------------------------------------------------------*/
541
542static inline void
543cs_quadrature_tet_1pt(const cs_real_3_t v1,
544 const cs_real_3_t v2,
545 const cs_real_3_t v3,
546 const cs_real_3_t v4,
547 double vol,
548 cs_real_3_t gpts[],
549 double weight[])
550{
551 gpts[0][0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
552 gpts[0][1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
553 gpts[0][2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
554 weight[0] = vol;
555}
556
557/*----------------------------------------------------------------------------*/
570/*----------------------------------------------------------------------------*/
571
572void
573cs_quadrature_tet_4pts(const cs_real_3_t v1,
574 const cs_real_3_t v2,
575 const cs_real_3_t v3,
576 const cs_real_3_t v4,
577 double vol,
578 cs_real_3_t gpts[],
579 double weights[]);
580
581/*----------------------------------------------------------------------------*/
594/*----------------------------------------------------------------------------*/
595
596void
597cs_quadrature_tet_5pts(const cs_real_3_t v1,
598 const cs_real_3_t v2,
599 const cs_real_3_t v3,
600 const cs_real_3_t v4,
601 double vol,
602 cs_real_3_t gpts[],
603 double weights[]);
604
605/*----------------------------------------------------------------------------*/
618/*----------------------------------------------------------------------------*/
619
620void
621cs_quadrature_tet_15pts(const cs_real_3_t v1,
622 const cs_real_3_t v2,
623 const cs_real_3_t v3,
624 const cs_real_3_t v4,
625 double vol,
626 cs_real_3_t gpts[],
627 double weights[]);
628
629/*----------------------------------------------------------------------------*/
643/*----------------------------------------------------------------------------*/
644
645static inline void
646cs_quadrature_edge_1pt_scal(double tcur,
647 const cs_real_3_t v1,
648 const cs_real_3_t v2,
649 double len,
650 cs_analytic_func_t *ana,
651 void *input,
652 double results[])
653{
654 cs_real_3_t xg;
655 double feval;
656
657 /* Copied from cs_quadrature_1pt */
658 xg[0] = .5 * (v1[0] + v2[0]);
659 xg[1] = .5 * (v1[1] + v2[1]);
660 xg[2] = .5 * (v1[2] + v2[2]);
661
662 /* Evaluate the function at the Gauss points */
663 ana(tcur, 1, NULL, xg, false, input, &feval);
664
665 /* Update the result with the quadrature rule */
666 *results += len * feval;
667}
668
669/*----------------------------------------------------------------------------*/
683/*----------------------------------------------------------------------------*/
684
685static inline void
686cs_quadrature_edge_2pts_scal(double tcur,
687 const cs_real_3_t v1,
688 const cs_real_3_t v2,
689 double len,
690 cs_analytic_func_t *ana,
691 void *input,
692 double results[])
693{
694 cs_real_3_t gauss_pts[2];
695 double feval[2], weights[2];
696
697 /* Compute Gauss points and its unique weight */
698 cs_quadrature_edge_2pts(v1, v2, len, gauss_pts, weights);
699
700 /* Evaluate the function at the Gauss points */
701 ana(tcur, 2, NULL, (const cs_real_t *)gauss_pts, false, input, feval);
702
703 /* Update the result with the quadrature rule */
704 *results += weights[0] * feval[0] + weights[1] * feval[1];
705}
706
707/*----------------------------------------------------------------------------*/
721/*----------------------------------------------------------------------------*/
722
723static inline void
724cs_quadrature_edge_3pts_scal(double tcur,
725 const cs_real_3_t v1,
726 const cs_real_3_t v2,
727 double len,
728 cs_analytic_func_t *ana,
729 void *input,
730 double results[])
731{
732 cs_real_3_t gauss_pts[3];
733 double feval[3], weights[3];
734
735 /* Compute Gauss points and its weights */
736 cs_quadrature_edge_3pts(v1, v2, len, gauss_pts, weights);
737
738 /* Evaluate the function at the Gauss points */
739 ana(tcur, 3, NULL, (const cs_real_t *)gauss_pts, false, input, feval);
740
741 /* Update the result with the quadrature rule */
742 *results += weights[0]*feval[0] + weights[1]*feval[1] + weights[2]*feval[2];
743}
744
745/*----------------------------------------------------------------------------*/
759/*----------------------------------------------------------------------------*/
760
761static inline void
762cs_quadrature_edge_1pt_vect(double tcur,
763 const cs_real_3_t v1,
764 const cs_real_3_t v2,
765 double len,
766 cs_analytic_func_t *ana,
767 void *input,
768 double results[])
769{
770 cs_real_3_t xg;
771 double feval[3];
772
773 /* Copied from cs_quadrature_1pt */
774 xg[0] = .5 * (v1[0] + v2[0]);
775 xg[1] = .5 * (v1[1] + v2[1]);
776 xg[2] = .5 * (v1[2] + v2[2]);
777
778 /* Evaluate the function at the Gauss points */
779 ana(tcur, 1, NULL, xg, false, input, feval);
780
781 /* Update the result with the quadrature rule */
782 results[0] += len * feval[0];
783 results[1] += len * feval[1];
784 results[2] += len * feval[2];
785}
786
787/*----------------------------------------------------------------------------*/
801/*----------------------------------------------------------------------------*/
802
803static inline void
804cs_quadrature_edge_2pts_vect(double tcur,
805 const cs_real_3_t v1,
806 const cs_real_3_t v2,
807 double len,
808 cs_analytic_func_t *ana,
809 void *input,
810 double results[])
811{
812 cs_real_3_t gauss_pts[2];
813 double feval[6], weights[2];
814
815 /* Compute Gauss points and its unique weight */
816 cs_quadrature_edge_2pts(v1, v2, len, gauss_pts, weights);
817
818 /* Evaluate the function at the Gauss points */
819 ana(tcur, 2, NULL, (const cs_real_t *)gauss_pts, false, input, feval);
820
821 /* Update the result with the quadrature rule */
822 results[0] += weights[0] * feval[0] + weights[1] * feval[3];
823 results[1] += weights[0] * feval[1] + weights[1] * feval[4];
824 results[2] += weights[0] * feval[2] + weights[1] * feval[5];
825}
826
827/*----------------------------------------------------------------------------*/
841/*----------------------------------------------------------------------------*/
842
843static inline void
844cs_quadrature_edge_3pts_vect(double tcur,
845 const cs_real_3_t v1,
846 const cs_real_3_t v2,
847 double len,
848 cs_analytic_func_t *ana,
849 void *input,
850 double results[])
851{
852 cs_real_3_t gauss_pts[3];
853 double feval[9], weights[3];
854
855 /* Compute Gauss points and its weights */
856 cs_quadrature_edge_3pts(v1, v2, len, gauss_pts, weights);
857
858 /* Evaluate the function at the Gauss points */
859 ana(tcur, 3, NULL, (const cs_real_t *)gauss_pts, false, input, feval);
860
861 /* Update the result with the quadrature rule */
862 for (int p = 0; p < 3; p++) {
863 results[0] += weights[p] * feval[3*p ];
864 results[1] += weights[p] * feval[3*p+1];
865 results[2] += weights[p] * feval[3*p+2];
866 }
867}
868
869/*----------------------------------------------------------------------------*/
884/*----------------------------------------------------------------------------*/
885
886static inline void
887cs_quadrature_tria_1pt_scal(double tcur,
888 const cs_real_3_t v1,
889 const cs_real_3_t v2,
890 const cs_real_3_t v3,
891 double area,
892 cs_analytic_func_t *ana,
893 void *input,
894 double results[])
895{
896 cs_real_3_t xg;
897 double evaluation;
898
899 /* Copied from cs_quadrature_1pt */
900#ifdef __cplusplus
901 constexpr cs_real_t c_1ov3 = 1./3.;
902 xg[0] = c_1ov3 * (v1[0] + v2[0] + v3[0]);
903 xg[1] = c_1ov3 * (v1[1] + v2[1] + v3[1]);
904 xg[2] = c_1ov3 * (v1[2] + v2[2] + v3[2]);
905#else
906 xg[0] = cs_math_1ov3 * (v1[0] + v2[0] + v3[0]);
907 xg[1] = cs_math_1ov3 * (v1[1] + v2[1] + v3[1]);
908 xg[2] = cs_math_1ov3 * (v1[2] + v2[2] + v3[2]);
909#endif
910
911 ana(tcur, 1, NULL, xg, false, input, &evaluation);
912
913 *results += area * evaluation;
914}
915
916/*----------------------------------------------------------------------------*/
931/*----------------------------------------------------------------------------*/
932
933static inline void
934cs_quadrature_tria_3pts_scal(double tcur,
935 const cs_real_3_t v1,
936 const cs_real_3_t v2,
937 const cs_real_3_t v3,
938 double area,
939 cs_analytic_func_t *ana,
940 void *input,
941 double results[])
942{
943 cs_real_3_t gauss_pts[3];
944 double evaluation[3], weights[3];
945
946 /* Compute Gauss points and its unique weight */
947 cs_quadrature_tria_3pts(v1, v2, v3, area, gauss_pts, weights);
948
949 ana(tcur, 3, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
950
951 /* Return results */
952 *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
953 weights[2] * evaluation[2];
954}
955
956/*----------------------------------------------------------------------------*/
971/*----------------------------------------------------------------------------*/
972
973static inline void
974cs_quadrature_tria_4pts_scal(double tcur,
975 const cs_real_3_t v1,
976 const cs_real_3_t v2,
977 const cs_real_3_t v3,
978 double area,
979 cs_analytic_func_t *ana,
980 void *input,
981 double results[])
982{
983 cs_real_3_t gauss_pts[4];
984 double evaluation[4], weights[4];
985
986 /* Compute Gauss points and its weights */
987 cs_quadrature_tria_4pts(v1, v2, v3, area, gauss_pts, weights);
988
989 ana(tcur, 4, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
990
991 *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
992 weights[2] * evaluation[2] + weights[3] * evaluation[3];
993}
994
995/*----------------------------------------------------------------------------*/
1010/*----------------------------------------------------------------------------*/
1011
1012static inline void
1013cs_quadrature_tria_7pts_scal(double tcur,
1014 const cs_real_3_t v1,
1015 const cs_real_3_t v2,
1016 const cs_real_3_t v3,
1017 double area,
1018 cs_analytic_func_t *ana,
1019 void *input,
1020 double results[])
1021{
1022 cs_real_3_t gauss_pts[7];
1023 double evaluation[7], weights[7];
1024
1025 /* Compute Gauss points and its weights */
1026 cs_quadrature_tria_7pts(v1, v2, v3, area, gauss_pts, weights);
1027
1028 ana(tcur, 7, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1029
1030 *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
1031 weights[2] * evaluation[2] + weights[3] * evaluation[3] +
1032 weights[4] * evaluation[4] + weights[5] * evaluation[5] +
1033 weights[6] * evaluation[6] ;
1034}
1035
1036/*----------------------------------------------------------------------------*/
1051/*----------------------------------------------------------------------------*/
1052
1053static inline void
1054cs_quadrature_tria_1pt_vect(double tcur,
1055 const cs_real_3_t v1,
1056 const cs_real_3_t v2,
1057 const cs_real_3_t v3,
1058 double area,
1059 cs_analytic_func_t *ana,
1060 void *input,
1061 double results[])
1062{
1063 cs_real_3_t xg;
1064 double evaluation[3];
1065
1066 /* Copied from cs_quadrature_1pt */
1067#ifdef __cplusplus
1068 constexpr cs_real_t c_1ov3 = 1./3.;
1069 xg[0] = c_1ov3 * (v1[0] + v2[0] + v3[0]);
1070 xg[1] = c_1ov3 * (v1[1] + v2[1] + v3[1]);
1071 xg[2] = c_1ov3 * (v1[2] + v2[2] + v3[2]);
1072#else
1073 xg[0] = cs_math_1ov3 * (v1[0] + v2[0] + v3[0]);
1074 xg[1] = cs_math_1ov3 * (v1[1] + v2[1] + v3[1]);
1075 xg[2] = cs_math_1ov3 * (v1[2] + v2[2] + v3[2]);
1076#endif
1077
1078 ana(tcur, 1, NULL, xg, false, input, evaluation);
1079
1080 results[0] += area * evaluation[0];
1081 results[1] += area * evaluation[1];
1082 results[2] += area * evaluation[2];
1083}
1084
1085/*----------------------------------------------------------------------------*/
1100/*----------------------------------------------------------------------------*/
1101
1102static inline void
1103cs_quadrature_tria_3pts_vect(double tcur,
1104 const cs_real_3_t v1,
1105 const cs_real_3_t v2,
1106 const cs_real_3_t v3,
1107 double area,
1108 cs_analytic_func_t *ana,
1109 void *input,
1110 double results[])
1111{
1112 cs_real_3_t gauss_pts[3];
1113 double evaluation[3*3], weights[3];
1114
1115 /* Compute Gauss points and its unique weight */
1116 cs_quadrature_tria_3pts(v1, v2, v3, area, gauss_pts, weights);
1117
1118 ana(tcur, 3, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1119
1120 for (int p = 0; p < 3; p++) {
1121 results[0] += weights[p] * evaluation[3*p ];
1122 results[1] += weights[p] * evaluation[3*p+1];
1123 results[2] += weights[p] * evaluation[3*p+2];
1124 }
1125}
1126
1127/*----------------------------------------------------------------------------*/
1142/*----------------------------------------------------------------------------*/
1143
1144static inline void
1145cs_quadrature_tria_4pts_vect(double tcur,
1146 const cs_real_3_t v1,
1147 const cs_real_3_t v2,
1148 const cs_real_3_t v3,
1149 double area,
1150 cs_analytic_func_t *ana,
1151 void *input,
1152 double results[])
1153{
1154 cs_real_3_t gauss_pts[4];
1155 double evaluation[3*4], weights[4];
1156
1157 /* Compute Gauss points and its weights */
1158 cs_quadrature_tria_4pts(v1, v2, v3, area, gauss_pts, weights);
1159
1160 ana(tcur, 4, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1161
1162 for (int p = 0; p < 4; p++) {
1163 results[0] += weights[p] * evaluation[3*p ];
1164 results[1] += weights[p] * evaluation[3*p+1];
1165 results[2] += weights[p] * evaluation[3*p+2];
1166 }
1167}
1168
1169/*----------------------------------------------------------------------------*/
1184/*----------------------------------------------------------------------------*/
1185
1186static inline void
1187cs_quadrature_tria_7pts_vect(double tcur,
1188 const cs_real_3_t v1,
1189 const cs_real_3_t v2,
1190 const cs_real_3_t v3,
1191 double area,
1192 cs_analytic_func_t *ana,
1193 void *input,
1194 double results[])
1195{
1196 cs_real_3_t gauss_pts[7];
1197 double evaluation[3*7], weights[7];
1198
1199 /* Compute Gauss points and its weights */
1200 cs_quadrature_tria_7pts(v1, v2, v3, area, gauss_pts, weights);
1201
1202 ana(tcur, 7, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1203
1204 for (int p = 0; p < 7; p++) {
1205 results[0] += weights[p] * evaluation[3*p ];
1206 results[1] += weights[p] * evaluation[3*p+1];
1207 results[2] += weights[p] * evaluation[3*p+2];
1208 }
1209}
1210
1211/*----------------------------------------------------------------------------*/
1226/*----------------------------------------------------------------------------*/
1227
1228static inline void
1229cs_quadrature_tria_1pt_tens(double tcur,
1230 const cs_real_3_t v1,
1231 const cs_real_3_t v2,
1232 const cs_real_3_t v3,
1233 double area,
1234 cs_analytic_func_t *ana,
1235 void *input,
1236 double results[])
1237{
1238 cs_real_3_t xg;
1239 double evaluation[9];
1240
1241 /* Copied from cs_quadrature_1pt */
1242#ifdef __cplusplus
1243 constexpr cs_real_t c_1ov3 = 1./3.;
1244 xg[0] = c_1ov3 * (v1[0] + v2[0] + v3[0]);
1245 xg[1] = c_1ov3 * (v1[1] + v2[1] + v3[1]);
1246 xg[2] = c_1ov3 * (v1[2] + v2[2] + v3[2]);
1247#else
1248 xg[0] = cs_math_1ov3 * (v1[0] + v2[0] + v3[0]);
1249 xg[1] = cs_math_1ov3 * (v1[1] + v2[1] + v3[1]);
1250 xg[2] = cs_math_1ov3 * (v1[2] + v2[2] + v3[2]);
1251#endif
1252
1253 ana(tcur, 1, NULL, xg, false, input, evaluation);
1254
1255 for (short int ij = 0; ij < 9; ij++)
1256 results[ij] += area * evaluation[ij];
1257}
1258
1259/*----------------------------------------------------------------------------*/
1274/*----------------------------------------------------------------------------*/
1275
1276static inline void
1277cs_quadrature_tria_3pts_tens(double tcur,
1278 const cs_real_3_t v1,
1279 const cs_real_3_t v2,
1280 const cs_real_3_t v3,
1281 double area,
1282 cs_analytic_func_t *ana,
1283 void *input,
1284 double results[])
1285{
1286 cs_real_3_t gauss_pts[3];
1287 double evaluation[9*3], weights[3];
1288
1289 /* Compute Gauss points and its unique weight */
1290 cs_quadrature_tria_3pts(v1, v2, v3, area, gauss_pts, weights);
1291
1292 ana(tcur, 3, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1293
1294 for (int p = 0; p < 3; p++) {
1295 const double wp = weights[p];
1296 double *eval_p = evaluation + 9*p;
1297 for (short int ij = 0; ij < 9; ij++)
1298 results[ij] += wp * eval_p[ij];
1299 }
1300}
1301
1302/*----------------------------------------------------------------------------*/
1317/*----------------------------------------------------------------------------*/
1318
1319static inline void
1320cs_quadrature_tria_4pts_tens(double tcur,
1321 const cs_real_3_t v1,
1322 const cs_real_3_t v2,
1323 const cs_real_3_t v3,
1324 double area,
1325 cs_analytic_func_t *ana,
1326 void *input,
1327 double results[])
1328{
1329 cs_real_3_t gauss_pts[4];
1330 double evaluation[9*4], weights[4];
1331
1332 /* Compute Gauss points and its weights */
1333 cs_quadrature_tria_4pts(v1, v2, v3, area, gauss_pts, weights);
1334
1335 ana(tcur, 4, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1336
1337 for (int p = 0; p < 4; p++) {
1338 const double wp = weights[p];
1339 double *eval_p = evaluation + 9*p;
1340 for (short int ij = 0; ij < 9; ij++)
1341 results[ij] += wp * eval_p[ij];
1342 }
1343}
1344
1345/*----------------------------------------------------------------------------*/
1360/*----------------------------------------------------------------------------*/
1361
1362static inline void
1363cs_quadrature_tria_7pts_tens(double tcur,
1364 const cs_real_3_t v1,
1365 const cs_real_3_t v2,
1366 const cs_real_3_t v3,
1367 double area,
1368 cs_analytic_func_t *ana,
1369 void *input,
1370 double results[])
1371{
1372 cs_real_3_t gauss_pts[7];
1373 double evaluation[9*7], weights[7];
1374
1375 /* Compute Gauss points and its weights */
1376 cs_quadrature_tria_7pts(v1, v2, v3, area, gauss_pts, weights);
1377
1378 ana(tcur, 7, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1379
1380 for (int p = 0; p < 7; p++) {
1381 const double wp = weights[p];
1382 double *eval_p = evaluation + 9*p;
1383 for (short int ij = 0; ij < 9; ij++)
1384 results[ij] += wp * eval_p[ij];
1385 }
1386}
1387
1388/*----------------------------------------------------------------------------*/
1404/*----------------------------------------------------------------------------*/
1405
1406static inline void
1407cs_quadrature_tet_1pt_scal(double tcur,
1408 const cs_real_3_t v1,
1409 const cs_real_3_t v2,
1410 const cs_real_3_t v3,
1411 const cs_real_3_t v4,
1412 double vol,
1413 cs_analytic_func_t *ana,
1414 void *input,
1415 double results[])
1416{
1417 cs_real_3_t xg;
1418 double evaluation;
1419
1420 /* Copied from cs_quadrature_tet_1pt */
1421 xg[0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
1422 xg[1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
1423 xg[2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
1424
1425 ana(tcur, 1, NULL, xg, false, input, &evaluation);
1426
1427 *results += vol * evaluation;
1428}
1429
1430/*----------------------------------------------------------------------------*/
1446/*----------------------------------------------------------------------------*/
1447
1448static inline void
1449cs_quadrature_tet_4pts_scal(double tcur,
1450 const cs_real_3_t v1,
1451 const cs_real_3_t v2,
1452 const cs_real_3_t v3,
1453 const cs_real_3_t v4,
1454 double vol,
1455 cs_analytic_func_t *ana,
1456 void *input,
1457 double results[])
1458{
1459 cs_real_3_t gauss_pts[4];
1460 double evaluation[4], weights[4];
1461
1462 /* Compute Gauss points and its unique weight */
1463 cs_quadrature_tet_4pts(v1, v2, v3, v4, vol, gauss_pts, weights);
1464
1465 ana(tcur, 4, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1466
1467 /* Return results */
1468 *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
1469 weights[2] * evaluation[2] + weights[3] * evaluation[3];
1470}
1471
1472/*----------------------------------------------------------------------------*/
1488/*----------------------------------------------------------------------------*/
1489
1490static inline void
1491cs_quadrature_tet_5pts_scal(double tcur,
1492 const cs_real_3_t v1,
1493 const cs_real_3_t v2,
1494 const cs_real_3_t v3,
1495 const cs_real_3_t v4,
1496 double vol,
1497 cs_analytic_func_t *ana,
1498 void *input,
1499 double results[])
1500{
1501 cs_real_3_t gauss_pts[5];
1502 double evaluation[5], weights[5];
1503
1504 /* Compute Gauss points and its weights */
1505 cs_quadrature_tet_5pts(v1, v2, v3, v4, vol, gauss_pts, weights);
1506
1507 ana(tcur, 5, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1508
1509 *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
1510 weights[2] * evaluation[2] + weights[3] * evaluation[3] +
1511 weights[4] * evaluation[4];
1512}
1513
1514/*----------------------------------------------------------------------------*/
1530/*----------------------------------------------------------------------------*/
1531
1532static inline void
1533cs_quadrature_tet_1pt_vect(double tcur,
1534 const cs_real_3_t v1,
1535 const cs_real_3_t v2,
1536 const cs_real_3_t v3,
1537 const cs_real_3_t v4,
1538 double vol,
1539 cs_analytic_func_t *ana,
1540 void *input,
1541 double results[])
1542{
1543 cs_real_3_t xg;
1544 double evaluation[3];
1545
1546 /* Copied from cs_quadrature_tet_1pt */
1547 xg[0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
1548 xg[1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
1549 xg[2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
1550
1551 ana(tcur, 1, NULL, xg, false, input, evaluation);
1552
1553 results[0] += vol * evaluation[0];
1554 results[1] += vol * evaluation[1];
1555 results[2] += vol * evaluation[2];
1556}
1557
1558/*----------------------------------------------------------------------------*/
1574/*----------------------------------------------------------------------------*/
1575
1576static inline void
1577cs_quadrature_tet_4pts_vect(double tcur,
1578 const cs_real_3_t v1,
1579 const cs_real_3_t v2,
1580 const cs_real_3_t v3,
1581 const cs_real_3_t v4,
1582 double vol,
1583 cs_analytic_func_t *ana,
1584 void *input,
1585 double results[])
1586{
1587 cs_real_3_t gauss_pts[4];
1588 double evaluation[3*4], weights[4];
1589
1590 /* Compute Gauss points and its unique weight */
1591 cs_quadrature_tet_4pts(v1, v2, v3, v4, vol, gauss_pts, weights);
1592
1593 ana(tcur, 4, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1594
1595 for (int p = 0; p < 4; p++) {
1596 results[0] += weights[p] * evaluation[3*p ];
1597 results[1] += weights[p] * evaluation[3*p+1];
1598 results[2] += weights[p] * evaluation[3*p+2];
1599 }
1600}
1601
1602/*----------------------------------------------------------------------------*/
1618/*----------------------------------------------------------------------------*/
1619
1620static inline void
1621cs_quadrature_tet_5pts_vect(double tcur,
1622 const cs_real_3_t v1,
1623 const cs_real_3_t v2,
1624 const cs_real_3_t v3,
1625 const cs_real_3_t v4,
1626 double vol,
1627 cs_analytic_func_t *ana,
1628 void *input,
1629 double results[])
1630{
1631 cs_real_3_t gauss_pts[5];
1632 double evaluation[3*5], weights[5];
1633
1634 /* Compute Gauss points and its weights */
1635 cs_quadrature_tet_5pts(v1, v2, v3, v4, vol, gauss_pts, weights);
1636
1637 ana(tcur, 5, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1638
1639 for (int p = 0; p < 5; p++) {
1640 results[0] += weights[p] * evaluation[3*p ];
1641 results[1] += weights[p] * evaluation[3*p+1];
1642 results[2] += weights[p] * evaluation[3*p+2];
1643 }
1644}
1645
1646/*----------------------------------------------------------------------------*/
1662/*----------------------------------------------------------------------------*/
1663
1664static inline void
1665cs_quadrature_tet_1pt_tens(double tcur,
1666 const cs_real_3_t v1,
1667 const cs_real_3_t v2,
1668 const cs_real_3_t v3,
1669 const cs_real_3_t v4,
1670 double vol,
1671 cs_analytic_func_t *ana,
1672 void *input,
1673 double results[])
1674{
1675 cs_real_3_t xg;
1676 double evaluation[9];
1677
1678 /* Copied from cs_quadrature_tet_1pt */
1679 xg[0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
1680 xg[1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
1681 xg[2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
1682
1683 ana(tcur, 1, NULL, xg, false, input, evaluation);
1684
1685 for (short int ij = 0; ij < 9; ij++)
1686 results[ij] += vol * evaluation[ij];
1687}
1688
1689/*----------------------------------------------------------------------------*/
1705/*----------------------------------------------------------------------------*/
1706
1707static inline void
1708cs_quadrature_tet_4pts_tens(double tcur,
1709 const cs_real_3_t v1,
1710 const cs_real_3_t v2,
1711 const cs_real_3_t v3,
1712 const cs_real_3_t v4,
1713 double vol,
1714 cs_analytic_func_t *ana,
1715 void *input,
1716 double results[])
1717{
1718 cs_real_3_t gauss_pts[4];
1719 double evaluation[9*4], weights[4];
1720
1721 /* Compute Gauss points and its unique weight */
1722 cs_quadrature_tet_4pts(v1, v2, v3, v4, vol, gauss_pts, weights);
1723
1724 ana(tcur, 4, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1725
1726 for (int p = 0; p < 4; p++) {
1727 const double wp = weights[p];
1728 double *eval_p = evaluation + 9*p;
1729 for (short int ij = 0; ij < 9; ij++)
1730 results[ij] += wp * eval_p[ij];
1731 }
1732}
1733
1734/*----------------------------------------------------------------------------*/
1750/*----------------------------------------------------------------------------*/
1751
1752static inline void
1753cs_quadrature_tet_5pts_tens(double tcur,
1754 const cs_real_3_t v1,
1755 const cs_real_3_t v2,
1756 const cs_real_3_t v3,
1757 const cs_real_3_t v4,
1758 double vol,
1759 cs_analytic_func_t *ana,
1760 void *input,
1761 double results[])
1762{
1763 cs_real_3_t gauss_pts[5];
1764 double evaluation[9*5], weights[5];
1765
1766 /* Compute Gauss points and its weights */
1767 cs_quadrature_tet_5pts(v1, v2, v3, v4, vol, gauss_pts, weights);
1768
1769 ana(tcur, 5, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1770
1771 for (int p = 0; p < 5; p++) {
1772 const double wp = weights[p];
1773 double *eval_p = evaluation + 9*p;
1774 for (short int ij = 0; ij < 9; ij++)
1775 results[ij] += wp * eval_p[ij];
1776 }
1777}
1778
1779/*----------------------------------------------------------------------------*/
1794/*----------------------------------------------------------------------------*/
1795
1796static inline void
1797cs_quadrature_hex_1pt_vect(double tcur,
1798 const cs_real_3_t vb,
1799 const double hx,
1800 const double hy,
1801 const double hz,
1802 cs_analytic_func_t *ana,
1803 void *input,
1804 double results[])
1805{
1806 cs_real_3_t xg[1];
1807 double weight[1], evaluation[3];
1808
1809 /* Copied from cs_quadrature_hew_1pt */
1810 cs_quadrature_hex_1pt(vb, hx, hy, hz, xg, weight);
1811
1812 ana(tcur, 1, NULL, (const cs_real_t *)xg, false, input, evaluation);
1813
1814 results[0] += weight[0] * evaluation[0];
1815 results[1] += weight[0] * evaluation[1];
1816 results[2] += weight[0] * evaluation[2];
1817}
1818
1819/*----------------------------------------------------------------------------*/
1834/*----------------------------------------------------------------------------*/
1835
1836static inline void
1837cs_quadrature_hex_8pts_vect(double tcur,
1838 const cs_real_3_t vb,
1839 const double hx,
1840 const double hy,
1841 const double hz,
1842 cs_analytic_func_t *ana,
1843 void *input,
1844 double results[])
1845{
1846 cs_real_3_t gauss_pts[8];
1847 double evaluation[3 * 8], weights[8];
1848
1849 /* Compute Gauss points and its unique weight */
1850 cs_quadrature_hex_8pts(vb, hx, hy, hz, gauss_pts, weights);
1851
1852 ana(tcur, 8, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1853
1854 for (int p = 0; p < 8; p++) {
1855 results[0] += weights[p] * evaluation[3 * p];
1856 results[1] += weights[p] * evaluation[3 * p + 1];
1857 results[2] += weights[p] * evaluation[3 * p + 2];
1858 }
1859}
1860
1861/*----------------------------------------------------------------------------*/
1876/*----------------------------------------------------------------------------*/
1877
1878static inline void
1879cs_quadrature_hex_27pts_vect(double tcur,
1880 const cs_real_3_t vb,
1881 const double hx,
1882 const double hy,
1883 const double hz,
1884 cs_analytic_func_t *ana,
1885 void *input,
1886 double results[])
1887{
1888 cs_real_3_t gauss_pts[27];
1889 double evaluation[3 * 27], weights[27];
1890
1891 /* Compute Gauss points and its weights */
1892 cs_quadrature_hex_27pts(vb, hx, hy, hz, gauss_pts, weights);
1893
1894 ana(tcur, 27, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1895
1896 for (int p = 0; p < 27; p++) {
1897 results[0] += weights[p] * evaluation[3 * p];
1898 results[1] += weights[p] * evaluation[3 * p + 1];
1899 results[2] += weights[p] * evaluation[3 * p + 2];
1900 }
1901}
1902
1903/*----------------------------------------------------------------------------*/
1914/*----------------------------------------------------------------------------*/
1915
1916static inline cs_quadrature_edge_integral_t *
1917cs_quadrature_get_edge_integral(int dim,
1918 cs_quadrature_type_t qtype)
1919{
1920 switch (dim) {
1921
1922 case 1: /* Scalar-valued integral */
1923
1924 switch (qtype) {
1925
1926 case CS_QUADRATURE_BARY:
1927 case CS_QUADRATURE_BARY_SUBDIV:
1928 return cs_quadrature_edge_1pt_scal;
1929 case CS_QUADRATURE_HIGHER:
1930 return cs_quadrature_edge_2pts_scal;
1931 case CS_QUADRATURE_HIGHEST:
1932 return cs_quadrature_edge_3pts_scal;
1933
1934 default:
1935 bft_error(__FILE__, __LINE__, 0,
1936 " %s: Invalid quadrature type\n", __func__);
1937 }
1938 break;
1939
1940 case 3: /* Vector-valued case */
1941
1942 switch (qtype) {
1943
1944 case CS_QUADRATURE_BARY:
1945 case CS_QUADRATURE_BARY_SUBDIV:
1946 return cs_quadrature_edge_1pt_vect;
1947 case CS_QUADRATURE_HIGHER:
1948 return cs_quadrature_edge_2pts_vect;
1949 case CS_QUADRATURE_HIGHEST:
1950 return cs_quadrature_edge_3pts_vect;
1951
1952 default:
1953 bft_error(__FILE__, __LINE__, 0,
1954 " %s: Invalid quadrature type\n", __func__);
1955 }
1956 break;
1957
1958 default:
1959 bft_error(__FILE__, __LINE__, 0,
1960 " %s: Invalid dimension value %d. Only 1 and 3 are valid.\n",
1961 __func__, dim);
1962
1963 } /* switch on dim */
1964
1965 return NULL; /* Should not go to this stage */
1966}
1967
1968/*----------------------------------------------------------------------------*/
1978/*----------------------------------------------------------------------------*/
1979
1980static inline cs_quadrature_tria_integral_t *
1981cs_quadrature_get_tria_integral(int dim,
1982 cs_quadrature_type_t qtype)
1983{
1984 switch (dim) {
1985
1986 case 1: /* Scalar-valued integral */
1987
1988 switch (qtype) {
1989
1990 case CS_QUADRATURE_BARY:
1991 case CS_QUADRATURE_BARY_SUBDIV:
1992 return cs_quadrature_tria_1pt_scal;
1993 case CS_QUADRATURE_HIGHER:
1994 return cs_quadrature_tria_4pts_scal;
1995 case CS_QUADRATURE_HIGHEST:
1996 return cs_quadrature_tria_7pts_scal;
1997
1998 default:
1999 bft_error(__FILE__, __LINE__, 0,
2000 " %s: Invalid quadrature type\n", __func__);
2001 }
2002 break;
2003
2004 case 3: /* Vector-valued case */
2005
2006 switch (qtype) {
2007
2008 case CS_QUADRATURE_BARY:
2009 case CS_QUADRATURE_BARY_SUBDIV:
2010 return cs_quadrature_tria_1pt_vect;
2011 case CS_QUADRATURE_HIGHER:
2012 return cs_quadrature_tria_4pts_vect;
2013 case CS_QUADRATURE_HIGHEST:
2014 return cs_quadrature_tria_7pts_vect;
2015
2016 default:
2017 bft_error(__FILE__, __LINE__, 0,
2018 " %s: Invalid quadrature type\n", __func__);
2019 }
2020 break;
2021
2022 case 9: /* Tensor-valued case */
2023
2024 switch (qtype) {
2025
2026 case CS_QUADRATURE_BARY:
2027 case CS_QUADRATURE_BARY_SUBDIV:
2028 return cs_quadrature_tria_1pt_tens;
2029 case CS_QUADRATURE_HIGHER:
2030 return cs_quadrature_tria_4pts_tens;
2031 case CS_QUADRATURE_HIGHEST:
2032 return cs_quadrature_tria_7pts_tens;
2033
2034 default:
2035 bft_error(__FILE__, __LINE__, 0,
2036 " %s: Invalid quadrature type\n", __func__);
2037 }
2038 break;
2039
2040 default:
2041 bft_error(__FILE__, __LINE__, 0,
2042 " %s: Invalid dimension value %d. Only 1, 3 and 9 are valid.\n",
2043 __func__, dim);
2044
2045 } /* switch on dim */
2046
2047 return NULL; /* Should not go to this stage */
2048}
2049
2050/*----------------------------------------------------------------------------*/
2060/*----------------------------------------------------------------------------*/
2061
2062static inline cs_quadrature_tetra_integral_t *
2063cs_quadrature_get_tetra_integral(int dim,
2064 cs_quadrature_type_t qtype)
2065{
2066 switch (dim) {
2067
2068 case 1: /* Scalar-valued case */
2069
2070 switch (qtype) {
2071
2072 case CS_QUADRATURE_BARY:
2073 case CS_QUADRATURE_BARY_SUBDIV:
2074 return cs_quadrature_tet_1pt_scal;
2075 case CS_QUADRATURE_HIGHER:
2076 return cs_quadrature_tet_4pts_scal;
2077 case CS_QUADRATURE_HIGHEST:
2078 return cs_quadrature_tet_5pts_scal;
2079
2080 default:
2081 bft_error(__FILE__, __LINE__, 0,
2082 " %s: Invalid quadrature type\n", __func__);
2083 }
2084 break;
2085
2086 case 3: /* Vector-valued case */
2087
2088 switch (qtype) {
2089
2090 case CS_QUADRATURE_BARY:
2091 case CS_QUADRATURE_BARY_SUBDIV:
2092 return cs_quadrature_tet_1pt_vect;
2093 case CS_QUADRATURE_HIGHER:
2094 return cs_quadrature_tet_4pts_vect;
2095 case CS_QUADRATURE_HIGHEST:
2096 return cs_quadrature_tet_5pts_vect;
2097
2098 default:
2099 bft_error(__FILE__, __LINE__, 0,
2100 " %s: Invalid quadrature type\n", __func__);
2101 }
2102 break;
2103
2104 case 9: /* Tensor-valued case */
2105
2106 switch (qtype) {
2107
2108 case CS_QUADRATURE_BARY:
2109 case CS_QUADRATURE_BARY_SUBDIV:
2110 return cs_quadrature_tet_1pt_tens;
2111 case CS_QUADRATURE_HIGHER:
2112 return cs_quadrature_tet_4pts_tens;
2113 case CS_QUADRATURE_HIGHEST:
2114 return cs_quadrature_tet_5pts_tens;
2115
2116 default:
2117 bft_error(__FILE__, __LINE__, 0,
2118 " %s: Invalid quadrature type\n", __func__);
2119 }
2120 break;
2121
2122 default:
2123 bft_error(__FILE__, __LINE__, 0,
2124 " %s: Invalid dimension value %d. Only 1, 3 and 9 are valid.\n",
2125 __func__, dim);
2126
2127 } /* Switch on dim */
2128
2129 /* Avoid no return warning */
2130 return NULL;
2131}
2132
2133/*----------------------------------------------------------------------------*/
2143/*----------------------------------------------------------------------------*/
2144
2145static inline cs_quadrature_hexa_integral_t *
2146cs_quadrature_get_hexa_integral(int dim, cs_quadrature_type_t qtype)
2147{
2148 switch (dim) {
2149
2150 case 3: /* Vector-valued case */
2151
2152 switch (qtype) {
2153
2154 case CS_QUADRATURE_BARY:
2155 case CS_QUADRATURE_BARY_SUBDIV:
2156 return cs_quadrature_hex_1pt_vect;
2157 case CS_QUADRATURE_HIGHER:
2158 return cs_quadrature_hex_8pts_vect;
2159 case CS_QUADRATURE_HIGHEST:
2160 return cs_quadrature_hex_27pts_vect;
2161
2162 default:
2163 bft_error(
2164 __FILE__, __LINE__, 0, " %s: Invalid quadrature type\n", __func__);
2165 }
2166 break;
2167
2168 default:
2169 bft_error(__FILE__,
2170 __LINE__,
2171 0,
2172 " %s: Invalid dimension value %d. Only 3 is valid.\n",
2173 __func__,
2174 dim);
2175
2176 } /* Switch on dim */
2177
2178 /* Avoid no return warning */
2179 return NULL;
2180}
2181
2182/*----------------------------------------------------------------------------*/
2194/*----------------------------------------------------------------------------*/
2195
2196cs_eflag_t
2197cs_quadrature_get_flag(const cs_quadrature_type_t qtype,
2198 const cs_flag_t loc);
2199
2200/*----------------------------------------------------------------------------*/
2201
2202END_C_DECLS
2203
2204#endif /* __CS_QUADRATURE_H__ */
bft_error
void bft_error(const char *const file_name, const int line_num, const int sys_error_code, const char *const format,...)
Calls the error handler (set by bft_error_handler_set() or default).
Definition: bft_error.cpp:193
bft_error.h
cs_base.h
cs_defs.h
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:542
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:342
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:359
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:543
cs_flag_t
unsigned short int cs_flag_t
Definition: cs_defs.h:344
p
@ p
Definition: cs_field_pointer.h:67
cs_flag.h
cs_eflag_t
unsigned int cs_eflag_t
Definition: cs_flag.h:192
cs_math.h
cs_math_1ov3
const cs_real_t cs_math_1ov3
cs_param_types.h
cs_analytic_func_t
void() cs_analytic_func_t(cs_real_t time, cs_lnum_t n_elts, const cs_lnum_t *elt_ids, const cs_real_t *coords, bool dense_output, void *input, cs_real_t *retval)
Generic function pointer for an evaluation relying on an analytic function.
Definition: cs_param_types.h:127
cs_quadrature_tet_5pts
void cs_quadrature_tet_5pts(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 weights[])
Compute the quadrature in a tetrehedra. Exact for 3rd order polynomials (order 4).
Definition: cs_quadrature.cpp:507
cs_quadrature_edge_1pt_scal
static 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[])
Compute the integral over an edge with the mid-point rule and add it to results Case of a scalar-valu...
Definition: cs_quadrature.h:646
cs_quadrature_hexa_integral_t
void() cs_quadrature_hexa_integral_t(double tcur, const cs_real_3_t vb, const double hx, const double hy, const double hz, cs_analytic_func_t *ana, void *input, double results[])
Compute the integral over an orthogonal hexahedron based on a specified quadrature rule and add it to...
Definition: cs_quadrature.h:254
cs_quadrature_tet_4pts_vect
static 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[])
Compute the integral over a tetrahedron with a quadrature rule using 4 Gauss points and a unique weig...
Definition: cs_quadrature.h:1577
cs_quadrature_tet_1pt
static 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[])
Compute the quadrature in a tetrehedra. Exact for 1st order polynomials (order 2).
Definition: cs_quadrature.h:543
cs_quadrature_tet_5pts_scal
static 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[])
Compute the integral over a tetrahedron with a quadrature rule using 5 Gauss points and 5 weights and...
Definition: cs_quadrature.h:1491
cs_quadrature_edge_3pts_scal
static 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[])
Compute the integral over an edge with a quadrature rule using 3 Gauss points and weights and add it ...
Definition: cs_quadrature.h:724
cs_quadrature_tet_5pts_tens
static 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[])
Compute the integral over a tetrahedron with a quadrature rule using 5 Gauss points and 5 weights and...
Definition: cs_quadrature.h:1753
cs_quadrature_tria_4pts_scal
static 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[])
Compute the integral over a triangle with a quadrature rule using 4 Gauss points and 4 weights and ad...
Definition: cs_quadrature.h:974
cs_quadrature_tria_7pts_scal
static 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[])
Compute the integral over a triangle with a quadrature rule using 7 Gauss points and 7 weights and ad...
Definition: cs_quadrature.h:1013
cs_quadrature_tria_7pts
void cs_quadrature_tria_7pts(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[])
Compute quadrature points for a triangle (7 points) Exact for polynomial function up to order 5.
Definition: cs_quadrature.cpp:423
cs_quadrature_get_edge_integral
static cs_quadrature_edge_integral_t * cs_quadrature_get_edge_integral(int dim, cs_quadrature_type_t qtype)
Retrieve the integral function according to the quadrature type and the stride provided Case of integ...
Definition: cs_quadrature.h:1917
cs_quadrature_hex_1pt
static void cs_quadrature_hex_1pt(const cs_real_3_t vb, const double hx, const double hy, const double hz, cs_real_3_t gpts[], double *w)
Compute quadrature points for an orthogonal hexahedron Exact for polynomial function up to order 1.
Definition: cs_quadrature.h:370
cs_quadrature_type_t
cs_quadrature_type_t
Type of quadrature to use when computing an integral quantity. This rationale is used for integrals a...
Definition: cs_quadrature.h:84
CS_QUADRATURE_HIGHEST
@ CS_QUADRATURE_HIGHEST
Definition: cs_quadrature.h:90
CS_QUADRATURE_NONE
@ CS_QUADRATURE_NONE
Definition: cs_quadrature.h:86
CS_QUADRATURE_BARY_SUBDIV
@ CS_QUADRATURE_BARY_SUBDIV
Definition: cs_quadrature.h:88
CS_QUADRATURE_N_TYPES
@ CS_QUADRATURE_N_TYPES
Definition: cs_quadrature.h:91
CS_QUADRATURE_HIGHER
@ CS_QUADRATURE_HIGHER
Definition: cs_quadrature.h:89
CS_QUADRATURE_BARY
@ CS_QUADRATURE_BARY
Definition: cs_quadrature.h:87
cs_quadrature_tria_integral_t
void() cs_quadrature_tria_integral_t(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[])
Compute the integral over a triangle based on a specified quadrature rule and add it to results.
Definition: cs_quadrature.h:201
cs_quadrature_tet_1pt_tens
static 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[])
Compute the integral over a tetrahedron using a barycentric quadrature rule and add it to results....
Definition: cs_quadrature.h:1665
cs_quadrature_edge_2pts_vect
static 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[])
Compute the integral over an edge with a quadrature rule using 2 Gauss points and a unique weight and...
Definition: cs_quadrature.h:804
cs_quadrature_tria_4pts_tens
static 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[])
Compute the integral over a triangle with a quadrature rule using 4 Gauss points and 4 weights and ad...
Definition: cs_quadrature.h:1320
cs_quadrature_tria_t
void() cs_quadrature_tria_t(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 *weights)
Generic functoin pointer to compute the quadrature points for a triangle.
Definition: cs_quadrature.h:130
cs_quadrature_edge_1pt
static 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)
Compute quadrature points for an edge from v1 -> v2 (2 points) Exact for polynomial function up to or...
Definition: cs_quadrature.h:303
cs_quadrature_tria_4pts_vect
static 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[])
Compute the integral over a triangle with a quadrature rule using 4 Gauss points and 4 weights and ad...
Definition: cs_quadrature.h:1145
cs_quadrature_tet_1pt_vect
static 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[])
Compute the integral over a tetrahedron using a barycentric quadrature rule and add it to results....
Definition: cs_quadrature.h:1533
cs_quadrature_tet_5pts_vect
static 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[])
Compute the integral over a tetrahedron with a quadrature rule using 5 Gauss points and 5 weights and...
Definition: cs_quadrature.h:1621
cs_quadrature_edge_3pts
void cs_quadrature_edge_3pts(const cs_real_3_t v1, const cs_real_3_t v2, double len, cs_real_3_t gpts[], double w[])
Compute quadrature points for an edge from v1 -> v2 (3 points) Exact for polynomial function up to or...
Definition: cs_quadrature.cpp:206
cs_quadrature_setup
void cs_quadrature_setup(void)
Compute constant weights for all quadratures used.
Definition: cs_quadrature.cpp:108
cs_quadrature_tria_7pts_vect
static 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[])
Compute the integral over a triangle with a quadrature rule using 7 Gauss points and 7 weights and ad...
Definition: cs_quadrature.h:1187
cs_quadrature_edge_2pts
void cs_quadrature_edge_2pts(const cs_real_3_t v1, const cs_real_3_t v2, double len, cs_real_3_t gpts[], double *w)
Compute quadrature points for an edge from v1 -> v2 (2 points) Exact for polynomial function up to or...
cs_quadrature_hex_8pts
void cs_quadrature_hex_8pts(const cs_real_3_t vb, const double hx, const double hy, const double hz, cs_real_3_t gpts[], double *w)
Compute quadrature points for an orthogonal hexahedron Exact for polynomial function up to order 3.
Definition: cs_quadrature.cpp:242
cs_quadrature_edge_2pts_scal
static 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[])
Compute the integral over an edge with a quadrature rule using 2 Gauss points and a unique weight and...
Definition: cs_quadrature.h:686
cs_quadrature_tet_1pt_scal
static 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[])
Compute the integral over a tetrahedron using a barycentric quadrature rule and add it to results....
Definition: cs_quadrature.h:1407
cs_quadrature_get_tetra_integral
static cs_quadrature_tetra_integral_t * cs_quadrature_get_tetra_integral(int dim, cs_quadrature_type_t qtype)
Retrieve the integral function according to the quadrature type and the stride provided.
Definition: cs_quadrature.h:2063
cs_quadrature_tria_3pts_vect
static 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[])
Compute the integral over a triangle with a quadrature rule using 3 Gauss points and a unique weight ...
Definition: cs_quadrature.h:1103
cs_quadrature_tria_1pt_scal
static 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[])
Compute the integral over a triangle using a barycentric quadrature rule and add it to results Case o...
Definition: cs_quadrature.h:887
cs_quadrature_tet_4pts_scal
static 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[])
Compute the integral over a tetrahedron with a quadrature rule using 4 Gauss points and a unique weig...
Definition: cs_quadrature.h:1449
cs_quadrature_hex_27pts
void cs_quadrature_hex_27pts(const cs_real_3_t vb, const double hx, const double hy, const double hz, cs_real_3_t gpts[], double w[])
Compute quadrature points for an orthogonal hexahedron Exact for polynomial function up to order 5.
Definition: cs_quadrature.cpp:295
cs_quadrature_tria_1pt
static 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)
Compute quadrature points for a triangle (1 point) Exact for polynomial function up to order 1 (baryc...
Definition: cs_quadrature.h:440
cs_quadrature_tria_4pts
void cs_quadrature_tria_4pts(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[])
Compute quadrature points for a triangle (4 points) Exact for polynomial function up to order 3.
Definition: cs_quadrature.cpp:386
cs_quadrature_tetra_integral_t
void() cs_quadrature_tetra_integral_t(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[])
Compute the integral over a tetrahedron based on a specified quadrature rule and add it to results.
Definition: cs_quadrature.h:228
cs_quadrature_tria_3pts
void cs_quadrature_tria_3pts(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)
Compute quadrature points for a triangle (3 points) Exact for polynomial function up to order 2.
cs_quadrature_tet_15pts
void cs_quadrature_tet_15pts(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 weights[])
Compute the quadrature in a tetrehedra. Exact for 5th order polynomials (order 6).
Definition: cs_quadrature.cpp:553
cs_quadrature_tria_1pt_vect
static 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[])
Compute the integral over a triangle using a barycentric quadrature rule and add it to results Case o...
Definition: cs_quadrature.h:1054
cs_quadrature_hex_1pt_vect
static void cs_quadrature_hex_1pt_vect(double tcur, const cs_real_3_t vb, const double hx, const double hy, const double hz, cs_analytic_func_t *ana, void *input, double results[])
Compute the integral over a orthoganal hexahedron using a barycentric quadrature rule and add it to r...
Definition: cs_quadrature.h:1797
cs_quadrature_get_type_name
const char * cs_quadrature_get_type_name(const cs_quadrature_type_t type)
Return th name associated to a type of quadrature.
Definition: cs_quadrature.cpp:155
cs_quadrature_tria_3pts_scal
static 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[])
Compute the integral over a triangle with a quadrature rule using 3 Gauss points and a unique weight ...
Definition: cs_quadrature.h:934
cs_quadrature_edge_t
void() cs_quadrature_edge_t(const cs_real_3_t v1, const cs_real_3_t v2, double len, cs_real_3_t gpts[], double *weights)
Generic function pointer to compute the quadrature points for an edge from v1 -> v2.
Definition: cs_quadrature.h:109
cs_quadrature_edge_1pt_vect
static 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[])
Compute the integral over an edge using a barycentric quadrature rule and add it to results Case of a...
Definition: cs_quadrature.h:762
cs_quadrature_hex_27pts_vect
static void cs_quadrature_hex_27pts_vect(double tcur, const cs_real_3_t vb, const double hx, const double hy, const double hz, cs_analytic_func_t *ana, void *input, double results[])
Compute the integral over an orthogonal hexahedron with a quadrature rule using 27 Gauss points add i...
Definition: cs_quadrature.h:1879
cs_quadrature_tet_t
void() cs_quadrature_tet_t(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 weights[])
Generic function to compute the quadrature points in a tetrehedra.
Definition: cs_quadrature.h:152
cs_quadrature_edge_integral_t
void() cs_quadrature_edge_integral_t(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[])
Compute the integral over an edge based on a specified quadrature rule and add it to results.
Definition: cs_quadrature.h:176
cs_quadrature_tria_3pts_tens
static 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[])
Compute the integral over a triangle with a quadrature rule using 3 Gauss points and a unique weight ...
Definition: cs_quadrature.h:1277
cs_quadrature_tria_7pts_tens
static 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[])
Compute the integral over a triangle with a quadrature rule using 7 Gauss points and 7 weights and ad...
Definition: cs_quadrature.h:1363
cs_quadrature_tet_4pts
void cs_quadrature_tet_4pts(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 weights[])
Compute the quadrature in a tetrehedra. Exact for 2nd order polynomials (order 3).
Definition: cs_quadrature.cpp:465
cs_quadrature_get_flag
cs_eflag_t cs_quadrature_get_flag(const cs_quadrature_type_t qtype, const cs_flag_t loc)
Get the flags adapted to the given quadrature type qtype and the location on which the quadrature sho...
Definition: cs_quadrature.cpp:616
cs_quadrature_edge_3pts_vect
static 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[])
Compute the integral over an edge with a quadrature rule using 3 Gauss points and weights and add it ...
Definition: cs_quadrature.h:844
cs_quadrature_tria_1pt_tens
static 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[])
Compute the integral over a triangle using a barycentric quadrature rule and add it to results....
Definition: cs_quadrature.h:1229
cs_quadrature_get_hexa_integral
static cs_quadrature_hexa_integral_t * cs_quadrature_get_hexa_integral(int dim, cs_quadrature_type_t qtype)
Retrieve the integral function according to the quadrature type and the stride provided for an hexahe...
Definition: cs_quadrature.h:2146
cs_quadrature_get_tria_integral
static cs_quadrature_tria_integral_t * cs_quadrature_get_tria_integral(int dim, cs_quadrature_type_t qtype)
Retrieve the integral function according to the quadrature type and the stride provided.
Definition: cs_quadrature.h:1981
cs_quadrature_hex_8pts_vect
static void cs_quadrature_hex_8pts_vect(double tcur, const cs_real_3_t vb, const double hx, const double hy, const double hz, cs_analytic_func_t *ana, void *input, double results[])
Compute the integral over an orthoganal hexahedron with a quadrature rule using 8 Gauss points and a ...
Definition: cs_quadrature.h:1837
cs_quadrature_tet_4pts_tens
static 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[])
Compute the integral over a tetrahedron with a quadrature rule using 4 Gauss points and a unique weig...
Definition: cs_quadrature.h:1708