function [X,W] = inttet(p)
% INTTET
% The function gives the quadrature points X as a (3,n)-matrix and
% the weights W as an (n,1)-matrix for integration over the unit
% tetrahedra whose nodes are (0,0,0), (1,0,0), (0,1,0) and (0,0,1).
%
% The input argument p indicates the maximal degree of the polynomials
% of two variables which will be integrated exactly (p<15). By a linear
% transformation, X and W were obtained from the paper:
%
% Lin-bo Zhang, Tao Cui and Hui Liu,
% A set of symmetric quadrature rules on triangles and tetrahedra,
% J. Comp. Math., 26(3), 2008.
%
% SYNTAX: [X,W] = inttet(p)
%
% IN: p degree of polynomial which is integrated exactly
%
% OUT: X points as a (3,n)-matrix
% W weighs as a (n,1)-matrix
%
%
% Author: Immanuel Anjam, University of Jyvaskyla
% Contact: immanuel.anjam@gmail.com
% Date: 01.06.2010
% Version: 1.0
%
X = [];
W = [];
%the quadrature for p=4 and p=5 contains the same amount of points,
%so let us use the better one for p=4 also
if ( p == 4 )
p = 5;
end
switch p
case 1
% 1 point
w = 1.0;
[X,W] = s4(w);
case 2
% 4 points
w = 1/4;
a = 0.1381966011250105151795413165634361;
[X,W] = s31(a,w);
case 3
% 8 points
w = 0.1385279665118621423236176983756412;
a = 0.3280546967114266473358058199811974;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = 0.1114720334881378576763823016243588;
a = 0.1069522739329306827717020415706165;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
%{
case 4
% 14 points
w = 0.0734930431163619493435869458636788;
a = 0.0927352503108912262865589206603214;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = 0.1126879257180158503650149284763889;
a = 0.3108859192633006097581474949404033;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = 0.0425460207770814668609320837732882;
a = 0.0455037041256496500000000000000000;
[points,weighs] = s22(a,w);
X = [X points];
W = [W; weighs];
%}
case 5
% 14 points
w = .1126879257180158507991856523332863;
a = .3108859192633006097973457337634578;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0734930431163619495437102054863275;
a = .0927352503108912264023239137370306;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0425460207770814664380694281202574;
a = .0455037041256496494918805262793394;
[points,weighs] = s22(a,w);
X = [X points];
W = [W; weighs];
case 6
% 24 points
w = .0399227502581674920996906275574800;
a = .2146028712591520292888392193862850;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0100772110553206429480132374459369;
a = .0406739585346113531155794489564101;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0553571815436547220951532778537260;
a = .3223378901422755103439944707624921;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0482142857142857142857142857142857;
a = .0636610018750175252992355276057270;
b = .6030056647916491413674311390609397;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
case 7
% 36 points
w = .0061834158394585176827283275896323;
a = .0406107071929452723515318677627212;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0785146502738723588282424885149114;
a = .1787522026964984761546314943983834;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0447395776143224792777362432057442;
a = .3249495905373373335715573286644841;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0121651445922912935604045907106762;
a = .1340777379721611918326213913378565;
b = .7270125070093171000000000000000000;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0280223074984909211766930561858945;
a = .0560275404597284769655799958528421;
b = .3265740998664049580757011076659178;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
case 8
% 46 points
w = .0063971477799023213214514203351730;
a = .0396754230703899012650713295393895;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0401904480209661724881611584798178;
a = .3144878006980963137841605626971483;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0243079755047703211748691087719226;
a = .1019866930627033000000000000000000;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0548588924136974404669241239903914;
a = .1842036969491915122759464173489092;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0357196122340991824649509689966176;
a = .0634362877545398924051412387018983;
[points,weighs] = s22(a,w);
X = [X points];
W = [W; weighs];
w = .0071831906978525394094511052198038;
a = .0216901620677280048026624826249302;
b = .7199319220394659358894349533527348;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0163721819453191175409381397561191;
a = .2044800806367957142413355748727453;
b = .5805771901288092241753981713906204;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
case 9
% 61 points
w = .0564266931795062065887150432761254;
[points,weighs] = s4(w);
X = [X points];
W = [W; weighs];
w = .0033410950747134804029997443047177;
a = .0340221770010448664654037088787676;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0301137547687737639073142384315749;
a = .3227703335338005253913766832549640;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0064909609200615346357621168945686;
a = .0604570774257749300000000000000000;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0098092858682545864319687425925550;
a = .4553629909472082118003081504416430;
b = .0056831773653301799061001601457447;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0281191538233654725516326174252926;
a = .1195022553938258009779737046961144;
b = .4631168324784899409762244936577296;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0078945869083315007683414920096088;
a = .0280219557834011581550575066541237;
b = .7252060768398674887385659542848099;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0194928120472399967169721944892460;
a = .1748330320115746157853246459722452;
b = .6166825717812564045706830909795407;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
case 10
% 81 points
w = .0473997735560207383847388211780511;
[points,weighs] = s4(w);
X = [X points];
W = [W; weighs];
w = .0269370599922686998027641610048821;
a = .3122500686951886477298083186868275;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0098691597167933832345577354301731;
a = .1143096538573461505873711976536504;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0003619443443392536242398783848085;
a = .0061380088247907478475937132484154;
b = .9429887673452048661976305869182508;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0101358716797557927885164701150168;
a = .0327794682164426707747210203323242;
b = .3401847940871076327889879249496713;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0113938812201952316236209348807143;
a = .4104307392189654942878978442515117;
b = .1654860256196110516044901244445264;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0065761472770359041674557402004507;
a = .0324852815648230478355149399784262;
b = .1338521522120095130978284359645666;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0257397319804560712790360122596547;
a = .1210501811455894259938950015950505;
b = .4771903799042803505441064082969072;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0129070357988619906392954302494990;
a = .1749793421839390242849492265283104;
b = .6280718454753660106932760722179097;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
case 11
% 109 points
w = .0394321080286588635073303344912044;
[points,weighs] = s4(w);
X = [X points];
W = [W; weighs];
w = .0156621262272791131500885627687651;
a = .1214913677765337944977023099080722;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0033321723749014081444092361540149;
a = .0323162591510728963539544520895810;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0140260774074897474374913609976924;
a = .3249261497886067978128419024144220;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0010859075293324663068220983772355;
a = .0041483569716600120000000000000100;
b = .5982659967901863502054538427761778;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0202359604306631789111165731654084;
a = .2246246106763771414144751511649864;
b = .4736622878323495714083696692020524;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0117902148721258635368493804677018;
a = .0519050877725656967442272164426589;
b = .5631447779082798987371019763030571;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0076903149825212959011315780207389;
a = .1349301312162402042237591723429930;
b = .7083588307858189538569950051271300;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0044373057034592039047307260214396;
a = .0251911921082524729200511850653055;
b = .7837195073400773754305740342999090;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0114295484671840404107705525985940;
a = .3653187797817336139693319800988672;
b = .1346039083168658000000000000000100;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0061856401712178114128192550838953;
a = .5229075395099384729652169275860292;
b = .1407536305436959018425391394912785;
c = .0097624381964526155082922803899778;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
case 12
% 140 points
w = .0127676377009707415020377859651250;
a = .1152997443514801453045572073891591;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0161211042379092682185815448957576;
a = .2023362822405909000000000000000100;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0003716126985784422000425581898608;
a = .0117175979576199515124790675483140;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0197174417866854576395533090381887;
a = .3133064413678010672776027996445893;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0025713909308627183621823475944855;
a = .2500057301155837000000000000000100;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0038172478705105759057531841278333;
a = .0209954743507580066902018252705902;
[points,weighs] = s22(a,w);
X = [X points];
W = [W; weighs];
w = .0120872270776631131786031841931461;
a = .1517740182474501000000000000000100;
[points,weighs] = s22(a,w);
X = [X points];
W = [W; weighs];
w = .0031058611584347334343168814992962;
a = .0244197787434353647831400090476166;
b = .8483292846978728506452088674348157;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0054595313364710306691274212676944;
a = .2562070985320183089638201070856221;
b = .4824873738738488478028928967297354;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0021428997484969975066685209365595;
a = .0167903209796029906147179602885794;
b = .6947719423657559269594985098841772;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0055246714672578296224493009816508;
a = .1261608211398720423997070384689592;
b = .7254104893029481189748595052126338;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0085369566944991804298517783667220;
a = .4314351745263798472167069506637196;
b = .1127219398928524152095997721100754;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0115101778483233069733364412340329;
a = .5016700624625056974751550716847613;
b = .2724718028695223917835104675306045;
c = .0720743288072989146501594845633582;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
w = .0052038786528856136039679242125245;
a = .2616448545378187456694550500639680;
b = .0862922919470617319174235194435249;
c = .0205654106558761383006248976271090;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
case 13
% 171 points
w = .0150136877730831467506297063161598;
[points,weighs] = s4(w);
X = [X points];
W = [W; weighs];
w = .0182252092801734253237906894149010;
a = .1552160935190895031411578433570474;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0070061092177414642403851869392631;
a = .3301226633396736002443319259519678;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0164235497439495482954057310790553;
a = .1668064038938624992893778260114423;
[points,weighs] = s22(a,w);
X = [X points];
W = [W; weighs];
w = .0051206100963605970726259694970217;
a = .0249237885477736177970140037486009;
[points,weighs] = s22(a,w);
X = [X points];
W = [W; weighs];
w = .0111966986529049163438203208635196;
a = .0971976299157510014307224371624082;
[points,weighs] = s22(a,w);
X = [X points];
W = [W; weighs];
w = .0156191497333799540095381130243197;
a = .2478592901573625669274691062082793;
b = .4336532423568514471872606143476738;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0024844230133164744190405677633847;
a = .0222315960818670029087952186089293;
b = .8369003204037340051450948659569859;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0016385985348182389384452530944075;
a = .1072786933130534104915045963958480;
b = .7749803059750018075658787727417929;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0059030304401249219717191465553587;
a = .1981768438839898114233184058214276;
b = .5875693057822053025917201790359592;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0110220824582180524044509798920153;
a = .0691792434773793164773253434746550;
b = .6042000666600664470793526487111530;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0004064518399641782258515551275585;
a = .0231147194719331600000000000000100;
b = .9308757927924442486492022888288831;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0026879699729685420974578192665173;
a = .1178892875101960892229011747064425;
b = .1165153642254072000000000000000100;
c = .0420240011255154209567663430372000;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
w = .0019795048055267119053189467551074;
a = .6770327986022842635503222132674659;
b = .0461653760246197108345804112217608;
c = .0008443403189050397572989969213590;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
w = .0054463191814257912094318704010867;
a = .4848900886736331220108009415479083;
b = .3588829429552020157242364690942109;
c = .1381828349176287299695508090791236;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
case 14
% 236 points
w = .0040651136652707670436208836835636;
a = .3272533625238485639093096692685289;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0022145385334455781437599569500071;
a = .0447613044666850808837942096478842;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0058134382678884505495373338821455;
a = .0861403311024363536537208740298857;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0196255433858357215975623333961715;
a = .2087626425004322968265357083976176;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0003875737905908214364538721248394;
a = .0141049738029209600635879152102928;
[points,weighs] = s31(a,w);
X = [X points];
W = [W; weighs];
w = .0116429719721770369855213401005552;
a = .1021653241807768123476692526982584;
b = .5739463675943338202814002893460107;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0052890429882817131317736883052856;
a = .4075700516600107157213295651301783;
b = .0922278701390201300000000000000000;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0018310854163600559376697823488069;
a = .0156640007402803585557586709578084;
b = .7012810959589440327139967673208426;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0082496473772146452067449669173660;
a = .2254963562525029053780724154201103;
b = .4769063974420887115860583354107011;
[points,weighs] = s211(a,b,w);
X = [X points];
W = [W; weighs];
w = .0030099245347082451376888748208987;
a = .3905984281281458000000000000000000;
b = .2013590544123922168123077327235092;
c = .0161122880710300298578026931548371;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
w = .0008047165617367534636261808760312;
a = .1061350679989021455556139029848079;
b = .0327358186817269284944004077912660;
c = .0035979076537271666907971523385925;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
w = .0029850412588493071187655692883922;
a = .5636383731697743896896816630648502;
b = .2302920722300657454502526874135652;
c = .1907199341743551862712487790637898;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
w = .0056896002418760766963361477811973;
a = .3676255095325860844092206775991167;
b = .2078851380230044950717102125250735;
c = .3312104885193449000000000000000000;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
w = .0041590865878545715670013980182613;
a = .7192323689817295295023401840796991;
b = .1763279118019329762157993033636973;
c = .0207602362571310090754973440611644;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
w = .0007282389204572724356136429745654;
a = .5278249952152987298409240075817276;
b = .4372890892203418165526238760841918;
c = .0092201651856641949463177554949220;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
w = .0054326500769958248216242340651926;
a = .5483674544948190728994910505607746;
b = .3447815506171641228703671870920331;
c = .0867217283322215394629438740085828;
[points,weighs] = s1111(a,b,c,w);
X = [X points];
W = [W; weighs];
end
%multiply the weighs by the volume of unit tetrahedron
W = 1/6 .* W;
end
%#########################################################################
function [X,W] = s4(w)
%the first star contains only this one point
X = [1/4 1/4 1/4]';
W = w;
end
%#########################################################################
function [X,W] = s31(a,w)
%for this star, (and the other remaining 3 stars) we first generate the
%barycentric co-ordinates, which contain 4 dimensions. the points are
%obtained by taking all the unique 3 dimensional permutations from the
%barycentric co-ordinates
baryc = [a a a (1-3*a)];
X = unique( perms(baryc), 'rows' );
X = X(:,1:3)';
W = w * ones( size(X,2), 1 );
end
%#########################################################################
function [X,W] = s22(a,w)
baryc = [a a (1/2 - a) (1/2 - a)];
X = unique( perms(baryc), 'rows' );
X = X(:,1:3)';
W = w * ones( size(X,2), 1 );
end
%#########################################################################
function [X,W] = s211(a,b,w)
baryc = [a a b (1 - 2*a - b)];
X = unique( perms(baryc), 'rows' );
X = X(:,1:3)';
W = w * ones( size(X,2), 1 );
end
%#########################################################################
function [X,W] = s1111(a,b,c,w)
baryc = [a b c (1 - a - b - c)];
X = unique( perms(baryc), 'rows' );
X = X(:,1:3)';
W = w * ones( size(X,2), 1 );
end
