Skip to content
Snippets Groups Projects
Commit 6817b6a5 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

Numeric : fix rounding in Quad and Hex quadrature order

parent 82dbf85b
No related branches found
No related tags found
No related merge requests found
...@@ -69,7 +69,7 @@ IntPt *getGQHPts(int order) ...@@ -69,7 +69,7 @@ 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 + 1) / (float)2 + 0.5;
int index = n-2 + 4; int index = n-2 + 4;
if(index >= (int)(sizeof(GQH) / sizeof(IntPt*))){ if(index >= (int)(sizeof(GQH) / sizeof(IntPt*))){
Msg::Error("Increase size of GQH in gauss quadrature hex"); Msg::Error("Increase size of GQH in gauss quadrature hex");
...@@ -100,5 +100,6 @@ int getNGQHPts(int order) ...@@ -100,5 +100,6 @@ int getNGQHPts(int order)
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); int n = (order + 1) / (float)2 + 0.5;
return n * n * n;
} }
...@@ -110,7 +110,7 @@ IntPt *getGQQPts(int order) ...@@ -110,7 +110,7 @@ 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 + 1) / (float)2 + 0.5;
int index = n-2 + 7; int index = n-2 + 7;
if(index >= (int)(sizeof(GQQ) / sizeof(IntPt*))){ if(index >= (int)(sizeof(GQQ) / sizeof(IntPt*))){
Msg::Error("Increase size of GQQ in gauss quadrature quad"); Msg::Error("Increase size of GQQ in gauss quadrature quad");
...@@ -139,5 +139,6 @@ int getNGQQPts(int order) ...@@ -139,5 +139,6 @@ int getNGQQPts(int order)
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); int n = (order + 1) / (float)2 + 0.5;
return n * n;
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment