SUBROUTINE TWODQ(F,N,X,Y,TOL,ICLOSE,MAXTRI,MEVALS,RESULT, * ERROR,NU,ND,NEVALS,IFLAG,DATA,IWORK) C***BEGIN PROLOGUE TWODQ C***DATE WRITTEN 840518 (YYMMDD) C***REVISION DATE 840518 (YYMMDD) C***CATEGORY NO. D1I C***KEYWORDS QUADRATURE,TWO DIMENSIONAL,ADAPTIVE,CUBATURE C***AUTHOR KAHANER,D.K.,N.B.S. C RECHARD,O.W.,N.B.S. C BARNHILL,ROBERT,UNIV. OF UTAH C***PURPOSE To compute the two-dimensional integral of a function C F over a region consisting of N triangles. C***DESCRIPTION C C This subroutine computes the two-dimensional integral of a C function F over a region consisting of N triangles. C A total error estimate is obtained and compared with a C tolerance - TOL - that is provided as input to the subroutine. C The error tolerance is treated as either relative or absolute C depending on the input value of IFLAG. A 'Local Quadrature C Module' is applied to each input triangle and estimates of the C total integral and the total error are computed. The local C quadrature module is either subroutine LQM0 or subroutine C LQM1 and the choice between them is determined by the C value of the input variable ICLOSE. C C If the total error estimate exceeds the tolerance, the triangle C with the largest absolute error is divided into two triangles C by a median to its longest side. The local quadrature module C is then applied to each of the subtriangles to obtain new C estimates of the integral and the error. This process is C repeated until either (1) the error tolerance is satisfied, C (2) the number of triangles generated exceeds the input C parameter MAXTRI, (3) the number of integrand evaluations C exceeds the input parameter MEVALS, or (4) the subroutine C senses that roundoff error is beginning to contaminate C the result. C C The user must specify MAXTRI, the maximum number of triangles C in the final triangulation of the region, and provide two C storage arrays - DATA and IWORK - whose sizes are at least C 9*MAXTRI and 2*MAXTRI respectively. The user must also C specify MEVALS, the maximum number of function evaluations C to be allowed. This number will be effective in limiting C the computation only if it is less than 94*MAXTRI when LQM1 C is specified or 56*MAXTRI when LQM0 is specified. C C After the subroutine has returned to the calling program C with output values, it can be called again with a smaller C value of TOL and/or a different value of MEVALS. The tolerance C can also be changed from relative to absolute C or vice-versa by changing IFLAG. Unless C the parameters NU and ND are reset to zero the subroutine C will restart with the final set of triangles and output C values from the previous call. C C ARGUMENTS: C C F function subprogram defining the integrand F(u,v); C the actual name for F needs to be declared EXTERNAL C by the calling program. C C N the number of input triangles. C C X a 3 by N array containing the abscissae of the vertices C of the N triangles. C C Y a 3 by N array containing the ordinates of the vertices C of the N triangles C C TOL the desired bound on the error. If IFLAG=0 on input, C TOL is interpreted as a bound on the relative error; C if IFLAG=1, the bound is on the absolute error. C C ICLOSE an integer parameter that determines the selection C of LQM0 or LQM1. If ICLOSE=1 then LQM1 is used. C Any other value of ICLOSE causes LQM0 to be used. C LQM0 uses function values only at interior points of C the triangle. LQM1 is usually more accurate than LQM0 C but involves evaluating the integrand at more points C including some on the boundary of the triangle. It C will usually be better to use LQM1 unless the integrand C has singularities on the boundary of the triangle. C C MAXTRI The maximum number of triangles that are allowed C to be generated by the computation. C C MEVALS The maximum number of function evaluations allowed. C C RESULT output of the estimate of the integral. C C ERROR output of the estimate of the absolute value of the C total error. C C NU an integer variable used for both input and output. Must C be set to 0 on first call of the subroutine. Subsequent C calls to restart the subroutine should use the previous C output value. C C ND an integer variable used for both input and output. Must C be set to 0 on first call of the subroutine. Subsequent C calls to restart the subroutine should use the previous C output value. C C NEVALS The actual number of function evaluations performed. C C IFLAG on input: C IFLAG=0 means TOL is a bound on the relative error; C IFLAG=1 means TOL is a bound on the absolute error; C any other input value for IFLAG causes the subroutine C to return immediately with IFLAG set equal to 9. C C on output: C IFLAG=0 means normal termination; C IFLAG=1 means termination for lack of space to divide C another triangle; C IFLAG=2 means termination because of roundoff noise C IFLAG=3 means termination with relative error <= C 5.0* machine epsilon; C IFLAG=4 means termination because the number of function C evaluations has exceeded MEVALS. C IFLAG=9 means termination because of error in input flag C C DATA a one dimensional array of length >= 9*MAXTRI C passed to the subroutine by the calling program. It is C used by the subroutine to store information C about triangles used in the quadrature. C C IWORK a one dimensional integer array of length >= 2*MAXTRI C passed to the subroutine by the calling program. C It is used by the subroutine to store pointers C to the information in the DATA array. C C C The information for each triangle is contained in a nine word C record consisting of the error estimate, the estimate of the C integral, the coordinates of the three vertices, and the area. C These records are stored in the DATA array C that is passed to the subroutine. The storage is organized C into two heaps of length NU and ND respectively. The first heap C contains those triangles for which the error exceeds C epsabs*a/ATOT where epsabs is a bound on the absolute error C derived from the input tolerance (which may refer to relative C or absolute error), a is the area of the triangle, and ATOT C is the total area of all triangles. The second heap contains C those triangles for which the error is less than or equal to C epsabs*a/ATOT. At the top of each heap is the triangle with C the largest absolute error. C C Pointers into the heaps are contained in the array IWORK. C Pointers to the first heap are contained C between IWORK(1) and IWORK(NU). Pointers to the second C heap are contained between IWORK(MAXTRI+1) and C IWORK(MAXTRI+ND). The user thus has access to the records C stored in the DATA array through the pointers in IWORK. C For example, the following two DO loops will print out C the records for each triangle in the two heaps: C C DO 10 I=1,NU C PRINT*,(DATA(IWORK(I)+J),J=0,8) C 10 CONTINUE C DO 20 I=1,ND C PRINT*,(DATA(IWORK(MAXTRI+I)+J),J=0,8 C 20 CONTINUE C C When the total number of triangles is equal to C MAXTRI, the program attempts to remove a triangle from the C bottom of the second heap and continue. If the second heap C is empty, the program returns with the current estimates of C the integral and the error and with IFLAG set equal to 1. C Note that in this case the actual number of triangles C processed may exceed MAXTRI and the triangles stored in C the DATA array may not constitute a complete triangulation C of the region. C C The following sample program will calculate the integral of C cos(x+y) over the square (0.,0.),(1.,0.),(1.,1.),(0.,1.) and C print out the values of the estimated integral, the estimated C error, the number of function evaluations, and IFLAG. C C REAL X(3,2),Y(3,2),DATA(450),RES,ERR C INTEGER IWORK(100),NU,ND,NEVALS,IFLAG C EXTERNAL F C X(1,1)=0. C Y(1,1)=0. C X(2,1)=1. C Y(2,1)=0. C X(3,1)=1. C Y(3,1)=1. C X(1,2)=0. C Y(1,2)=0. C X(2,2)=1. C Y(2,2)=1. C X(3,2)=0. C Y(3,2)=1. C NU=0 C ND=0 C IFLAG=1 C CALL TWODQ(F,2,X,Y,1.E-04,1,50,4000,RES,ERR,NU,ND, C * NEVALS,IFLAG,DATA,IWORK) C PRINT*,RES,ERR,NEVALS,IFLAG C END C FUNCTION F(X,Y) C F=COS(X+Y) C RETURN C END C C***REFERENCES (NONE) C C***ROUTINES CALLED HINITD,HINITU,HPACC,HPDEL,HPINS,LQM0,LQM1, C TRIDV,R1MACH C***END PROLOGUE TWODQ