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 <vector>
7 #include "GmshMessage.h"
8 #include "GaussIntegration.h"
9 #include "GaussLegendre1D.h"
10
11 IntPt GQQ1[1] = {{{0.0, 0.0, 0.0}, 4.0}};
12
13 const double xq3[3] = {0.816496580928, -0.408248290464, -0.408248290464};
14 const double yq3[3] = {0.0, 0.840896415255, -0.840896415255};
15 const double pq3[3] = {1.3333333333333, 1.3333333333333, 1.3333333333333};
16 IntPt GQQ3[3] = {{{xq3[0], yq3[0], 0}, pq3[0]},
17 {{xq3[1], yq3[1], 0}, pq3[1]},
18 {{xq3[2], yq3[2], 0}, pq3[2]}};
19
20 const double xq4[4] = {0.577350269189626, -0.577350269189626, 0.577350269189626,
21 -0.577350269189626};
22 const double yq4[4] = {0.577350269189626, 0.577350269189626, -0.577350269189626,
23 -0.577350269189626};
24 const double pq4[4] = {1., 1., 1., 1.};
25 IntPt GQQ4[4] = {
26 {{xq4[0], yq4[0], 0}, pq4[0]},
27 {{xq4[1], yq4[1], 0}, pq4[1]},
28 {{xq4[2], yq4[2], 0}, pq4[2]},
29 {{xq4[3], yq4[3], 0}, pq4[3]},
30 };
31
32 const double xq7[7] = {0.,
33 0.,
34 0.,
35 0.7745966692414834,
36 0.7745966692414834,
37 -0.7745966692414834,
38 -0.7745966692414834};
39 const double yq7[7] = {0.,
40 0.9660917830792959,
41 -0.9660917830792959,
42 0.7745966692414834,
43 -0.7745966692414834,
44 0.7745966692414834,
45 -0.7745966692414834};
46 const double pq7[7] = {1.1428571428571428, 0.31746031746031744,
47 0.31746031746031744, 0.5555555555555556,
48 0.5555555555555556, 0.5555555555555556,
49 0.5555555555555556};
50 IntPt GQQ7[7] = {{{xq7[0], yq7[0], 0}, pq7[0]}, {{xq7[1], yq7[1], 0}, pq7[1]},
51 {{xq7[2], yq7[2], 0}, pq7[2]}, {{xq7[3], yq7[3], 0}, pq7[3]},
52 {{xq7[4], yq7[4], 0}, pq7[4]}, {{xq7[5], yq7[5], 0}, pq7[5]},
53 {{xq7[6], yq7[6], 0}, pq7[6]}};
54
55 const double a9 = 0.774596669241483;
56 const double z = 0.0;
57 const double xq9[9] = {-a9, z, a9, -a9, z, a9, -a9, z, a9};
58 const double yq9[9] = {-a9, -a9, -a9, z, z, z, a9, a9, a9};
59 const double pb2 = 0.888888888888889 * 0.888888888888889;
60 const double pc2 = 0.555555555555556 * 0.555555555555556;
61 const double pbc = 0.555555555555556 * 0.888888888888889;
62 const double pq9[9] = {pc2, pbc, pc2, pbc, pb2, pbc, pc2, pbc, pc2};
63
64 IntPt GQQ9[9] = {{{xq9[0], yq9[0], 0}, pq9[0]}, {{xq9[1], yq9[1], 0}, pq9[1]},
65 {{xq9[2], yq9[2], 0}, pq9[2]}, {{xq9[3], yq9[3], 0}, pq9[3]},
66 {{xq9[4], yq9[4], 0}, pq9[4]}, {{xq9[5], yq9[5], 0}, pq9[5]},
67 {{xq9[6], yq9[6], 0}, pq9[6]}, {{xq9[7], yq9[7], 0}, pq9[7]},
68 {{xq9[8], yq9[8], 0}, pq9[8]}};
69
70 const double a16 = 0.861136311594053;
71 const double b16 = 0.339981043584856;
72 const double xq16[16] = {-a16, -b16, b16, a16, -a16, -b16, b16, a16,
73 -a16, -b16, b16, a16, -a16, -b16, b16, a16};
74 const double yq16[16] = {-a16, -a16, -a16, -a16, -b16, -b16, -b16, -b16,
75 b16, b16, b16, b16, a16, a16, a16, a16};
76 const double pe2 = 0.347854845137454 * 0.347854845137454;
77 const double pf2 = 0.652145154862546 * 0.652145154862546;
78 const double pef = 0.347854845137454 * 0.652145154862546;
79
80 double pq16[16] = {pe2, pef, pef, pe2, pef, pf2, pf2, pef,
81 pef, pf2, pf2, pef, pe2, pef, pef, pe2};
82
83 IntPt GQQ16[16] = {
84 {{xq16[0], yq16[0], 0}, pq16[0]}, {{xq16[1], yq16[1], 0}, pq16[1]},
85 {{xq16[2], yq16[2], 0}, pq16[2]}, {{xq16[3], yq16[3], 0}, pq16[3]},
86 {{xq16[4], yq16[4], 0}, pq16[4]}, {{xq16[5], yq16[5], 0}, pq16[5]},
87 {{xq16[6], yq16[6], 0}, pq16[6]}, {{xq16[7], yq16[7], 0}, pq16[7]},
88 {{xq16[8], yq16[8], 0}, pq16[8]}, {{xq16[9], yq16[9], 0}, pq16[9]},
89 {{xq16[10], yq16[10], 0}, pq16[10]}, {{xq16[11], yq16[11], 0}, pq16[11]},
90 {{xq16[12], yq16[12], 0}, pq16[12]}, {{xq16[13], yq16[13], 0}, pq16[13]},
91 {{xq16[14], yq16[14], 0}, pq16[14]}, {{xq16[15], yq16[15], 0}, pq16[15]}};
92
93 static IntPt *GQQ[3] = {GQQ1, GQQ3, GQQ7};
94 static int GQQnPt[3] = {1, 3, 7};
95 static std::vector<IntPt *> GQQGL(40, nullptr);
96
getGQQPts(int order,bool forceTensorRule)97 IntPt *getGQQPts(int order, bool forceTensorRule)
98 {
99 if(!forceTensorRule && order <= 2) return GQQ[order];
100
101 if(static_cast<int>(GQQGL.size()) < order + 1)
102 GQQGL.resize(order + 1, nullptr);
103 if(!GQQGL[order]) {
104 int n = (order + 1) / (float)2 + 0.5;
105 double *pt, *wt;
106 gmshGaussLegendre1D(n, &pt, &wt);
107 IntPt *intpt = new IntPt[n * n];
108 int k = 0;
109 for(int i = 0; i < n; i++) {
110 for(int j = 0; j < n; j++) {
111 intpt[k].pt[0] = pt[i];
112 intpt[k].pt[1] = pt[j];
113 intpt[k].pt[2] = 0.0;
114 intpt[k++].weight = wt[i] * wt[j];
115 }
116 }
117 GQQGL[order] = intpt;
118 }
119 return GQQGL[order];
120 }
121
getNGQQPts(int order,bool forceTensorRule)122 int getNGQQPts(int order, bool forceTensorRule)
123 {
124 if(!forceTensorRule && order <= 2) return GQQnPt[order];
125
126 int n = (order + 1) / (float)2 + 0.5;
127 return n * n;
128 }
129