My Project
programmer's documentation
cs_quadrature.h
Go to the documentation of this file.
1 #ifndef __CS_QUADRATURE_H__
2 #define __CS_QUADRATURE_H__
3 
4 /*============================================================================
5  * Routines 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-2019 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_math.h"
37 #include "cs_param.h"
38 
39 /*----------------------------------------------------------------------------*/
40 
42 
43 /*============================================================================
44  * Macro definitions
45  *============================================================================*/
46 
47 /*============================================================================
48  * Type definitions
49  *============================================================================*/
50 
51 typedef enum {
52 
54  CS_QUADRATURE_BARY, /* Value at the barycenter * meas */
55  CS_QUADRATURE_BARY_SUBDIV, /* Value at the barycenter * meas on a sub-mesh */
56  CS_QUADRATURE_HIGHER, /* Unique weight but several Gauss points */
57  CS_QUADRATURE_HIGHEST, /* Specific weight for each Gauss points */
59 
61 
62 /*----------------------------------------------------------------------------*/
73 /*----------------------------------------------------------------------------*/
74 
75 typedef void
77  const cs_real_3_t v2,
78  double len,
79  cs_real_3_t gpts[],
80  double *weights);
81 
82 /*----------------------------------------------------------------------------*/
94 /*----------------------------------------------------------------------------*/
95 
96 typedef void
98  const cs_real_3_t v2,
99  const cs_real_3_t v3,
100  double area,
101  cs_real_3_t gpts[],
102  double *weights);
103 
104 /*----------------------------------------------------------------------------*/
116 /*----------------------------------------------------------------------------*/
117 
118 typedef void
120  const cs_real_3_t v2,
121  const cs_real_3_t v3,
122  const cs_real_3_t v4,
123  double vol,
124  cs_real_3_t gpts[],
125  double weights[]);
126 
127 /*----------------------------------------------------------------------------*/
140 /*----------------------------------------------------------------------------*/
141 
142 typedef void
144  const cs_real_3_t v1,
145  const cs_real_3_t v2,
146  double len,
147  cs_analytic_func_t *ana,
148  void *input,
149  double results[]);
150 
151 /*----------------------------------------------------------------------------*/
165 /*----------------------------------------------------------------------------*/
166 
167 typedef void
169  const cs_real_3_t v1,
170  const cs_real_3_t v2,
171  const cs_real_3_t v3,
172  double area,
173  cs_analytic_func_t *ana,
174  void *input,
175  double results[]);
176 
177 /*----------------------------------------------------------------------------*/
192 /*----------------------------------------------------------------------------*/
193 
194 typedef void
196  const cs_real_3_t v1,
197  const cs_real_3_t v2,
198  const cs_real_3_t v3,
199  const cs_real_3_t v4,
200  double vol,
201  cs_analytic_func_t *ana,
202  void *input,
203  double results[]);
204 
205 
206 /*============================================================================
207  * Public function prototypes
208  *============================================================================*/
209 
210 /*----------------------------------------------------------------------------*/
214 /*----------------------------------------------------------------------------*/
215 
216 void
217 cs_quadrature_setup(void);
218 
219 /*----------------------------------------------------------------------------*/
227 /*----------------------------------------------------------------------------*/
228 
229 const char *
231 
232 /*----------------------------------------------------------------------------*/
243 /*----------------------------------------------------------------------------*/
244 
245 static inline void
247  const cs_real_3_t v2,
248  double len,
249  cs_real_3_t gpts[],
250  double *w)
251 {
252  gpts[0][0] = 0.5*(v1[0] + v2[0]);
253  gpts[0][1] = 0.5*(v1[1] + v2[1]);
254  gpts[0][2] = 0.5*(v1[2] + v2[2]);
255  w[0] = len;
256 }
257 
258 /*----------------------------------------------------------------------------*/
269 /*----------------------------------------------------------------------------*/
270 
271 void
273  const cs_real_3_t v2,
274  double len,
275  cs_real_3_t gpts[],
276  double *w);
277 
278 /*----------------------------------------------------------------------------*/
289 /*----------------------------------------------------------------------------*/
290 
291 void
293  const cs_real_3_t v2,
294  double len,
295  cs_real_3_t gpts[],
296  double w[]);
297 
298 /*----------------------------------------------------------------------------*/
310 /*----------------------------------------------------------------------------*/
311 
312 static inline void
314  const cs_real_3_t v2,
315  const cs_real_3_t v3,
316  double area,
317  cs_real_3_t gpts[],
318  double *w)
319 {
320  gpts[0][0] = cs_math_1ov3 * (v1[0] + v2[0] + v3[0]);
321  gpts[0][1] = cs_math_1ov3 * (v1[1] + v2[1] + v3[1]);
322  gpts[0][2] = cs_math_1ov3 * (v1[2] + v2[2] + v3[2]);
323  w[0] = area;
324 }
325 
326 /*----------------------------------------------------------------------------*/
338 /*----------------------------------------------------------------------------*/
339 
340 void
342  const cs_real_3_t v2,
343  const cs_real_3_t v3,
344  double area,
345  cs_real_3_t gpts[],
346  double *w);
347 
348 /*----------------------------------------------------------------------------*/
360 /*----------------------------------------------------------------------------*/
361 
362 void
364  const cs_real_3_t v2,
365  const cs_real_3_t v3,
366  double area,
367  cs_real_3_t gpts[],
368  double w[]);
369 
370 /*----------------------------------------------------------------------------*/
382 /*----------------------------------------------------------------------------*/
383 
384 void
386  const cs_real_3_t v2,
387  const cs_real_3_t v3,
388  double area,
389  cs_real_3_t gpts[],
390  double w[]);
391 
392 /*----------------------------------------------------------------------------*/
405 /*----------------------------------------------------------------------------*/
406 
407 static inline void
409  const cs_real_3_t v2,
410  const cs_real_3_t v3,
411  const cs_real_3_t v4,
412  double vol,
413  cs_real_3_t gpts[],
414  double weight[])
415 {
416  gpts[0][0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
417  gpts[0][1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
418  gpts[0][2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
419  weight[0] = vol;
420 }
421 
422 /*----------------------------------------------------------------------------*/
435 /*----------------------------------------------------------------------------*/
436 
437 void
439  const cs_real_3_t v2,
440  const cs_real_3_t v3,
441  const cs_real_3_t v4,
442  double vol,
443  cs_real_3_t gpts[],
444  double weights[]);
445 
446 /*----------------------------------------------------------------------------*/
459 /*----------------------------------------------------------------------------*/
460 
461 void
463  const cs_real_3_t v2,
464  const cs_real_3_t v3,
465  const cs_real_3_t v4,
466  double vol,
467  cs_real_3_t gpts[],
468  double weights[]);
469 
470 /*----------------------------------------------------------------------------*/
483 /*----------------------------------------------------------------------------*/
484 
485 void
487  const cs_real_3_t v2,
488  const cs_real_3_t v3,
489  const cs_real_3_t v4,
490  double vol,
491  cs_real_3_t gpts[],
492  double weights[]);
493 
494 /*----------------------------------------------------------------------------*/
508 /*----------------------------------------------------------------------------*/
509 
510 static inline void
512  const cs_real_3_t v1,
513  const cs_real_3_t v2,
514  double len,
515  cs_analytic_func_t *ana,
516  void *input,
517  double results[])
518 {
519  cs_real_3_t xg;
520  double evaluation;
521 
522  // Copied from cs_quadrature_1pt
523  xg[0] = .5 * (v1[0] + v2[0]);
524  xg[1] = .5 * (v1[1] + v2[1]);
525  xg[2] = .5 * (v1[2] + v2[2]);
526 
527  ana(tcur, 1, NULL, xg, false, input, &evaluation);
528 
529  *results += len * evaluation;
530 }
531 
532 /*----------------------------------------------------------------------------*/
546 /*----------------------------------------------------------------------------*/
547 
548 static inline void
550  const cs_real_3_t v1,
551  const cs_real_3_t v2,
552  double len,
553  cs_analytic_func_t *ana,
554  void *input,
555  double results[])
556 {
557  cs_real_3_t gauss_pts[2];
558  double evaluation[2], weights[2];
559 
560  /* Compute Gauss points and its unique weight */
561  cs_quadrature_edge_2pts(v1, v2, len, gauss_pts, weights);
562 
563  ana(tcur, 2, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
564 
565  /* Return results */
566  *results += weights[0] * evaluation[0] + weights[1] * evaluation[1];
567 }
568 
569 /*----------------------------------------------------------------------------*/
583 /*----------------------------------------------------------------------------*/
584 
585 static inline void
587  const cs_real_3_t v1,
588  const cs_real_3_t v2,
589  double len,
590  cs_analytic_func_t *ana,
591  void *input,
592  double results[])
593 {
594  cs_real_3_t gauss_pts[3];
595  double evaluation[3], weights[3];
596 
597  /* Compute Gauss points and its weights */
598  cs_quadrature_edge_3pts(v1, v2, len, gauss_pts, weights);
599 
600  ana(tcur, 3, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
601 
602  *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
603  weights[2] * evaluation[2];
604 }
605 
606 
607 /*----------------------------------------------------------------------------*/
622 /*----------------------------------------------------------------------------*/
623 
624 static inline void
626  const cs_real_3_t v1,
627  const cs_real_3_t v2,
628  const cs_real_3_t v3,
629  double area,
630  cs_analytic_func_t *ana,
631  void *input,
632  double results[])
633 {
634  cs_real_3_t xg;
635  double evaluation;
636 
637  /* Copied from cs_quadrature_1pt */
638  xg[0] = cs_math_1ov3 * (v1[0] + v2[0] + v3[0]);
639  xg[1] = cs_math_1ov3 * (v1[1] + v2[1] + v3[1]);
640  xg[2] = cs_math_1ov3 * (v1[2] + v2[2] + v3[2]);
641 
642  ana(tcur, 1, NULL, xg, false, input, &evaluation);
643 
644  *results += area * evaluation;
645 }
646 
647 /*----------------------------------------------------------------------------*/
662 /*----------------------------------------------------------------------------*/
663 
664 static inline void
666  const cs_real_3_t v1,
667  const cs_real_3_t v2,
668  const cs_real_3_t v3,
669  double area,
670  cs_analytic_func_t *ana,
671  void *input,
672  double results[])
673 {
674  cs_real_3_t gauss_pts[3];
675  double evaluation[3], weights[3];
676 
677  /* Compute Gauss points and its unique weight */
678  cs_quadrature_tria_3pts(v1, v2, v3, area, gauss_pts, weights);
679 
680  ana(tcur, 3, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
681 
682  /* Return results */
683  *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
684  weights[2] * evaluation[2];
685 }
686 
687 /*----------------------------------------------------------------------------*/
702 /*----------------------------------------------------------------------------*/
703 
704 static inline void
706  const cs_real_3_t v1,
707  const cs_real_3_t v2,
708  const cs_real_3_t v3,
709  double area,
710  cs_analytic_func_t *ana,
711  void *input,
712  double results[])
713 {
714  cs_real_3_t gauss_pts[4];
715  double evaluation[4], weights[4];
716 
717  /* Compute Gauss points and its weights */
718  cs_quadrature_tria_4pts(v1, v2, v3, area, gauss_pts, weights);
719 
720  ana(tcur, 4, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
721 
722  *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
723  weights[2] * evaluation[2] + weights[3] * evaluation[3];
724 }
725 
726 /*----------------------------------------------------------------------------*/
741 /*----------------------------------------------------------------------------*/
742 
743 static inline void
745  const cs_real_3_t v1,
746  const cs_real_3_t v2,
747  const cs_real_3_t v3,
748  double area,
749  cs_analytic_func_t *ana,
750  void *input,
751  double results[])
752 {
753  cs_real_3_t xg;
754  double evaluation[3];
755 
756  /* Copied from cs_quadrature_1pt */
757  xg[0] = cs_math_1ov3 * (v1[0] + v2[0] + v3[0]);
758  xg[1] = cs_math_1ov3 * (v1[1] + v2[1] + v3[1]);
759  xg[2] = cs_math_1ov3 * (v1[2] + v2[2] + v3[2]);
760 
761  ana(tcur, 1, NULL, xg, false, input, evaluation);
762 
763  results[0] += area * evaluation[0];
764  results[1] += area * evaluation[1];
765  results[2] += area * evaluation[2];
766 }
767 
768 /*----------------------------------------------------------------------------*/
783 /*----------------------------------------------------------------------------*/
784 
785 static inline void
787  const cs_real_3_t v1,
788  const cs_real_3_t v2,
789  const cs_real_3_t v3,
790  double area,
791  cs_analytic_func_t *ana,
792  void *input,
793  double results[])
794 {
795  cs_real_3_t gauss_pts[3];
796  double evaluation[3*3], weights[3];
797 
798  /* Compute Gauss points and its unique weight */
799  cs_quadrature_tria_3pts(v1, v2, v3, area, gauss_pts, weights);
800 
801  ana(tcur, 3, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
802 
803  for (int p = 0; p < 3; p++) {
804  results[0] += weights[p] * evaluation[3*p ];
805  results[1] += weights[p] * evaluation[3*p+1];
806  results[2] += weights[p] * evaluation[3*p+2];
807  }
808 }
809 
810 /*----------------------------------------------------------------------------*/
825 /*----------------------------------------------------------------------------*/
826 
827 static inline void
829  const cs_real_3_t v1,
830  const cs_real_3_t v2,
831  const cs_real_3_t v3,
832  double area,
833  cs_analytic_func_t *ana,
834  void *input,
835  double results[])
836 {
837  cs_real_3_t gauss_pts[4];
838  double evaluation[3*4], weights[4];
839 
840  /* Compute Gauss points and its weights */
841  cs_quadrature_tria_4pts(v1, v2, v3, area, gauss_pts, weights);
842 
843  ana(tcur, 4, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
844 
845  for (int p = 0; p < 4; p++) {
846  results[0] += weights[p] * evaluation[3*p ];
847  results[1] += weights[p] * evaluation[3*p+1];
848  results[2] += weights[p] * evaluation[3*p+2];
849  }
850 }
851 
852 /*----------------------------------------------------------------------------*/
867 /*----------------------------------------------------------------------------*/
868 
869 static inline void
871  const cs_real_3_t v1,
872  const cs_real_3_t v2,
873  const cs_real_3_t v3,
874  double area,
875  cs_analytic_func_t *ana,
876  void *input,
877  double results[])
878 {
879  cs_real_3_t xg;
880  double evaluation[9];
881 
882  /* Copied from cs_quadrature_1pt */
883  xg[0] = cs_math_1ov3 * (v1[0] + v2[0] + v3[0]);
884  xg[1] = cs_math_1ov3 * (v1[1] + v2[1] + v3[1]);
885  xg[2] = cs_math_1ov3 * (v1[2] + v2[2] + v3[2]);
886 
887  ana(tcur, 1, NULL, xg, false, input, evaluation);
888 
889  for (short int ij = 0; ij < 9; ij++)
890  results[ij] += area * evaluation[ij];
891 }
892 
893 /*----------------------------------------------------------------------------*/
908 /*----------------------------------------------------------------------------*/
909 
910 static inline void
912  const cs_real_3_t v1,
913  const cs_real_3_t v2,
914  const cs_real_3_t v3,
915  double area,
916  cs_analytic_func_t *ana,
917  void *input,
918  double results[])
919 {
920  cs_real_3_t gauss_pts[3];
921  double evaluation[9*3], weights[3];
922 
923  /* Compute Gauss points and its unique weight */
924  cs_quadrature_tria_3pts(v1, v2, v3, area, gauss_pts, weights);
925 
926  ana(tcur, 3, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
927 
928  for (int p = 0; p < 3; p++) {
929  const double wp = weights[p];
930  double *eval_p = evaluation + 9*p;
931  for (short int ij = 0; ij < 9; ij++)
932  results[ij] += wp * eval_p[ij];
933  }
934 }
935 
936 /*----------------------------------------------------------------------------*/
951 /*----------------------------------------------------------------------------*/
952 
953 static inline void
955  const cs_real_3_t v1,
956  const cs_real_3_t v2,
957  const cs_real_3_t v3,
958  double area,
959  cs_analytic_func_t *ana,
960  void *input,
961  double results[])
962 {
963  cs_real_3_t gauss_pts[4];
964  double evaluation[9*4], weights[4];
965 
966  /* Compute Gauss points and its weights */
967  cs_quadrature_tria_4pts(v1, v2, v3, area, gauss_pts, weights);
968 
969  ana(tcur, 4, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
970 
971  for (int p = 0; p < 4; p++) {
972  const double wp = weights[p];
973  double *eval_p = evaluation + 9*p;
974  for (short int ij = 0; ij < 9; ij++)
975  results[ij] += wp * eval_p[ij];
976  }
977 }
978 
979 /*----------------------------------------------------------------------------*/
995 /*----------------------------------------------------------------------------*/
996 
997 static inline void
999  const cs_real_3_t v1,
1000  const cs_real_3_t v2,
1001  const cs_real_3_t v3,
1002  const cs_real_3_t v4,
1003  double vol,
1004  cs_analytic_func_t *ana,
1005  void *input,
1006  double results[])
1007 {
1008  cs_real_3_t xg;
1009  double evaluation;
1010 
1011  /* Copied from cs_quadrature_tet_1pt */
1012  xg[0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
1013  xg[1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
1014  xg[2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
1015 
1016  ana(tcur, 1, NULL, xg, false, input, &evaluation);
1017 
1018  *results += vol * evaluation;
1019 }
1020 
1021 /*----------------------------------------------------------------------------*/
1037 /*----------------------------------------------------------------------------*/
1038 
1039 static inline void
1041  const cs_real_3_t v1,
1042  const cs_real_3_t v2,
1043  const cs_real_3_t v3,
1044  const cs_real_3_t v4,
1045  double vol,
1046  cs_analytic_func_t *ana,
1047  void *input,
1048  double results[])
1049 {
1050  cs_real_3_t gauss_pts[4];
1051  double evaluation[4], weights[4];
1052 
1053  /* Compute Gauss points and its unique weight */
1054  cs_quadrature_tet_4pts(v1, v2, v3, v4, vol, gauss_pts, weights);
1055 
1056  ana(tcur, 4, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1057 
1058  /* Return results */
1059  *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
1060  weights[2] * evaluation[2] + weights[3] * evaluation[3];
1061 }
1062 
1063 /*----------------------------------------------------------------------------*/
1079 /*----------------------------------------------------------------------------*/
1080 
1081 static inline void
1083  const cs_real_3_t v1,
1084  const cs_real_3_t v2,
1085  const cs_real_3_t v3,
1086  const cs_real_3_t v4,
1087  double vol,
1088  cs_analytic_func_t *ana,
1089  void *input,
1090  double results[])
1091 {
1092  cs_real_3_t gauss_pts[5];
1093  double evaluation[5], weights[5];
1094 
1095  /* Compute Gauss points and its weights */
1096  cs_quadrature_tet_5pts(v1, v2, v3, v4, vol, gauss_pts, weights);
1097 
1098  ana(tcur, 5, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1099 
1100  *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
1101  weights[2] * evaluation[2] + weights[3] * evaluation[3] +
1102  weights[4] * evaluation[4];
1103 }
1104 
1105 /*----------------------------------------------------------------------------*/
1121 /*----------------------------------------------------------------------------*/
1122 
1123 static inline void
1125  const cs_real_3_t v1,
1126  const cs_real_3_t v2,
1127  const cs_real_3_t v3,
1128  const cs_real_3_t v4,
1129  double vol,
1130  cs_analytic_func_t *ana,
1131  void *input,
1132  double results[])
1133 {
1134  cs_real_3_t xg;
1135  double evaluation[3];
1136 
1137  /* Copied from cs_quadrature_tet_1pt */
1138  xg[0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
1139  xg[1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
1140  xg[2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
1141 
1142  ana(tcur, 1, NULL, xg, false, input, evaluation);
1143 
1144  results[0] += vol * evaluation[0];
1145  results[1] += vol * evaluation[1];
1146  results[2] += vol * evaluation[2];
1147 }
1148 
1149 /*----------------------------------------------------------------------------*/
1165 /*----------------------------------------------------------------------------*/
1166 
1167 static inline void
1169  const cs_real_3_t v1,
1170  const cs_real_3_t v2,
1171  const cs_real_3_t v3,
1172  const cs_real_3_t v4,
1173  double vol,
1174  cs_analytic_func_t *ana,
1175  void *input,
1176  double results[])
1177 {
1178  cs_real_3_t gauss_pts[4];
1179  double evaluation[3*4], weights[4];
1180 
1181  /* Compute Gauss points and its unique weight */
1182  cs_quadrature_tet_4pts(v1, v2, v3, v4, vol, gauss_pts, weights);
1183 
1184  ana(tcur, 4, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1185 
1186  for (int p = 0; p < 4; p++) {
1187  results[0] += weights[p] * evaluation[3*p ];
1188  results[1] += weights[p] * evaluation[3*p+1];
1189  results[2] += weights[p] * evaluation[3*p+2];
1190  }
1191 }
1192 
1193 /*----------------------------------------------------------------------------*/
1209 /*----------------------------------------------------------------------------*/
1210 
1211 static inline void
1213  const cs_real_3_t v1,
1214  const cs_real_3_t v2,
1215  const cs_real_3_t v3,
1216  const cs_real_3_t v4,
1217  double vol,
1218  cs_analytic_func_t *ana,
1219  void *input,
1220  double results[])
1221 {
1222  cs_real_3_t gauss_pts[5];
1223  double evaluation[3*5], weights[5];
1224 
1225  /* Compute Gauss points and its weights */
1226  cs_quadrature_tet_5pts(v1, v2, v3, v4, vol, gauss_pts, weights);
1227 
1228  ana(tcur, 5, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1229 
1230  for (int p = 0; p < 5; p++) {
1231  results[0] += weights[p] * evaluation[3*p ];
1232  results[1] += weights[p] * evaluation[3*p+1];
1233  results[2] += weights[p] * evaluation[3*p+2];
1234  }
1235 }
1236 
1237 /*----------------------------------------------------------------------------*/
1253 /*----------------------------------------------------------------------------*/
1254 
1255 static inline void
1257  const cs_real_3_t v1,
1258  const cs_real_3_t v2,
1259  const cs_real_3_t v3,
1260  const cs_real_3_t v4,
1261  double vol,
1262  cs_analytic_func_t *ana,
1263  void *input,
1264  double results[])
1265 {
1266  cs_real_3_t xg;
1267  double evaluation[9];
1268 
1269  /* Copied from cs_quadrature_tet_1pt */
1270  xg[0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
1271  xg[1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
1272  xg[2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
1273 
1274  ana(tcur, 1, NULL, xg, false, input, evaluation);
1275 
1276  for (short int ij = 0; ij < 9; ij++)
1277  results[ij] += vol * evaluation[ij];
1278 }
1279 
1280 /*----------------------------------------------------------------------------*/
1296 /*----------------------------------------------------------------------------*/
1297 
1298 static inline void
1300  const cs_real_3_t v1,
1301  const cs_real_3_t v2,
1302  const cs_real_3_t v3,
1303  const cs_real_3_t v4,
1304  double vol,
1305  cs_analytic_func_t *ana,
1306  void *input,
1307  double results[])
1308 {
1309  cs_real_3_t gauss_pts[4];
1310  double evaluation[9*4], weights[4];
1311 
1312  /* Compute Gauss points and its unique weight */
1313  cs_quadrature_tet_4pts(v1, v2, v3, v4, vol, gauss_pts, weights);
1314 
1315  ana(tcur, 4, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1316 
1317  for (int p = 0; p < 4; p++) {
1318  const double wp = weights[p];
1319  double *eval_p = evaluation + 9*p;
1320  for (short int ij = 0; ij < 9; ij++)
1321  results[ij] += wp * eval_p[ij];
1322  }
1323 }
1324 
1325 /*----------------------------------------------------------------------------*/
1341 /*----------------------------------------------------------------------------*/
1342 
1343 static inline void
1345  const cs_real_3_t v1,
1346  const cs_real_3_t v2,
1347  const cs_real_3_t v3,
1348  const cs_real_3_t v4,
1349  double vol,
1350  cs_analytic_func_t *ana,
1351  void *input,
1352  double results[])
1353 {
1354  cs_real_3_t gauss_pts[5];
1355  double evaluation[9*5], weights[5];
1356 
1357  /* Compute Gauss points and its weights */
1358  cs_quadrature_tet_5pts(v1, v2, v3, v4, vol, gauss_pts, weights);
1359 
1360  ana(tcur, 5, NULL, (const cs_real_t *)gauss_pts, false, input, evaluation);
1361 
1362  for (int p = 0; p < 5; p++) {
1363  const double wp = weights[p];
1364  double *eval_p = evaluation + 9*p;
1365  for (short int ij = 0; ij < 9; ij++)
1366  results[ij] += wp * eval_p[ij];
1367  }
1368 }
1369 
1370 /*----------------------------------------------------------------------------*/
1380 /*----------------------------------------------------------------------------*/
1381 
1382 static inline cs_quadrature_tria_integral_t *
1384  cs_quadrature_type_t qtype)
1385 {
1386  switch (dim) {
1387 
1388  case 1: /* Scalar-valued integral */
1389 
1390  switch (qtype) {
1391 
1392  case CS_QUADRATURE_BARY:
1395  case CS_QUADRATURE_HIGHER:
1397  case CS_QUADRATURE_HIGHEST:
1399 
1400  default:
1401  bft_error(__FILE__, __LINE__, 0,
1402  " %s: Invalid quadrature type\n", __func__);
1403  }
1404  break;
1405 
1406  case 3: /* Vector-valued case */
1407 
1408  switch (qtype) {
1409 
1410  case CS_QUADRATURE_BARY:
1413  case CS_QUADRATURE_HIGHER:
1415  case CS_QUADRATURE_HIGHEST:
1417 
1418  default:
1419  bft_error(__FILE__, __LINE__, 0,
1420  " %s: Invalid quadrature type\n", __func__);
1421  }
1422  break;
1423 
1424  case 9: /* Tensor-valued case */
1425 
1426  switch (qtype) {
1427 
1428  case CS_QUADRATURE_BARY:
1431  case CS_QUADRATURE_HIGHER:
1433  case CS_QUADRATURE_HIGHEST:
1435 
1436  default:
1437  bft_error(__FILE__, __LINE__, 0,
1438  " %s: Invalid quadrature type\n", __func__);
1439  }
1440  break;
1441 
1442  default:
1443  bft_error(__FILE__, __LINE__, 0,
1444  " %s: Invalid dimension value %d. Only 1, 3 and 9 are valid.\n",
1445  __func__, dim);
1446 
1447  } /* switch on dim */
1448 
1449  return NULL; /* Should not go to this stage */
1450 }
1451 
1452 /*----------------------------------------------------------------------------*/
1462 /*----------------------------------------------------------------------------*/
1463 
1464 static inline cs_quadrature_tetra_integral_t *
1466  cs_quadrature_type_t qtype)
1467 {
1468  switch (dim) {
1469 
1470  case 1: /* Scalar-valued case */
1471 
1472  switch (qtype) {
1473 
1474  case CS_QUADRATURE_BARY:
1477  case CS_QUADRATURE_HIGHER:
1479  case CS_QUADRATURE_HIGHEST:
1481 
1482  default:
1483  bft_error(__FILE__, __LINE__, 0,
1484  " %s: Invalid quadrature type\n", __func__);
1485  }
1486  break;
1487 
1488  case 3: /* Vector-valued case */
1489 
1490  switch (qtype) {
1491 
1492  case CS_QUADRATURE_BARY:
1495  case CS_QUADRATURE_HIGHER:
1497  case CS_QUADRATURE_HIGHEST:
1499 
1500  default:
1501  bft_error(__FILE__, __LINE__, 0,
1502  " %s: Invalid quadrature type\n", __func__);
1503  }
1504  break;
1505 
1506  case 9: /* Tensor-valued case */
1507 
1508  switch (qtype) {
1509 
1510  case CS_QUADRATURE_BARY:
1513  case CS_QUADRATURE_HIGHER:
1515  case CS_QUADRATURE_HIGHEST:
1517 
1518  default:
1519  bft_error(__FILE__, __LINE__, 0,
1520  " %s: Invalid quadrature type\n", __func__);
1521  }
1522  break;
1523 
1524  default:
1525  bft_error(__FILE__, __LINE__, 0,
1526  " %s: Invalid dimension value %d. Only 1, 3 and 9 are valid.\n",
1527  __func__, dim);
1528 
1529  } /* Switch on dim */
1530 
1531  /* Avoid no return warning */
1532  return NULL;
1533 }
1534 
1535 /*----------------------------------------------------------------------------*/
1547 /*----------------------------------------------------------------------------*/
1548 
1549 cs_flag_t
1551  const cs_flag_t loc);
1552 
1553 /*----------------------------------------------------------------------------*/
1554 
1556 
1557 #endif /* __CS_QUADRATURE_H__ */
CS_QUADRATURE_HIGHER
Definition: cs_quadrature.h:56
CS_QUADRATURE_N_TYPES
Definition: cs_quadrature.h:58
input
static int input(void)
cs_defs.h
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:1212
bft_error.h
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:549
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:911
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.c:200
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:246
len
size_t len
Definition: mei_scanner.c:560
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:1383
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:1465
cs_real_3_t
cs_real_t cs_real_3_t[3]
vector of 3 floating-point values
Definition: cs_defs.h:315
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:744
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:1256
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.c:308
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:870
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:1124
END_C_DECLS
#define END_C_DECLS
Definition: cs_defs.h:468
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:302
BEGIN_C_DECLS
#define BEGIN_C_DECLS
Definition: cs_defs.h:467
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:954
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:828
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:76
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.c:193
CS_QUADRATURE_NONE
Definition: cs_quadrature.h:53
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:1299
CS_QUADRATURE_BARY
Definition: cs_quadrature.h:54
cs_math.h
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:625
cs_quadrature_get_flag
cs_flag_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.c:498
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.c:351
CS_QUADRATURE_HIGHEST
Definition: cs_quadrature.h:57
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:408
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:705
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:168
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:1344
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:195
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_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:1082
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.c:435
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:511
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:1040
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:97
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 compact, void *input, cs_real_t *retval)
Generic function pointer for an analytic function elt_ids is optional. If not NULL,...
Definition: cs_param.h:66
cs_math_1ov3
const cs_real_t cs_math_1ov3
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:586
cs_quadrature_setup
void cs_quadrature_setup(void)
Compute constant weights for all quadratures used.
Definition: cs_quadrature.c:107
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:143
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_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:665
cs_flag_t
unsigned short int cs_flag_t
Definition: cs_defs.h:304
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:998
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:1168
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:119
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:786
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.c:391
cs_quadrature_type_t
cs_quadrature_type_t
Definition: cs_quadrature.h:51
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:313
CS_QUADRATURE_BARY_SUBDIV
Definition: cs_quadrature.h:55
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.c:149
type
void const cs_int_t * type
Definition: cs_measures_util.h:425
cs_param.h
cs_base.h
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.c:271
p
Definition: cs_field_pointer.h:67