program main !*****************************************************************************80 ! !! CUBE_ARBQ_RULE_TEST() tests CUBE_ARBQ_RULE(). ! ! Licensing: ! ! This code is distributed under the GNU GPL license. ! ! Modified: ! ! 07 July 2014 ! ! Author: ! ! Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas. ! This FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Hong Xiao, Zydrunas Gimbutas, ! A numerical algorithm for the construction of efficient quadrature ! rules in two and higher dimensions, ! Computers and Mathematics with Applications, ! Volume 59, 2010, pages 663-676. ! implicit none integer cube_arbq_size integer degree character * ( 255 ) header integer n call timestamp ( ) write ( *, '(a)' ) '' write ( *, '(a)' ) 'CUBE_ARBQ_RULE_TEST' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' Test CUBE_ARBQ_RULE().' degree = 8 n = cube_arbq_size ( degree ) header = 'cube08' call test01 ( degree, n ) call test02 ( degree, n, header ) call test03 ( degree, n, header ) call test04 ( degree, n ) ! ! Terminate. ! write ( *, '(a)' ) '' write ( *, '(a)' ) 'CUBE_ARBQ_RULE_TEST' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) '' call timestamp ( ) stop 0 end subroutine test01 ( degree, n ) !*****************************************************************************80 ! !! TEST01 calls CUBE_ARBQ for a quadrature rule of given order. ! ! Licensing: ! ! This code is distributed under the GNU GPL license. ! ! Modified: ! ! 07 July 2014 ! ! Author: ! ! Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas. ! This FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Hong Xiao, Zydrunas Gimbutas, ! A numerical algorithm for the construction of efficient quadrature ! rules in two and higher dimensions, ! Computers and Mathematics with Applications, ! Volume 59, 2010, pages 663-676. ! ! Parameters: ! ! Input, integer DEGREE, the desired total polynomial degree exactness ! of the quadrature rule. ! ! Input, integer N, the number of nodes. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) d integer degree integer j real ( kind = rk ) volume real ( kind = rk ) w(n) real ( kind = rk ) x(3,n) write ( *, '(a)' ) '' write ( *, '(a)' ) 'TEST01' write ( *, '(a)' ) ' Quadrature rule for the symmetric cube.' write ( *, '(a,i4)' ) & ' Polynomial exactness degree DEGREE = ', degree volume = 8.0D+00 ! ! Retrieve and print a quadrature rule. ! call cube_arbq ( degree, n, x, w ) write ( *, '(a)' ) '' write ( *, '(a,i6)' ) ' Number of nodes N = ', n write ( *, '(a)' ) '' write ( *, '(a)' ) & ' J W X Y Z' write ( *, '(a)' ) '' do j = 1, n write ( *, '(2x,i4,2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6)' ) & j, w(j), x(1,j), x(2,j), x(3,j) end do d = sum ( w(1:n) ) write ( *, '(a,2x,g14.6)' ) ' Sum ', d write ( *, '(a,2x,g14.6)' ) ' Volume', volume return end subroutine test02 ( degree, n, header ) !*****************************************************************************80 ! !! TEST02 gets a rule and writes it to a file. ! ! Licensing: ! ! This code is distributed under the GNU GPL license. ! ! Modified: ! ! 07 July 2014 ! ! Author: ! ! Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas. ! This FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Hong Xiao, Zydrunas Gimbutas, ! A numerical algorithm for the construction of efficient quadrature ! rules in two and higher dimensions, ! Computers and Mathematics with Applications, ! Volume 59, 2010, pages 663-676. ! ! Parameters: ! ! Input, integer DEGREE, the desired total polynomial degree exactness ! of the quadrature rule. 0 <= DEGREE <= 15. ! ! Input, integer N, the number of nodes to be used by the rule. ! ! Input, character * ( * ) HEADER, an identifier for the filenames. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer degree character * ( * ) header integer i integer rule_unit character * ( 255 ) rule_filename real ( kind = rk ) w(n) real ( kind = rk ) x(3,n) write ( *, '(a)' ) '' write ( *, '(a)' ) 'TEST02' write ( *, '(a)' ) & ' Get a quadrature rule for the symmetric cube.' write ( *, '(a)' ) ' Then write it to a file.' write ( *, '(a,i4)' ) & ' Polynomial exactness degree DEGREE = ', degree ! ! Retrieve a quadrature rule. ! call cube_arbq ( degree, n, x, w ) ! ! Write the points and weights to a file. ! call get_unit ( rule_unit ) rule_filename = trim ( header ) // '.txt' open ( unit = rule_unit, file = rule_filename, & status = 'replace' ) do i = 1, n write ( rule_unit, '(4(e21.15,2x))' ) & x(1,i), x(2,i), x(3,i), w(i) end do close ( unit = rule_unit ) write ( *, '(a)' ) '' write ( *, '(a)' ) ' Quadrature rule written to file "' & // trim ( rule_filename ) // '".' return end subroutine test03 ( degree, n, header ) !*****************************************************************************80 ! !! TEST03 gets a rule and creates GNUPLOT input files. ! ! Licensing: ! ! This code is distributed under the GNU GPL license. ! ! Modified: ! ! 07 July 2014 ! ! Author: ! ! Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas. ! This FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Hong Xiao, Zydrunas Gimbutas, ! A numerical algorithm for the construction of efficient quadrature ! rules in two and higher dimensions, ! Computers and Mathematics with Applications, ! Volume 59, 2010, pages 663-676. ! ! Parameters: ! ! Input, integer DEGREE, the desired total polynomial degree exactness ! of the quadrature rule. 0 <= DEGREE <= 15. ! ! Input, integer N, the number of nodes to be used by the rule. ! ! Input, character * ( * ) HEADER, an identifier for the filenames. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer degree character * ( * ) header real ( kind = rk ) w(n) real ( kind = rk ) x(3,n) write ( *, '(a)' ) '' write ( *, '(a)' ) 'TEST03' write ( *, '(a)' ) & ' Get a quadrature rule for the symmetric cube.' write ( *, '(a)' ) ' Set up GNUPLOT graphics input.' write ( *, '(a,i4)' ) & ' Polynomial exactness degree DEGREE = ', degree ! ! Retrieve a quadrature rule. ! call cube_arbq ( degree, n, x, w ) ! ! Create files for input to GNUPLOT. ! call cube_arbq_gnuplot ( n, x, header ) return end subroutine test04 ( degree, n ) !*****************************************************************************80 ! !! TEST04 gets a rule and tests its accuracy. ! ! Licensing: ! ! This code is distributed under the GNU GPL license. ! ! Modified: ! ! 08 July 2014 ! ! Author: ! ! Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas. ! This FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Hong Xiao, Zydrunas Gimbutas, ! A numerical algorithm for the construction of efficient quadrature ! rules in two and higher dimensions, ! Computers and Mathematics with Applications, ! Volume 59, 2010, pages 663-676. ! ! Parameters: ! ! Input, integer DEGREE, the desired total polynomial degree exactness ! of the quadrature rule. 0 <= DEGREE <= 15. ! ! Input, integer N, the number of nodes to be used by the rule. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) d integer degree integer i integer j integer npols real ( kind = rk ) pols(((degree+1)*(degree+2)*(degree+3))/6) real ( kind = rk ) rints(((degree+1)*(degree+2)*(degree+3))/6) real ( kind = rk ) volume real ( kind = rk ) w(n) real ( kind = rk ) x(3,n) real ( kind = rk ) z(3) write ( *, '(a)' ) '' write ( *, '(a)' ) 'TEST04' write ( *, '(a)' ) & ' Get a quadrature rule for the symmetric cube.' write ( *, '(a)' ) ' Test its accuracy.' write ( *, '(a,i4)' ) & ' Polynomial exactness degree DEGREE = ', degree ! ! Retrieve a quadrature rule. ! call cube_arbq ( degree, n, x, w ) npols = ( ( degree + 1 ) * ( degree + 2 ) * ( degree + 3 ) ) / 6 rints(1:npols) = 0.0D+00 do i = 1, n z(1) = x(1,i) z(2) = x(2,i) z(3) = x(3,i) call lege3eva ( degree, z, pols ) do j = 1, npols rints(j) = rints(j) + w(i) * pols(j) end do end do volume = 8.0D+00 d = 0.0D+00 d = ( rints(1) - sqrt ( volume ) )**2 do i = 2, npols d = d + rints(i)**2 end do d = sqrt ( d ) / real ( npols, kind = rk ) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' RMS error = ', d return end