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 "GaussIntegration.h"
8 #include "GaussLegendre1D.h"
9 #include "GaussJacobi1D.h"
10
11 static std::vector<IntPt *> GQPyr(40, nullptr);
12
getGQPyrPts(int order)13 IntPt *getGQPyrPts(int order)
14 {
15 if(static_cast<int>(GQPyr.size()) < order + 1)
16 GQPyr.resize(order + 1, nullptr);
17 if(!GQPyr[order]) {
18 int nbPtUV = order / 2 + 1;
19 int nbPtW = order / 2 + 1;
20 int nbPtUV2 = nbPtUV * nbPtUV;
21
22 double *linPt, *linWt;
23 gmshGaussLegendre1D(nbPtUV, &linPt, &linWt);
24
25 double *GJ20Pt, *GJ20Wt;
26 getGaussJacobiQuadrature(2, 0, nbPtW, &GJ20Pt, &GJ20Wt);
27
28 IntPt *intpt = new IntPt[getNGQPyrPts(order)];
29
30 int l = 0;
31 for(int i = 0; i < getNGQPyrPts(order); i++) {
32 // compose an integration rule for (1-w)^2 f(u,v,w) on the standard
33 // hexahedron
34
35 int iW = i / (nbPtUV2);
36 int iU = (i - iW * nbPtUV2) / nbPtUV;
37 int iV = (i - iW * nbPtUV2 - iU * nbPtUV);
38
39 double wt = linWt[iU] * linWt[iV] * GJ20Wt[iW];
40
41 double up = linPt[iU];
42 double vp = linPt[iV];
43 double wp = GJ20Pt[iW];
44
45 // now incorporate the Duffy transformation from pyramid to hexahedron
46
47 intpt[l].pt[0] = 0.5 * (1 - wp) * up;
48 intpt[l].pt[1] = 0.5 * (1 - wp) * vp;
49 intpt[l].pt[2] = 0.5 * (1 + wp);
50
51 wt *= 0.125;
52 intpt[l++].weight = wt * 4. / 3.;
53 }
54 GQPyr[order] = intpt;
55 }
56 return GQPyr[order];
57 }
58
getNGQPyrPts(int order)59 int getNGQPyrPts(int order)
60 {
61 int nbPtUV = order / 2 + 1;
62 int nbPtW = order / 2 + 1;
63 return nbPtUV * nbPtUV * nbPtW;
64 }
65