Code covered by the BSD License  

Highlights from
inttet

image thumbnail

inttet

by

 

11 Jul 2012 (Updated )

3D integration quadrature for unit tetrahedron

inttet(p)
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

Contact us