'************************************************************************* '* Test programm for cubature over triangles via summed Gaussian n-point * '* formula * '* --------------------------------------------------------------------- * '* SAMPLE RUN: * '* (Integrate function EXP(-(x*x + y*y)) over triangle defined by three * '* points (0,0), (10,0) and (0,10). ) * '* * '* # Method Number Value Error Code * '* Subtriangles Integral (must be 0) * '* ------------------------------------------------- * '* 1 1 0.0000000112 0 * '* 1 4 0.0483240045 0 * '* 1 9 0.4706075316 0 * '* 1 16 0.7913533084 0 * '* 1 25 0.8814598010 0 * '* 2 1 0.0000000000 0 * '* 2 4 0.0120496973 0 * '* 2 9 0.2538493244 0 * '* 2 16 0.5593099250 0 * '* 2 25 0.7075499085 0 * '* 3 1 0.0644320023 0 * '* 3 4 1.0390297430 0 * '* 3 9 1.0188794444 0 * '* 3 16 0.8590974900 0 * '* 3 25 0.7974378946 0 * '* 7 1 0.8091876107 0 * '* 7 4 0.9654995491 0 * '* 7 9 0.7922397459 0 * '* 7 16 0.7849640612 0 * '* 7 25 0.7862848602 0 * '* ------------------------------------------------- * '* * '* --------------------------------------------------------------------- * '* Reference: * '* "Numerical Algorithms with C, By Gisela Engeln-Muellges * '* and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * '* * '* Basic Release 1.0 By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '************************************************************************* 'Program TK3GAN DEFDBL A-H, O-Z DEFINT I-N OPTION BASE 0 DIM Verfahren AS INTEGER 'define summits of triangle Px = 0#: Py = 0# Qx = 10#: Qy = 0# Rx = 0#: Ry = 10# F$ = " #######.##########" Init = 1 'Initialize constants on first time 'print header CLS PRINT " -------------------------------------------------" PRINT " # Method Number Value Error Code " PRINT " Subtriangles Integral (must be 0) " PRINT " -------------------------------------------------" 'main integration loop FOR Verfahren = 1 TO 7 'only methods 1,2,3,7 are implemented IF Verfahren < 4 OR Verfahren > 6 THEN FOR KantenTeilung = 1 TO 5 'call appropriate Gaussian cubature method GOSUB 1000 'call Kub3GauN (Px, Py, Qx, Qy, Rx, Ry, Verfahren,Kantenteilung, ' Wert,icode) 'print one line of results PRINT " "; Verfahren; PRINT USING " ##"; KantenTeilung ^ 2; PRINT USING F$; Wert; PRINT " "; icode NEXT KantenTeilung END IF NEXT Verfahren PRINT " -------------------------------------------------" END 'of main program 'user defined function(x,y) 500 func = EXP(-(X * X + Y * Y)) RETURN 1000 'Subroutine Kub3GauN (Px, Py, Qx, Qy, Rx, Ry, n, m, Wert, icode) '************************************************************************ '* Cubature over triangular regions using the n-point Gauss formula * '* * '* Integrate the function f (x,y) over the triangle PQR using the summed* '* n-point Gauss cubature formula on m x m subtriangles. * '* The subtriangles have edges of length 1/m of the original edge * '* lengths. * '* * '* Input parameters: * '* Double Px,Py coordinates of P * '* Double Qx,Ry coordinates of Q * '* Double Rx,Ry coordinates of R * '* integer n order of method (= number of points in each sub- * '* triangle) * '* integer m number of subtriangles along one edge * '* * '* Output parameter: * '* Double Wert value of integral * '* * '* Return value : * '* 0: o.k. * '* 1: m improper * '* 2: the corners P, Q and R are collinear * '* 3: nth order method not implemented * '* * '* Author: Uli Eggermann, 8.1.1990 (C version) * '************************************************************************ DIM Gau1Konst(3) DIM Gau2Konst(2, 3) DIM Gau3Konst(3, 3) DIM Gau7Konst(7, 3) DIM GauNmax AS INTEGER DIM d AS INTEGER ' GauNEpsilon serves as a check of collinearity; if the area of P, Q ' and R has area less than Epsilon/2, we judge the three points to be ' collinear. GauNEpsilon = .0001# ' Initialize Gauss constants on first time IF Init <> 0 THEN Gau1Konst(0) = .5#: Gau1Konst(1) = 1# / 3#: Gau1Konst(2) = 1# / 3# Gau2Konst(0, 0) = .25#: Gau2Konst(0, 1) = 1# / 6#: Gau2Konst(0, 2) = .5# Gau2Konst(1, 0) = .25#: Gau2Konst(1, 1) = .5#: Gau2Konst(1, 2) = 1# / 6# Gau3Konst(0, 0) = 1# / 6#: Gau3Konst(0, 1) = 1# / 6#: Gau3Konst(0, 2) = 1# / 6# Gau3Konst(1, 0) = 1# / 6#: Gau3Konst(1, 1) = 2# / 3#: Gau3Konst(1, 2) = 1# / 6# Gau3Konst(2, 0) = 1# / 6#: Gau3Konst(2, 1) = 1# / 6#: Gau3Konst(2, 2) = 2# / 3# Gau7Konst(0, 0) = .1125#: Gau7Konst(0, 1) = 1# / 3#: Gau7Konst(0, 2) = 1# / 3# Gau7Konst(1, 0) = .0661970763942531#: Gau7Konst(1, 1) = .4701420641051151# Gau7Konst(1, 2) = .4701420641051151# Gau7Konst(2, 0) = .0661970763942531#: Gau7Konst(2, 1) = 5.971587178976981D-02 Gau7Konst(2, 2) = .4701420641051151# Gau7Konst(3, 0) = .0661970763942531#: Gau7Konst(3, 1) = .4701420641051151# Gau7Konst(3, 2) = 5.971587178976981D-02 Gau7Konst(4, 0) = 6.296959027241357D-02: Gau7Konst(4, 1) = .1012865073234563# Gau7Konst(4, 2) = .1012865073234563# Gau7Konst(5, 0) = 6.296959027241357D-02: Gau7Konst(5, 1) = .7974269853530873# Gau7Konst(5, 2) = .1012865073234563# Gau7Konst(6, 0) = 6.296959027241357D-02: Gau7Konst(6, 1) = .1012865073234563# Gau7Konst(6, 2) = .7974269853530873# Init = 0 END IF GauNmax = 13: m = KantenTeilung: n = Verfahren IF m < 1 THEN ' m o.k. ? icode = 1 RETURN END IF 'Test collinearity Area = Px * (Qy - Ry) + Qx * (Ry - Py) + Rx * (Py - Qy) IF ABS(Area) < GauNEpsilon THEN icode = 2 RETURN END IF IF n < 1 AND n > 3 AND n <> 7 THEN ' desired method implemented ? icode = 3 RETURN END IF Wert = 0# ' initialize Area = Area / (m * m) ' double triangle area hPQx = (Qx - Px) / m hPRx = (Rx - Px) / m ' edge vectors for the hPQy = (Qy - Py) / m ' m * m hPRy = (Ry - Py) / m ' subtriangles FOR d = 0 TO 1 ' types of triangles d IF d <> 0 THEN Fak = -1# ELSE Fak = 1# ' d = 1: reflected END IF FOR j = d TO m - 1 ' j: along PR FOR i = d TO m - j + d - 1 ' i: along PQ Dx = Px + i * hPQx + j * hPRx ' (Dx,Dy) ist top Dy = Py + i * hPQy + j * hPRy ' corner of subtriangle FOR k = 0 TO n - 1 ' Sum of weighted IF n = 1 THEN Gw = Gau1Konst(0) Gx = Gau1Konst(1) Gy = Gau1Konst(2) ELSEIF n = 2 THEN Gw = Gau2Konst(k, 0) Gx = Gau2Konst(k, 1) Gy = Gau2Konst(k, 2) ELSEIF n = 3 THEN Gw = Gau3Konst(k, 0) Gx = Gau3Konst(k, 1) Gy = Gau3Konst(k, 2) ELSEIF n = 7 THEN Gw = Gau7Konst(k, 0) Gx = Gau7Konst(k, 1) Gy = Gau7Konst(k, 2) END IF X = Dx + Fak * (Gx * hPQx + Gy * hPRx) Y = Dy + Fak * (Gx * hPQy + Gy * hPRy) GOSUB 500 'func(X,Y) Wert = Wert + Gw * func NEXT k NEXT i NEXT j NEXT d Wert = Wert * Area 'multiply by double the subarea icode = 0 RETURN ' end of file tk3gan.bas