subroutine comp_next ( n, k, a, more, h, t ) !*****************************************************************************80 ! !! comp_next() computes the compositions of the integer N into K parts. ! ! Discussion: ! ! A composition of the integer N into K parts is an ordered sequence ! of K nonnegative integers which sum to N. The compositions (1,2,1) ! and (1,1,2) are considered to be distinct. ! ! The routine computes one composition on each call until there are no more. ! For instance, one composition of 6 into 3 parts is ! 3+2+1, another would be 6+0+0. ! ! On the first call to this routine, set MORE = FALSE. The routine ! will compute the first element in the sequence of compositions, and ! return it, as well as setting MORE = TRUE. If more compositions ! are desired, call again, and again. Each time, the routine will ! return with a new composition. ! ! However, when the LAST composition in the sequence is computed ! and returned, the routine will reset MORE to FALSE, signaling that ! the end of the sequence has been reached. ! ! There are 28 compositions of 6 into three parts. This routine will ! produce those compositions in the following order: ! ! I A ! - --------- ! 1 6 0 0 ! 2 5 1 0 ! 3 4 2 0 ! 4 3 3 0 ! 5 2 4 0 ! 6 1 5 0 ! 7 0 6 0 ! 8 5 0 1 ! 9 4 1 1 ! 10 3 2 1 ! 11 2 3 1 ! 12 1 4 1 ! 13 0 5 1 ! 14 4 0 2 ! 15 3 1 2 ! 16 2 2 2 ! 17 1 3 2 ! 18 0 4 2 ! 19 3 0 3 ! 20 2 1 3 ! 21 1 2 3 ! 22 0 3 3 ! 23 2 0 4 ! 24 1 1 4 ! 25 0 2 4 ! 26 1 0 5 ! 27 0 1 5 ! 28 0 0 6 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 July 2008 ! ! Author: ! ! Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf. ! This version by John Burkardt. ! ! Reference: ! ! Albert Nijenhuis, Herbert Wilf, ! Combinatorial Algorithms for Computers and Calculators, ! Second Edition, ! Academic Press, 1978, ! ISBN: 0-12-519260-6, ! LC: QA164.N54. ! ! Parameters: ! ! Input, integer N, the integer whose compositions are desired. ! ! Input, integer K, the number of parts in the composition. ! ! Input/output, integer A(K), the parts of the composition. ! ! Input/output, logical MORE, set by the user to start the computation, ! and by the routine to terminate it. ! ! Input/output, integer H, T, two internal parameters needed ! for the computation. The user should allocate space for these in the ! calling program, include them in the calling sequence, but never alter ! them! ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer k integer a(k) integer h logical more integer n integer t ! ! The first computation. ! if ( .not. more ) then t = n h = 0 a(1) = n a(2:k) = 0 ! ! The next computation. ! else ! ! If the first entry A(1) is positive, then set H to zero, ! so that when we increment H, it points to A(1); we will decrement A(1) by 1 ! and increment A(2). ! if ( 1 < t ) then h = 0 end if ! ! Otherwise, A(1) is 0. Then by H + 1 is the entry we incremented last time. ! Set H = H + 1, zero A(H), adding all but one of its value to A(1), ! and incrementing A(H+1) by 1. ! h = h + 1 t = a(h) a(h) = 0 a(1) = t - 1 a(h+1) = a(h+1) + 1 end if ! ! This is the last element of the sequence if all the ! items are in the last slot. ! more = ( a(k) /= n ) return end subroutine monomial_value ( d, n, e, x, v ) !*****************************************************************************80 ! !! monomial_value() evaluates a monomial. ! ! Discussion: ! ! This routine evaluates a monomial of the form ! ! product ( 1 <= i <= d ) x(i)^e(i) ! ! The combination 0.0^0, if encountered, is treated as 1.0. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 May 2023 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer D, the spatial dimension. ! ! integer N, the number of evaluation points. ! ! integer E(D), the exponents. ! ! real ( kind = rk ) X(N,D), the point coordinates. ! ! Output: ! ! real ( kind = rk ) V(N), the monomial values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer d integer n integer e(d) integer j real ( kind = rk ) v(n) real ( kind = rk ) x(n,d) v(1:n) = 1.0D+00 do j = 1, d if ( 0 /= e(j) ) then v(1:n) = v(1:n) * x(1:n,j) ** e(j) end if end do return end function quadrilateral_unit_area ( ) !*****************************************************************************80 ! !! quadrilateral_unit_area() returns the area of a unit quadrilateral. ! ! Discussion: ! ! The unit quadrilateral has vertices (0,0), (1,0), (1,1), (0,1). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Output: ! ! real quadrilateral_unit_area: the area. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) quadrilateral_unit_area quadrilateral_unit_area = 1.0D+00 return end subroutine quadrilateral_unit_monomial_integral ( expon, value ) !*****************************************************************************80 ! !! quadrilateral_unit_monomial_integral(): monomial integral in a unit quadrilateral. ! ! Discussion: ! ! This routine returns the integral of ! ! product ( 1 <= I <= 2 ) X(I)^EXPON(I) ! ! over the unit quadrilateral. ! ! The unit quadrilateral Q is bounded by the vertices ! (0,0), (1,0), (1,1), (0,1)). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 May 2023 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer EXPON(2): the exponents. ! ! Output: ! ! real VALUE: the integral of the monomial. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer expon(2) real ( kind = rk ) value value = 1.0D+00 / ( expon(1) + 1 ) & * 1.0D+00 / ( expon(2) + 1 ) return end subroutine quadrilateral_witherden_rule ( p, n, x, y, w ) !*****************************************************************************80 ! !! quadrilateral_witherden_rule() returns a quadrilateral quadrature rule of given precision. ! ! Discussion: ! ! The unit quadrilateral Q is bounded by the vertices ! (0,0), (1,0), (1,1), (0,1). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! int p: the precision, 0 <= p <= 21. ! ! int n: the order of the rule. ! ! Output: ! ! double x(n), y(n): the coordinates of quadrature points. ! ! double w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer p real ( kind = rk ) w(n) real ( kind = rk ) x(n) real ( kind = rk ) y(n) if ( p < 0 ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'quadrilateral_witherden_rule(): Fatal error!' write ( *, '(a)' ) ' Input p < 0.' stop ( 1 ) end if if ( 21 < p ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'quadrilateral_witherden_rule(): Fatal error!' write ( *, '(a)' ) ' Input 21 < p.' stop ( 1 ) end if if ( p <= 1 ) then call rule01 ( n, x, y, w ) else if ( p <= 3 ) then call rule03 ( n, x, y, w ) else if ( p <= 5 ) then call rule05 ( n, x, y, w ) else if ( p <= 7 ) then call rule07 ( n, x, y, w ) else if ( p <= 9 ) then call rule09 ( n, x, y, w ) else if ( p <= 11 ) then call rule11 ( n, x, y, w ) else if ( p <= 13 ) then call rule13 ( n, x, y, w ) else if ( p <= 15 ) then call rule15 ( n, x, y, w ) else if ( p <= 17 ) then call rule17 ( n, x, y, w ) else if ( p <= 19 ) then call rule19 ( n, x, y, w ) else if ( p <= 21 ) then call rule21 ( n, x, y, w ) end if return end subroutine rule_order ( p, order ) !*****************************************************************************80 ! !! rule_order() returns the order of a quadrilateral quadrature rule of given precision. ! ! Discussion: ! ! The "order" of a quadrature rule is the number of points involved. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer p: the precision, 0 <= p <= 21. ! ! Output: ! ! integer order: the order of the rule. ! implicit none integer order integer, dimension ( 0 : 21 ) :: order_vec = (/ & 1, & 1, 4, 4, 8, 8, 12, 12, 20, 20, 28, & 28, 37, 37, 48, 48, 60, 60, 72, 72, 85, & 85 /) integer p if ( p < 0 ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'rule_order(): Fatal error!' write ( *, '(a)' ) ' Input p < 0.' stop ( 1 ) end if if ( 21 < p ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'rule_order(): Fatal error!' write ( *, '(a)' ) ' Input 21 < p.' stop ( 1 ) end if order = order_vec(p) return end subroutine rule01 ( n, x, y, w ) !*****************************************************************************80 ! !! rule01() returns the prism rule of precision 1. ! ! Discussion: ! ! We suppose we are given a quadrilateral Quad with vertices ! (0,0), (1,0), (1,1), (0,1), ! We call a rule with n points, returning coordinates ! x, y, and weights w. Then the integral I of f(x,y,z) over Quad is ! approximated by Q as follows: ! ! Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 1 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.5000000000000000D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 1.0000000000000000D+00 /) x = x_save(1:n_save) y = y_save(1:n_save) w = w_save(1:n_save) return end subroutine rule03 ( n, x, y, w ) !*****************************************************************************80 ! !! rule03() returns the quadrilateral rule of precision 3. ! ! Discussion: ! ! We suppose we are given a quadrilateral Quad with vertices ! (0,0), (1,0), (1,1), (0,1), ! We call a rule with n points, returning coordinates ! x, y, and weights w. Then the integral I of f(x,y,z) over Quad is ! approximated by Q as follows: ! ! Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 29 April 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 4 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.7886751345948129D+00, & 0.7886751345948129D+00, & 0.2113248654051871D+00, & 0.2113248654051871D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.7886751345948129D+00, & 0.2113248654051871D+00, & 0.7886751345948129D+00, & 0.2113248654051871D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.2500000000000000D+00, & 0.2500000000000000D+00, & 0.2500000000000000D+00, & 0.2500000000000000D+00 /) x = x_save(1:n_save) y = y_save(1:n_save) w = w_save(1:n_save) return end subroutine rule05 ( n, x, y, w ) !*****************************************************************************80 ! !! rule05() returns the prism rule of precision 5. ! ! Discussion: ! ! We suppose we are given a quadrilateral Quad with vertices ! (0,0), (1,0), (1,1), (0,1), ! We call a rule with n points, returning coordinates ! x, y, and weights w. Then the integral I of f(x,y,z) over Quad is ! approximated by Q as follows: ! ! Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 8 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.8415650255319866D+00, & 0.5000000000000000D+00, & 0.1584349744680134D+00, & 0.5000000000000000D+00, & 0.9409585518440984D+00, & 0.9409585518440984D+00, & 0.0590414481559016D+00, & 0.0590414481559016D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.8415650255319866D+00, & 0.5000000000000000D+00, & 0.1584349744680134D+00, & 0.9409585518440984D+00, & 0.0590414481559016D+00, & 0.9409585518440984D+00, & 0.0590414481559016D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.2040816326530612D+00, & 0.2040816326530612D+00, & 0.2040816326530612D+00, & 0.2040816326530612D+00, & 0.0459183673469388D+00, & 0.0459183673469388D+00, & 0.0459183673469388D+00, & 0.0459183673469388D+00 /) x = x_save(1:n_save) y = y_save(1:n_save) w = w_save(1:n_save) return end subroutine rule07 ( n, x, y, w ) !*****************************************************************************80 ! !! rule07() returns the prism rule of precision 7. ! ! Discussion: ! ! We suppose we are given a quadrilateral Quad with vertices ! (0,0), (1,0), (1,1), (0,1), ! We call a rule with n points, returning coordinates ! x, y, and weights w. Then the integral I of f(x,y,z) over Quad is ! approximated by Q as follows: ! ! Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 12 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.9629100498862757D+00, & 0.5000000000000000D+00, & 0.0370899501137243D+00, & 0.5000000000000000D+00, & 0.9029898914592993D+00, & 0.9029898914592993D+00, & 0.0970101085407006D+00, & 0.0970101085407006D+00, & 0.6902772166041579D+00, & 0.6902772166041579D+00, & 0.3097227833958421D+00, & 0.3097227833958421D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.9629100498862757D+00, & 0.5000000000000000D+00, & 0.0370899501137243D+00, & 0.9029898914592993D+00, & 0.0970101085407006D+00, & 0.9029898914592993D+00, & 0.0970101085407006D+00, & 0.6902772166041579D+00, & 0.3097227833958421D+00, & 0.6902772166041579D+00, & 0.3097227833958421D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0604938271604938D+00, & 0.0604938271604938D+00, & 0.0604938271604938D+00, & 0.0604938271604938D+00, & 0.0593579436726576D+00, & 0.0593579436726576D+00, & 0.0593579436726576D+00, & 0.0593579436726576D+00, & 0.1301482291668486D+00, & 0.1301482291668486D+00, & 0.1301482291668486D+00, & 0.1301482291668486D+00 /) x = x_save(1:n_save) y = y_save(1:n_save) w = w_save(1:n_save) return end subroutine rule09 ( n, x, y, w ) !*****************************************************************************80 ! !! rule09() returns the prism rule of precision 9. ! ! Discussion: ! ! We suppose we are given a quadrilateral Quad with vertices ! (0,0), (1,0), (1,1), (0,1), ! We call a rule with n points, returning coordinates ! x, y, and weights w. Then the integral I of f(x,y,z) over Quad is ! approximated by Q as follows: ! ! Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 20 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.7444634284871845D+00, & 0.5000000000000000D+00, & 0.2555365715128155D+00, & 0.5000000000000000D+00, & 0.9698276290484189D+00, & 0.9698276290484189D+00, & 0.0301723709515812D+00, & 0.0301723709515812D+00, & 0.8454402752431720D+00, & 0.8454402752431720D+00, & 0.1545597247568281D+00, & 0.1545597247568281D+00, & 0.9593102205283611D+00, & 0.6724360126822018D+00, & 0.9593102205283611D+00, & 0.3275639873177982D+00, & 0.0406897794716389D+00, & 0.6724360126822018D+00, & 0.0406897794716389D+00, & 0.3275639873177982D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.7444634284871845D+00, & 0.5000000000000000D+00, & 0.2555365715128155D+00, & 0.9698276290484189D+00, & 0.0301723709515812D+00, & 0.9698276290484189D+00, & 0.0301723709515812D+00, & 0.8454402752431720D+00, & 0.1545597247568281D+00, & 0.8454402752431720D+00, & 0.1545597247568281D+00, & 0.6724360126822018D+00, & 0.9593102205283611D+00, & 0.3275639873177982D+00, & 0.9593102205283611D+00, & 0.6724360126822018D+00, & 0.0406897794716389D+00, & 0.3275639873177982D+00, & 0.0406897794716389D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.1135409901716873D+00, & 0.1135409901716873D+00, & 0.1135409901716873D+00, & 0.1135409901716873D+00, & 0.0106828079664439D+00, & 0.0106828079664439D+00, & 0.0106828079664439D+00, & 0.0106828079664439D+00, & 0.0535500902317154D+00, & 0.0535500902317154D+00, & 0.0535500902317154D+00, & 0.0535500902317154D+00, & 0.0361130558150767D+00, & 0.0361130558150767D+00, & 0.0361130558150767D+00, & 0.0361130558150767D+00, & 0.0361130558150767D+00, & 0.0361130558150767D+00, & 0.0361130558150767D+00, & 0.0361130558150767D+00 /) x = x_save(1:n_save) y = y_save(1:n_save) w = w_save(1:n_save) return end subroutine rule11 ( n, x, y, w ) !*****************************************************************************80 ! !! rule11() returns the prism rule of precision 11. ! ! Discussion: ! ! We suppose we are given a quadrilateral Quad with vertices ! (0,0), (1,0), (1,1), (0,1), ! We call a rule with n points, returning coordinates ! x, y, and weights w. Then the integral I of f(x,y,z) over Quad is ! approximated by Q as follows: ! ! Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 28 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.8573089148323030D+00, & 0.5000000000000000D+00, & 0.1426910851676970D+00, & 0.5000000000000000D+00, & 0.6368286050857298D+00, & 0.6368286050857298D+00, & 0.3631713949142702D+00, & 0.3631713949142702D+00, & 0.8183019661061506D+00, & 0.8183019661061506D+00, & 0.1816980338938495D+00, & 0.1816980338938495D+00, & 0.9758151943920168D+00, & 0.9077827168448191D+00, & 0.9758151943920168D+00, & 0.0922172831551808D+00, & 0.0241848056079833D+00, & 0.9077827168448191D+00, & 0.0241848056079833D+00, & 0.0922172831551808D+00, & 0.6731036000238227D+00, & 0.9677839357437956D+00, & 0.6731036000238227D+00, & 0.0322160642562044D+00, & 0.3268963999761773D+00, & 0.9677839357437956D+00, & 0.3268963999761773D+00, & 0.0322160642562044D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.8573089148323030D+00, & 0.5000000000000000D+00, & 0.1426910851676970D+00, & 0.6368286050857298D+00, & 0.3631713949142702D+00, & 0.6368286050857298D+00, & 0.3631713949142702D+00, & 0.8183019661061506D+00, & 0.1816980338938495D+00, & 0.8183019661061506D+00, & 0.1816980338938495D+00, & 0.9077827168448191D+00, & 0.9758151943920168D+00, & 0.0922172831551808D+00, & 0.9758151943920168D+00, & 0.9077827168448191D+00, & 0.0241848056079833D+00, & 0.0922172831551808D+00, & 0.0241848056079833D+00, & 0.9677839357437956D+00, & 0.6731036000238227D+00, & 0.0322160642562044D+00, & 0.6731036000238227D+00, & 0.9677839357437956D+00, & 0.3268963999761773D+00, & 0.0322160642562044D+00, & 0.3268963999761773D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0543501099671780D+00, & 0.0543501099671780D+00, & 0.0543501099671780D+00, & 0.0543501099671780D+00, & 0.0693185257459628D+00, & 0.0693185257459628D+00, & 0.0693185257459628D+00, & 0.0693185257459628D+00, & 0.0534834094695620D+00, & 0.0534834094695620D+00, & 0.0534834094695620D+00, & 0.0534834094695620D+00, & 0.0110186422787458D+00, & 0.0110186422787458D+00, & 0.0110186422787458D+00, & 0.0110186422787458D+00, & 0.0110186422787458D+00, & 0.0110186422787458D+00, & 0.0110186422787458D+00, & 0.0110186422787458D+00, & 0.0254053351299028D+00, & 0.0254053351299028D+00, & 0.0254053351299028D+00, & 0.0254053351299028D+00, & 0.0254053351299028D+00, & 0.0254053351299028D+00, & 0.0254053351299028D+00, & 0.0254053351299028D+00 /) x = x_save(1:n_save) y = y_save(1:n_save) w = w_save(1:n_save) return end subroutine rule13 ( n, x, y, w ) !*****************************************************************************80 ! !! rule13() returns the prism rule of precision 13. ! ! Discussion: ! ! We suppose we are given a quadrilateral Quad with vertices ! (0,0), (1,0), (1,1), (0,1), ! We call a rule with n points, returning coordinates ! x, y, and weights w. Then the integral I of f(x,y,z) over Quad is ! approximated by Q as follows: ! ! Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 37 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.5000000000000000D+00, & 0.9917306663066228D+00, & 0.5000000000000000D+00, & 0.0082693336933772D+00, & 0.5000000000000000D+00, & 0.8199307091835548D+00, & 0.5000000000000000D+00, & 0.1800692908164451D+00, & 0.5000000000000000D+00, & 0.9593892440398557D+00, & 0.9593892440398557D+00, & 0.0406107559601443D+00, & 0.0406107559601443D+00, & 0.6898142525933719D+00, & 0.6898142525933719D+00, & 0.3101857474066281D+00, & 0.3101857474066281D+00, & 0.8497771082566756D+00, & 0.8497771082566756D+00, & 0.1502228917433244D+00, & 0.1502228917433244D+00, & 0.8217986874983090D+00, & 0.9875344419529918D+00, & 0.8217986874983090D+00, & 0.0124655580470082D+00, & 0.1782013125016910D+00, & 0.9875344419529918D+00, & 0.1782013125016910D+00, & 0.0124655580470082D+00, & 0.6667699405823916D+00, & 0.9322138046335309D+00, & 0.6667699405823916D+00, & 0.0677861953664691D+00, & 0.3332300594176084D+00, & 0.9322138046335309D+00, & 0.3332300594176084D+00, & 0.0677861953664691D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.9917306663066228D+00, & 0.5000000000000000D+00, & 0.0082693336933772D+00, & 0.5000000000000000D+00, & 0.8199307091835548D+00, & 0.5000000000000000D+00, & 0.1800692908164451D+00, & 0.9593892440398557D+00, & 0.0406107559601443D+00, & 0.9593892440398557D+00, & 0.0406107559601443D+00, & 0.6898142525933719D+00, & 0.3101857474066281D+00, & 0.6898142525933719D+00, & 0.3101857474066281D+00, & 0.8497771082566756D+00, & 0.1502228917433244D+00, & 0.8497771082566756D+00, & 0.1502228917433244D+00, & 0.9875344419529918D+00, & 0.8217986874983090D+00, & 0.0124655580470082D+00, & 0.8217986874983090D+00, & 0.9875344419529918D+00, & 0.1782013125016910D+00, & 0.0124655580470082D+00, & 0.1782013125016910D+00, & 0.9322138046335309D+00, & 0.6667699405823916D+00, & 0.0677861953664691D+00, & 0.6667699405823916D+00, & 0.9322138046335309D+00, & 0.3332300594176084D+00, & 0.0677861953664691D+00, & 0.3332300594176084D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0749983360897274D+00, & 0.0095321563374632D+00, & 0.0095321563374632D+00, & 0.0095321563374632D+00, & 0.0095321563374632D+00, & 0.0461338672424518D+00, & 0.0461338672424518D+00, & 0.0461338672424518D+00, & 0.0461338672424518D+00, & 0.0098767851618631D+00, & 0.0098767851618631D+00, & 0.0098767851618631D+00, & 0.0098767851618631D+00, & 0.0578496298501713D+00, & 0.0578496298501713D+00, & 0.0578496298501713D+00, & 0.0578496298501713D+00, & 0.0343105241782509D+00, & 0.0343105241782509D+00, & 0.0343105241782509D+00, & 0.0343105241782509D+00, & 0.0083799450125954D+00, & 0.0083799450125954D+00, & 0.0083799450125954D+00, & 0.0083799450125954D+00, & 0.0083799450125954D+00, & 0.0083799450125954D+00, & 0.0083799450125954D+00, & 0.0083799450125954D+00, & 0.0283937815910886D+00, & 0.0283937815910886D+00, & 0.0283937815910886D+00, & 0.0283937815910886D+00, & 0.0283937815910886D+00, & 0.0283937815910886D+00, & 0.0283937815910886D+00, & 0.0283937815910886D+00 /) x = x_save(1:n_save) y = y_save(1:n_save) w = w_save(1:n_save) return end subroutine rule15 ( n, x, y, w ) !*****************************************************************************80 ! !! rule15() returns the prism rule of precision 15. ! ! Discussion: ! ! We suppose we are given a quadrilateral Quad with vertices ! (0,0), (1,0), (1,1), (0,1), ! We call a rule with n points, returning coordinates ! x, y, and weights w. Then the integral I of f(x,y,z) over Quad is ! approximated by Q as follows: ! ! Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 48 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.8990494939485035D+00, & 0.5000000000000000D+00, & 0.1009505060514965D+00, & 0.5000000000000000D+00, & 0.6519265131992299D+00, & 0.5000000000000000D+00, & 0.3480734868007701D+00, & 0.5000000000000000D+00, & 0.9412111350534427D+00, & 0.9412111350534427D+00, & 0.0587888649465573D+00, & 0.0587888649465573D+00, & 0.9888948997699514D+00, & 0.9888948997699514D+00, & 0.0111051002300486D+00, & 0.0111051002300486D+00, & 0.9043560679097518D+00, & 0.7836067866982954D+00, & 0.9043560679097518D+00, & 0.2163932133017046D+00, & 0.0956439320902482D+00, & 0.7836067866982954D+00, & 0.0956439320902482D+00, & 0.2163932133017046D+00, & 0.6523273995185263D+00, & 0.7893680970179033D+00, & 0.6523273995185263D+00, & 0.2106319029820967D+00, & 0.3476726004814737D+00, & 0.7893680970179033D+00, & 0.3476726004814737D+00, & 0.2106319029820967D+00, & 0.9902524218622659D+00, & 0.8487318159954835D+00, & 0.9902524218622659D+00, & 0.1512681840045165D+00, & 0.0097475781377340D+00, & 0.8487318159954835D+00, & 0.0097475781377340D+00, & 0.1512681840045165D+00, & 0.6324220779361581D+00, & 0.9778712599047559D+00, & 0.6324220779361581D+00, & 0.0221287400952441D+00, & 0.3675779220638419D+00, & 0.9778712599047559D+00, & 0.3675779220638419D+00, & 0.0221287400952441D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.8990494939485035D+00, & 0.5000000000000000D+00, & 0.1009505060514965D+00, & 0.5000000000000000D+00, & 0.6519265131992299D+00, & 0.5000000000000000D+00, & 0.3480734868007701D+00, & 0.9412111350534427D+00, & 0.0587888649465573D+00, & 0.9412111350534427D+00, & 0.0587888649465573D+00, & 0.9888948997699514D+00, & 0.0111051002300486D+00, & 0.9888948997699514D+00, & 0.0111051002300486D+00, & 0.7836067866982954D+00, & 0.9043560679097518D+00, & 0.2163932133017046D+00, & 0.9043560679097518D+00, & 0.7836067866982954D+00, & 0.0956439320902482D+00, & 0.2163932133017046D+00, & 0.0956439320902482D+00, & 0.7893680970179033D+00, & 0.6523273995185263D+00, & 0.2106319029820967D+00, & 0.6523273995185263D+00, & 0.7893680970179033D+00, & 0.3476726004814737D+00, & 0.2106319029820967D+00, & 0.3476726004814737D+00, & 0.8487318159954835D+00, & 0.9902524218622659D+00, & 0.1512681840045165D+00, & 0.9902524218622659D+00, & 0.8487318159954835D+00, & 0.0097475781377340D+00, & 0.1512681840045165D+00, & 0.0097475781377340D+00, & 0.9778712599047559D+00, & 0.6324220779361581D+00, & 0.0221287400952441D+00, & 0.6324220779361581D+00, & 0.9778712599047559D+00, & 0.3675779220638419D+00, & 0.0221287400952441D+00, & 0.3675779220638419D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0287332318158813D+00, & 0.0287332318158813D+00, & 0.0287332318158813D+00, & 0.0287332318158813D+00, & 0.0454214397418180D+00, & 0.0454214397418180D+00, & 0.0454214397418180D+00, & 0.0454214397418180D+00, & 0.0103096209471915D+00, & 0.0103096209471915D+00, & 0.0103096209471915D+00, & 0.0103096209471915D+00, & 0.0014834366214598D+00, & 0.0014834366214598D+00, & 0.0014834366214598D+00, & 0.0014834366214598D+00, & 0.0252872858693581D+00, & 0.0252872858693581D+00, & 0.0252872858693581D+00, & 0.0252872858693581D+00, & 0.0252872858693581D+00, & 0.0252872858693581D+00, & 0.0252872858693581D+00, & 0.0252872858693581D+00, & 0.0369591804082203D+00, & 0.0369591804082203D+00, & 0.0369591804082203D+00, & 0.0369591804082203D+00, & 0.0369591804082203D+00, & 0.0369591804082203D+00, & 0.0369591804082203D+00, & 0.0369591804082203D+00, & 0.0058182985361683D+00, & 0.0058182985361683D+00, & 0.0058182985361683D+00, & 0.0058182985361683D+00, & 0.0058182985361683D+00, & 0.0058182985361683D+00, & 0.0058182985361683D+00, & 0.0058182985361683D+00, & 0.0139613706230780D+00, & 0.0139613706230780D+00, & 0.0139613706230780D+00, & 0.0139613706230780D+00, & 0.0139613706230780D+00, & 0.0139613706230780D+00, & 0.0139613706230780D+00, & 0.0139613706230780D+00 /) x = x_save(1:n_save) y = y_save(1:n_save) w = w_save(1:n_save) return end subroutine rule17 ( n, x, y, w ) !*****************************************************************************80 ! !! rule17() returns the prism rule of precision 17. ! ! Discussion: ! ! We suppose we are given a quadrilateral Quad with vertices ! (0,0), (1,0), (1,1), (0,1), ! We call a rule with n points, returning coordinates ! x, y, and weights w. Then the integral I of f(x,y,z) over Quad is ! approximated by Q as follows: ! ! Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 60 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.7725214936045741D+00, & 0.5000000000000000D+00, & 0.2274785063954258D+00, & 0.5000000000000000D+00, & 0.9587485941095766D+00, & 0.5000000000000000D+00, & 0.0412514058904234D+00, & 0.5000000000000000D+00, & 0.5962771213070235D+00, & 0.5962771213070235D+00, & 0.4037228786929766D+00, & 0.4037228786929766D+00, & 0.9356643423477216D+00, & 0.9356643423477216D+00, & 0.0643356576522784D+00, & 0.0643356576522784D+00, & 0.8471940450447754D+00, & 0.8471940450447754D+00, & 0.1528059549552246D+00, & 0.1528059549552246D+00, & 0.7288414964415777D+00, & 0.7288414964415777D+00, & 0.2711585035584222D+00, & 0.2711585035584222D+00, & 0.9846021280457544D+00, & 0.9846021280457544D+00, & 0.0153978719542456D+00, & 0.0153978719542456D+00, & 0.8810441380400991D+00, & 0.6397388833575961D+00, & 0.8810441380400991D+00, & 0.3602611166424039D+00, & 0.1189558619599009D+00, & 0.6397388833575961D+00, & 0.1189558619599009D+00, & 0.3602611166424039D+00, & 0.9536971760991118D+00, & 0.7734493352071550D+00, & 0.9536971760991118D+00, & 0.2265506647928450D+00, & 0.0463028239008882D+00, & 0.7734493352071550D+00, & 0.0463028239008882D+00, & 0.2265506647928450D+00, & 0.8806243332195414D+00, & 0.9918939857507747D+00, & 0.8806243332195414D+00, & 0.0081060142492252D+00, & 0.1193756667804586D+00, & 0.9918939857507747D+00, & 0.1193756667804586D+00, & 0.0081060142492252D+00, & 0.9935473535223840D+00, & 0.6514037408944306D+00, & 0.9935473535223840D+00, & 0.3485962591055693D+00, & 0.0064526464776161D+00, & 0.6514037408944306D+00, & 0.0064526464776161D+00, & 0.3485962591055693D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.7725214936045741D+00, & 0.5000000000000000D+00, & 0.2274785063954258D+00, & 0.5000000000000000D+00, & 0.9587485941095766D+00, & 0.5000000000000000D+00, & 0.0412514058904234D+00, & 0.5962771213070235D+00, & 0.4037228786929766D+00, & 0.5962771213070235D+00, & 0.4037228786929766D+00, & 0.9356643423477216D+00, & 0.0643356576522784D+00, & 0.9356643423477216D+00, & 0.0643356576522784D+00, & 0.8471940450447754D+00, & 0.1528059549552246D+00, & 0.8471940450447754D+00, & 0.1528059549552246D+00, & 0.7288414964415777D+00, & 0.2711585035584222D+00, & 0.7288414964415777D+00, & 0.2711585035584222D+00, & 0.9846021280457544D+00, & 0.0153978719542456D+00, & 0.9846021280457544D+00, & 0.0153978719542456D+00, & 0.6397388833575961D+00, & 0.8810441380400991D+00, & 0.3602611166424039D+00, & 0.8810441380400991D+00, & 0.6397388833575961D+00, & 0.1189558619599009D+00, & 0.3602611166424039D+00, & 0.1189558619599009D+00, & 0.7734493352071550D+00, & 0.9536971760991118D+00, & 0.2265506647928450D+00, & 0.9536971760991118D+00, & 0.7734493352071550D+00, & 0.0463028239008882D+00, & 0.2265506647928450D+00, & 0.0463028239008882D+00, & 0.9918939857507747D+00, & 0.8806243332195414D+00, & 0.0081060142492252D+00, & 0.8806243332195414D+00, & 0.9918939857507747D+00, & 0.1193756667804586D+00, & 0.0081060142492252D+00, & 0.1193756667804586D+00, & 0.6514037408944306D+00, & 0.9935473535223840D+00, & 0.3485962591055693D+00, & 0.9935473535223840D+00, & 0.6514037408944306D+00, & 0.0064526464776161D+00, & 0.3485962591055693D+00, & 0.0064526464776161D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0347010944218050D+00, & 0.0347010944218050D+00, & 0.0347010944218050D+00, & 0.0347010944218050D+00, & 0.0163369559518662D+00, & 0.0163369559518662D+00, & 0.0163369559518662D+00, & 0.0163369559518662D+00, & 0.0357141049055499D+00, & 0.0357141049055499D+00, & 0.0357141049055499D+00, & 0.0357141049055499D+00, & 0.0091789771213129D+00, & 0.0091789771213129D+00, & 0.0091789771213129D+00, & 0.0091789771213129D+00, & 0.0215639249095360D+00, & 0.0215639249095360D+00, & 0.0215639249095360D+00, & 0.0215639249095360D+00, & 0.0336390683029792D+00, & 0.0336390683029792D+00, & 0.0336390683029792D+00, & 0.0336390683029792D+00, & 0.0019139083537275D+00, & 0.0019139083537275D+00, & 0.0019139083537275D+00, & 0.0019139083537275D+00, & 0.0257941811962338D+00, & 0.0257941811962338D+00, & 0.0257941811962338D+00, & 0.0257941811962338D+00, & 0.0257941811962338D+00, & 0.0257941811962338D+00, & 0.0257941811962338D+00, & 0.0257941811962338D+00, & 0.0137066000593229D+00, & 0.0137066000593229D+00, & 0.0137066000593229D+00, & 0.0137066000593229D+00, & 0.0137066000593229D+00, & 0.0137066000593229D+00, & 0.0137066000593229D+00, & 0.0137066000593229D+00, & 0.0037799525970642D+00, & 0.0037799525970642D+00, & 0.0037799525970642D+00, & 0.0037799525970642D+00, & 0.0037799525970642D+00, & 0.0037799525970642D+00, & 0.0037799525970642D+00, & 0.0037799525970642D+00, & 0.0051952491639908D+00, & 0.0051952491639908D+00, & 0.0051952491639908D+00, & 0.0051952491639908D+00, & 0.0051952491639908D+00, & 0.0051952491639908D+00, & 0.0051952491639908D+00, & 0.0051952491639908D+00 /) x = x_save(1:n_save) y = y_save(1:n_save) w = w_save(1:n_save) return end subroutine rule19 ( n, x, y, w ) !*****************************************************************************80 ! !! rule19() returns the prism rule of precision 19. ! ! Discussion: ! ! We suppose we are given a quadrilateral Quad with vertices ! (0,0), (1,0), (1,1), (0,1), ! We call a rule with n points, returning coordinates ! x, y, and weights w. Then the integral I of f(x,y,z) over Quad is ! approximated by Q as follows: ! ! Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 72 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.8571721297863971D+00, & 0.5000000000000000D+00, & 0.1428278702136029D+00, & 0.5000000000000000D+00, & 0.6328360260604818D+00, & 0.5000000000000000D+00, & 0.3671639739395182D+00, & 0.5000000000000000D+00, & 0.9822171346289836D+00, & 0.5000000000000000D+00, & 0.0177828653710164D+00, & 0.5000000000000000D+00, & 0.9016898647338425D+00, & 0.9016898647338425D+00, & 0.0983101352661575D+00, & 0.0983101352661575D+00, & 0.9960827029485674D+00, & 0.9960827029485674D+00, & 0.0039172970514326D+00, & 0.0039172970514326D+00, & 0.9647248013998047D+00, & 0.9647248013998047D+00, & 0.0352751986001953D+00, & 0.0352751986001953D+00, & 0.7551391286846717D+00, & 0.6333201572972811D+00, & 0.7551391286846717D+00, & 0.3666798427027189D+00, & 0.2448608713153283D+00, & 0.6333201572972811D+00, & 0.2448608713153283D+00, & 0.3666798427027189D+00, & 0.6953671028876249D+00, & 0.9939464665708677D+00, & 0.6953671028876249D+00, & 0.0060535334291323D+00, & 0.3046328971123751D+00, & 0.9939464665708677D+00, & 0.3046328971123751D+00, & 0.0060535334291323D+00, & 0.8585839606548726D+00, & 0.7562459386080489D+00, & 0.8585839606548726D+00, & 0.2437540613919512D+00, & 0.1414160393451274D+00, & 0.7562459386080489D+00, & 0.1414160393451274D+00, & 0.2437540613919512D+00, & 0.6327200039056480D+00, & 0.9344512170772521D+00, & 0.6327200039056480D+00, & 0.0655487829227479D+00, & 0.3672799960943520D+00, & 0.9344512170772521D+00, & 0.3672799960943520D+00, & 0.0655487829227479D+00, & 0.8100176993466282D+00, & 0.9631514779035646D+00, & 0.8100176993466282D+00, & 0.0368485220964354D+00, & 0.1899823006533718D+00, & 0.9631514779035646D+00, & 0.1899823006533718D+00, & 0.0368485220964354D+00, & 0.9008357923592984D+00, & 0.9942232653419869D+00, & 0.9008357923592984D+00, & 0.0057767346580131D+00, & 0.0991642076407016D+00, & 0.9942232653419869D+00, & 0.0991642076407016D+00, & 0.0057767346580131D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.8571721297863971D+00, & 0.5000000000000000D+00, & 0.1428278702136029D+00, & 0.5000000000000000D+00, & 0.6328360260604818D+00, & 0.5000000000000000D+00, & 0.3671639739395182D+00, & 0.5000000000000000D+00, & 0.9822171346289836D+00, & 0.5000000000000000D+00, & 0.0177828653710164D+00, & 0.9016898647338425D+00, & 0.0983101352661575D+00, & 0.9016898647338425D+00, & 0.0983101352661575D+00, & 0.9960827029485674D+00, & 0.0039172970514326D+00, & 0.9960827029485674D+00, & 0.0039172970514326D+00, & 0.9647248013998047D+00, & 0.0352751986001953D+00, & 0.9647248013998047D+00, & 0.0352751986001953D+00, & 0.6333201572972811D+00, & 0.7551391286846717D+00, & 0.3666798427027189D+00, & 0.7551391286846717D+00, & 0.6333201572972811D+00, & 0.2448608713153283D+00, & 0.3666798427027189D+00, & 0.2448608713153283D+00, & 0.9939464665708677D+00, & 0.6953671028876249D+00, & 0.0060535334291323D+00, & 0.6953671028876249D+00, & 0.9939464665708677D+00, & 0.3046328971123751D+00, & 0.0060535334291323D+00, & 0.3046328971123751D+00, & 0.7562459386080489D+00, & 0.8585839606548726D+00, & 0.2437540613919512D+00, & 0.8585839606548726D+00, & 0.7562459386080489D+00, & 0.1414160393451274D+00, & 0.2437540613919512D+00, & 0.1414160393451274D+00, & 0.9344512170772521D+00, & 0.6327200039056480D+00, & 0.0655487829227479D+00, & 0.6327200039056480D+00, & 0.9344512170772521D+00, & 0.3672799960943520D+00, & 0.0655487829227479D+00, & 0.3672799960943520D+00, & 0.9631514779035646D+00, & 0.8100176993466282D+00, & 0.0368485220964354D+00, & 0.8100176993466282D+00, & 0.9631514779035646D+00, & 0.1899823006533718D+00, & 0.0368485220964354D+00, & 0.1899823006533718D+00, & 0.9942232653419869D+00, & 0.9008357923592984D+00, & 0.0057767346580131D+00, & 0.9008357923592984D+00, & 0.9942232653419869D+00, & 0.0991642076407016D+00, & 0.0057767346580131D+00, & 0.0991642076407016D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0244343411732072D+00, & 0.0244343411732072D+00, & 0.0244343411732072D+00, & 0.0244343411732072D+00, & 0.0348269032306206D+00, & 0.0348269032306206D+00, & 0.0348269032306206D+00, & 0.0348269032306206D+00, & 0.0087173958729720D+00, & 0.0087173958729720D+00, & 0.0087173958729720D+00, & 0.0087173958729720D+00, & 0.0119511367906989D+00, & 0.0119511367906989D+00, & 0.0119511367906989D+00, & 0.0119511367906989D+00, & 0.0004044287944779D+00, & 0.0004044287944779D+00, & 0.0004044287944779D+00, & 0.0004044287944779D+00, & 0.0043602620085889D+00, & 0.0043602620085889D+00, & 0.0043602620085889D+00, & 0.0043602620085889D+00, & 0.0294314582100140D+00, & 0.0294314582100140D+00, & 0.0294314582100140D+00, & 0.0294314582100140D+00, & 0.0294314582100140D+00, & 0.0294314582100140D+00, & 0.0294314582100140D+00, & 0.0294314582100140D+00, & 0.0040448944029145D+00, & 0.0040448944029145D+00, & 0.0040448944029145D+00, & 0.0040448944029145D+00, & 0.0040448944029145D+00, & 0.0040448944029145D+00, & 0.0040448944029145D+00, & 0.0040448944029145D+00, & 0.0207116547458512D+00, & 0.0207116547458512D+00, & 0.0207116547458512D+00, & 0.0207116547458512D+00, & 0.0207116547458512D+00, & 0.0207116547458512D+00, & 0.0207116547458512D+00, & 0.0207116547458512D+00, & 0.0162472703731490D+00, & 0.0162472703731490D+00, & 0.0162472703731490D+00, & 0.0162472703731490D+00, & 0.0162472703731490D+00, & 0.0162472703731490D+00, & 0.0162472703731490D+00, & 0.0162472703731490D+00, & 0.0097451380991026D+00, & 0.0097451380991026D+00, & 0.0097451380991026D+00, & 0.0097451380991026D+00, & 0.0097451380991026D+00, & 0.0097451380991026D+00, & 0.0097451380991026D+00, & 0.0097451380991026D+00, & 0.0024723502336859D+00, & 0.0024723502336859D+00, & 0.0024723502336859D+00, & 0.0024723502336859D+00, & 0.0024723502336859D+00, & 0.0024723502336859D+00, & 0.0024723502336859D+00, & 0.0024723502336859D+00 /) x = x_save(1:n_save) y = y_save(1:n_save) w = w_save(1:n_save) return end subroutine rule21 ( n, x, y, w ) !*****************************************************************************80 ! !! rule21() returns the quadrilateral rule of precision 21. ! ! Discussion: ! ! We suppose we are given a quadrilateral Quad with vertices ! (0,0), (1,0), (1,1), (0,1), ! We call a rule with n points, returning coordinates ! x, y, and weights w. Then the integral I of f(x,y,z) over Quad is ! approximated by Q as follows: ! ! Q = area(Quad) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 29 April 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 85 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.5000000000000000D+00, & 0.7366755371791120D+00, & 0.5000000000000000D+00, & 0.2633244628208880D+00, & 0.5000000000000000D+00, & 0.9176392148779342D+00, & 0.5000000000000000D+00, & 0.0823607851220658D+00, & 0.5000000000000000D+00, & 0.6286859503036145D+00, & 0.6286859503036145D+00, & 0.3713140496963855D+00, & 0.3713140496963855D+00, & 0.9816689160578117D+00, & 0.9816689160578117D+00, & 0.0183310839421883D+00, & 0.0183310839421883D+00, & 0.9312259626898258D+00, & 0.9312259626898258D+00, & 0.0687740373101742D+00, & 0.0687740373101742D+00, & 0.7484489812596729D+00, & 0.7484489812596729D+00, & 0.2515510187403271D+00, & 0.2515510187403271D+00, & 0.8521660875977003D+00, & 0.8521660875977003D+00, & 0.1478339124022997D+00, & 0.1478339124022997D+00, & 0.6209390927383510D+00, & 0.8370731099776589D+00, & 0.6209390927383510D+00, & 0.1629268900223411D+00, & 0.3790609072616490D+00, & 0.8370731099776589D+00, & 0.3790609072616490D+00, & 0.1629268900223411D+00, & 0.7400784831563976D+00, & 0.9123236876354603D+00, & 0.7400784831563976D+00, & 0.0876763123645397D+00, & 0.2599215168436024D+00, & 0.9123236876354603D+00, & 0.2599215168436024D+00, & 0.0876763123645397D+00, & 0.9661209679608770D+00, & 0.6353495920508325D+00, & 0.9661209679608770D+00, & 0.3646504079491675D+00, & 0.0338790320391230D+00, & 0.6353495920508325D+00, & 0.0338790320391230D+00, & 0.3646504079491675D+00, & 0.9674511629120053D+00, & 0.8375197837185376D+00, & 0.9674511629120053D+00, & 0.1624802162814624D+00, & 0.0325488370879947D+00, & 0.8375197837185376D+00, & 0.0325488370879947D+00, & 0.1624802162814624D+00, & 0.9952277337647621D+00, & 0.7444581255885834D+00, & 0.9952277337647621D+00, & 0.2555418744114166D+00, & 0.0047722662352379D+00, & 0.7444581255885834D+00, & 0.0047722662352379D+00, & 0.2555418744114166D+00, & 0.9155676229879507D+00, & 0.9944803056932863D+00, & 0.9155676229879507D+00, & 0.0055196943067138D+00, & 0.0844323770120493D+00, & 0.9944803056932863D+00, & 0.0844323770120493D+00, & 0.0055196943067138D+00, & 0.5457348004611456D+00, & 0.9920085742447995D+00, & 0.5457348004611456D+00, & 0.0079914257552005D+00, & 0.4542651995388544D+00, & 0.9920085742447995D+00, & 0.4542651995388544D+00, & 0.0079914257552005D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.7366755371791120D+00, & 0.5000000000000000D+00, & 0.2633244628208880D+00, & 0.5000000000000000D+00, & 0.9176392148779342D+00, & 0.5000000000000000D+00, & 0.0823607851220658D+00, & 0.6286859503036145D+00, & 0.3713140496963855D+00, & 0.6286859503036145D+00, & 0.3713140496963855D+00, & 0.9816689160578117D+00, & 0.0183310839421883D+00, & 0.9816689160578117D+00, & 0.0183310839421883D+00, & 0.9312259626898258D+00, & 0.0687740373101742D+00, & 0.9312259626898258D+00, & 0.0687740373101742D+00, & 0.7484489812596729D+00, & 0.2515510187403271D+00, & 0.7484489812596729D+00, & 0.2515510187403271D+00, & 0.8521660875977003D+00, & 0.1478339124022997D+00, & 0.8521660875977003D+00, & 0.1478339124022997D+00, & 0.8370731099776589D+00, & 0.6209390927383510D+00, & 0.1629268900223411D+00, & 0.6209390927383510D+00, & 0.8370731099776589D+00, & 0.3790609072616490D+00, & 0.1629268900223411D+00, & 0.3790609072616490D+00, & 0.9123236876354603D+00, & 0.7400784831563976D+00, & 0.0876763123645397D+00, & 0.7400784831563976D+00, & 0.9123236876354603D+00, & 0.2599215168436024D+00, & 0.0876763123645397D+00, & 0.2599215168436024D+00, & 0.6353495920508325D+00, & 0.9661209679608770D+00, & 0.3646504079491675D+00, & 0.9661209679608770D+00, & 0.6353495920508325D+00, & 0.0338790320391230D+00, & 0.3646504079491675D+00, & 0.0338790320391230D+00, & 0.8375197837185376D+00, & 0.9674511629120053D+00, & 0.1624802162814624D+00, & 0.9674511629120053D+00, & 0.8375197837185376D+00, & 0.0325488370879947D+00, & 0.1624802162814624D+00, & 0.0325488370879947D+00, & 0.7444581255885834D+00, & 0.9952277337647621D+00, & 0.2555418744114166D+00, & 0.9952277337647621D+00, & 0.7444581255885834D+00, & 0.0047722662352379D+00, & 0.2555418744114166D+00, & 0.0047722662352379D+00, & 0.9944803056932863D+00, & 0.9155676229879507D+00, & 0.0055196943067138D+00, & 0.9155676229879507D+00, & 0.9944803056932863D+00, & 0.0844323770120493D+00, & 0.0055196943067138D+00, & 0.0844323770120493D+00, & 0.9920085742447995D+00, & 0.5457348004611456D+00, & 0.0079914257552005D+00, & 0.5457348004611456D+00, & 0.9920085742447995D+00, & 0.4542651995388544D+00, & 0.0079914257552005D+00, & 0.4542651995388544D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0337982100890357D+00, & 0.0266985396964851D+00, & 0.0266985396964851D+00, & 0.0266985396964851D+00, & 0.0266985396964851D+00, & 0.0165643020012361D+00, & 0.0165643020012361D+00, & 0.0165643020012361D+00, & 0.0165643020012361D+00, & 0.0299711209115851D+00, & 0.0299711209115851D+00, & 0.0299711209115851D+00, & 0.0299711209115851D+00, & 0.0020573653223050D+00, & 0.0020573653223050D+00, & 0.0020573653223050D+00, & 0.0020573653223050D+00, & 0.0076509102321939D+00, & 0.0076509102321939D+00, & 0.0076509102321939D+00, & 0.0076509102321939D+00, & 0.0241979479733988D+00, & 0.0241979479733988D+00, & 0.0241979479733988D+00, & 0.0241979479733988D+00, & 0.0153790853703291D+00, & 0.0153790853703291D+00, & 0.0153790853703291D+00, & 0.0153790853703291D+00, & 0.0216821248783021D+00, & 0.0216821248783021D+00, & 0.0216821248783021D+00, & 0.0216821248783021D+00, & 0.0216821248783021D+00, & 0.0216821248783021D+00, & 0.0216821248783021D+00, & 0.0216821248783021D+00, & 0.0139697793510318D+00, & 0.0139697793510318D+00, & 0.0139697793510318D+00, & 0.0139697793510318D+00, & 0.0139697793510318D+00, & 0.0139697793510318D+00, & 0.0139697793510318D+00, & 0.0139697793510318D+00, & 0.0092822068602285D+00, & 0.0092822068602285D+00, & 0.0092822068602285D+00, & 0.0092822068602285D+00, & 0.0092822068602285D+00, & 0.0092822068602285D+00, & 0.0092822068602285D+00, & 0.0092822068602285D+00, & 0.0071939002127551D+00, & 0.0071939002127551D+00, & 0.0071939002127551D+00, & 0.0071939002127551D+00, & 0.0071939002127551D+00, & 0.0071939002127551D+00, & 0.0071939002127551D+00, & 0.0071939002127551D+00, & 0.0028773801106233D+00, & 0.0028773801106233D+00, & 0.0028773801106233D+00, & 0.0028773801106233D+00, & 0.0028773801106233D+00, & 0.0028773801106233D+00, & 0.0028773801106233D+00, & 0.0028773801106233D+00, & 0.0018821205326004D+00, & 0.0018821205326004D+00, & 0.0018821205326004D+00, & 0.0018821205326004D+00, & 0.0018821205326004D+00, & 0.0018821205326004D+00, & 0.0018821205326004D+00, & 0.0018821205326004D+00, & 0.0026280760395628D+00, & 0.0026280760395628D+00, & 0.0026280760395628D+00, & 0.0026280760395628D+00, & 0.0026280760395628D+00, & 0.0026280760395628D+00, & 0.0026280760395628D+00, & 0.0026280760395628D+00 /) x = x_save(1:n_save) y = y_save(1:n_save) w = w_save(1:n_save) return end