13

Usually one would consult a paper or book to find integration points and weights for unit triangle and tetrahedra. I am looking for a method to automatically compute such points and weights. The following Mathematica code example computes integration weights and points for unit line (quad/hexahedron) elements:

unitGaussianQuadraturePoints[points_] := 
  Sort[x /. 
    Solve[Evaluate[LegendreP[points, x] == 0], {x}], ! 
     OrderedQ[N[{#1, #2}]] &];

unitGaussianQuadratureWeights[points_] := 
  Module[{gps, f, int, integr, vars, eqns}, 
   gps = unitGaussianQuadraturePoints[points];
   f[0, 0] := 1;
   f[0., 0] := 1.;
   f[x_, n_] := x^n;
   int = Integrate[f[x, #], x] & /@ Range[0, points - 1];
   integr = Subtract @@@ (int /. x :> {1, -1});
   vars = Table[Unique[c], {Length[gps]}];
   eqns = 
    Table[Plus @@ Thread[Times[vars, f[#, i - 1] & /@ gps]] == 
      integr[[i]], {i, points}];
   Return[(vars /. Solve[eqns, vars])];];


unitGaussianQuadratureWeights[2]

{{1, 1}}

unitGaussianQuadraturePoints[2]

{1/Sqrt[3], -(1/Sqrt[3])}

I am looking for a paper/book that describes algorithmically how this is done for triangles and/or tetrahedra. Can someone point me to some information about this. Thanks.

David Ketcheson
  • 16,522
  • 4
  • 54
  • 105
  • 1
    There's an easier way to do your Gauss-Legendre quadrature rules in Mathematica: {points, weights} = MapThread[Map, {{2 # - 1 &, 2 # &}, Most[NIntegrate\GaussRuleData[n, prec]]}]`. – J. M. Jan 22 '12 at 16:31
  • In any event: have you seen this? – J. M. Jan 22 '12 at 16:35
  • @J.M, your above proposed method does, unfortunately, not work for prec=Infinity; but thanks for that too. –  Jan 22 '12 at 17:47
  • 2
    In that case, here's a method that works, due to Golub and Welsch: Transpose[MapAt[2(First /@ #)^2 &, Eigensystem[SparseArray[{Band[{2, 1}] -> #, Band[{1, 2}] -> #}, {n, n}]], {2}]] &[Table[k/Sqrt[(2 k - 1)(2 k + 1)], {k, n - 1}]]. – J. M. Jan 22 '12 at 17:59
  • @J.M., nice to know, since it is faster for higher n, do you perhaps happen to have an exact reference? However, I am more interested in something like this for triangles and tets. You don't happen to have something like this for those too? That would be fantastic. –  Jan 22 '12 at 18:17
  • 1
    Here is the paper by Golub and Welsch. I'll dig through my papers and see if there's something for simplices... – J. M. Jan 22 '12 at 23:14
  • Have a look at this as well. – J. M. Jan 23 '12 at 00:41
  • See also this paper, this paper, and SIMPACK. – J. M. Jan 23 '12 at 03:15
  • For future reference, Alg. 824 mentions (among others) DCUTRI (Berntsen and Espelid [1990], Alg 706) and DCUTET (Berntsen et al. [1990,1993], Alg, 720) for triangles and tets. –  Jan 24 '12 at 13:18
  • For future reference, Alg. 706 and 720 mention only one high degree rule, at least I am not able to deduce a general algorithm from that. –  Jan 25 '12 at 08:38
  • The symbolic version of Eigensystem needs to be normalized, to give the same results as in a numerical case: MapAt[2 (First /@ Normalize /@ #)^2 &, Eigensystem[ SparseArray[{Band[{2, 1}] -> #, Band[{1, 2}] -> #}, {n, n}]], {2}] &[Table[k/Sqrt[(2 k - 1) (2 k + 1)], {k, n - 1}]]; –  Jan 25 '12 at 15:25

2 Answers2

8

The Encyclopaedia of Cubature Formulas has an extensive list of techniques for this purpose and associated references.

Peter Brune
  • 1,675
  • 12
  • 20
3

Here is a paper http://journal.library.iisc.ernet.in/vol200405/paper6/rathod.pdf that describes how to map the unit triangle to the standard 2-square in order to calculate the weights and sampling points for the triangle in terms of Gauss-Legendre points for the standard 2-square.

James Custer
  • 397
  • 2
  • 10
  • That's an interesting idea, it looks like for n=2 this needs 4 points, for the typical literature reference for triangles for n=2, 3 points are given. Do you know anything about this? –  Jan 22 '12 at 19:31
  • This stems from the fact that they are using a mapping from the triangle to the square. I can't say anything beyond that since I don't work with triangles (I use quadrilaterals), so I don't know what is usually done in practice. I just found the paper and thought it seemed like a pretty straightforward thing to do. – James Custer Jan 22 '12 at 19:36
  • indeed it is quite straightforward and I'll see that the other papers suggest, but the simplicity for this one and the elegance of using something I already have are a plus for this one. The downside is then the additional function evaluation. Thanks in any case. –  Jan 22 '12 at 19:42
  • another down side is that the points are not symmetrical. –  Jan 24 '12 at 13:19