From: David Wilkinson Subject: Re: Gaussian integration formula on triangle Date: Thu, 2 Sep 1999 20:26:44 +0100 Newsgroups: sci.math.num-analysis Have you tried Abramowitz & Stegun, Handbook of Mathematical Functions, page 893. They give Gaussian Integration for a triangle with 7 points and errors of order h^6 In article , zt writes >I am looking for a high order Gaussian integration formula on a triangle. >The formula should be exact for all x^i y^j (0<=i+j<=n), and n is at least >8. I am interested in a real 2D formula with symmetric abscissas and >weights, not tensor product of 1D Gaussian integrals. Thanks. > >Z. Tan >zt@cfdrc.com > > > -- David Wilkinson ============================================================================== From: Olaf Michelsson Subject: Re: Gaussian integration formula on triangle Date: Mon, 06 Sep 1999 16:17:59 +0200 Newsgroups: sci.math.num-analysis zt wrote: > > I am looking for a high order Gaussian integration formula on a triangle. > The formula should be exact for all x^i y^j (0<=i+j<=n), and n is at least > 8. I am interested in a real 2D formula with symmetric abscissas and > weights, not tensor product of 1D Gaussian integrals. Thanks. > > Z. Tan > zt@cfdrc.com search for DCUTRI in toms at netlib. It includes an adaptive integrator for triangles, where a 37 point Gauss-qudrature scheme is used. This formula is fully symmetric and has only positive weights inside the triangle. -- Olaf Michelsson mailto:michelsson@e-technik.tu-ilmenau.de Faculty E/I Dept. of Theoretical Electrical Engineering Technical University of Ilmenau, Germany ============================================================================== From: smithj@smith-james-p.jsc.nasa.gov (James P. Smith) Subject: Re: Gaussian integration formula on triangle Date: 2 Sep 1999 20:24:47 GMT Newsgroups: sci.math.num-analysis This is what I used in the programming of triangular p-elements with points (0,0), (1,0), and (0,1) in the parametric space. Went up to p=8. SUBROUTINE GAUSSPTRI(ZETA,ETA,WEIGHT,PORDER,NUMGAU) C C Copyright (c) 1997 National Aeronautics and Space Administration. C No copyright claimed in USA under Title 17, U.S. Code. All other C rights reserved. C C James P. Smith C NASA/Johnson Space Center C Structures and Dynamics Branch C james.p.smith1@jsc.nasa.gov C C THIS ROUTINE ESTABLISHES THE GAUSS POINTS AND WEIGHTS FOR A TRIANGULAR C REGION. A DIFFERENCE STORAGE SCHEME IS REQUIRED FOR THESE POINTS C SINCE IT'S NOT AN N-INTEGRATION POINT BY N-INTEGRATION POINT SCHEME. C INCLUDE 'PARAMS.INC' IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION ZETA(52),ETA(52),WEIGHT(52),T(5) INTEGER PORDER C C ZETA - GAUSS POINT IN ZETA DIRECTION FOR A TRIANGULAR REGION C C ETA - GAUSS POINT IN ETA DIRECTION FOR A TRIANGULAR REGION C C WEIGHT - WEIGHTING FACTOR FOR THE GAUSS POINT C C T - TEMPORARY STORAGE C IF (PORDER .EQ. 1) NUMGAU=16 IF (PORDER .EQ. 2 .OR. PORDER .EQ. 3) NUMGAU=16 IF (PORDER .EQ. 4) NUMGAU=16 IF (PORDER .EQ. 5) NUMGAU=25 IF (PORDER .EQ. 6) NUMGAU=42 IF (PORDER .EQ. 7) NUMGAU=42 IF (PORDER .EQ. 8) NUMGAU=52 IF (NUMGAU .EQ. 16) THEN ZETA(1)=1.D0/3.D0 ETA(1)=1.D0/3.D0 WEIGHT(1)=.144315607677787D0 ZETA(2)=.081414823414554D0 ETA(2)=.459292588292723D0 WEIGHT(2)=.095091634267285D0 ZETA(5)=.658861384496480D0 ETA(5)=.170569307751760D0 WEIGHT(5)=.103217370534718D0 ZETA(8)=.898905543365938D0 ETA(8)=.050547228317031D0 WEIGHT(8)=.032458497623198D0 DO 10 I=3,9,3 ZETA(I)=ETA(I-1) ETA(I)=ZETA(I-1) WEIGHT(I)=WEIGHT(I-1) ZETA(I+1)=ETA(I-1) ETA(I+1)=ETA(I-1) WEIGHT(I+1)=WEIGHT(I-1) 10 CONTINUE ZETA(11)=.008394777409958D0 ETA(11)=.263112829634638D0 WEIGHT(11)=.027230314174435D0 T(1)=.728492392955404D0 ZETA(12)=ZETA(11) ETA(12)=T(1) ZETA(13)=ETA(11) ETA(13)=ZETA(11) ZETA(14)=ETA(11) ETA(14)=T(1) ZETA(15)=T(1) ETA(15)=ZETA(11) ZETA(16)=T(1) ETA(16)=ETA(11) DO 20 I=2,6 WEIGHT(10+I)=WEIGHT(11) 20 CONTINUE ELSEIF (NUMGAU .EQ. 25) THEN ZETA(1)=1.D0/3.D0 ETA(1)=1.D0/3.D0 WEIGHT(1)=.090817990382754D0 ZETA(2)=.028844733232685D0 ETA(2)=.485577633383657D0 WEIGHT(2)=.036725957756467D0 ZETA(5)=.781036849029926D0 ETA(5)=.109481575485037D0 WEIGHT(5)=.045321059435528D0 DO 30 I=3,6,3 ZETA(I)=ETA(I-1) ETA(I)=ZETA(I-1) WEIGHT(I)=WEIGHT(I-1) ZETA(I+1)=ETA(I-1) ETA(I+1)=ETA(I-1) WEIGHT(I+1)=WEIGHT(I-1) 30 CONTINUE ZETA(8)=.141707219414880D0 ETA(8)=.307939838764121D0 WEIGHT(8)=.072757916845420D0 T(1)=.550352941820999D0 ZETA(14)=.025003534762686D0 ETA(14)=.246672560639903D0 WEIGHT(14)=.028327242531057D0 T(2)=.728323904597411D0 ZETA(20)=.009540815400299D0 ETA(20)=.066803251012200D0 WEIGHT(20)=.009421666963733D0 T(3)=.923655933587500D0 IT=1 DO 40 I=9,21,6 ZETA(I)=ZETA(I-1) ETA(I)=T(IT) ZETA(I+1)=ETA(I-1) ETA(I+1)=ZETA(I-1) ZETA(I+2)=ETA(I-1) ETA(I+2)=T(IT) ZETA(I+3)=T(IT) ETA(I+3)=ZETA(I-1) ZETA(I+4)=T(IT) ETA(I+4)=ETA(I-1) DO 50 J=0,4 WEIGHT(I+J)=WEIGHT(I-1) 50 CONTINUE IT=IT+1 40 CONTINUE ELSEIF (NUMGAU .EQ. 33) THEN ZETA(1)=.023565220452390D0 ETA(1)=.488217389773805D0 WEIGHT(1)=.025731066440455D0 ZETA(4)=.120551215411079D0 ETA(4)=.439724392294460D0 WEIGHT(4)=.043692544538038D0 ZETA(7)=.457579229975768D0 ETA(7)=.271210385012116D0 WEIGHT(7)=.062858224217885D0 ZETA(10)=.744847708916828D0 ETA(10)=.127576145541586D0 WEIGHT(10)=.034796112930709D0 ZETA(13)=.957365299093579D0 ETA(13)=.021317350453210D0 WEIGHT(13)=.006166261051559D0 DO 60 I=2,14,3 ZETA(I)=ETA(I-1) ETA(I)=ZETA(I-1) WEIGHT(I)=WEIGHT(I-1) ZETA(I+1)=ETA(I-1) ETA(I+1)=ETA(I-1) WEIGHT(I+1)=WEIGHT(I-1) 60 CONTINUE ZETA(16)=.225343494534698D0 ETA(16)=.275713269685514D0 WEIGHT(16)=.040371557766381D0 T(1)=.608943235779788D0 ZETA(22)=.022838332222257D0 ETA(22)=.281325580989940D0 WEIGHT(22)=.022356773202303D0 T(2)=.695836086787803D0 ZETA(28)=.025734050548330D0 ETA(28)=.116251915907597D0 WEIGHT(28)=.017316231108659D0 T(3)=.858014033544073D0 IT=1 DO 70 I=17,29,6 ZETA(I)=ZETA(I-1) ETA(I)=T(IT) ZETA(I+1)=ETA(I-1) ETA(I+1)=ZETA(I-1) ZETA(I+2)=ETA(I-1) ETA(I+2)=T(IT) ZETA(I+3)=T(IT) ETA(I+3)=ZETA(I-1) ZETA(I+4)=T(IT) ETA(I+4)=ETA(I-1) DO 80 J=0,4 WEIGHT(I+J)=WEIGHT(I-1) 80 CONTINUE IT=IT+1 70 CONTINUE ELSEIF (NUMGAU .EQ. 42) THEN ZETA(1)=.022072179275643D0 ETA(1)=.488963910362179D0 WEIGHT(1)=.021883581368429D0 ZETA(4)=.164710561319092D0 ETA(4)=.417644719340454D0 WEIGHT(4)=.032788353544125D0 ZETA(7)=.453044943382323D0 ETA(7)=.273477528308839D0 WEIGHT(7)=.051774104507292D0 ZETA(10)=.645588935174913D0 ETA(10)=.177205532412543D0 WEIGHT(10)=.042162588736993D0 ZETA(13)=.876400233818255D0 ETA(13)=.061799883090873D0 WEIGHT(13)=.014433699669777D0 ZETA(16)=.961218077502598D0 ETA(16)=.019390961248701D0 WEIGHT(16)=.004923403602400D0 DO 90 I=2,17,3 ZETA(I)=ETA(I-1) ETA(I)=ZETA(I-1) WEIGHT(I)=WEIGHT(I-1) ZETA(I+1)=ETA(I-1) ETA(I+1)=ETA(I-1) WEIGHT(I+1)=WEIGHT(I-1) 90 CONTINUE ZETA(19)=.057124757403648D0 ETA(19)=.172266687821356D0 WEIGHT(19)=.024665753212564D0 T(1)=.770608554774996D0 ZETA(25)=.092916249356972D0 ETA(25)=.336861459796345D0 WEIGHT(25)=.038571510787061D0 T(2)=.570222290846683D0 ZETA(31)=.014646950055654D0 ETA(31)=.298372882136258D0 WEIGHT(31)=.014436308113534D0 T(3)=.686980167808088D0 ZETA(37)=.001268330932872D0 ETA(37)=.118974497696957D0 WEIGHT(37)=.005010228838501D0 T(4)=.879757171370171D0 IT=1 DO 100 I=20,38,6 ZETA(I)=ZETA(I-1) ETA(I)=T(IT) ZETA(I+1)=ETA(I-1) ETA(I+1)=ZETA(I-1) ZETA(I+2)=ETA(I-1) ETA(I+2)=T(IT) ZETA(I+3)=T(IT) ETA(I+3)=ZETA(I-1) ZETA(I+4)=T(IT) ETA(I+4)=ETA(I-1) DO 110 J=0,4 WEIGHT(I+J)=WEIGHT(I-1) 110 CONTINUE IT=IT+1 100 CONTINUE ELSEIF (NUMGAU .EQ. 52) THEN ZETA(1)=1.D0/3.D0 ETA(1)=1.D0/3.D0 WEIGHT(1)=.046875697427642D0 ZETA(2)=.005238916103123D0 ETA(2)=.497380541948438D0 WEIGHT(2)=.006405878578585D0 ZETA(5)=.173061122901295D0 ETA(5)=.413469438549352D0 WEIGHT(5)=.041710296739387D0 ZETA(8)=.059082801866017D0 ETA(8)=.470458599066991D0 WEIGHT(8)=.026891484250064D0 ZETA(11)=.518892500060958D0 ETA(11)=.240553749969521D0 WEIGHT(11)=.042132522761650D0 ZETA(14)=.704068411554854D0 ETA(14)=.147965794222573D0 WEIGHT(14)=.030000266842773D0 ZETA(17)=.849069624685052D0 ETA(17)=.075465187657474D0 WEIGHT(17)=.014200098925024D0 ZETA(20)=.966807194753950D0 ETA(20)=.016596402623025D0 WEIGHT(20)=.003582462351273D0 DO 120 I=3,21,3 ZETA(I)=ETA(I-1) ETA(I)=ZETA(I-1) WEIGHT(I)=WEIGHT(I-1) ZETA(I+1)=ETA(I-1) ETA(I+1)=ETA(I-1) WEIGHT(I+1)=WEIGHT(I-1) 120 CONTINUE ZETA(23)=.103575692245252D0 ETA(23)=.296555596579887D0 WEIGHT(23)=.032773147460627D0 T(1)=.599868711174861D0 ZETA(29)=.020083411655416D0 ETA(29)=.337723063403079D0 WEIGHT(29)=.015298306248441D0 T(2)=.642183524941505D0 ZETA(35)=-.004341002614139D0 ETA(35)=.204748281642812D0 WEIGHT(35)=.002386244192839D0 T(3)=.799592720971327D0 ZETA(41)=.041941786468010D0 ETA(41)=.189358492130623D0 WEIGHT(41)=.019084792755899D0 T(4)=.768699721401368D0 ZETA(47)=.014317320230681D0 ETA(47)=.085283615682657D0 WEIGHT(47)=.006850054546542D0 T(5)=.900399064086661D0 IT=1 DO 130 I=24,48,6 ZETA(I)=ZETA(I-1) ETA(I)=T(IT) ZETA(I+1)=ETA(I-1) ETA(I+1)=ZETA(I-1) ZETA(I+2)=ETA(I-1) ETA(I+2)=T(IT) ZETA(I+3)=T(IT) ETA(I+3)=ZETA(I-1) ZETA(I+4)=T(IT) ETA(I+4)=ETA(I-1) DO 140 J=0,4 WEIGHT(I+J)=WEIGHT(I-1) 140 CONTINUE IT=IT+1 130 CONTINUE ENDIF RETURN END zt (zt@cfdrc.com) wrote: : I am looking for a high order Gaussian integration formula on a triangle. : The formula should be exact for all x^i y^j (0<=i+j<=n), and n is at least : 8. I am interested in a real 2D formula with symmetric abscissas and : weights, not tensor product of 1D Gaussian integrals. Thanks. : Z. Tan : zt@cfdrc.com -- ---------------------------------------------------------------------- James P. Smith NASA/JSC, Mail Code ES2 Houston, Texas 77058 james.p.smith1@jsc.nasa.gov smithj@smith-james-p.jsc.nasa.gov