1 #ifndef __CS_QUADRATURE_H__
2 #define __CS_QUADRATURE_H__
253 gpts[0][0] = 0.5*(v1[0] + v2[0]);
254 gpts[0][1] = 0.5*(v1[1] + v2[1]);
255 gpts[0][2] = 0.5*(v1[2] + v2[2]);
417 gpts[0][0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
418 gpts[0][1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
419 gpts[0][2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
512 cs_quadrature_edge_1pt_scal(
double tcur,
524 xg[0] = .5 * (v1[0] + v2[0]);
525 xg[1] = .5 * (v1[1] + v2[1]);
526 xg[2] = .5 * (v1[2] + v2[2]);
529 ana(tcur, 1, NULL, xg,
false, input, &feval);
532 *results += len * feval;
552 cs_quadrature_edge_2pts_scal(
double tcur,
561 double feval[2], weights[2];
567 ana(tcur, 2, NULL, (
const cs_real_t *)gauss_pts,
false, input, feval);
570 *results += weights[0] * feval[0] + weights[1] * feval[1];
590 cs_quadrature_edge_3pts_scal(
double tcur,
599 double feval[3], weights[3];
605 ana(tcur, 3, NULL, (
const cs_real_t *)gauss_pts,
false, input, feval);
608 *results += weights[0]*feval[0] + weights[1]*feval[1] + weights[2]*feval[2];
628 cs_quadrature_edge_1pt_vect(
double tcur,
640 xg[0] = .5 * (v1[0] + v2[0]);
641 xg[1] = .5 * (v1[1] + v2[1]);
642 xg[2] = .5 * (v1[2] + v2[2]);
645 ana(tcur, 1, NULL, xg,
false, input, feval);
648 results[0] += len * feval[0];
649 results[1] += len * feval[1];
650 results[2] += len * feval[2];
670 cs_quadrature_edge_2pts_vect(
double tcur,
679 double feval[6], weights[2];
685 ana(tcur, 2, NULL, (
const cs_real_t *)gauss_pts,
false, input, feval);
688 results[0] += weights[0] * feval[0] + weights[1] * feval[3];
689 results[1] += weights[0] * feval[1] + weights[1] * feval[4];
690 results[2] += weights[0] * feval[2] + weights[1] * feval[5];
710 cs_quadrature_edge_3pts_vect(
double tcur,
719 double feval[9], weights[3];
725 ana(tcur, 3, NULL, (
const cs_real_t *)gauss_pts,
false, input, feval);
728 for (
int p = 0;
p < 3;
p++) {
729 results[0] += weights[
p] * feval[3*
p ];
730 results[1] += weights[
p] * feval[3*
p+1];
731 results[2] += weights[
p] * feval[3*
p+2];
753 cs_quadrature_tria_1pt_scal(
double tcur,
770 ana(tcur, 1, NULL, xg,
false, input, &evaluation);
772 *results += area * evaluation;
793 cs_quadrature_tria_3pts_scal(
double tcur,
803 double evaluation[3], weights[3];
808 ana(tcur, 3, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
811 *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
812 weights[2] * evaluation[2];
833 cs_quadrature_tria_4pts_scal(
double tcur,
843 double evaluation[4], weights[4];
848 ana(tcur, 4, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
850 *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
851 weights[2] * evaluation[2] + weights[3] * evaluation[3];
872 cs_quadrature_tria_7pts_scal(
double tcur,
882 double evaluation[7], weights[7];
887 ana(tcur, 7, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
889 *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
890 weights[2] * evaluation[2] + weights[3] * evaluation[3] +
891 weights[4] * evaluation[4] + weights[5] * evaluation[5] +
892 weights[6] * evaluation[6] ;
913 cs_quadrature_tria_1pt_vect(
double tcur,
923 double evaluation[3];
930 ana(tcur, 1, NULL, xg,
false, input, evaluation);
932 results[0] += area * evaluation[0];
933 results[1] += area * evaluation[1];
934 results[2] += area * evaluation[2];
955 cs_quadrature_tria_3pts_vect(
double tcur,
965 double evaluation[3*3], weights[3];
970 ana(tcur, 3, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
972 for (
int p = 0;
p < 3;
p++) {
973 results[0] += weights[
p] * evaluation[3*
p ];
974 results[1] += weights[
p] * evaluation[3*
p+1];
975 results[2] += weights[
p] * evaluation[3*
p+2];
997 cs_quadrature_tria_4pts_vect(
double tcur,
1007 double evaluation[3*4], weights[4];
1012 ana(tcur, 4, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
1014 for (
int p = 0;
p < 4;
p++) {
1015 results[0] += weights[
p] * evaluation[3*
p ];
1016 results[1] += weights[
p] * evaluation[3*
p+1];
1017 results[2] += weights[
p] * evaluation[3*
p+2];
1039 cs_quadrature_tria_7pts_vect(
double tcur,
1049 double evaluation[3*7], weights[7];
1054 ana(tcur, 7, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
1056 for (
int p = 0;
p < 7;
p++) {
1057 results[0] += weights[
p] * evaluation[3*
p ];
1058 results[1] += weights[
p] * evaluation[3*
p+1];
1059 results[2] += weights[
p] * evaluation[3*
p+2];
1081 cs_quadrature_tria_1pt_tens(
double tcur,
1091 double evaluation[9];
1098 ana(tcur, 1, NULL, xg,
false, input, evaluation);
1100 for (
short int ij = 0; ij < 9; ij++)
1101 results[ij] += area * evaluation[ij];
1122 cs_quadrature_tria_3pts_tens(
double tcur,
1132 double evaluation[9*3], weights[3];
1137 ana(tcur, 3, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
1139 for (
int p = 0;
p < 3;
p++) {
1140 const double wp = weights[
p];
1141 double *eval_p = evaluation + 9*
p;
1142 for (
short int ij = 0; ij < 9; ij++)
1143 results[ij] += wp * eval_p[ij];
1165 cs_quadrature_tria_4pts_tens(
double tcur,
1175 double evaluation[9*4], weights[4];
1180 ana(tcur, 4, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
1182 for (
int p = 0;
p < 4;
p++) {
1183 const double wp = weights[
p];
1184 double *eval_p = evaluation + 9*
p;
1185 for (
short int ij = 0; ij < 9; ij++)
1186 results[ij] += wp * eval_p[ij];
1208 cs_quadrature_tria_7pts_tens(
double tcur,
1218 double evaluation[9*7], weights[7];
1223 ana(tcur, 7, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
1225 for (
int p = 0;
p < 7;
p++) {
1226 const double wp = weights[
p];
1227 double *eval_p = evaluation + 9*
p;
1228 for (
short int ij = 0; ij < 9; ij++)
1229 results[ij] += wp * eval_p[ij];
1252 cs_quadrature_tet_1pt_scal(
double tcur,
1266 xg[0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
1267 xg[1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
1268 xg[2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
1270 ana(tcur, 1, NULL, xg,
false, input, &evaluation);
1272 *results += vol * evaluation;
1294 cs_quadrature_tet_4pts_scal(
double tcur,
1305 double evaluation[4], weights[4];
1310 ana(tcur, 4, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
1313 *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
1314 weights[2] * evaluation[2] + weights[3] * evaluation[3];
1336 cs_quadrature_tet_5pts_scal(
double tcur,
1347 double evaluation[5], weights[5];
1352 ana(tcur, 5, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
1354 *results += weights[0] * evaluation[0] + weights[1] * evaluation[1] +
1355 weights[2] * evaluation[2] + weights[3] * evaluation[3] +
1356 weights[4] * evaluation[4];
1378 cs_quadrature_tet_1pt_vect(
double tcur,
1389 double evaluation[3];
1392 xg[0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
1393 xg[1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
1394 xg[2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
1396 ana(tcur, 1, NULL, xg,
false, input, evaluation);
1398 results[0] += vol * evaluation[0];
1399 results[1] += vol * evaluation[1];
1400 results[2] += vol * evaluation[2];
1422 cs_quadrature_tet_4pts_vect(
double tcur,
1433 double evaluation[3*4], weights[4];
1438 ana(tcur, 4, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
1440 for (
int p = 0;
p < 4;
p++) {
1441 results[0] += weights[
p] * evaluation[3*
p ];
1442 results[1] += weights[
p] * evaluation[3*
p+1];
1443 results[2] += weights[
p] * evaluation[3*
p+2];
1466 cs_quadrature_tet_5pts_vect(
double tcur,
1477 double evaluation[3*5], weights[5];
1482 ana(tcur, 5, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
1484 for (
int p = 0;
p < 5;
p++) {
1485 results[0] += weights[
p] * evaluation[3*
p ];
1486 results[1] += weights[
p] * evaluation[3*
p+1];
1487 results[2] += weights[
p] * evaluation[3*
p+2];
1510 cs_quadrature_tet_1pt_tens(
double tcur,
1521 double evaluation[9];
1524 xg[0] = 0.25 * (v1[0] + v2[0] + v3[0] + v4[0]);
1525 xg[1] = 0.25 * (v1[1] + v2[1] + v3[1] + v4[1]);
1526 xg[2] = 0.25 * (v1[2] + v2[2] + v3[2] + v4[2]);
1528 ana(tcur, 1, NULL, xg,
false, input, evaluation);
1530 for (
short int ij = 0; ij < 9; ij++)
1531 results[ij] += vol * evaluation[ij];
1553 cs_quadrature_tet_4pts_tens(
double tcur,
1564 double evaluation[9*4], weights[4];
1569 ana(tcur, 4, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
1571 for (
int p = 0;
p < 4;
p++) {
1572 const double wp = weights[
p];
1573 double *eval_p = evaluation + 9*
p;
1574 for (
short int ij = 0; ij < 9; ij++)
1575 results[ij] += wp * eval_p[ij];
1598 cs_quadrature_tet_5pts_tens(
double tcur,
1609 double evaluation[9*5], weights[5];
1614 ana(tcur, 5, NULL, (
const cs_real_t *)gauss_pts,
false, input, evaluation);
1616 for (
int p = 0;
p < 5;
p++) {
1617 const double wp = weights[
p];
1618 double *eval_p = evaluation + 9*
p;
1619 for (
short int ij = 0; ij < 9; ij++)
1620 results[ij] += wp * eval_p[ij];
1638 cs_quadrature_get_edge_integral(
int dim,
1649 return cs_quadrature_edge_1pt_scal;
1651 return cs_quadrature_edge_2pts_scal;
1653 return cs_quadrature_edge_3pts_scal;
1657 " %s: Invalid quadrature type\n", __func__);
1667 return cs_quadrature_edge_1pt_vect;
1669 return cs_quadrature_edge_2pts_vect;
1671 return cs_quadrature_edge_3pts_vect;
1675 " %s: Invalid quadrature type\n", __func__);
1681 " %s: Invalid dimension value %d. Only 1 and 3 are valid.\n",
1702 cs_quadrature_get_tria_integral(
int dim,
1713 return cs_quadrature_tria_1pt_scal;
1715 return cs_quadrature_tria_4pts_scal;
1717 return cs_quadrature_tria_7pts_scal;
1721 " %s: Invalid quadrature type\n", __func__);
1731 return cs_quadrature_tria_1pt_vect;
1733 return cs_quadrature_tria_4pts_vect;
1735 return cs_quadrature_tria_7pts_vect;
1739 " %s: Invalid quadrature type\n", __func__);
1749 return cs_quadrature_tria_1pt_tens;
1751 return cs_quadrature_tria_4pts_tens;
1753 return cs_quadrature_tria_7pts_tens;
1757 " %s: Invalid quadrature type\n", __func__);
1763 " %s: Invalid dimension value %d. Only 1, 3 and 9 are valid.\n",
1784 cs_quadrature_get_tetra_integral(
int dim,
1795 return cs_quadrature_tet_1pt_scal;
1797 return cs_quadrature_tet_4pts_scal;
1799 return cs_quadrature_tet_5pts_scal;
1803 " %s: Invalid quadrature type\n", __func__);
1813 return cs_quadrature_tet_1pt_vect;
1815 return cs_quadrature_tet_4pts_vect;
1817 return cs_quadrature_tet_5pts_vect;
1821 " %s: Invalid quadrature type\n", __func__);
1831 return cs_quadrature_tet_1pt_tens;
1833 return cs_quadrature_tet_4pts_tens;
1835 return cs_quadrature_tet_5pts_tens;
1839 " %s: Invalid quadrature type\n", __func__);
1845 " %s: Invalid dimension value %d. Only 1, 3 and 9 are valid.\n",