gauss_points

GAUSS_POINTS
Weights and coordinate for Gauss integration over reference element as a
function of the element type. Specifically, we have rules that integrate
products of gradients of shape functions.
[w,x,y,z] = gauss_points(DIM,ORDER)

NOTE: Integrals of f on specific reference element
       I = int(f(x,y,z)) dV =  w_{i}*f(x_{i},y_{i},z_{i}
Reference triangle dependent on dimension
   DIM=1 - Reference = [-1,1]
   DIM=2 - Reference = [0,1]x[0,1]     
   DIM=3 - Reference = [0,1]x[0,1]x[0,1]

INPUT:
1. DIM - Dimension of integration rule - 1,2,3
2. ORDER - Order integration exact for - rules below are written up
   DIM = 1, ORDER=2,4,6,9
   DIM = 2, ORDER=0,2,4,5,7
   DIM = 3, ORDER=0,2,4
OUTPUT
1. w - the weights matrix of shape function derivatives size(ndim,nshape)
2. x,y,z - evaluation points of the function f

M Crabb - 29.06.2012

0001 function [w,x,y,z]=gauss_points(dim,order)
0002 %GAUSS_POINTS
0003 %Weights and coordinate for Gauss integration over reference element as a
0004 %function of the element type. Specifically, we have rules that integrate
0005 %products of gradients of shape functions.
0006 %[w,x,y,z] = gauss_points(DIM,ORDER)
0007 %
0008 %NOTE: Integrals of f on specific reference element
0009 %       I = int(f(x,y,z)) dV =  w_{i}*f(x_{i},y_{i},z_{i}
0010 %Reference triangle dependent on dimension
0011 %   DIM=1 - Reference = [-1,1]
0012 %   DIM=2 - Reference = [0,1]x[0,1]
0013 %   DIM=3 - Reference = [0,1]x[0,1]x[0,1]
0014 %
0015 %INPUT:
0016 %1. DIM - Dimension of integration rule - 1,2,3
0017 %2. ORDER - Order integration exact for - rules below are written up
0018 %   DIM = 1, ORDER=2,4,6,9
0019 %   DIM = 2, ORDER=0,2,4,5,7
0020 %   DIM = 3, ORDER=0,2,4
0021 %OUTPUT
0022 %1. w - the weights matrix of shape function derivatives size(ndim,nshape)
0023 %2. x,y,z - evaluation points of the function f
0024 %
0025 %M Crabb - 29.06.2012
0026 
0027 
0028 %Switch between cases of dimension and polynomial order
0029 if(dim==1) 
0030     if(order==2)
0031         [w,x,y,z]=dim_1_order_2;
0032     elseif(order==4)
0033         [w,x,y,z]=dim_1_order_4;
0034     elseif(order==6)
0035         [w,x,y,z]=dim_1_order_6;        
0036     elseif(order==9)
0037         [w,x,y,z]=dim_1_order_9;
0038     else
0039        error('Integration order not recognised for dimension'); 
0040     end
0041 elseif(dim==2)
0042     if(order==0)
0043         [w,x,y,z]=dim_2_order_0;       
0044     elseif(order==2)
0045         [w,x,y,z]=dim_2_order_2;        
0046     elseif(order==4)
0047         [w,x,y,z]=dim_2_order_4;        
0048     elseif(order==5)
0049         [w,x,y,z]=dim_2_order_5;        
0050     elseif(order==7)
0051         [w,x,y,z]=dim_2_order_7;
0052     else
0053         error('Integration order not recognised for dimension');
0054     end            
0055 elseif(dim==3)
0056     if(order==0)
0057         [w,x,y,z]=dim_3_order_0;        
0058     elseif(order==2)
0059         [w,x,y,z]=dim_3_order_2;        
0060     elseif(order==4)
0061         [w,x,y,z]=dim_3_order_4;    
0062     else
0063         error('Integration order not recognised for dimension');
0064     end
0065     
0066 else
0067     error('Dimension not recognised');
0068 end
0069 
0070 %1 DIMENSIONAL INTEGRATION RULES
0071 %Order 2 polynomial on triangle edge (need 2 point Gauss rule )
0072 function [w,x,y,z] = dim_1_order_2
0073     %Integration scheme integrates quadratic exactly on [-1,1]
0074     w(1)=1.0; x(1)=1.0/sqrt(3.0);     y(1)=0.0; z(1)=0.0;
0075     w(2)=1.0; x(2)=-1.0/sqrt(3.0);    y(2)=0.0; z(2)=0.0;
0076 end
0077 
0078 %Order 4 polynomial on triangle edge (need 3 point Gauss ruler)
0079 function [w,x,y,z] = dim_1_order_4
0080     %Integration scheme integrates quartic exactly on [-1,1]
0081     w1=sqrt(3.0/5.0);
0082     w(1)=8.0/9.0; x(1)=0.0; y(1)=0.0; z(1)=0.0;
0083     w(2)=5.0/9.0; x(2)=w1;  y(2)=0.0; z(2)=0.0;
0084     w(3)=5.0/9.0; x(3)=-w1; y(3)=0.0; z(3)=0.0;
0085 end
0086 
0087 %Order 6 polynomial on triangle edge (need 4 point Gauss rule)
0088 function [w,x,y,z] = dim_1_order_6
0089     %Integration scheme integrates hexic exactly on [-1,1]
0090     w1=(18.0+sqrt(30.0))/36.0; w2=(18.0-sqrt(30.0))/36.0;
0091     p1=sqrt((3.0-2.0*sqrt(6.0/5.0))/7.0);
0092     p2=sqrt((3.0+2.0*sqrt(6.0/5.0))/7.0);
0093     w(1)=w1;   x(1)=p1;    y(1)=0.0; z(1)=0.0;
0094     w(2)=w1;   x(2)=-p1;   y(2)=0.0; z(2)=0.0;
0095     w(3)=w2;   x(3)=p2;    y(3)=0.0; z(3)=0.0;
0096     w(4)=w2;   x(4)=-p2;   y(4)=0.0; z(4)=0.0;
0097 end
0098 
0099 %Order 9 polynomial on triangle edge (need 5 point Gauss rule)
0100 function [w,x,y,z] = dim_1_order_9
0101     %Integration scheme integrates nonic exactly on [-1,1]
0102     w1=128/225; 
0103     w2=(322.0 + 13.0*sqrt(70.0))/900.0;
0104     w3=(322.0 - 13.0*sqrt(70.0))/900.0;
0105     p1=sqrt( 5.0 - 2.0*sqrt(10.0/7.0) )/3.0;
0106     p2=sqrt( 5.0 + 2.0*sqrt(10.0/7.0) )/3.0;
0107     w(1)=w1;   x(1)=0;     y(1)=0.0; z(1)=0.0;
0108     w(2)=w2;   x(2)=p1;    y(2)=0.0; z(2)=0.0;
0109     w(3)=w2;   x(3)=-p1;   y(3)=0.0; z(3)=0.0;
0110     w(4)=w3;   x(4)=p2;    y(4)=0.0; z(4)=0.0;
0111     w(5)=w3;   x(5)=-p2;   y(5)=0.0; z(5)=0.0;
0112 end
0113 
0114 %2 DIMENSIONAL INTEGRATION RULES
0115 %Order 0 polynomial
0116 function [w,x,y,z] = dim_2_order_0
0117     %Integration scheme integrates linear exactly on unit triangle
0118     w(1)=1.0/2.0; x(1)=1.0/3.0; y(1)=1.0/3.0; z(1)=0.0;
0119 end
0120 
0121 %Order 2 polynomial
0122 function [w,x,y,z] = dim_2_order_2
0123     %Integration scheme integrate quadratic exactly on unit triangle
0124     w(1)=1.0/6.0; x(1)=2.0/3.0; y(1)=1.0/6.0; z(1)=0.0;
0125     w(2)=1.0/6.0; x(2)=1.0/6.0; y(2)=1.0/6.0; z(2)=0.0;
0126     w(3)=1.0/6.0; x(3)=1.0/6.0; y(3)=2.0/3.0; z(3)=0.0;   
0127 end
0128 
0129 %Order 4 polynomial
0130 function [w,x,y,z] = dim_2_order_4
0131 %Integration scheme integrates quartic exactly on unit triangle
0132     w(1)=0.5*0.2233815896780; x(1)=0.1081030181681; y(1)=0.4459484909160; z(1)=0.0;
0133     w(2)=0.5*0.2233815896780; x(2)=0.4459484909160; y(2)=0.1081030181681; z(2)=0.0;
0134     w(3)=0.5*0.2233815896780; x(3)=0.4459484909160; y(3)=0.4459484909160; z(3)=0.0;
0135     w(4)=0.5*0.1099517436553; x(4)=0.8168475729805; y(4)=0.0915762135098; z(4)=0.0;
0136     w(5)=0.5*0.1099517436553; x(5)=0.0915762135098; y(5)=0.0915762135098; z(5)=0.0;
0137     w(6)=0.5*0.1099517436553; x(6)=0.0915762135098; y(6)=0.8168475729805; z(6)=0.0;
0138 end
0139 
0140 %Order 5 polynomial
0141 function [w,x,y,z]=dim_2_order_5
0142     %Integration scheme integrates quintic exactly on unit triangle
0143     w(1)=0.5*0.1259391805448; x(1)=0.1012865073235; y(1)=0.1012865073235; z(1)=0.0;
0144     w(2)=0.5*0.1259391805448; x(2)=0.7974269853531; y(2)=0.1012865073235; z(2)=0.0;
0145     w(3)=0.5*0.1259391805448; x(3)=0.1012865073235; y(3)=0.7974269853531; z(3)=0.0;
0146     w(4)=0.5*0.1323941527885; x(4)=0.4701420641051; y(4)=0.0597158717898; z(4)=0.0;
0147     w(5)=0.5*0.1323941527885; x(5)=0.4701420641051; y(5)=0.4701420641051; z(5)=0.0;
0148     w(6)=0.5*0.1323941527885; x(6)=0.0597158717898; y(6)=0.4701420641051; z(6)=0.0;
0149     w(7)=0.5*0.225;           x(7)=0.3333333333333; y(7)=0.3333333333333; z(7)=0.0;
0150 end
0151 
0152 %Order 7 polynomial
0153 function [w,x,y,z]=dim_2_order_7
0154     %Integration scheme integrate heptic exactly on unit triangle
0155     w(1)=0.5*0.0533472356088;x(1)=0.0651301029022; y(1)=0.0651301029022; z(1)=0.0;
0156     w(2)=0.5*0.0533472356088;x(2)=0.8697397941956; y(2)=0.0651301029022; z(2)=0.0;
0157     w(3)=0.5*0.0533472356088;x(3)=0.0651301029022; y(3)=0.8697397941956; z(3)=0.0;
0158     w(4)=0.5*0.0771137608903;x(4)=0.3128654960049; y(4)=0.0486903154253; z(4)=0.0;
0159     w(5)=0.5*0.0771137608903;x(5)=0.6384441885698; y(5)=0.3128654960049; z(5)=0.0;
0160     w(6)=0.5*0.0771137608903;x(6)=0.0486903154253; y(6)=0.6384441885698; z(6)=0.0;
0161     w(7)=0.5*0.0771137608903;x(7)=0.6384441885698; y(7)=0.0486903154253; z(7)=0.0;
0162     w(8)=0.5*0.0771137608903;x(8)=0.3128654960049; y(8)=0.6384441885698; z(8)=0.0;
0163     w(9)=0.5*0.0771137608903;x(9)=0.0486903154253; y(9)=0.3128654960049; z(9)=0.0;
0164     w(10)=0.5*0.1756152574332;x(10)=0.2603459660790; y(10)=0.2603459660790; z(10)=0.0;
0165     w(11)=0.5*0.1756152574332;x(11)=0.4793080678419; y(11)=0.2603459660790; z(11)=0.0;
0166     w(12)=0.5*0.1756152574332;x(12)=0.2603459660790; y(12)=0.4793080678419; z(12)=0.0;
0167     w(13)=-0.5*0.1495700444677;x(13)=0.3333333333333; y(13)=0.3333333333333; z(13)=0.0;
0168 end
0169 
0170 %3 DIMENSIONAL INTEGRATION RULES
0171 %Order 0 polynomial
0172 function [w,x,y,z] = dim_3_order_0
0173     %Integration scheme integrates exactly linear on unit tetrahedron
0174     w(1)=1.0/6.0; x(1)=1.0/4.0; y(1)=1.0/4.0; z(1)=1.0/4.0;
0175 end
0176 
0177 %Order 2 polynomial
0178 function [w,x,y,z] = dim_3_order_2
0179     %Integration scheme integrates exactly quadratic on unit tetrahedron
0180     a=(5+3*sqrt(5))/20; b=(5-sqrt(5))/20;
0181     w(1)=1.0/24.0; x(1)=a; y(1)=b; z(1)=b;
0182     w(2)=1.0/24.0; x(2)=b; y(2)=a; z(2)=b;
0183     w(3)=1.0/24.0; x(3)=b; y(3)=b; z(3)=a;
0184     w(4)=1.0/24.0; x(4)=b; y(4)=b; z(4)=b;
0185 end
0186 
0187 %Order 4 polynomial
0188 function [w,x,y,z] = dim_3_order_4
0189     %Integration scheme integrates exactly quartic on unit tetrahedron
0190     a=0.25;
0191     w(1)=-0.0131555555555555556;
0192     x(1)=a; y(1)=a; z(1)=a;
0193     
0194     a=0.785714285714285714; b=0.0714285714285714285;
0195     w(2)=0.00762222222222222222;
0196     x(2)=a; y(2)=b; z(2)=b;
0197     w(3)=0.00762222222222222222;
0198     x(3)=b; y(3)=a; z(3)=b;
0199     w(4)=0.00762222222222222222;
0200     x(4)=b; y(4)=b; z(4)=a;
0201     w(5)=0.00762222222222222222;
0202     x(5)=b; y(5)=b; z(5)=b;
0203     
0204     a=0.399403576166799219; b=0.100596423833200785;
0205     w(6)=0.0248888888888888889;
0206     x(6)=a; y(6)=a; z(6)=b;
0207     w(7)=0.0248888888888888889;
0208     x(7)=a; y(7)=b; z(7)=b;
0209     w(8)=0.0248888888888888889;
0210     x(8)=b; y(8)=b; z(8)=a;
0211     w(9)=0.0248888888888888889;
0212     x(9)=b; y(9)=a; z(9)=b;
0213     w(10)=0.0248888888888888889;
0214     x(10)=b; y(10)=a; z(10)=a;
0215     w(11)=0.0248888888888888889;
0216     x(11)=a; y(11)=b; z(11)=a;
0217 end
0218 
0219 end