diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index f16e3bd60a47b1c79bf73f51facc9c179e78b20c..eb4267ecf0c6eb7cb4127d40dfb8b83008714893 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -267,14 +267,10 @@ bool MLineChild::isInside(double u, double v, double w)
 
 void MLineChild::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 {
-
-
   *npts = 0;
   double jac[3][3];
   if(_intpt) delete [] _intpt;
-
-  int nbP = pOrder / 2 + 1;  // MLine::getIntegrationPoints()
-  _intpt = new IntPt[nbP];
+  _intpt = new IntPt[getNGQLPts(pOrder)];
 
   int nptsi;
   IntPt *ptsi;
@@ -291,30 +287,30 @@ void MLineChild::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 //  -----------mich mach---------------------//
 
 
-		MVertex *vo = _orig->getVertex(0);
-    MVertex *vf = _orig->getVertex(1);
+  MVertex *vo = _orig->getVertex(0);
+  MVertex *vf = _orig->getVertex(1);
+
+  SPoint3 P(vo->x(), vo->y(), vo->z());
+  SPoint3 Q(vf->x(), vf->y(), vf->z());
+
+  SPoint3 R;
+
+  R = P + Q;
+  R/=2;
 
-		SPoint3 P(vo->x(), vo->y(), vo->z());
-    SPoint3 Q(vf->x(), vf->y(), vf->z());
+  vo = getVertex(0);
+  vf = getVertex(1);
 
-		SPoint3 R;
-  
-    R = P + Q;
-    R/=2;
+  SPoint3 A(vo->x(), vo->y(), vo->z());
+  SPoint3 B(vf->x(), vf->y(), vf->z());
 
-		vo = getVertex(0);
-    vf = getVertex(1);
+  double lengthDemi = R.distance(Q);
 
-    SPoint3 A(vo->x(), vo->y(), vo->z());
-    SPoint3 B(vf->x(), vf->y(), vf->z());
+  if (P.distance(A) < Q.distance(A)) {v_uvw[0][0] = - R.distance(A)/lengthDemi;v_uvw[0][1]=0;v_uvw[0][2]=0;}
+  else {v_uvw[0][0] = R.distance(A)/lengthDemi;v_uvw[0][1]=0;v_uvw[0][2]=0;}
 
-    double lengthDemi = R.distance(Q);
-	    
-    if (P.distance(A) < Q.distance(A)) {v_uvw[0][0] = - R.distance(A)/lengthDemi;v_uvw[0][1]=0;v_uvw[0][2]=0;}
-    else {v_uvw[0][0] = R.distance(A)/lengthDemi;v_uvw[0][1]=0;v_uvw[0][2]=0;}
-		
-    if (P.distance(B) < Q.distance(B)) {v_uvw[1][0] = - R.distance(B)/lengthDemi;v_uvw[1][1]=0;v_uvw[1][2]=0;}
-    else {v_uvw[1][0] = R.distance(B)/lengthDemi;v_uvw[1][1]=0;v_uvw[1][2]=0;}
+  if (P.distance(B) < Q.distance(B)) {v_uvw[1][0] = - R.distance(B)/lengthDemi;v_uvw[1][1]=0;v_uvw[1][2]=0;}
+  else {v_uvw[1][0] = R.distance(B)/lengthDemi;v_uvw[1][1]=0;v_uvw[1][2]=0;}
 
 //  ---------------------------------------//
 
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index b6cdc2fc7f795d6bdc6f854f789e6c500cab18f7..eb79c54a6c2dc9babbd02c2b998af1e3b2a85415 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -31,18 +31,8 @@ const polynomialBasis* MLine::getFunctionSpace(int o) const
 
 void MLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 {
-  static IntPt GQL[100]; 
-  double *t, *w;
-  int nbP = pOrder / 2 + 1;
-  gmshGaussLegendre1D(nbP, &t, &w);
-  for (int i = 0; i < nbP; i++){
-    GQL[i].pt[0] = t[i];
-    GQL[i].pt[1] = 0;
-    GQL[i].pt[2] = 0;
-    GQL[i].weight = w[i];
-  }
-  *npts = nbP;
-  *pts = GQL;
+  *npts = getNGQLPts(pOrder);
+  *pts = getGQLPts(pOrder);
 }
 
 double MLine::getInnerRadius()
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 4dd2e58936940d5d474ac7aab9475286630a6793..48b7bbf9b166523d40d6e2c914c1b3d2137d7ee1 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -3,7 +3,7 @@
 
 #include "Integration3D.h"
 #include "recurCut.h"
-#include "../../Numeric/Gauss.h"
+#include "Gauss.h"
 
 #define ZERO_LS_TOL  1.e-7
 #define EQUALITY_TOL 1.e-15