PFV quadature rules
From CFD-Wiki
Revision as of 16:51, 15 March 2013 by Jonas Holdeman (Talk | contribs)
function [ xa,ya,wt,nq ] = GQuad2( ord ) %GQuad2 % Returns arrays of integration points and weights for specified degree % Usage: % [xa,ya,wt,nq] = GQuad2( ord ) % where % ord - highest degree of term integrated exactly on interval [-1,1] % if ord<1, error('integration order < 1'); elseif ord>15, error('max quad data order = 15'); end nt = fix(ord/2)+1; switch nt; case {1, 2} % degree 3 Aq=[-0.5773502691896257, 0.5773502691896257]; Hq=[ 1.0, 1.0]; case 3 % degree 5 Aq=[-0.7745966692414833, 0.000000000000000, 0.7745966692414833]; Hq=[ 0.5555555555555558, 0.8888888888888888,0.5555555555555558]; case 4 % degree 7 Aq=[-0.8611363115940526,-0.3399810435848563, 0.3399810435848563, ... 0.8611363115940526]; Hq=[ 0.3478548451374538, 0.6521451548625461, 0.6521451548625461, ... 0.3478548451374538]; case 5 % degree 9 Aq=[-0.9061798459386640,-0.5384693101056831, 0.0000000000000000, ... 0.5384693101056831, 0.9061798459386640]; Hq=[ .23692688505618909, .47862867049936647, .56888888888888889, ... .47862867049936647, .23692688505618909]; case 6 % degree 11 Aq=[-.932469514203152,-.661209386466265,-.238619186083197, ... .238619186083197, .661209386466265, .932469514203152]; Hq=[ .171324492379170, .360761573048139, .467913934572691, ... .467913934572691, .360761573048139, .171324492379170]; case 7 % degree 13 Aq=[-0.9491079123427585,-0.7415311855993945,-0.4058451513773972, ... 0.0000000000000000, 0.4058451513773972, 0.7415311855993945, ... 0.9491079123427585]; Hq=[ 0.1294849661688699, 0.2797053914892767, 0.3818300505051190, ... 0.4179591836734694, 0.3818300505051190, 0.2797053914892767, ... 0.1294849661688699]; case 8 % degree 15 Aq=[-0.9602898564975363,-0.7966664774136267,-0.5255324099163290, ... -0.1834346424956498, 0.1834346424956498, 0.5255324099163290, ... 0.7966664774136267, 0.9602898564975363]; Hq=[ 0.1012285362903761, 0.2223810344533745, 0.3137066458778873, ... 0.3626837833783620, 0.3626837833783620, 0.3137066458778873, ... 0.2223810344533745, 0.1012285362903761]; end nt=size(Aq,2); xa = repmat(Aq,nt,1); ya = xa'; wt = repmat(Hq,nt,1); wt = wt.*wt'; nq = nt*nt; return; end