Code: Quadrature on Tetrahedra

From CFD-Wiki

Some Gauss quadrature rules for the right-tetrahedron with corners at (0,0,0), (1,0,0), (0,1,0), (0,0,1) (Matlab).

function [xa,ya,za,wt]=TetQuadDat(n)
% Quadrature data for tetrahedron
% Refs
%  P Keast, Moderate degree tetrahedral quadrature formulas, CMAME 55: 339-348 (1986)
%  O. C. Zienkiewicz, The Finite Element Method,  Sixth Edition,

switch n ;
  
   case 1   %  Z4, K1,  N=4
 xa= [0.5854101966249685, 0.1381966011250105, 0.1381966011250105, 0.1381966011250105]; 
 ya= [0.1381966011250105, 0.1381966011250105, 0.1381966011250105, 0.5854101966249685];
 za= [0.1381966011250105, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105];
 wt= [0.2500000000000000, 0.2500000000000000, 0.2500000000000000, 0.2500000000000000]/6;

   case 2   %  Z5, K2  N=5
 xa= [ 0.2500000000000000, 0.5000000000000000, 0.1666666666666667, 0.1666666666666667, 0.1666666666666667];
 ya= [ 0.2500000000000000, 0.1666666666666667, 0.1666666666666667, 0.1666666666666667, 0.5000000000000000];
 za= [ 0.2500000000000000, 0.1666666666666667, 0.1666666666666667, 0.5000000000000000, 0.1666666666666667];
 wt= [-0.8000000000000000, 0.4500000000000000, 0.4500000000000000, 0.4500000000000000, 0.4500000000000000]/6;

   case 3   %  K4  N=11
 xa= [0.2500000000000000, 0.7857142857142857, 0.0714285714285714, 0.0714285714285714, 0.0714285714285714, ...
     0.1005964238332008, 0.3994035761667992, 0.3994035761667992, 0.3994035761667992, 0.1005964238332008, 0.1005964238332008];
 ya= [0.2500000000000000, 0.0714285714285714, 0.0714285714285714, 0.0714285714285714, 0.7857142857142857, ...
     0.3994035761667992, 0.1005964238332008, 0.3994035761667992, 0.1005964238332008, 0.3994035761667992, 0.1005964238332008];
 za= [0.2500000000000000, 0.0714285714285714, 0.0714285714285714, 0.7857142857142857, 0.0714285714285714, ...
     0.3994035761667992, 0.3994035761667992, 0.1005964238332008, 0.1005964238332008, 0.1005964238332008, 0.3994035761667992];
 wt=[-0.0789333333333333, 0.0457333333333333, 0.0457333333333333, 0.0457333333333333, 0.0457333333333333, ...
     0.1493333333333333, 0.1493333333333333, 0.1493333333333333, 0.1493333333333333, 0.1493333333333333, 0.1493333333333333]/6;

   case 4  %K6  N=15
 xa=[0.2500000000000000, 0.0000000000000000, 0.3333333333333333, 0.3333333333333333, 0.3333333333333333, ...
     0.7272727272727273, 0.0909090909090909, 0.0909090909090909, 0.0909090909090909, 0.4334498464263357, ...
     0.0665501535736643, 0.0665501535736643, 0.0665501535736643, 0.4334498464263357, 0.4334498464263357];
 ya=[0.2500000000000000, 0.3333333333333333, 0.3333333333333333, 0.3333333333333333, 0.0000000000000000, ...
     0.0909090909090909, 0.0909090909090909, 0.0909090909090909, 0.7272727272727273, 0.0665501535736643, ...
     0.4334498464263357, 0.0665501535736643, 0.4334498464263357, 0.0665501535736643, 0.4334498464263357];
 za=[0.2500000000000000, 0.3333333333333333, 0.3333333333333333, 0.0000000000000000, 0.3333333333333333, ...
     0.0909090909090909, 0.0909090909090909, 0.7272727272727273, 0.0909090909090909, 0.0665501535736643, ...
     0.0665501535736643, 0.4334498464263357, 0.4334498464263357, 0.4334498464263357, 0.0665501535736643];
 wt=[0.1817020685825351, 0.0361607142857143, 0.0361607142857143, 0.0361607142857143, 0.0361607142857143, ...
     0.0698714945161738, 0.0698714945161738, 0.0698714945161738, 0.0698714945161738, 0.0656948493683187, ...
     0.0656948493683187, 0.0656948493683187, 0.0656948493683187, 0.0656948493683187, 0.0656948493683187]/6;
end   % switch n

For Python, those above and many more can be found in quadrature[1]. For example,

import quadrature

scheme = quadrature.tetrahedron.ShunnHam(4)

print(scheme.points)
print(scheme.weights)

gives

[[ 0.03235259  0.03235259  0.90294222]
 [ 0.03235259  0.90294222  0.03235259]
 [ 0.90294222  0.03235259  0.03235259]
 [ 0.03235259  0.03235259  0.03235259]
 [ 0.06036044  0.26268258  0.61659653]
 [ 0.26268258  0.06036044  0.61659653]
 [ 0.06036044  0.06036044  0.61659653]
 [ 0.26268258  0.61659653  0.06036044]
 [ 0.06036044  0.61659653  0.06036044]
 [ 0.61659653  0.06036044  0.06036044]
 [ 0.06036044  0.61659653  0.26268258]
 [ 0.61659653  0.06036044  0.26268258]
 [ 0.06036044  0.06036044  0.26268258]
 [ 0.61659653  0.26268258  0.06036044]
 [ 0.06036044  0.26268258  0.06036044]
 [ 0.26268258  0.06036044  0.06036044]
 [ 0.3097693   0.3097693   0.07069209]
 [ 0.3097693   0.07069209  0.3097693 ]
 [ 0.07069209  0.3097693   0.3097693 ]
 [ 0.3097693   0.3097693   0.3097693 ]]
[ 0.00706707  0.00706707  0.00706707  0.00706707  0.04699867  0.04699867
  0.04699867  0.04699867  0.04699867  0.04699867  0.04699867  0.04699867
  0.04699867  0.04699867  0.04699867  0.04699867  0.10193692  0.10193692
  0.10193692  0.10193692]