Generación automática de puntos de integración y pesos para triángulos y tetraedros.

12

Por lo general, se consultaría un documento o libro para encontrar puntos de integración y pesos para triángulos unitarios y tetraedros. Estoy buscando un método para calcular automáticamente dichos puntos y pesos. El siguiente ejemplo de código de Mathematica calcula los pesos y puntos de integración para elementos de línea de unidad (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])}

Estoy buscando un documento / libro que describa algorítmicamente cómo se hace esto para triángulos y / o tetraedros. ¿Alguien puede señalarme alguna información sobre esto? Gracias.

symbolic-computation reference-request David Ketcheson
fuente
1
Hay una manera más fácil de hacer sus reglas de cuadratura de Gauss-Legendre en Mathematica : {points, weights} = MapThread[Map, {{2 # - 1 &, 2 # &}, Most[NIntegrate`GaussRuleData[n, prec]]}].
JM
En cualquier caso: ¿has visto esto ?
JM
@JM, su método propuesto anteriormente, desafortunadamente, no funciona para prec = Infinity; pero gracias por eso también.
2
En ese caso, aquí hay un método que funciona, debido a Golub y 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
Aquí está el artículo de Golub y Welsch. Examinaré mis papeles y veré si hay algo para simplificar ...
JM

Respuestas:

3

Aquí hay un documento http://journal.library.iisc.ernet.in/vol200405/paper6/rathod.pdf que describe cómo mapear el triángulo unitario al cuadrado estándar de 2 para calcular los pesos y los puntos de muestreo para triángulo en términos de puntos de Gauss-Legendre para el estándar de 2 cuadrados.

James Custer
fuente
Esa es una idea interesante, parece que para n = 2 esto necesita 4 puntos, para la referencia de literatura típica para triángulos para n = 2, se dan 3 puntos. ¿Sabes algo sobre esto?
Esto se debe al hecho de que están utilizando una asignación del triángulo al cuadrado. No puedo decir nada más allá de eso, ya que no trabajo con triángulos (uso cuadriláteros), así que no sé qué se hace generalmente en la práctica. Acabo de encontrar el papel y pensé que parecía algo bastante sencillo de hacer.
James Custer
de hecho, es bastante sencillo y veré que los otros documentos sugieren, pero la simplicidad de este y la elegancia de usar algo que ya tengo son una ventaja para este. La desventaja es entonces la evaluación de la función adicional. Gracias de todas formas.
Otra desventaja es que los puntos no son simétricos.