CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > Code: Quadrature on Tetrahedra

Code: Quadrature on Tetrahedra

From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
(Created page with "Some Gauss quadrature rules for the right-tetrahedron with corners at (0,0,0), (1,0,0), (0,1,0), (0,0,1). <pre> 4-point quadrature, k=2 ---------------------...")
 
(3 intermediate revisions not shown)
Line 1: Line 1:
-
Some Gauss quadrature rules for the right-tetrahedron with corners at (0,0,0), (1,0,0), (0,1,0), (0,0,1).  
+
Some Gauss quadrature rules for the right-tetrahedron with corners at (0,0,0), (1,0,0), (0,1,0), (0,0,1) (Matlab).  
<pre>
<pre>
-
              4-point quadrature,  k=2
+
function [xa,ya,za,wt]=TetQuadDat(n)
-
   ------------------------------------------------------
+
% Quadrature data for tetrahedron
-
            a=0.58541020,   b=0.13819660
+
% Refs
-
  ------------------------------------------------------
+
%  P Keast, Moderate degree tetrahedral quadrature formulas, CMAME 55: 339-348 (1986)
-
   x-coordinate y-coordinate z-coordinate    6 x Weight
+
%  O. C. Zienkiewicz, The Finite Element Method,  Sixth Edition,
-
   ------------------------------------------------------
+
 
-
        a            b            b          0.25
+
switch n ;
-
        b            a            b          0.25
+
    
-
        b            b            a          0.25
+
  case 1   %  Z4, K1,  N=4
-
        a            a            a          0.25
+
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
</pre>
</pre>
 +
 +
 +
For Python, those above and many more can be found in quadrature[https://github.com/nschloe/quadrature]. For example,
<pre>
<pre>
-
              4-point quadrature,  k=3
+
import quadrature
-
  ------------------------------------------------------
+
 
-
  x-coordinate  y-coordinate  z-coordinate    6 x Weight
+
scheme = quadrature.tetrahedron.ShunnHam(4)
-
  ------------------------------------------------------
+
 
-
      1/4           1/4          1/4          -8/10
+
print(scheme.points)
-
      1/2          1/6          1/6          9/20
+
print(scheme.weights)
-
      1/6          1/2          1/6          9/20
+
-
      1/6          1/6          1/2          9/20
+
-
      1/6          1/6          1/6          9/20
+
-
  ------------------------------------------------------
+
</pre>
</pre>
 +
gives
 +
<pre>
 +
[[ 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]
</pre>
</pre>

Latest revision as of 21:13, 7 November 2016

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]
Retrieved from "https://openfoam.cfd-online.com/Wiki/Code:_Quadrature_on_Tetrahedra"
My wiki