1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5
6 #include <math.h>
7 #include "GaussIntegration.h"
8 #include "GaussLegendre1D.h"
9
brickToTet(double xi,double eta,double zeta,double * r,double * s,double * t,double * J)10 static void brickToTet(double xi, double eta, double zeta, double *r, double *s,
11 double *t, double *J)
12 {
13 double r1, rs1;
14 *r = 0.5e0 * (1.0e0 + xi);
15 r1 = 1.0e0 - (*r);
16 *s = 0.5e0 * (1.0e0 + eta) * r1;
17 rs1 = 1.0e0 - (*r) - (*s);
18 *t = 0.5e0 * (1.0e0 + zeta) * rs1;
19 *J = 0.125e0 * r1 * rs1;
20 }
21
quadToTri(double xi,double eta,double * r,double * s,double * J)22 void quadToTri(double xi, double eta, double *r, double *s, double *J)
23 {
24 double r1;
25 *r = 0.5e0 * (1.0e0 + xi);
26 r1 = 1.0e0 - (*r);
27 *s = 0.5e0 * (1.0e0 + eta) * r1;
28 *J = 0.25e0 * r1;
29 }
30
quadToTriJac(double,double eta)31 double quadToTriJac(double, double eta) { return 0.125e0 * (1.0e0 - eta); }
32
GaussLegendreTet(int n1,int n2,int n3,IntPt * pts)33 int GaussLegendreTet(int n1, int n2, int n3, IntPt *pts)
34 {
35 /* get degenerate n1Xn2Xn3 Gauss-Legendre scheme to integrate over a tet */
36 int i, j, k, index = 0;
37 double *pt1, *pt2, *pt3, *wt1, *wt2, *wt3, dJ;
38
39 gmshGaussLegendre1D(n1, &pt1, &wt1);
40 gmshGaussLegendre1D(n2, &pt2, &wt2);
41 gmshGaussLegendre1D(n3, &pt3, &wt3);
42 for(i = 0; i < n1; i++) {
43 for(j = 0; j < n2; j++) {
44 for(k = 0; k < n3; k++) {
45 brickToTet(pt1[i], pt2[j], pt3[k], &(pts[index].pt[0]),
46 &(pts[index].pt[1]), &(pts[index].pt[2]), &dJ);
47 pts[index++].weight = dJ * wt1[i] * wt2[j] * wt3[k];
48 }
49 }
50 }
51 return index;
52 }
53
GaussLegendreTri(int n1,int n2,IntPt * pts)54 int GaussLegendreTri(int n1, int n2, IntPt *pts)
55 {
56 /* get degenerate n1Xn2 Gauss-Legendre scheme to integrate over a tri */
57 int i, j, index = 0;
58 double *pt1, *pt2, *wt1, *wt2, dJ;
59 // const double two = 2.0000000000000000;
60
61 gmshGaussLegendre1D(n1, &pt1, &wt1);
62 gmshGaussLegendre1D(n2, &pt2, &wt2);
63 for(i = 0; i < n1; i++) {
64 for(j = 0; j < n2; j++) {
65 quadToTri(pt1[i], pt2[j], &(pts[index].pt[0]), &(pts[index].pt[1]), &dJ);
66 pts[index].pt[2] = 0;
67 pts[index++].weight = dJ * wt1[i] * wt2[j];
68 }
69 }
70 return index;
71 }
72