Skip to content
Snippets Groups Projects
Commit 66ce8f1d authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

prevent out-of-bound access

parent 74d40d3b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
+ 2
− 2
View file @ 66ce8f1d
...@@ -91,8 +91,8 @@ opt(WRAP_JAVA "Enable generation of Java wrappers (experimental)" OFF) ...@@ -91,8 +91,8 @@ opt(WRAP_JAVA "Enable generation of Java wrappers (experimental)" OFF)
opt(WRAP_PYTHON "Enable generation of Python wrappers" OFF) opt(WRAP_PYTHON "Enable generation of Python wrappers" OFF)
set(GMSH_MAJOR_VERSION 2) set(GMSH_MAJOR_VERSION 2)
set(GMSH_MINOR_VERSION 8) set(GMSH_MINOR_VERSION 9)
set(GMSH_PATCH_VERSION 5) set(GMSH_PATCH_VERSION 0)
set(GMSH_EXTRA_VERSION "" CACHE STRING "Gmsh extra version string") set(GMSH_EXTRA_VERSION "" CACHE STRING "Gmsh extra version string")
set(GMSH_VERSION "${GMSH_MAJOR_VERSION}.${GMSH_MINOR_VERSION}") set(GMSH_VERSION "${GMSH_MAJOR_VERSION}.${GMSH_MINOR_VERSION}")
......
Numeric/GaussQuadratureHex.cpp
+ 35
− 40
View file @ 66ce8f1d
...@@ -20,12 +20,12 @@ const double yh6[6] = { b1, mb1, b1, mb1, 0., 0.}; ...@@ -20,12 +20,12 @@ const double yh6[6] = { b1, mb1, b1, mb1, 0., 0.};
const double zh6[6] = {mc1, mc1, c1, c1, mc1, c1}; const double zh6[6] = {mc1, mc1, c1, c1, mc1, c1};
const double ph6[6] = { w1, w1, w1, w1, w1, w1}; const double ph6[6] = { w1, w1, w1, w1, w1, w1};
IntPt GQH1[1] = IntPt GQH1[1] =
{ {
{ {0.0,0.0,0.0},8.0 } { {0.0,0.0,0.0},8.0 }
}; };
IntPt GQH6[6] = IntPt GQH6[6] =
{ {
{ {xh6[0],yh6[0],zh6[0]},ph6[0] }, { {xh6[0],yh6[0],zh6[0]},ph6[0] },
{ {xh6[1],yh6[1],zh6[1]},ph6[1] }, { {xh6[1],yh6[1],zh6[1]},ph6[1] },
...@@ -35,16 +35,18 @@ IntPt GQH6[6] = ...@@ -35,16 +35,18 @@ IntPt GQH6[6] =
{ {xh6[5],yh6[5],zh6[5]},ph6[5] } { {xh6[5],yh6[5],zh6[5]},ph6[5] }
}; };
const double xh8[8] = {0.577350269189626,-0.577350269189626,0.577350269189626,-0.577350269189626, const double xh8[8] =
0.577350269189626,-0.577350269189626,0.577350269189626,-0.577350269189626}; {0.577350269189626,-0.577350269189626,0.577350269189626,-0.577350269189626,
const double yh8[8] = {0.577350269189626,0.577350269189626,-0.577350269189626,-0.577350269189626, 0.577350269189626,-0.577350269189626,0.577350269189626,-0.577350269189626};
0.577350269189626,0.577350269189626,-0.577350269189626,-0.577350269189626}; const double yh8[8] =
{0.577350269189626,0.577350269189626,-0.577350269189626,-0.577350269189626,
const double zh8[8] = {-0.577350269189626,-0.577350269189626,-0.577350269189626,-0.577350269189626, 0.577350269189626,0.577350269189626,-0.577350269189626,-0.577350269189626};
0.577350269189626,0.577350269189626,0.577350269189626,0.577350269189626}; const double zh8[8] =
{-0.577350269189626,-0.577350269189626,-0.577350269189626,-0.577350269189626,
0.577350269189626,0.577350269189626,0.577350269189626,0.577350269189626};
const double ph8[8] = {1.,1.,1.,1.,1.,1.,1.,1.}; const double ph8[8] = {1.,1.,1.,1.,1.,1.,1.,1.};
IntPt GQH8[8] = IntPt GQH8[8] =
{ {
{ {xh8[0],yh8[0],zh8[0]},ph8[0] }, { {xh8[0],yh8[0],zh8[0]},ph8[0] },
{ {xh8[1],yh8[1],zh8[1]},ph8[1] }, { {xh8[1],yh8[1],zh8[1]},ph8[1] },
...@@ -56,7 +58,6 @@ IntPt GQH8[8] = ...@@ -56,7 +58,6 @@ IntPt GQH8[8] =
{ {xh8[7],yh8[7],zh8[7]},ph8[7] } { {xh8[7],yh8[7],zh8[7]},ph8[7] }
}; };
IntPt *getGQHPts(int order); IntPt *getGQHPts(int order);
int getNGQHPts(int order); int getNGQHPts(int order);
...@@ -64,46 +65,40 @@ IntPt * GQH[17] = {GQH1,GQH1,GQH6,GQH8,0,0,0,0,0,0,0,0,0,0,0,0,0}; ...@@ -64,46 +65,40 @@ IntPt * GQH[17] = {GQH1,GQH1,GQH6,GQH8,0,0,0,0,0,0,0,0,0,0,0,0,0};
int GQHnPt[4] = {1,1,6,8}; int GQHnPt[4] = {1,1,6,8};
IntPt *getGQHPts(int order) IntPt *getGQHPts(int order)
{ {
if(order < 2) return GQH[order];
if(order<2)return GQH[order]; if(order == 2) return GQH[3];
if(order == 2)return GQH[3]; if(order == 3) return GQH[3];
if(order == 3)return GQH[3];
int n = (order+3)/2; int n = (order+3)/2;
int index = n-2 + 4; int index = n-2 + 4;
if(!GQH[index]) if(index >= (int)(sizeof(GQH) / sizeof(IntPt*))){
{ Msg::Error("Increase size of GQH in gauss quadrature hex");
double *pt,*wt; index = 0;
// printf("o = %d n = %d i = %d\n",order,n*n*n,index); }
gmshGaussLegendre1D(n,&pt,&wt); if(!GQH[index]){
GQH[index] = new IntPt[n*n*n]; double *pt,*wt;
int l = 0; gmshGaussLegendre1D(n,&pt,&wt);
for(int i=0; i < n; i++) { GQH[index] = new IntPt[n*n*n];
for(int j=0; j < n; j++) { int l = 0;
for(int k=0; k < n; k++) { for(int i=0; i < n; i++) {
GQH[index][l].pt[0] = pt[i]; for(int j=0; j < n; j++) {
GQH[index][l].pt[1] = pt[j]; for(int k=0; k < n; k++) {
GQH[index][l].pt[2] = pt[k]; GQH[index][l].pt[0] = pt[i];
GQH[index][l++].weight = wt[i]*wt[j]*wt[k]; GQH[index][l].pt[1] = pt[j];
// printf ("%f %f %f %f\n",pt[i],pt[j],pt[k],wt[i]*wt[j]*wt[k]); GQH[index][l].pt[2] = pt[k];
} GQH[index][l++].weight = wt[i]*wt[j]*wt[k];
} }
} }
} }
}
return GQH[index]; return GQH[index];
} }
int getNGQHPts(int order) int getNGQHPts(int order)
{ {
if(order == 3)return 8; if(order == 3)return 8;
if(order == 2)return 8; if(order == 2)return 8;
if(order < 2) if(order < 2)
return GQHnPt[order]; return GQHnPt[order];
return ((order+3)/2)*((order+3)/2)*((order+3)/2); return ((order+3)/2)*((order+3)/2)*((order+3)/2);
} }
Numeric/GaussQuadratureLin.cpp
+ 7
− 4
View file @ 66ce8f1d
...@@ -9,12 +9,15 @@ ...@@ -9,12 +9,15 @@
IntPt * GQL[20] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; IntPt * GQL[20] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
IntPt *getGQLPts(int order) IntPt *getGQLPts(int order)
{ {
// Number of Gauss Point: // Number of Gauss Point:
// (order + 1) / 2 *ROUNDED UP* // (order + 1) / 2 *ROUNDED UP*
int n = (order + 1) / (double)2 + 0.5;
int n = (order + 1) / (double)2 + 0.5;
int index = n; int index = n;
if(index >= (int)(sizeof(GQL) / sizeof(IntPt*))){
Msg::Error("Increase size of GQL in gauss quadrature line");
index = 0;
}
if(!GQL[index]) { if(!GQL[index]) {
double *pt,*wt; double *pt,*wt;
gmshGaussLegendre1D(n,&pt,&wt); gmshGaussLegendre1D(n,&pt,&wt);
...@@ -30,6 +33,6 @@ IntPt *getGQLPts(int order) ...@@ -30,6 +33,6 @@ IntPt *getGQLPts(int order)
} }
int getNGQLPts(int order) int getNGQLPts(int order)
{ {
return (order + 1) / (double)2 + 0.5; return (order + 1) / (double)2 + 0.5;
} }
Numeric/GaussQuadraturePri.cpp
+ 4
− 7
View file @ 66ce8f1d
...@@ -19,8 +19,10 @@ IntPt *getGQPriPts(int order) ...@@ -19,8 +19,10 @@ IntPt *getGQPriPts(int order)
int nTri = getNGQTPts(order); int nTri = getNGQTPts(order);
int n = nLin*nTri; int n = nLin*nTri;
int index = order; int index = order;
if (order >= (int)(sizeof(GQP) / sizeof(IntPt*))) if (index >= (int)(sizeof(GQP) / sizeof(IntPt*))){
Msg::Fatal("Increase size of GQP in gauss quadrature prism"); Msg::Error("Increase size of GQP in gauss quadrature prism");
index = 0;
}
if(!GQP[index]){ if(!GQP[index]){
double *linPt,*linWt; double *linPt,*linWt;
IntPt *triPts = getGQTPts(order); IntPt *triPts = getGQTPts(order);
...@@ -51,8 +53,3 @@ int getNGQPriPts(int order) ...@@ -51,8 +53,3 @@ int getNGQPriPts(int order)
// return GQPnPt[order]; // return GQPnPt[order];
// return ((order+3)/2)*((order+3)/2)*((order+3)/2); // return ((order+3)/2)*((order+3)/2)*((order+3)/2);
} }
Numeric/GaussQuadraturePyr.cpp
+ 8
− 6
View file @ 66ce8f1d
...@@ -12,13 +12,17 @@ IntPt *getGQPyrPts(int order); ...@@ -12,13 +12,17 @@ IntPt *getGQPyrPts(int order);
int getNGQPyrPts(int order); int getNGQPyrPts(int order);
IntPt * GQPyr[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, IntPt * GQPyr[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
IntPt *getGQPyrPts(int order) IntPt *getGQPyrPts(int order)
{ {
int index = order; int index = order;
if(index >= (int)(sizeof(GQPyr) / sizeof(IntPt*))){
Msg::Error("Increase size of GQPyr in gauss quadrature pyr");
index = 0;
}
if(!GQPyr[index]) { if(!GQPyr[index]) {
int nbPtUV = order/2 + 1; int nbPtUV = order/2 + 1;
...@@ -32,8 +36,6 @@ IntPt *getGQPyrPts(int order) ...@@ -32,8 +36,6 @@ IntPt *getGQPyrPts(int order)
getGaussJacobiQuadrature(2, 0, nbPtW, &GJ20Pt, &GJ20Wt); getGaussJacobiQuadrature(2, 0, nbPtW, &GJ20Pt, &GJ20Wt);
GQPyr[index] = new IntPt[getNGQPyrPts(order)]; GQPyr[index] = new IntPt[getNGQPyrPts(order)];
if (order >= (int)(sizeof(GQPyr) / sizeof(IntPt*)))
Msg::Fatal("Increase size of GQPyr in gauss quadrature prism");
int l = 0; int l = 0;
for (int i = 0; i < getNGQPyrPts(order); i++) { for (int i = 0; i < getNGQPyrPts(order); i++) {
...@@ -51,7 +53,7 @@ IntPt *getGQPyrPts(int order) ...@@ -51,7 +53,7 @@ IntPt *getGQPyrPts(int order)
double up = linPt[iU]; double up = linPt[iU];
double vp = linPt[iV]; double vp = linPt[iV];
double wp = GJ20Pt[iW]; double wp = GJ20Pt[iW];
// now incorporate the Duffy transformation from pyramid to hexahedron // now incorporate the Duffy transformation from pyramid to hexahedron
GQPyr[index][l].pt[0] = 0.5*(1-wp)*up; GQPyr[index][l].pt[0] = 0.5*(1-wp)*up;
...@@ -71,7 +73,7 @@ IntPt *getGQPyrPts(int order) ...@@ -71,7 +73,7 @@ IntPt *getGQPyrPts(int order)
int getNGQPyrPts(int order) int getNGQPyrPts(int order)
{ {
int nbPtUV = order/2 + 1; int nbPtUV = order/2 + 1;
int nbPtW = order/2 + 1; int nbPtW = order/2 + 1;
return nbPtUV*nbPtUV*nbPtW; return nbPtUV*nbPtUV*nbPtW;
} }
Numeric/GaussQuadratureQuad.cpp
+ 27
− 34
View file @ 66ce8f1d
...@@ -50,7 +50,7 @@ const double a9 = 0.774596669241483; ...@@ -50,7 +50,7 @@ const double a9 = 0.774596669241483;
const double z = 0.0; const double z = 0.0;
const double xq9[9] = {-a9, z, a9, -a9, z, a9, -a9, z, a9}; const double xq9[9] = {-a9, z, a9, -a9, z, a9, -a9, z, a9};
const double yq9[9] = {-a9, -a9, -a9, z, z, z, a9, a9, a9}; const double yq9[9] = {-a9, -a9, -a9, z, z, z, a9, a9, a9};
const double pb2 = 0.888888888888889*0.888888888888889; const double pb2 = 0.888888888888889*0.888888888888889;
const double pc2 = 0.555555555555556*0.555555555555556; const double pc2 = 0.555555555555556*0.555555555555556;
const double pbc = 0.555555555555556*0.888888888888889; const double pbc = 0.555555555555556*0.888888888888889;
const double pq9[9] = {pc2, pbc, pc2, pbc, pb2, pbc, pc2, pbc , pc2}; const double pq9[9] = {pc2, pbc, pc2, pbc, pb2, pbc, pc2, pbc , pc2};
...@@ -69,7 +69,7 @@ IntPt GQQ9[9] = { ...@@ -69,7 +69,7 @@ IntPt GQQ9[9] = {
const double a16 = 0.861136311594053; const double a16 = 0.861136311594053;
const double b16 = 0.339981043584856; const double b16 = 0.339981043584856;
const double xq16[16] = {-a16, -b16, b16, a16, -a16, -b16, b16, a16, const double xq16[16] = {-a16, -b16, b16, a16, -a16, -b16, b16, a16,
-a16, -b16, b16, a16, -a16, -b16, b16, a16}; -a16, -b16, b16, a16, -a16, -b16, b16, a16};
const double yq16[16] = {-a16, -a16, -a16, -a16, -b16, -b16, -b16, -b16, const double yq16[16] = {-a16, -a16, -a16, -a16, -b16, -b16, -b16, -b16,
b16, b16, b16, b16, a16, a16, a16, a16}; b16, b16, b16, b16, a16, a16, a16, a16};
...@@ -77,7 +77,7 @@ const double pe2 = 0.347854845137454 * 0.347854845137454; ...@@ -77,7 +77,7 @@ const double pe2 = 0.347854845137454 * 0.347854845137454;
const double pf2 = 0.652145154862546 * 0.652145154862546; const double pf2 = 0.652145154862546 * 0.652145154862546;
const double pef = 0.347854845137454 * 0.652145154862546; const double pef = 0.347854845137454 * 0.652145154862546;
double pq16[16] = { pe2, pef, pef, pe2, pef, pf2, pf2, pef, double pq16[16] = { pe2, pef, pef, pe2, pef, pf2, pf2, pef,
pef, pf2, pf2, pef, pe2, pef, pef, pe2}; pef, pf2, pf2, pef, pe2, pef, pef, pe2};
IntPt GQQ16[16] = { IntPt GQQ16[16] = {
...@@ -106,45 +106,38 @@ IntPt * GQQ[27] = {GQQ1,GQQ1,GQQ3,GQQ4,GQQ7,GQQ9,GQQ16,0,0,0,0,0,0,0,0,0,0,0,0,0 ...@@ -106,45 +106,38 @@ IntPt * GQQ[27] = {GQQ1,GQQ1,GQQ3,GQQ4,GQQ7,GQQ9,GQQ16,0,0,0,0,0,0,0,0,0,0,0,0,0
int GQQnPt[7] = {1,1,3,4,7,9,16}; int GQQnPt[7] = {1,1,3,4,7,9,16};
IntPt *getGQQPts(int order) IntPt *getGQQPts(int order)
{ {
if(order < 2) return GQQ[order];
if(order<2)return GQQ[order]; if(order == 2) return GQQ[3];
if(order==2)return GQQ[3]; if(order == 3) return GQQ[3];
if(order==3)return GQQ[3];
int n = (order+3)/2; int n = (order+3)/2;
int index = n-2 + 7; int index = n-2 + 7;
if(!GQQ[index]) if(index >= (int)(sizeof(GQQ) / sizeof(IntPt*))){
{ Msg::Error("Increase size of GQQ in gauss quadrature quad");
double *pt,*wt; index = 0;
gmshGaussLegendre1D(n,&pt,&wt); }
GQQ[index] = new IntPt[n*n]; if(!GQQ[index]){
int k = 0; double *pt,*wt;
for(int i=0; i < n; i++) { gmshGaussLegendre1D(n,&pt,&wt);
for(int j=0; j < n; j++) { GQQ[index] = new IntPt[n*n];
GQQ[index][k].pt[0] = pt[i]; int k = 0;
GQQ[index][k].pt[1] = pt[j]; for(int i=0; i < n; i++) {
GQQ[index][k].pt[2] = 0.0; for(int j=0; j < n; j++) {
GQQ[index][k++].weight = wt[i]*wt[j]; GQQ[index][k].pt[0] = pt[i];
// printf ("%f %f %f\n",pt[i],pt[j],wt[i]*wt[j]); GQQ[index][k].pt[1] = pt[j];
} GQQ[index][k].pt[2] = 0.0;
GQQ[index][k++].weight = wt[i]*wt[j];
} }
} }
}
return GQQ[index]; return GQQ[index];
} }
int getNGQQPts(int order) int getNGQQPts(int order)
{ {
if(order == 3)return 4; if(order == 3) return 4;
if(order == 2)return 4; if(order == 2) return 4;
if(order < 2) if(order < 2)
return GQQnPt[order]; return GQQnPt[order];
return ((order+3)/2)*((order+3)/2); return ((order+3)/2)*((order+3)/2);
} }
Numeric/GaussQuadratureTri.cpp
+ 11
− 12
View file @ 66ce8f1d
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
#include "GaussIntegration.h" #include "GaussIntegration.h"
#include "GaussLegendre1D.h" #include "GaussLegendre1D.h"
IntPt GQT1[1] = { IntPt GQT1[1] = {
{{0.333333333333333, 0.333333333333333, 0.}, 0.500000000000000} {{0.333333333333333, 0.333333333333333, 0.}, 0.500000000000000}
}; };
IntPt GQT2[3] = { IntPt GQT2[3] = {
...@@ -83,7 +83,7 @@ IntPt GQT8[16] = { ...@@ -83,7 +83,7 @@ IntPt GQT8[16] = {
{{0.170569307751760, 0.170569307751760, 0.}, 0.051608685267359}, {{0.170569307751760, 0.170569307751760, 0.}, 0.051608685267359},
{{0.898905543365938, 0.050547228317031, 0.}, 0.016229248811599}, {{0.898905543365938, 0.050547228317031, 0.}, 0.016229248811599},
{{0.050547228317031, 0.898905543365938, 0.}, 0.016229248811599}, {{0.050547228317031, 0.898905543365938, 0.}, 0.016229248811599},
{{0.050547228317031, 0.050547228317031, 0.}, 0.016229248811599}, {{0.050547228317031, 0.050547228317031, 0.}, 0.016229248811599},
{{0.008394777409958, 0.728492392955404, 0.}, 0.013615157087217}, {{0.008394777409958, 0.728492392955404, 0.}, 0.013615157087217},
{{0.728492392955404, 0.008394777409958, 0.}, 0.013615157087217}, {{0.728492392955404, 0.008394777409958, 0.}, 0.013615157087217},
{{0.263112829634638, 0.008394777409958, 0.}, 0.013615157087217}, {{0.263112829634638, 0.008394777409958, 0.}, 0.013615157087217},
...@@ -98,9 +98,6 @@ IntPt * GQT[9] = {GQT1,GQT1,GQT2,GQT3,GQT4,GQT5,GQT6,GQT7,GQT8}; ...@@ -98,9 +98,6 @@ IntPt * GQT[9] = {GQT1,GQT1,GQT2,GQT3,GQT4,GQT5,GQT6,GQT7,GQT8};
IntPt * GQTdegen[17] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; IntPt * GQTdegen[17] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
int GQTnPt[9] = {1,1,3,4,6,7,12,13,16}; int GQTnPt[9] = {1,1,3,4,6,7,12,13,16};
// ----------------------------------------------------------------------------- // -----------------------------------------------------------------------------
/*! Quadrature rule for an interpolation of order 1 on the triangle */ /*! Quadrature rule for an interpolation of order 1 on the triangle */
/* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */ /* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */
...@@ -912,7 +909,7 @@ IntPt * GQTSolin[21] = { ...@@ -912,7 +909,7 @@ IntPt * GQTSolin[21] = {
triP19Solin, triP19Solin,
triP20Solin triP20Solin
}; };
int GQTnPtSolin[21] = { int GQTnPtSolin[21] = {
1, 1,
1, 1,
...@@ -940,22 +937,24 @@ IntPt *getGQTPts(int order); ...@@ -940,22 +937,24 @@ IntPt *getGQTPts(int order);
int getNGQTPts(int order); int getNGQTPts(int order);
IntPt *getGQTPts(int order) IntPt *getGQTPts(int order)
{ {
if (order < 21) return GQTSolin[order]; if (order < 21) return GQTSolin[order];
int n = (order+3)/2; int n = (order+3)/2;
int index = n-4; int index = n-4;
if(index >= (int)(sizeof(GQTdegen) / sizeof(IntPt*))){
Msg::Error("Increase size of GQTdegen in gauss quadrature tri");
index = 0;
}
if(!GQTdegen[index]){ if(!GQTdegen[index]){
int npts = n*n; int npts = n*n;
GQTdegen[index] = new IntPt[npts]; GQTdegen[index] = new IntPt[npts];
GaussLegendreTri(n,n,GQTdegen[index]); GaussLegendreTri(n,n,GQTdegen[index]);
} }
return GQTdegen[index]; return GQTdegen[index];
} }
int getNGQTPts(int order) int getNGQTPts(int order)
{ {
if (order < 21) return GQTnPtSolin[order]; if (order < 21) return GQTnPtSolin[order];
return ((order+3)/2)*((order+3)/2); return ((order+3)/2)*((order+3)/2);
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment