From 66ce8f1dfbd78a008fb7e4dc92655caa9f6c27c0 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Thu, 10 Jul 2014 05:40:09 +0000
Subject: [PATCH] prevent out-of-bound access

---
 CMakeLists.txt                  |  4 +-
 Numeric/GaussQuadratureHex.cpp  | 75 +++++++++++++++------------------
 Numeric/GaussQuadratureLin.cpp  | 11 +++--
 Numeric/GaussQuadraturePri.cpp  | 11 ++---
 Numeric/GaussQuadraturePyr.cpp  | 14 +++---
 Numeric/GaussQuadratureQuad.cpp | 61 ++++++++++++---------------
 Numeric/GaussQuadratureTri.cpp  | 23 +++++-----
 7 files changed, 94 insertions(+), 105 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0c00f7e01f..ddcc217f08 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -91,8 +91,8 @@ opt(WRAP_JAVA "Enable generation of Java wrappers (experimental)" OFF)
 opt(WRAP_PYTHON "Enable generation of Python wrappers" OFF)
 
 set(GMSH_MAJOR_VERSION 2)
-set(GMSH_MINOR_VERSION 8)
-set(GMSH_PATCH_VERSION 5)
+set(GMSH_MINOR_VERSION 9)
+set(GMSH_PATCH_VERSION 0)
 set(GMSH_EXTRA_VERSION "" CACHE STRING "Gmsh extra version string")
 
 set(GMSH_VERSION "${GMSH_MAJOR_VERSION}.${GMSH_MINOR_VERSION}")
diff --git a/Numeric/GaussQuadratureHex.cpp b/Numeric/GaussQuadratureHex.cpp
index c40f83f836..7f84c8e871 100644
--- a/Numeric/GaussQuadratureHex.cpp
+++ b/Numeric/GaussQuadratureHex.cpp
@@ -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 ph6[6] = { w1,  w1,  w1,  w1,  w1,  w1};
 
-IntPt GQH1[1] = 
+IntPt GQH1[1] =
 {
   { {0.0,0.0,0.0},8.0 }
 };
 
-IntPt GQH6[6] = 
+IntPt GQH6[6] =
 {
   { {xh6[0],yh6[0],zh6[0]},ph6[0] },
   { {xh6[1],yh6[1],zh6[1]},ph6[1] },
@@ -35,16 +35,18 @@ IntPt GQH6[6] =
   { {xh6[5],yh6[5],zh6[5]},ph6[5] }
 };
 
-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};
-
-const double zh8[8] = {-0.577350269189626,-0.577350269189626,-0.577350269189626,-0.577350269189626,
-                   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};
+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.};
 
-IntPt GQH8[8] = 
+IntPt GQH8[8] =
 {
   { {xh8[0],yh8[0],zh8[0]},ph8[0] },
   { {xh8[1],yh8[1],zh8[1]},ph8[1] },
@@ -56,7 +58,6 @@ IntPt GQH8[8] =
   { {xh8[7],yh8[7],zh8[7]},ph8[7] }
 };
 
-
 IntPt *getGQHPts(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};
 int GQHnPt[4] = {1,1,6,8};
 
 IntPt *getGQHPts(int order)
