LineGaussRuleInfo[ {rule_, numer_ }, point_ ] := Module [
{ g2 = {-1,1}/Sqrt[3],
w3={5/9,8/9,5/9},
g3={-Sqrt[3/5],0,Sqrt[3/5]},
w4={(1/2)-Sqrt[5/6]/6,
(1/2)+Sqrt[5/6]/6,
(1/2)+Sqrt[5/6]/6,
(1/2)-Sqrt[5/6]/6},
g4={-Sqrt[(3+2*Sqrt[6/5])/7],
-Sqrt[(3-2*Sqrt[6/5])/7],
Sqrt[(3-2*Sqrt[6/5])/7],
Sqrt[(3+2*Sqrt[6/5])/7]},
g5={-Sqrt[5+2*Sqrt[10/7]],
-Sqrt[5-2*Sqrt[10/7]],
0,
Sqrt[5-2*Sqrt[10/7]],
Sqrt[5+2*Sqrt[10/7]]}/3,
w5={322-13*Sqrt[70],
322+13*Sqrt[70],
512,
322+13*Sqrt[70],
322-13*Sqrt[70]}/900,
i=point,p=rule, info={Null,0}},
If [ p == 1,info = {0,2}];
If [ p == 2, info={g2[[i]],1}];
If [ p == 3, info={g3[[i]],w3[[i]]}];
If [ p == 4, info ={g4[[i]],w4[[i]]}];
If [ p == 5, info = {g5[[i]],w5[[i]]}];
If [ numer, Return [ N[info,20]], Return [ Simplify[info]]];
];

HexaGaussRuleInfo[{rule_,numer_},point_] := Module [
  {x,h,m,p1,p2,p3,p12,i,j,jj,k,mm,w1,w2,w3,
  info={{Null,Null,Null},0}},
  If [ Length[rule] == 3, {p1,p2,p3} = rule, p1=p2=p3=rule];
  If [ Length[point] == 3, {i,j,k}=point, mm= point;
      p12=p1*p2; k=Floor[(mm-1)/p12]+1; jj=mm-p12*(k-1);
      j=Floor[(jj-1)/p1]+1; i=jj-p1*(j-1)];
      {x,w1} = LineGaussRuleInfo[{p1,numer},i];
      {h,w2} = LineGaussRuleInfo[{p2,numer},j];
      {m,w3} = LineGaussRuleInfo[{p3,numer},k];
      info = {{x,h,m},w1*w2*w3};
      If [ numer, Return[N[info,20]], Return[Simplify[info]]];
];

Print[" "];
Print["HexaGaussRuleInfoTest.nb"];
Print["  Mathematica version"];
Print[" "];
Print["  Test the quadrature rules for a unit hexahedron,"];
Print["  returned as {{x,y,z},w} = HexaGaussRuleInfo[{rule,numer},i]"];
Print[" "];
Print["  INPUT:"];
Print[" "];
Print["  RULE is 1, 2, 3, 4, 5, the order in both dimensions,"];
Print["  or RULE = { order1, order2, order3 }, allowing the user to"];
Print["  specify different orders in each dimension."];
Print[" "];
Print["  NUMER is True to return numerical values for the output, or"];
Print["  False to return symbolic values."];
Print[" "];
Print["  I is the index, between 1 and RULE*RULE*RULE, of the point."]
Print[" The index can also be specified as {I,J,K}."]
Print[" "];
Print["  OUTPUT:"];
Print[" "];
Print["  {X,Y,Z} are the coordinates of the abscissas,"]
Print["  For the unit hexahedron, -1 <= X, Y, Z <= 1."]
Print[" "];
Print["  W is the weight of the quadrature point."]
Print[" "]
Print["  The following cells print out each rule, and the sum"]
Print["  of the weights."]

For [ rule = 1, rule ≤ 5, rule++,
numer = True;
wsum = 0;
Print[" "];
Print["Hexahedron rule number ", rule];
Print[" "];
Print["W"];
Print[" "];
For [ i=1, i ≤ rule * rule * rule, i++,
{ { x, y, z},w} = HexaGaussRuleInfo[{rule,numer},i];
Print [ w ];
wsum = wsum +w;
];
Print["sum(w) = ", wsum ];
Print[" "];
Print["X, Y, Z"];
Print[" "];
For [ i=1, i ≤ rule * rule * rule, i++,
{ { x, y, z },w} = HexaGaussRuleInfo[{rule,numer},i];
Print [ x,"  ", y, "  ", z ];
];
];