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 ---------------------...")
(Matlab code)
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];
-
</pre>
+
za= [0.1381966011250105, 0.1381966011250105, 0.5854101966249685, 0.1381966011250105];
-
<pre>
+
wt= [0.2500000000000000, 0.2500000000000000, 0.2500000000000000, 0.2500000000000000];
-
              4-point quadrature,  k=3
+
 
-
  ------------------------------------------------------
+
  case 2   % Z5, K2 N=5
-
  x-coordinate y-coordinate  z-coordinate    6 x Weight
+
xa= [0.2500000000000000, 0.5000000000000000, 0.1666666666666667, 0.1666666666666667, 0.1666666666666667];
-
  ------------------------------------------------------
+
ya= [0.2500000000000000, 0.1666666666666667, 0.1666666666666667, 0.1666666666666667, 0.5000000000000000];
-
      1/4          1/4          1/4          -8/10
+
za= [0.2500000000000000, 0.1666666666666667, 0.1666666666666667, 0.5000000000000000, 0.1666666666666667];
-
      1/2          1/6          1/6          9/20
+
wt=-[0.8000000000000000, 0.4500000000000000, 0.4500000000000000, 0.4500000000000000, 0.4500000000000000];
-
      1/6          1/2          1/6          9/20
+
 
-
      1/6          1/6          1/2          9/20
+
  case 3   %  K4  N=11
-
      1/6          1/6          1/6          9/20
+
xa= [0.2500000000000000, 0.7857142857142857, 0.0714285714285714, 0.0714285714285714, 0.0714285714285714, ...
-
   ------------------------------------------------------
+
    0.1005964238332008, 0.3994035761667992, 0.3994035761667992, 0.3994035761667992, 0.1005964238332008, 0.1005964238332008];
-
</pre>
+
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];
 +
 
 +
  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];
 +
end   % switch n
</pre>
</pre>

Revision as of 19:30, 9 July 2011

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];

   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];

   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];

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