-{ 
-
-  if(order<2)return GQH[order];
-  if(order == 2)return GQH[3]; 
-  if(order == 3)return GQH[3]; 
-  
+{
+  if(order < 2) return GQH[order];
+  if(order == 2) return GQH[3];
+  if(order == 3) return GQH[3];
   int n = (order+3)/2;
   int index = n-2 + 4;
-  if(!GQH[index])
-    {
-      double *pt,*wt;
-      //      printf("o = %d n = %d i = %d\n",order,n*n*n,index);
-      gmshGaussLegendre1D(n,&pt,&wt);
-      GQH[index] = new IntPt[n*n*n];
-      int l = 0;
-      for(int i=0; i < n; i++) {
-        for(int j=0; j < n; j++) {
-          for(int k=0; k < n; k++) {
-            GQH[index][l].pt[0] = pt[i];
-            GQH[index][l].pt[1] = pt[j];
-            GQH[index][l].pt[2] = pt[k];
-            GQH[index][l++].weight = wt[i]*wt[j]*wt[k];
-            //      printf ("%f %f %f %f\n",pt[i],pt[j],pt[k],wt[i]*wt[j]*wt[k]);
-          }
+  if(index >= (int)(sizeof(GQH) / sizeof(IntPt*))){
+    Msg::Error("Increase size of GQH in gauss quadrature hex");
+    index = 0;
+  }
+  if(!GQH[index]){
+    double *pt,*wt;
+    gmshGaussLegendre1D(n,&pt,&wt);
+    GQH[index] = new IntPt[n*n*n];
+    int l = 0;
+    for(int i=0; i < n; i++) {
+      for(int j=0; j < n; j++) {
+        for(int k=0; k < n; k++) {
+          GQH[index][l].pt[0] = pt[i];
+          GQH[index][l].pt[1] = pt[j];
+          GQH[index][l].pt[2] = pt[k];
+          GQH[index][l++].weight = wt[i]*wt[j]*wt[k];
         }
       }
     }
+  }
   return GQH[index];
 }
 
 int getNGQHPts(int order)
-{ 
+{
   if(order == 3)return 8;
   if(order == 2)return 8;
   if(order < 2)
-    return GQHnPt[order]; 
+    return GQHnPt[order];
   return ((order+3)/2)*((order+3)/2)*((order+3)/2);
 }
-
-
-
-
-
diff --git a/Numeric/GaussQuadratureLin.cpp b/Numeric/GaussQuadratureLin.cpp
index 62743a20e0..0b922188fc 100644
--- a/Numeric/GaussQuadratureLin.cpp
+++ b/Numeric/GaussQuadratureLin.cpp
@@ -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 *getGQLPts(int order)
-{  
+{
   // Number of Gauss Point:
   // (order + 1) / 2 *ROUNDED UP*
-
-  int n     = (order + 1) / (double)2 + 0.5;
+  int n = (order + 1) / (double)2 + 0.5;
   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]) {
     double *pt,*wt;
     gmshGaussLegendre1D(n,&pt,&wt);
@@ -30,6 +33,6 @@ IntPt *getGQLPts(int order)
 }
 
 int getNGQLPts(int order)
-{ 
+{
   return (order + 1) / (double)2 + 0.5;
 }
diff --git a/Numeric/GaussQuadraturePri.cpp b/Numeric/GaussQuadraturePri.cpp
index 227e4cebc6..355c40fdce 100644
--- a/Numeric/GaussQuadraturePri.cpp
+++ b/Numeric/GaussQuadraturePri.cpp
@@ -19,8 +19,10 @@ IntPt *getGQPriPts(int order)
   int nTri = getNGQTPts(order);
   int n = nLin*nTri;
   int index = order;
-  if (order >= (int)(sizeof(GQP) / sizeof(IntPt*)))
-    Msg::Fatal("Increase size of GQP in gauss quadrature prism");
+  if (index >= (int)(sizeof(GQP) / sizeof(IntPt*))){
+    Msg::Error("Increase size of GQP in gauss quadrature prism");
+    index = 0;
+  }
   if(!GQP[index]){
     double *linPt,*linWt;
     IntPt *triPts = getGQTPts(order);
@@ -51,8 +53,3 @@ int getNGQPriPts(int order)
 //     return GQPnPt[order];
 //   return ((order+3)/2)*((order+3)/2)*((order+3)/2);
 }
-
-
-
-
-
diff --git a/Numeric/GaussQuadraturePyr.cpp b/Numeric/GaussQuadraturePyr.cpp
index 6c1c5b9a79..a83864dc4b 100644
--- a/Numeric/GaussQuadraturePyr.cpp
+++ b/Numeric/GaussQuadraturePyr.cpp
@@ -12,13 +12,17 @@ IntPt *getGQPyrPts(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,
-                 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)
 {
-
   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]) {
 
     int nbPtUV = order/2 + 1;
@@ -32,8 +36,6 @@ IntPt *getGQPyrPts(int order)
     getGaussJacobiQuadrature(2, 0, nbPtW, &GJ20Pt, &GJ20Wt);
 
     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;
     for (int i = 0; i < getNGQPyrPts(order); i++) {
@@ -51,7 +53,7 @@ IntPt *getGQPyrPts(int order)
       double up = linPt[iU];
       double vp = linPt[iV];
       double wp = GJ20Pt[iW];
-      
+
       // now incorporate the Duffy transformation from pyramid to hexahedron
 
       GQPyr[index][l].pt[0] = 0.5*(1-wp)*up;
@@ -71,7 +73,7 @@ IntPt *getGQPyrPts(int order)
 int getNGQPyrPts(int order)
 {
   int nbPtUV = order/2 + 1;
-  int nbPtW  = order/2 + 1; 
+  int nbPtW  = order/2 + 1;
 
   return nbPtUV*nbPtUV*nbPtW;
 }
diff --git a/Numeric/GaussQuadratureQuad.cpp b/Numeric/GaussQuadratureQuad.cpp
index 09293f42fa..e581b2bfdb 100644
--- a/Numeric/GaussQuadratureQuad.cpp
+++ b/Numeric/GaussQuadratureQuad.cpp
@@ -50,7 +50,7 @@ const double a9 = 0.774596669241483;
 const double z = 0.0;
 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 pb2 = 0.888888888888889*0.888888888888889; 
+const double pb2 = 0.888888888888889*0.888888888888889;
 const double pc2 = 0.555555555555556*0.555555555555556;
 const double pbc = 0.555555555555556*0.888888888888889;
 const double pq9[9] = {pc2, pbc, pc2, pbc, pb2, pbc, pc2, pbc , pc2};
@@ -69,7 +69,7 @@ IntPt GQQ9[9] = {
 
 const double a16 = 0.861136311594053;
 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};
 const double yq16[16] = {-a16, -a16, -a16, -a16, -b16, -b16, -b16, -b16,
                          b16, b16, b16, b16, a16, a16, a16, a16};
@@ -77,7 +77,7 @@ const double pe2 = 0.347854845137454 * 0.347854845137454;
 const double pf2 = 0.652145154862546 * 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};
 
 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
 int GQQnPt[7] = {1,1,3,4,7,9,16};
 
 IntPt *getGQQPts(int order)
-{ 
-  
-  if(order<2)return GQQ[order];
-  if(order==2)return GQQ[3];
-  if(order==3)return GQQ[3];
-
+{
+  if(order < 2) return GQQ[order];
+  if(order == 2) return GQQ[3];
+  if(order == 3) return GQQ[3];
   int n = (order+3)/2;
   int index = n-2 + 7;
-  if(!GQQ[index])
-    {
-      double *pt,*wt;
-      gmshGaussLegendre1D(n,&pt,&wt);
-      GQQ[index] = new IntPt[n*n];
-      int k = 0;
-      for(int i=0; i < n; i++) {
-        for(int j=0; j < n; j++) {
-          GQQ[index][k].pt[0] = pt[i];
-          GQQ[index][k].pt[1] = pt[j];
-          GQQ[index][k].pt[2] = 0.0;
-          GQQ[index][k++].weight = wt[i]*wt[j];
-          //      printf ("%f %f %f\n",pt[i],pt[j],wt[i]*wt[j]);
-        }
+  if(index >= (int)(sizeof(GQQ) / sizeof(IntPt*))){
+    Msg::Error("Increase size of GQQ in gauss quadrature quad");
+    index = 0;
+  }
+  if(!GQQ[index]){
+    double *pt,*wt;
+    gmshGaussLegendre1D(n,&pt,&wt);
+    GQQ[index] = new IntPt[n*n];
+    int k = 0;
+    for(int i=0; i < n; i++) {
+      for(int j=0; j < n; j++) {
+        GQQ[index][k].pt[0] = pt[i];
+        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];
 }
 
 int getNGQQPts(int order)
-{ 
-  if(order == 3)return 4;
-  if(order == 2)return 4;
-
+{
+  if(order == 3) return 4;
+  if(order == 2) return 4;
   if(order < 2)
-    return GQQnPt[order]; 
+    return GQQnPt[order];
   return ((order+3)/2)*((order+3)/2);
 }
-
-
-
-
-
-
diff --git a/Numeric/GaussQuadratureTri.cpp b/Numeric/GaussQuadratureTri.cpp
index 06831a73e5..9486636cb5 100644
--- a/Numeric/GaussQuadratureTri.cpp
+++ b/Numeric/GaussQuadratureTri.cpp
@@ -6,7 +6,7 @@
 #include "GaussIntegration.h"
 #include "GaussLegendre1D.h"
 
-IntPt GQT1[1] = { 
+IntPt GQT1[1] = {
   {{0.333333333333333, 0.333333333333333, 0.}, 0.500000000000000}
 };
 IntPt GQT2[3] = {
@@ -83,7 +83,7 @@ IntPt GQT8[16] = {
   {{0.170569307751760, 0.170569307751760, 0.}, 0.051608685267359},
   {{0.898905543365938, 0.050547228317031, 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.728492392955404, 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};
 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};
 
-
-
-
 // -----------------------------------------------------------------------------
 /*! Quadrature rule for an interpolation of order 1 on the triangle */
 /* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel */
@@ -912,7 +909,7 @@ IntPt * GQTSolin[21] = {
   triP19Solin,
   triP20Solin
 };
-  
+
 int GQTnPtSolin[21] = {
   1,
   1,
@@ -940,22 +937,24 @@ IntPt *getGQTPts(int order);
 int getNGQTPts(int order);
 
 IntPt *getGQTPts(int order)
-{ 
+{
   if (order < 21) return GQTSolin[order];
   int n = (order+3)/2;
   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]){
     int npts = n*n;
     GQTdegen[index] = new IntPt[npts];
     GaussLegendreTri(n,n,GQTdegen[index]);
   }
-  return GQTdegen[index]; 
-
+  return GQTdegen[index];
 }
 
 int getNGQTPts(int order)
-{ 
-  if (order < 21) return GQTnPtSolin[order];  
+{
+  if (order < 21) return GQTnPtSolin[order];
   return ((order+3)/2)*((order+3)/2);
 }
-
-- 
GitLab