Generasi otomatis titik integrasi dan bobot untuk segitiga dan tetrahedra

12

Biasanya orang akan berkonsultasi dengan kertas atau buku untuk menemukan titik integrasi dan bobot untuk satuan segitiga dan tetrahedra. Saya mencari metode untuk secara otomatis menghitung poin dan bobot tersebut. Contoh kode Mathematica berikut menghitung bobot dan poin integrasi untuk elemen baris unit (quad / hexahedron):

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])}

Saya mencari kertas / buku yang menjelaskan secara algoritmik bagaimana hal ini dilakukan untuk segitiga dan / atau tetrahedra. Dapatkah seseorang mengarahkan saya ke beberapa informasi tentang ini. Terima kasih.

David Ketcheson
sumber
1
Ada cara yang lebih mudah untuk melakukan aturan quadrature Gauss-Legendre Anda di Mathematica : {points, weights} = MapThread[Map, {{2 # - 1 &, 2 # &}, Most[NIntegrate`GaussRuleData[n, prec]]}].
JM
Dalam hal apa pun: pernahkah Anda melihat ini ?
JM
@ JM, metode yang Anda usulkan di atas, sayangnya, tidak berfungsi untuk prec = Infinity; tapi terima kasih untuk itu juga.
2
Dalam hal ini, inilah metode yang bekerja, karena Golub dan 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}]].
JM
1
Ini kertas karya Golub dan Welsch. Saya akan menggali surat-surat saya dan melihat apakah ada sesuatu untuk kesederhanaan ...
JM

Jawaban:

3

Berikut ini adalah makalah http://journal.library.iisc.ernet.in/vol200405/paper6/rathod.pdf yang menjelaskan cara memetakan satuan segitiga ke standar 2-persegi untuk menghitung bobot dan titik pengambilan sampel untuk segitiga dalam hal poin Gauss-Legendre untuk 2-square standar.

James Custer
sumber
Itu ide yang menarik, sepertinya untuk n = 2 ini membutuhkan 4 poin, untuk referensi literatur tipikal untuk segitiga untuk n = 2, 3 poin diberikan. Apakah Anda tahu sesuatu tentang ini?
Ini berasal dari fakta bahwa mereka menggunakan pemetaan dari segitiga ke bujur sangkar. Saya tidak bisa mengatakan apa-apa selain itu karena saya tidak bekerja dengan segitiga (saya menggunakan segiempat), jadi saya tidak tahu apa yang biasanya dilakukan dalam praktek. Saya baru saja menemukan kertas itu dan berpikir itu seperti hal yang mudah dilakukan.
James Custer
memang itu cukup mudah dan saya akan melihat bahwa makalah lain menyarankan, tetapi kesederhanaan untuk yang satu ini dan keanggunan menggunakan sesuatu yang sudah saya miliki adalah nilai tambah untuk yang ini. Kelemahannya adalah evaluasi fungsi tambahan. Bagaimanapun, terima kasih.
sisi lain adalah bahwa poin tidak simetris.