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 subroutine rule_order ( p, order ) !*****************************************************************************80 ! !! rule_order() returns the order of a 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: ! ! 23 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 p: the precision, 0 <= p <= 20. ! ! Output: ! ! integer order: the order of the rule. ! implicit none integer order integer, dimension ( 0:20 ) :: order_vec = (/ & 1, & 1, 3, 6, 6, 7, & 12, 15, 16, 19, 25, & 28, 33, 37, 42, 49, & 55, 60, 67, 73, 79 /) integer p if ( p < 0 ) then write ( *, '(a)' ) 'rule_order(): Input p < 0.' end if if ( 20 < p ) then write ( *, '(a)' ) 'rule_order(): Input 20 < p.' end if order = order_vec(p) return end subroutine rule00 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule00() returns the rule of precision 0. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 1 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.5000000000000000D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule01 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule01() returns the rule of precision 1. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 1 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.5000000000000000D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule02 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule02() returns the rule of precision 2. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 3 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.1666666666666667D+00, & 0.6666666666666666D+00, & 0.1666666666666667D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.6666666666666666D+00, & 0.1666666666666667D+00, & 0.1666666666666667D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.1666666666666666D+00, & 0.1666666666666667D+00, & 0.6666666666666665D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.1666666666666667D+00, & 0.1666666666666667D+00, & 0.1666666666666667D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule03 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule03() returns the rule of precision 3. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 6 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.4459484909159649D+00, & 0.1081030181680702D+00, & 0.4459484909159649D+00, & 0.0915762135097707D+00, & 0.8168475729804585D+00, & 0.0915762135097707D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.1081030181680702D+00, & 0.4459484909159649D+00, & 0.4459484909159649D+00, & 0.8168475729804585D+00, & 0.0915762135097707D+00, & 0.0915762135097707D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.4459484909159649D+00, & 0.4459484909159649D+00, & 0.1081030181680702D+00, & 0.0915762135097707D+00, & 0.0915762135097707D+00, & 0.8168475729804585D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.1116907948390057D+00, & 0.1116907948390057D+00, & 0.1116907948390057D+00, & 0.0549758718276609D+00, & 0.0549758718276609D+00, & 0.0549758718276609D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule04 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule04() returns the rule of precision 4. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 6 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.4459484909159649D+00, & 0.1081030181680702D+00, & 0.4459484909159649D+00, & 0.0915762135097707D+00, & 0.8168475729804585D+00, & 0.0915762135097707D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.1081030181680702D+00, & 0.4459484909159649D+00, & 0.4459484909159649D+00, & 0.8168475729804585D+00, & 0.0915762135097707D+00, & 0.0915762135097707D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.4459484909159649D+00, & 0.4459484909159649D+00, & 0.1081030181680702D+00, & 0.0915762135097707D+00, & 0.0915762135097707D+00, & 0.8168475729804585D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.1116907948390057D+00, & 0.1116907948390057D+00, & 0.1116907948390057D+00, & 0.0549758718276609D+00, & 0.0549758718276609D+00, & 0.0549758718276609D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule05 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule05() returns the rule of precision 5. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 7 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00, & 0.1012865073234563D+00, & 0.7974269853530873D+00, & 0.1012865073234563D+00, & 0.4701420641051151D+00, & 0.0597158717897698D+00, & 0.4701420641051151D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00, & 0.7974269853530873D+00, & 0.1012865073234563D+00, & 0.1012865073234563D+00, & 0.0597158717897698D+00, & 0.4701420641051151D+00, & 0.4701420641051151D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00, & 0.1012865073234563D+00, & 0.1012865073234563D+00, & 0.7974269853530873D+00, & 0.4701420641051151D+00, & 0.4701420641051150D+00, & 0.0597158717897698D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.1125000000000000D+00, & 0.0629695902724136D+00, & 0.0629695902724136D+00, & 0.0629695902724136D+00, & 0.0661970763942531D+00, & 0.0661970763942531D+00, & 0.0661970763942531D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule06 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule06() returns the rule of precision 6. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 12 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.0630890144915022D+00, & 0.8738219710169955D+00, & 0.0630890144915022D+00, & 0.2492867451709104D+00, & 0.5014265096581791D+00, & 0.2492867451709104D+00, & 0.0531450498448169D+00, & 0.6365024991213987D+00, & 0.3103524510337844D+00, & 0.6365024991213987D+00, & 0.3103524510337844D+00, & 0.0531450498448169D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.8738219710169955D+00, & 0.0630890144915022D+00, & 0.0630890144915022D+00, & 0.5014265096581791D+00, & 0.2492867451709104D+00, & 0.2492867451709104D+00, & 0.6365024991213987D+00, & 0.0531450498448169D+00, & 0.6365024991213987D+00, & 0.3103524510337844D+00, & 0.0531450498448169D+00, & 0.3103524510337844D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.0630890144915023D+00, & 0.0630890144915022D+00, & 0.8738219710169957D+00, & 0.2492867451709104D+00, & 0.2492867451709104D+00, & 0.5014265096581791D+00, & 0.3103524510337844D+00, & 0.3103524510337844D+00, & 0.0531450498448169D+00, & 0.0531450498448169D+00, & 0.6365024991213987D+00, & 0.6365024991213987D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0254224531851034D+00, & 0.0254224531851034D+00, & 0.0254224531851034D+00, & 0.0583931378631897D+00, & 0.0583931378631897D+00, & 0.0583931378631897D+00, & 0.0414255378091868D+00, & 0.0414255378091868D+00, & 0.0414255378091868D+00, & 0.0414255378091868D+00, & 0.0414255378091868D+00, & 0.0414255378091868D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule07 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule07() returns the rule of precision 7. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 15 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.0337306485545878D+00, & 0.9325387028908243D+00, & 0.0337306485545878D+00, & 0.2415773825954036D+00, & 0.5168452348091929D+00, & 0.2415773825954036D+00, & 0.4743096925047182D+00, & 0.0513806149905635D+00, & 0.4743096925047182D+00, & 0.0470366446525952D+00, & 0.7542800405500532D+00, & 0.1986833147973516D+00, & 0.7542800405500532D+00, & 0.1986833147973516D+00, & 0.0470366446525952D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.9325387028908243D+00, & 0.0337306485545878D+00, & 0.0337306485545878D+00, & 0.5168452348091929D+00, & 0.2415773825954036D+00, & 0.2415773825954036D+00, & 0.0513806149905635D+00, & 0.4743096925047182D+00, & 0.4743096925047182D+00, & 0.7542800405500532D+00, & 0.0470366446525952D+00, & 0.7542800405500532D+00, & 0.1986833147973516D+00, & 0.0470366446525952D+00, & 0.1986833147973516D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.0337306485545878D+00, & 0.0337306485545878D+00, & 0.9325387028908243D+00, & 0.2415773825954036D+00, & 0.2415773825954036D+00, & 0.5168452348091930D+00, & 0.4743096925047183D+00, & 0.4743096925047182D+00, & 0.0513806149905636D+00, & 0.1986833147973517D+00, & 0.1986833147973516D+00, & 0.0470366446525953D+00, & 0.0470366446525953D+00, & 0.7542800405500532D+00, & 0.7542800405500533D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0082725250553961D+00, & 0.0082725250553961D+00, & 0.0082725250553961D+00, & 0.0639720856150778D+00, & 0.0639720856150778D+00, & 0.0639720856150778D+00, & 0.0385433230929930D+00, & 0.0385433230929930D+00, & 0.0385433230929930D+00, & 0.0279393664515999D+00, & 0.0279393664515999D+00, & 0.0279393664515999D+00, & 0.0279393664515999D+00, & 0.0279393664515999D+00, & 0.0279393664515999D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule08 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule08() returns the rule of precision 8. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 16 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00, & 0.4592925882927231D+00, & 0.0814148234145537D+00, & 0.4592925882927231D+00, & 0.1705693077517602D+00, & 0.6588613844964796D+00, & 0.1705693077517602D+00, & 0.0505472283170310D+00, & 0.8989055433659381D+00, & 0.0505472283170310D+00, & 0.0083947774099576D+00, & 0.7284923929554042D+00, & 0.2631128296346381D+00, & 0.7284923929554042D+00, & 0.2631128296346381D+00, & 0.0083947774099576D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00, & 0.0814148234145537D+00, & 0.4592925882927231D+00, & 0.4592925882927231D+00, & 0.6588613844964796D+00, & 0.1705693077517602D+00, & 0.1705693077517602D+00, & 0.8989055433659381D+00, & 0.0505472283170310D+00, & 0.0505472283170310D+00, & 0.7284923929554042D+00, & 0.0083947774099576D+00, & 0.7284923929554042D+00, & 0.2631128296346381D+00, & 0.0083947774099576D+00, & 0.2631128296346381D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00, & 0.4592925882927232D+00, & 0.4592925882927231D+00, & 0.0814148234145537D+00, & 0.1705693077517602D+00, & 0.1705693077517602D+00, & 0.6588613844964795D+00, & 0.0505472283170310D+00, & 0.0505472283170310D+00, & 0.8989055433659381D+00, & 0.2631128296346382D+00, & 0.2631128296346382D+00, & 0.0083947774099576D+00, & 0.0083947774099576D+00, & 0.7284923929554044D+00, & 0.7284923929554044D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0721578038388936D+00, & 0.0475458171336423D+00, & 0.0475458171336423D+00, & 0.0475458171336423D+00, & 0.0516086852673591D+00, & 0.0516086852673591D+00, & 0.0516086852673591D+00, & 0.0162292488115990D+00, & 0.0162292488115990D+00, & 0.0162292488115990D+00, & 0.0136151570872175D+00, & 0.0136151570872175D+00, & 0.0136151570872175D+00, & 0.0136151570872175D+00, & 0.0136151570872175D+00, & 0.0136151570872175D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule09 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule09() returns the rule of precision 9. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 19 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00, & 0.4370895914929366D+00, & 0.1258208170141267D+00, & 0.4370895914929366D+00, & 0.1882035356190327D+00, & 0.6235929287619345D+00, & 0.1882035356190327D+00, & 0.4896825191987376D+00, & 0.0206349616025248D+00, & 0.4896825191987376D+00, & 0.0447295133944527D+00, & 0.9105409732110945D+00, & 0.0447295133944527D+00, & 0.0368384120547363D+00, & 0.7411985987844980D+00, & 0.2219629891607657D+00, & 0.7411985987844980D+00, & 0.2219629891607657D+00, & 0.0368384120547363D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00, & 0.1258208170141267D+00, & 0.4370895914929366D+00, & 0.4370895914929366D+00, & 0.6235929287619345D+00, & 0.1882035356190327D+00, & 0.1882035356190327D+00, & 0.0206349616025248D+00, & 0.4896825191987376D+00, & 0.4896825191987376D+00, & 0.9105409732110945D+00, & 0.0447295133944527D+00, & 0.0447295133944527D+00, & 0.7411985987844980D+00, & 0.0368384120547363D+00, & 0.7411985987844980D+00, & 0.2219629891607657D+00, & 0.0368384120547363D+00, & 0.2219629891607657D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00, & 0.4370895914929366D+00, & 0.4370895914929366D+00, & 0.1258208170141267D+00, & 0.1882035356190327D+00, & 0.1882035356190327D+00, & 0.6235929287619344D+00, & 0.4896825191987376D+00, & 0.4896825191987376D+00, & 0.0206349616025248D+00, & 0.0447295133944527D+00, & 0.0447295133944528D+00, & 0.9105409732110945D+00, & 0.2219629891607657D+00, & 0.2219629891607657D+00, & 0.0368384120547363D+00, & 0.0368384120547363D+00, & 0.7411985987844980D+00, & 0.7411985987844980D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0485678981413994D+00, & 0.0389137705023871D+00, & 0.0389137705023871D+00, & 0.0389137705023871D+00, & 0.0398238694636051D+00, & 0.0398238694636051D+00, & 0.0398238694636051D+00, & 0.0156673501135695D+00, & 0.0156673501135695D+00, & 0.0156673501135695D+00, & 0.0127888378293490D+00, & 0.0127888378293490D+00, & 0.0127888378293490D+00, & 0.0216417696886447D+00, & 0.0216417696886447D+00, & 0.0216417696886447D+00, & 0.0216417696886447D+00, & 0.0216417696886447D+00, & 0.0216417696886447D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule10 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule10() returns the rule of precision 10. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 25 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00, & 0.0320553732169435D+00, & 0.9358892535661130D+00, & 0.0320553732169435D+00, & 0.1421611010565644D+00, & 0.7156777978868712D+00, & 0.1421611010565644D+00, & 0.3218129952888354D+00, & 0.5300541189273440D+00, & 0.1481328857838206D+00, & 0.5300541189273440D+00, & 0.1481328857838206D+00, & 0.3218129952888354D+00, & 0.0296198894887298D+00, & 0.6012333286834592D+00, & 0.3691467818278110D+00, & 0.6012333286834592D+00, & 0.3691467818278110D+00, & 0.0296198894887298D+00, & 0.0283676653399385D+00, & 0.8079306009228790D+00, & 0.1637017337371825D+00, & 0.8079306009228790D+00, & 0.1637017337371825D+00, & 0.0283676653399385D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00, & 0.9358892535661130D+00, & 0.0320553732169435D+00, & 0.0320553732169435D+00, & 0.7156777978868712D+00, & 0.1421611010565644D+00, & 0.1421611010565644D+00, & 0.5300541189273440D+00, & 0.3218129952888354D+00, & 0.5300541189273440D+00, & 0.1481328857838206D+00, & 0.3218129952888354D+00, & 0.1481328857838206D+00, & 0.6012333286834592D+00, & 0.0296198894887298D+00, & 0.6012333286834592D+00, & 0.3691467818278110D+00, & 0.0296198894887298D+00, & 0.3691467818278110D+00, & 0.8079306009228790D+00, & 0.0283676653399385D+00, & 0.8079306009228790D+00, & 0.1637017337371825D+00, & 0.0283676653399385D+00, & 0.1637017337371825D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00, & 0.0320553732169435D+00, & 0.0320553732169435D+00, & 0.9358892535661130D+00, & 0.1421611010565643D+00, & 0.1421611010565644D+00, & 0.7156777978868711D+00, & 0.1481328857838206D+00, & 0.1481328857838206D+00, & 0.3218129952888354D+00, & 0.3218129952888354D+00, & 0.5300541189273440D+00, & 0.5300541189273440D+00, & 0.3691467818278109D+00, & 0.3691467818278110D+00, & 0.0296198894887297D+00, & 0.0296198894887297D+00, & 0.6012333286834592D+00, & 0.6012333286834591D+00, & 0.1637017337371824D+00, & 0.1637017337371825D+00, & 0.0283676653399385D+00, & 0.0283676653399385D+00, & 0.8079306009228790D+00, & 0.8079306009228790D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0408716645731430D+00, & 0.0066764844065748D+00, & 0.0066764844065748D+00, & 0.0066764844065748D+00, & 0.0229789818023724D+00, & 0.0229789818023724D+00, & 0.0229789818023724D+00, & 0.0319524531982120D+00, & 0.0319524531982120D+00, & 0.0319524531982120D+00, & 0.0319524531982120D+00, & 0.0319524531982120D+00, & 0.0319524531982120D+00, & 0.0170923240814797D+00, & 0.0170923240814797D+00, & 0.0170923240814797D+00, & 0.0170923240814797D+00, & 0.0170923240814797D+00, & 0.0170923240814797D+00, & 0.0126488788536442D+00, & 0.0126488788536442D+00, & 0.0126488788536442D+00, & 0.0126488788536442D+00, & 0.0126488788536442D+00, & 0.0126488788536442D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule11 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule11() returns the rule of precision 11. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 28 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00, & 0.0284854176143719D+00, & 0.9430291647712562D+00, & 0.0284854176143719D+00, & 0.2102199567031783D+00, & 0.5795600865936434D+00, & 0.2102199567031783D+00, & 0.1026354827122464D+00, & 0.7947290345755071D+00, & 0.1026354827122464D+00, & 0.4958919009658909D+00, & 0.0082161980682182D+00, & 0.4958919009658909D+00, & 0.4384659267643522D+00, & 0.1230681464712956D+00, & 0.4384659267643522D+00, & 0.1493247886520824D+00, & 0.8433497836618531D+00, & 0.0073254276860645D+00, & 0.8433497836618531D+00, & 0.0073254276860645D+00, & 0.1493247886520824D+00, & 0.0460105001654300D+00, & 0.6644083741968642D+00, & 0.2895811256377059D+00, & 0.6644083741968642D+00, & 0.2895811256377059D+00, & 0.0460105001654300D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00, & 0.9430291647712562D+00, & 0.0284854176143719D+00, & 0.0284854176143719D+00, & 0.5795600865936434D+00, & 0.2102199567031783D+00, & 0.2102199567031783D+00, & 0.7947290345755071D+00, & 0.1026354827122464D+00, & 0.1026354827122464D+00, & 0.0082161980682182D+00, & 0.4958919009658909D+00, & 0.4958919009658909D+00, & 0.1230681464712956D+00, & 0.4384659267643522D+00, & 0.4384659267643522D+00, & 0.8433497836618531D+00, & 0.1493247886520824D+00, & 0.8433497836618531D+00, & 0.0073254276860645D+00, & 0.1493247886520824D+00, & 0.0073254276860645D+00, & 0.6644083741968642D+00, & 0.0460105001654300D+00, & 0.6644083741968642D+00, & 0.2895811256377059D+00, & 0.0460105001654300D+00, & 0.2895811256377059D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00, & 0.0284854176143720D+00, & 0.0284854176143719D+00, & 0.9430291647712563D+00, & 0.2102199567031783D+00, & 0.2102199567031783D+00, & 0.5795600865936434D+00, & 0.1026354827122464D+00, & 0.1026354827122464D+00, & 0.7947290345755071D+00, & 0.4958919009658910D+00, & 0.4958919009658909D+00, & 0.0082161980682182D+00, & 0.4384659267643521D+00, & 0.4384659267643522D+00, & 0.1230681464712955D+00, & 0.0073254276860646D+00, & 0.0073254276860645D+00, & 0.1493247886520824D+00, & 0.1493247886520824D+00, & 0.8433497836618531D+00, & 0.8433497836618532D+00, & 0.2895811256377059D+00, & 0.2895811256377059D+00, & 0.0460105001654300D+00, & 0.0460105001654300D+00, & 0.6644083741968642D+00, & 0.6644083741968642D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0428805898661121D+00, & 0.0052159352564473D+00, & 0.0052159352564473D+00, & 0.0052159352564473D+00, & 0.0352578420558583D+00, & 0.0352578420558583D+00, & 0.0352578420558583D+00, & 0.0193153796185097D+00, & 0.0193153796185097D+00, & 0.0193153796185097D+00, & 0.0083031365272927D+00, & 0.0083031365272927D+00, & 0.0083031365272927D+00, & 0.0336580770397341D+00, & 0.0336580770397341D+00, & 0.0336580770397341D+00, & 0.0051451447864766D+00, & 0.0051451447864766D+00, & 0.0051451447864766D+00, & 0.0051451447864766D+00, & 0.0051451447864766D+00, & 0.0051451447864766D+00, & 0.0201662383202503D+00, & 0.0201662383202503D+00, & 0.0201662383202503D+00, & 0.0201662383202503D+00, & 0.0201662383202503D+00, & 0.0201662383202503D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule12 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule12() returns the rule of precision 12. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 33 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.4882037509455415D+00, & 0.0235924981089169D+00, & 0.4882037509455415D+00, & 0.1092578276593543D+00, & 0.7814843446812914D+00, & 0.1092578276593543D+00, & 0.2714625070149261D+00, & 0.4570749859701478D+00, & 0.2714625070149261D+00, & 0.0246463634363356D+00, & 0.9507072731273288D+00, & 0.0246463634363356D+00, & 0.4401116486585931D+00, & 0.1197767026828138D+00, & 0.4401116486585931D+00, & 0.2916556797383409D+00, & 0.6853101639063919D+00, & 0.0230341563552671D+00, & 0.6853101639063919D+00, & 0.0230341563552671D+00, & 0.2916556797383409D+00, & 0.1162960196779266D+00, & 0.6282497516835561D+00, & 0.2554542286385174D+00, & 0.6282497516835561D+00, & 0.2554542286385174D+00, & 0.1162960196779266D+00, & 0.8513377925102401D+00, & 0.1272797172335894D+00, & 0.0213824902561706D+00, & 0.1272797172335894D+00, & 0.0213824902561706D+00, & 0.8513377925102401D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.0235924981089169D+00, & 0.4882037509455415D+00, & 0.4882037509455415D+00, & 0.7814843446812914D+00, & 0.1092578276593543D+00, & 0.1092578276593543D+00, & 0.4570749859701478D+00, & 0.2714625070149261D+00, & 0.2714625070149261D+00, & 0.9507072731273288D+00, & 0.0246463634363356D+00, & 0.0246463634363356D+00, & 0.1197767026828138D+00, & 0.4401116486585931D+00, & 0.4401116486585931D+00, & 0.6853101639063919D+00, & 0.2916556797383409D+00, & 0.6853101639063919D+00, & 0.0230341563552671D+00, & 0.2916556797383409D+00, & 0.0230341563552671D+00, & 0.6282497516835561D+00, & 0.1162960196779266D+00, & 0.6282497516835561D+00, & 0.2554542286385174D+00, & 0.1162960196779266D+00, & 0.2554542286385174D+00, & 0.1272797172335894D+00, & 0.8513377925102401D+00, & 0.1272797172335894D+00, & 0.0213824902561706D+00, & 0.8513377925102401D+00, & 0.0213824902561706D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.4882037509455416D+00, & 0.4882037509455416D+00, & 0.0235924981089170D+00, & 0.1092578276593543D+00, & 0.1092578276593544D+00, & 0.7814843446812914D+00, & 0.2714625070149260D+00, & 0.2714625070149261D+00, & 0.4570749859701478D+00, & 0.0246463634363355D+00, & 0.0246463634363356D+00, & 0.9507072731273287D+00, & 0.4401116486585931D+00, & 0.4401116486585931D+00, & 0.1197767026828138D+00, & 0.0230341563552672D+00, & 0.0230341563552672D+00, & 0.2916556797383409D+00, & 0.2916556797383410D+00, & 0.6853101639063919D+00, & 0.6853101639063919D+00, & 0.2554542286385173D+00, & 0.2554542286385173D+00, & 0.1162960196779266D+00, & 0.1162960196779265D+00, & 0.6282497516835561D+00, & 0.6282497516835561D+00, & 0.0213824902561705D+00, & 0.0213824902561706D+00, & 0.8513377925102401D+00, & 0.8513377925102401D+00, & 0.1272797172335893D+00, & 0.1272797172335893D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0121334190407260D+00, & 0.0121334190407260D+00, & 0.0121334190407260D+00, & 0.0142430260344388D+00, & 0.0142430260344388D+00, & 0.0142430260344388D+00, & 0.0312706065979514D+00, & 0.0312706065979514D+00, & 0.0312706065979514D+00, & 0.0039658212549868D+00, & 0.0039658212549868D+00, & 0.0039658212549868D+00, & 0.0249591674640305D+00, & 0.0249591674640305D+00, & 0.0249591674640305D+00, & 0.0108917925193038D+00, & 0.0108917925193038D+00, & 0.0108917925193038D+00, & 0.0108917925193038D+00, & 0.0108917925193038D+00, & 0.0108917925193038D+00, & 0.0216136818297071D+00, & 0.0216136818297071D+00, & 0.0216136818297071D+00, & 0.0216136818297071D+00, & 0.0216136818297071D+00, & 0.0216136818297071D+00, & 0.0075418387882557D+00, & 0.0075418387882557D+00, & 0.0075418387882557D+00, & 0.0075418387882557D+00, & 0.0075418387882557D+00, & 0.0075418387882557D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule13 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule13() returns the rule of precision 13. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 37 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00, & 0.4890769464525394D+00, & 0.0218461070949213D+00, & 0.4890769464525394D+00, & 0.2213722862918329D+00, & 0.5572554274163342D+00, & 0.2213722862918329D+00, & 0.4269414142598004D+00, & 0.1461171714803992D+00, & 0.4269414142598004D+00, & 0.0215096811088432D+00, & 0.9569806377823136D+00, & 0.0215096811088432D+00, & 0.0878954830321973D+00, & 0.7485071158999522D+00, & 0.1635974010678505D+00, & 0.7485071158999522D+00, & 0.1635974010678505D+00, & 0.0878954830321973D+00, & 0.1109220428034634D+00, & 0.8647077702954428D+00, & 0.0243701869010938D+00, & 0.8647077702954428D+00, & 0.0243701869010938D+00, & 0.1109220428034634D+00, & 0.3084417608921178D+00, & 0.6235459955536755D+00, & 0.0680122435542067D+00, & 0.6235459955536755D+00, & 0.0680122435542067D+00, & 0.3084417608921178D+00, & 0.0051263891023824D+00, & 0.7223577931241880D+00, & 0.2725158177734297D+00, & 0.7223577931241880D+00, & 0.2725158177734297D+00, & 0.0051263891023824D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00, & 0.0218461070949213D+00, & 0.4890769464525394D+00, & 0.4890769464525394D+00, & 0.5572554274163342D+00, & 0.2213722862918329D+00, & 0.2213722862918329D+00, & 0.1461171714803992D+00, & 0.4269414142598004D+00, & 0.4269414142598004D+00, & 0.9569806377823136D+00, & 0.0215096811088432D+00, & 0.0215096811088432D+00, & 0.7485071158999522D+00, & 0.0878954830321973D+00, & 0.7485071158999522D+00, & 0.1635974010678505D+00, & 0.0878954830321973D+00, & 0.1635974010678505D+00, & 0.8647077702954428D+00, & 0.1109220428034634D+00, & 0.8647077702954428D+00, & 0.0243701869010938D+00, & 0.1109220428034634D+00, & 0.0243701869010938D+00, & 0.6235459955536755D+00, & 0.3084417608921178D+00, & 0.6235459955536755D+00, & 0.0680122435542067D+00, & 0.3084417608921178D+00, & 0.0680122435542067D+00, & 0.7223577931241880D+00, & 0.0051263891023824D+00, & 0.7223577931241880D+00, & 0.2725158177734297D+00, & 0.0051263891023824D+00, & 0.2725158177734297D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00, & 0.4890769464525394D+00, & 0.4890769464525394D+00, & 0.0218461070949214D+00, & 0.2213722862918329D+00, & 0.2213722862918329D+00, & 0.5572554274163342D+00, & 0.4269414142598003D+00, & 0.4269414142598003D+00, & 0.1461171714803991D+00, & 0.0215096811088433D+00, & 0.0215096811088433D+00, & 0.9569806377823138D+00, & 0.1635974010678505D+00, & 0.1635974010678505D+00, & 0.0878954830321973D+00, & 0.0878954830321973D+00, & 0.7485071158999522D+00, & 0.7485071158999522D+00, & 0.0243701869010938D+00, & 0.0243701869010938D+00, & 0.1109220428034634D+00, & 0.1109220428034634D+00, & 0.8647077702954428D+00, & 0.8647077702954428D+00, & 0.0680122435542068D+00, & 0.0680122435542067D+00, & 0.3084417608921178D+00, & 0.3084417608921178D+00, & 0.6235459955536755D+00, & 0.6235459955536756D+00, & 0.2725158177734296D+00, & 0.2725158177734296D+00, & 0.0051263891023823D+00, & 0.0051263891023823D+00, & 0.7223577931241879D+00, & 0.7223577931241879D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0339800182934158D+00, & 0.0119972009644474D+00, & 0.0119972009644474D+00, & 0.0119972009644474D+00, & 0.0291392425596000D+00, & 0.0291392425596000D+00, & 0.0291392425596000D+00, & 0.0278009837652267D+00, & 0.0278009837652267D+00, & 0.0278009837652267D+00, & 0.0030261685517696D+00, & 0.0030261685517696D+00, & 0.0030261685517696D+00, & 0.0120895199057969D+00, & 0.0120895199057969D+00, & 0.0120895199057969D+00, & 0.0120895199057969D+00, & 0.0120895199057969D+00, & 0.0120895199057969D+00, & 0.0074827005525828D+00, & 0.0074827005525828D+00, & 0.0074827005525828D+00, & 0.0074827005525828D+00, & 0.0074827005525828D+00, & 0.0074827005525828D+00, & 0.0173206380704242D+00, & 0.0173206380704242D+00, & 0.0173206380704242D+00, & 0.0173206380704242D+00, & 0.0173206380704242D+00, & 0.0173206380704242D+00, & 0.0047953405017716D+00, & 0.0047953405017716D+00, & 0.0047953405017716D+00, & 0.0047953405017716D+00, & 0.0047953405017716D+00, & 0.0047953405017716D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule14 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule14() returns the rule of precision 14. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 42 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.1772055324125434D+00, & 0.6455889351749131D+00, & 0.1772055324125434D+00, & 0.4176447193404539D+00, & 0.1647105613190922D+00, & 0.4176447193404539D+00, & 0.0617998830908726D+00, & 0.8764002338182548D+00, & 0.0617998830908726D+00, & 0.4889639103621786D+00, & 0.0220721792756427D+00, & 0.4889639103621786D+00, & 0.2734775283088386D+00, & 0.4530449433823227D+00, & 0.2734775283088386D+00, & 0.0193909612487010D+00, & 0.9612180775025979D+00, & 0.0193909612487010D+00, & 0.2983728821362578D+00, & 0.6869801678080878D+00, & 0.0146469500556544D+00, & 0.6869801678080878D+00, & 0.0146469500556544D+00, & 0.2983728821362578D+00, & 0.0571247574036479D+00, & 0.7706085547749965D+00, & 0.1722666878213556D+00, & 0.7706085547749965D+00, & 0.1722666878213556D+00, & 0.0571247574036479D+00, & 0.3368614597963450D+00, & 0.5702222908466832D+00, & 0.0929162493569718D+00, & 0.5702222908466832D+00, & 0.0929162493569718D+00, & 0.3368614597963450D+00, & 0.0012683309328720D+00, & 0.8797571713701711D+00, & 0.1189744976969568D+00, & 0.8797571713701711D+00, & 0.1189744976969568D+00, & 0.0012683309328720D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.6455889351749131D+00, & 0.1772055324125434D+00, & 0.1772055324125434D+00, & 0.1647105613190922D+00, & 0.4176447193404539D+00, & 0.4176447193404539D+00, & 0.8764002338182548D+00, & 0.0617998830908726D+00, & 0.0617998830908726D+00, & 0.0220721792756427D+00, & 0.4889639103621786D+00, & 0.4889639103621786D+00, & 0.4530449433823227D+00, & 0.2734775283088386D+00, & 0.2734775283088386D+00, & 0.9612180775025979D+00, & 0.0193909612487010D+00, & 0.0193909612487010D+00, & 0.6869801678080878D+00, & 0.2983728821362578D+00, & 0.6869801678080878D+00, & 0.0146469500556544D+00, & 0.2983728821362578D+00, & 0.0146469500556544D+00, & 0.7706085547749965D+00, & 0.0571247574036479D+00, & 0.7706085547749965D+00, & 0.1722666878213556D+00, & 0.0571247574036479D+00, & 0.1722666878213556D+00, & 0.5702222908466832D+00, & 0.3368614597963450D+00, & 0.5702222908466832D+00, & 0.0929162493569718D+00, & 0.3368614597963450D+00, & 0.0929162493569718D+00, & 0.8797571713701711D+00, & 0.0012683309328720D+00, & 0.8797571713701711D+00, & 0.1189744976969568D+00, & 0.0012683309328720D+00, & 0.1189744976969568D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.1772055324125434D+00, & 0.1772055324125434D+00, & 0.6455889351749131D+00, & 0.4176447193404538D+00, & 0.4176447193404538D+00, & 0.1647105613190921D+00, & 0.0617998830908726D+00, & 0.0617998830908726D+00, & 0.8764002338182548D+00, & 0.4889639103621787D+00, & 0.4889639103621787D+00, & 0.0220721792756428D+00, & 0.2734775283088386D+00, & 0.2734775283088386D+00, & 0.4530449433823227D+00, & 0.0193909612487010D+00, & 0.0193909612487010D+00, & 0.9612180775025978D+00, & 0.0146469500556544D+00, & 0.0146469500556544D+00, & 0.2983728821362578D+00, & 0.2983728821362578D+00, & 0.6869801678080878D+00, & 0.6869801678080878D+00, & 0.1722666878213556D+00, & 0.1722666878213556D+00, & 0.0571247574036480D+00, & 0.0571247574036480D+00, & 0.7706085547749966D+00, & 0.7706085547749965D+00, & 0.0929162493569718D+00, & 0.0929162493569718D+00, & 0.3368614597963450D+00, & 0.3368614597963450D+00, & 0.5702222908466832D+00, & 0.5702222908466832D+00, & 0.1189744976969569D+00, & 0.1189744976969569D+00, & 0.0012683309328720D+00, & 0.0012683309328721D+00, & 0.8797571713701711D+00, & 0.8797571713701711D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0210812943684965D+00, & 0.0210812943684965D+00, & 0.0210812943684965D+00, & 0.0163941767720627D+00, & 0.0163941767720627D+00, & 0.0163941767720627D+00, & 0.0072168498348883D+00, & 0.0072168498348883D+00, & 0.0072168498348883D+00, & 0.0109417906847144D+00, & 0.0109417906847144D+00, & 0.0109417906847144D+00, & 0.0258870522536458D+00, & 0.0258870522536458D+00, & 0.0258870522536458D+00, & 0.0024617018012000D+00, & 0.0024617018012000D+00, & 0.0024617018012000D+00, & 0.0072181540567669D+00, & 0.0072181540567669D+00, & 0.0072181540567669D+00, & 0.0072181540567669D+00, & 0.0072181540567669D+00, & 0.0072181540567669D+00, & 0.0123328766062818D+00, & 0.0123328766062818D+00, & 0.0123328766062818D+00, & 0.0123328766062818D+00, & 0.0123328766062818D+00, & 0.0123328766062818D+00, & 0.0192857553935303D+00, & 0.0192857553935303D+00, & 0.0192857553935303D+00, & 0.0192857553935303D+00, & 0.0192857553935303D+00, & 0.0192857553935303D+00, & 0.0025051144192503D+00, & 0.0025051144192503D+00, & 0.0025051144192503D+00, & 0.0025051144192503D+00, & 0.0025051144192503D+00, & 0.0025051144192503D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule15 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule15() returns the rule of precision 15. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! We suppose we are given a triangle T with vertices A, B, C. ! We call a rule with n points, returning barycentric coordinates ! a, b, c, and weights w. Then the integral I of f(x,y) over T is ! approximated by Q as follows: ! ! (x,y) = a(1:n) * A + b(1:n) * B + c(1:n) * C ! Q = area(T) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 49 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00, & 0.4053622141339755D+00, & 0.1892755717320491D+00, & 0.4053622141339755D+00, & 0.0701735528999860D+00, & 0.8596528942000279D+00, & 0.0701735528999860D+00, & 0.4741706814380198D+00, & 0.0516586371239604D+00, & 0.4741706814380198D+00, & 0.2263787134203496D+00, & 0.5472425731593008D+00, & 0.2263787134203496D+00, & 0.4949969567691262D+00, & 0.0100060864617476D+00, & 0.4949969567691262D+00, & 0.0158117262509886D+00, & 0.9683765474980227D+00, & 0.0158117262509886D+00, & 0.0183761123856811D+00, & 0.6669756448018681D+00, & 0.3146482428124509D+00, & 0.6669756448018681D+00, & 0.3146482428124509D+00, & 0.0183761123856811D+00, & 0.0091392370373084D+00, & 0.9199121577262361D+00, & 0.0709486052364555D+00, & 0.9199121577262361D+00, & 0.0709486052364555D+00, & 0.0091392370373084D+00, & 0.1905355894763939D+00, & 0.7152223569314506D+00, & 0.0942420535921554D+00, & 0.7152223569314506D+00, & 0.0942420535921554D+00, & 0.1905355894763939D+00, & 0.1680686452224144D+00, & 0.8132926410494192D+00, & 0.0186387137281664D+00, & 0.8132926410494192D+00, & 0.0186387137281664D+00, & 0.1680686452224144D+00, & 0.3389506114752772D+00, & 0.5652526648771142D+00, & 0.0957967236476086D+00, & 0.5652526648771142D+00, & 0.0957967236476086D+00, & 0.3389506114752772D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00, & 0.1892755717320491D+00, & 0.4053622141339755D+00, & 0.4053622141339755D+00, & 0.8596528942000279D+00, & 0.0701735528999860D+00, & 0.0701735528999860D+00, & 0.0516586371239604D+00, & 0.4741706814380198D+00, & 0.4741706814380198D+00, & 0.5472425731593008D+00, & 0.2263787134203496D+00, & 0.2263787134203496D+00, & 0.0100060864617476D+00, & 0.4949969567691262D+00, & 0.4949969567691262D+00, & 0.9683765474980227D+00, & 0.0158117262509886D+00, & 0.0158117262509886D+00, & 0.6669756448018681D+00, & 0.0183761123856811D+00, & 0.6669756448018681D+00, & 0.3146482428124509D+00, & 0.0183761123856811D+00, & 0.3146482428124509D+00, & 0.9199121577262361D+00, & 0.0091392370373084D+00, & 0.9199121577262361D+00, & 0.0709486052364555D+00, & 0.0091392370373084D+00, & 0.0709486052364555D+00, & 0.7152223569314506D+00, & 0.1905355894763939D+00, & 0.7152223569314506D+00, & 0.0942420535921554D+00, & 0.1905355894763939D+00, & 0.0942420535921554D+00, & 0.8132926410494192D+00, & 0.1680686452224144D+00, & 0.8132926410494192D+00, & 0.0186387137281664D+00, & 0.1680686452224144D+00, & 0.0186387137281664D+00, & 0.5652526648771142D+00, & 0.3389506114752772D+00, & 0.5652526648771142D+00, & 0.0957967236476086D+00, & 0.3389506114752772D+00, & 0.0957967236476086D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00, & 0.4053622141339755D+00, & 0.4053622141339754D+00, & 0.1892755717320491D+00, & 0.0701735528999861D+00, & 0.0701735528999861D+00, & 0.8596528942000279D+00, & 0.4741706814380198D+00, & 0.4741706814380198D+00, & 0.0516586371239605D+00, & 0.2263787134203495D+00, & 0.2263787134203496D+00, & 0.5472425731593007D+00, & 0.4949969567691263D+00, & 0.4949969567691262D+00, & 0.0100060864617477D+00, & 0.0158117262509886D+00, & 0.0158117262509886D+00, & 0.9683765474980226D+00, & 0.3146482428124507D+00, & 0.3146482428124508D+00, & 0.0183761123856810D+00, & 0.0183761123856810D+00, & 0.6669756448018680D+00, & 0.6669756448018680D+00, & 0.0709486052364554D+00, & 0.0709486052364555D+00, & 0.0091392370373083D+00, & 0.0091392370373083D+00, & 0.9199121577262361D+00, & 0.9199121577262360D+00, & 0.0942420535921554D+00, & 0.0942420535921554D+00, & 0.1905355894763940D+00, & 0.1905355894763940D+00, & 0.7152223569314508D+00, & 0.7152223569314506D+00, & 0.0186387137281664D+00, & 0.0186387137281664D+00, & 0.1680686452224144D+00, & 0.1680686452224144D+00, & 0.8132926410494192D+00, & 0.8132926410494192D+00, & 0.0957967236476086D+00, & 0.0957967236476086D+00, & 0.3389506114752772D+00, & 0.3389506114752772D+00, & 0.5652526648771142D+00, & 0.5652526648771142D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0443353873821841D+00, & 0.0427137815714606D+00, & 0.0427137815714606D+00, & 0.0427137815714606D+00, & 0.0164447375626252D+00, & 0.0164447375626252D+00, & 0.0164447375626252D+00, & 0.0173961480007634D+00, & 0.0173961480007634D+00, & 0.0173961480007634D+00, & 0.0467833617287096D+00, & 0.0467833617287096D+00, & 0.0467833617287096D+00, & 0.0095738461824601D+00, & 0.0095738461824601D+00, & 0.0095738461824601D+00, & 0.0029607746379054D+00, & 0.0029607746379054D+00, & 0.0029607746379054D+00, & 0.0156025728305760D+00, & 0.0156025728305760D+00, & 0.0156025728305760D+00, & 0.0156025728305760D+00, & 0.0156025728305760D+00, & 0.0156025728305760D+00, & 0.0040298533720181D+00, & 0.0040298533720181D+00, & 0.0040298533720181D+00, & 0.0040298533720181D+00, & 0.0040298533720181D+00, & 0.0040298533720181D+00, & 0.0287205869252013D+00, & 0.0287205869252013D+00, & 0.0287205869252013D+00, & 0.0287205869252013D+00, & 0.0287205869252013D+00, & 0.0287205869252013D+00, & 0.0116726211815758D+00, & 0.0116726211815758D+00, & 0.0116726211815758D+00, & 0.0116726211815758D+00, & 0.0116726211815758D+00, & 0.0116726211815758D+00, & 0.0313154762849693D+00, & 0.0313154762849693D+00, & 0.0313154762849693D+00, & 0.0313154762849693D+00, & 0.0313154762849693D+00, & 0.0313154762849693D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights already sum to 1. ! return end subroutine rule16 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule16() returns the rule of precision 16. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 55 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00, & 0.2459900704671417D+00, & 0.5080198590657166D+00, & 0.2459900704671417D+00, & 0.4155848968854205D+00, & 0.1688302062291590D+00, & 0.4155848968854205D+00, & 0.0853555665867003D+00, & 0.8292888668265994D+00, & 0.0853555665867003D+00, & 0.1619186441912712D+00, & 0.6761627116174576D+00, & 0.1619186441912712D+00, & 0.5000000000000000D+00, & 0.0000000000000000D+00, & 0.5000000000000000D+00, & 0.4752807275459421D+00, & 0.0494385449081158D+00, & 0.4752807275459421D+00, & 0.0547551749147031D+00, & 0.7541700614447677D+00, & 0.1910747636405291D+00, & 0.7541700614447677D+00, & 0.1910747636405291D+00, & 0.0547551749147031D+00, & 0.0232034277688137D+00, & 0.9682443680309587D+00, & 0.0085522042002276D+00, & 0.9682443680309587D+00, & 0.0085522042002276D+00, & 0.0232034277688137D+00, & 0.0189317782804059D+00, & 0.6493036982454464D+00, & 0.3317645234741476D+00, & 0.6493036982454464D+00, & 0.3317645234741476D+00, & 0.0189317782804059D+00, & 0.0190301297436974D+00, & 0.9002737032704295D+00, & 0.0806961669858730D+00, & 0.9002737032704295D+00, & 0.0806961669858730D+00, & 0.0190301297436974D+00, & 0.1026061902393981D+00, & 0.5891488405642479D+00, & 0.3082449691963540D+00, & 0.5891488405642479D+00, & 0.3082449691963540D+00, & 0.1026061902393981D+00, & 0.0059363500168222D+00, & 0.8066218674993957D+00, & 0.1874417824837821D+00, & 0.8066218674993957D+00, & 0.1874417824837821D+00, & 0.0059363500168222D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00, & 0.5080198590657166D+00, & 0.2459900704671417D+00, & 0.2459900704671417D+00, & 0.1688302062291590D+00, & 0.4155848968854205D+00, & 0.4155848968854205D+00, & 0.8292888668265994D+00, & 0.0853555665867003D+00, & 0.0853555665867003D+00, & 0.6761627116174576D+00, & 0.1619186441912712D+00, & 0.1619186441912712D+00, & 0.0000000000000000D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.0494385449081158D+00, & 0.4752807275459421D+00, & 0.4752807275459421D+00, & 0.7541700614447677D+00, & 0.0547551749147031D+00, & 0.7541700614447677D+00, & 0.1910747636405291D+00, & 0.0547551749147031D+00, & 0.1910747636405291D+00, & 0.9682443680309587D+00, & 0.0232034277688137D+00, & 0.9682443680309587D+00, & 0.0085522042002276D+00, & 0.0232034277688137D+00, & 0.0085522042002276D+00, & 0.6493036982454464D+00, & 0.0189317782804059D+00, & 0.6493036982454464D+00, & 0.3317645234741476D+00, & 0.0189317782804059D+00, & 0.3317645234741476D+00, & 0.9002737032704295D+00, & 0.0190301297436974D+00, & 0.9002737032704295D+00, & 0.0806961669858730D+00, & 0.0190301297436974D+00, & 0.0806961669858730D+00, & 0.5891488405642479D+00, & 0.1026061902393981D+00, & 0.5891488405642479D+00, & 0.3082449691963540D+00, & 0.1026061902393981D+00, & 0.3082449691963540D+00, & 0.8066218674993957D+00, & 0.0059363500168222D+00, & 0.8066218674993957D+00, & 0.1874417824837821D+00, & 0.0059363500168222D+00, & 0.1874417824837821D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00, & 0.2459900704671417D+00, & 0.2459900704671417D+00, & 0.5080198590657166D+00, & 0.4155848968854205D+00, & 0.4155848968854206D+00, & 0.1688302062291590D+00, & 0.0853555665867003D+00, & 0.0853555665867002D+00, & 0.8292888668265994D+00, & 0.1619186441912712D+00, & 0.1619186441912712D+00, & 0.6761627116174576D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.0000000000000000D+00, & 0.4752807275459421D+00, & 0.4752807275459421D+00, & 0.0494385449081158D+00, & 0.1910747636405292D+00, & 0.1910747636405292D+00, & 0.0547551749147033D+00, & 0.0547551749147032D+00, & 0.7541700614447679D+00, & 0.7541700614447679D+00, & 0.0085522042002275D+00, & 0.0085522042002276D+00, & 0.0232034277688137D+00, & 0.0232034277688137D+00, & 0.9682443680309587D+00, & 0.9682443680309586D+00, & 0.3317645234741476D+00, & 0.3317645234741476D+00, & 0.0189317782804059D+00, & 0.0189317782804059D+00, & 0.6493036982454464D+00, & 0.6493036982454464D+00, & 0.0806961669858730D+00, & 0.0806961669858730D+00, & 0.0190301297436974D+00, & 0.0190301297436975D+00, & 0.9002737032704295D+00, & 0.9002737032704295D+00, & 0.3082449691963540D+00, & 0.3082449691963540D+00, & 0.1026061902393981D+00, & 0.1026061902393981D+00, & 0.5891488405642479D+00, & 0.5891488405642479D+00, & 0.1874417824837821D+00, & 0.1874417824837821D+00, & 0.0059363500168222D+00, & 0.0059363500168222D+00, & 0.8066218674993957D+00, & 0.8066218674993957D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0226322830369094D+00, & 0.0205464615718495D+00, & 0.0205464615718495D+00, & 0.0205464615718495D+00, & 0.0203559166562127D+00, & 0.0203559166562127D+00, & 0.0203559166562127D+00, & 0.0073908173451122D+00, & 0.0073908173451122D+00, & 0.0073908173451122D+00, & 0.0147092048494940D+00, & 0.0147092048494940D+00, & 0.0147092048494940D+00, & 0.0022092731560753D+00, & 0.0022092731560753D+00, & 0.0022092731560753D+00, & 0.0129871666491386D+00, & 0.0129871666491386D+00, & 0.0129871666491386D+00, & 0.0094691362322078D+00, & 0.0094691362322078D+00, & 0.0094691362322078D+00, & 0.0094691362322078D+00, & 0.0094691362322078D+00, & 0.0094691362322078D+00, & 0.0008272333574175D+00, & 0.0008272333574175D+00, & 0.0008272333574175D+00, & 0.0008272333574175D+00, & 0.0008272333574175D+00, & 0.0008272333574175D+00, & 0.0075043008921429D+00, & 0.0075043008921429D+00, & 0.0075043008921429D+00, & 0.0075043008921429D+00, & 0.0075043008921429D+00, & 0.0075043008921429D+00, & 0.0039737969666962D+00, & 0.0039737969666962D+00, & 0.0039737969666962D+00, & 0.0039737969666962D+00, & 0.0039737969666962D+00, & 0.0039737969666962D+00, & 0.0159918050396850D+00, & 0.0159918050396850D+00, & 0.0159918050396850D+00, & 0.0159918050396850D+00, & 0.0159918050396850D+00, & 0.0159918050396850D+00, & 0.0026955935584244D+00, & 0.0026955935584244D+00, & 0.0026955935584244D+00, & 0.0026955935584244D+00, & 0.0026955935584244D+00, & 0.0026955935584244D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule17 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule17() returns the rule of precision 17. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 60 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.4171034443615992D+00, & 0.1657931112768016D+00, & 0.4171034443615992D+00, & 0.0147554916607540D+00, & 0.9704890166784921D+00, & 0.0147554916607540D+00, & 0.4655978716188903D+00, & 0.0688042567622194D+00, & 0.4655978716188903D+00, & 0.1803581162663706D+00, & 0.6392837674672588D+00, & 0.1803581162663706D+00, & 0.0666540634795969D+00, & 0.8666918730408062D+00, & 0.0666540634795969D+00, & 0.2857065024365866D+00, & 0.4285869951268267D+00, & 0.2857065024365866D+00, & 0.0160176423621193D+00, & 0.8247900701650881D+00, & 0.1591922874727927D+00, & 0.8247900701650881D+00, & 0.1591922874727927D+00, & 0.0160176423621193D+00, & 0.3062815917461865D+00, & 0.6263690303864522D+00, & 0.0673493778673612D+00, & 0.6263690303864522D+00, & 0.0673493778673612D+00, & 0.3062815917461865D+00, & 0.0132296727600869D+00, & 0.5712948679446841D+00, & 0.4154754592952291D+00, & 0.5712948679446841D+00, & 0.4154754592952291D+00, & 0.0132296727600869D+00, & 0.0780423405682824D+00, & 0.7532351459364581D+00, & 0.1687225134952595D+00, & 0.7532351459364581D+00, & 0.1687225134952595D+00, & 0.0780423405682824D+00, & 0.0131358708340027D+00, & 0.7150722591106424D+00, & 0.2717918700553549D+00, & 0.7150722591106424D+00, & 0.2717918700553549D+00, & 0.0131358708340027D+00, & 0.0115751759031806D+00, & 0.9159193532978169D+00, & 0.0725054707990024D+00, & 0.9159193532978169D+00, & 0.0725054707990024D+00, & 0.0115751759031806D+00, & 0.1575054779268699D+00, & 0.5432755795961598D+00, & 0.2992189424769703D+00, & 0.5432755795961598D+00, & 0.2992189424769703D+00, & 0.1575054779268699D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.1657931112768016D+00, & 0.4171034443615992D+00, & 0.4171034443615992D+00, & 0.9704890166784921D+00, & 0.0147554916607540D+00, & 0.0147554916607540D+00, & 0.0688042567622194D+00, & 0.4655978716188903D+00, & 0.4655978716188903D+00, & 0.6392837674672588D+00, & 0.1803581162663706D+00, & 0.1803581162663706D+00, & 0.8666918730408062D+00, & 0.0666540634795969D+00, & 0.0666540634795969D+00, & 0.4285869951268267D+00, & 0.2857065024365866D+00, & 0.2857065024365866D+00, & 0.8247900701650881D+00, & 0.0160176423621193D+00, & 0.8247900701650881D+00, & 0.1591922874727927D+00, & 0.0160176423621193D+00, & 0.1591922874727927D+00, & 0.6263690303864522D+00, & 0.3062815917461865D+00, & 0.6263690303864522D+00, & 0.0673493778673612D+00, & 0.3062815917461865D+00, & 0.0673493778673612D+00, & 0.5712948679446841D+00, & 0.0132296727600869D+00, & 0.5712948679446841D+00, & 0.4154754592952291D+00, & 0.0132296727600869D+00, & 0.4154754592952291D+00, & 0.7532351459364581D+00, & 0.0780423405682824D+00, & 0.7532351459364581D+00, & 0.1687225134952595D+00, & 0.0780423405682824D+00, & 0.1687225134952595D+00, & 0.7150722591106424D+00, & 0.0131358708340027D+00, & 0.7150722591106424D+00, & 0.2717918700553549D+00, & 0.0131358708340027D+00, & 0.2717918700553549D+00, & 0.9159193532978169D+00, & 0.0115751759031806D+00, & 0.9159193532978169D+00, & 0.0725054707990024D+00, & 0.0115751759031806D+00, & 0.0725054707990024D+00, & 0.5432755795961598D+00, & 0.1575054779268699D+00, & 0.5432755795961598D+00, & 0.2992189424769703D+00, & 0.1575054779268699D+00, & 0.2992189424769703D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.4171034443615992D+00, & 0.4171034443615992D+00, & 0.1657931112768016D+00, & 0.0147554916607540D+00, & 0.0147554916607540D+00, & 0.9704890166784921D+00, & 0.4655978716188903D+00, & 0.4655978716188903D+00, & 0.0688042567622194D+00, & 0.1803581162663707D+00, & 0.1803581162663706D+00, & 0.6392837674672589D+00, & 0.0666540634795969D+00, & 0.0666540634795969D+00, & 0.8666918730408062D+00, & 0.2857065024365867D+00, & 0.2857065024365866D+00, & 0.4285869951268268D+00, & 0.1591922874727927D+00, & 0.1591922874727926D+00, & 0.0160176423621192D+00, & 0.0160176423621192D+00, & 0.8247900701650881D+00, & 0.8247900701650881D+00, & 0.0673493778673613D+00, & 0.0673493778673613D+00, & 0.3062815917461866D+00, & 0.3062815917461866D+00, & 0.6263690303864523D+00, & 0.6263690303864523D+00, & 0.4154754592952290D+00, & 0.4154754592952290D+00, & 0.0132296727600869D+00, & 0.0132296727600869D+00, & 0.5712948679446841D+00, & 0.5712948679446841D+00, & 0.1687225134952595D+00, & 0.1687225134952595D+00, & 0.0780423405682824D+00, & 0.0780423405682824D+00, & 0.7532351459364581D+00, & 0.7532351459364581D+00, & 0.2717918700553549D+00, & 0.2717918700553549D+00, & 0.0131358708340027D+00, & 0.0131358708340027D+00, & 0.7150722591106424D+00, & 0.7150722591106424D+00, & 0.0725054707990025D+00, & 0.0725054707990024D+00, & 0.0115751759031806D+00, & 0.0115751759031806D+00, & 0.9159193532978169D+00, & 0.9159193532978169D+00, & 0.2992189424769703D+00, & 0.2992189424769703D+00, & 0.1575054779268699D+00, & 0.1575054779268699D+00, & 0.5432755795961597D+00, & 0.5432755795961598D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0136554632640511D+00, & 0.0136554632640511D+00, & 0.0136554632640511D+00, & 0.0013869437888188D+00, & 0.0013869437888188D+00, & 0.0013869437888188D+00, & 0.0125097254752487D+00, & 0.0125097254752487D+00, & 0.0125097254752487D+00, & 0.0131563152940090D+00, & 0.0131563152940090D+00, & 0.0131563152940090D+00, & 0.0062295004011527D+00, & 0.0062295004011527D+00, & 0.0062295004011527D+00, & 0.0188581185763976D+00, & 0.0188581185763976D+00, & 0.0188581185763976D+00, & 0.0039891501029648D+00, & 0.0039891501029648D+00, & 0.0039891501029648D+00, & 0.0039891501029648D+00, & 0.0039891501029648D+00, & 0.0039891501029648D+00, & 0.0112438862733455D+00, & 0.0112438862733455D+00, & 0.0112438862733455D+00, & 0.0112438862733455D+00, & 0.0112438862733455D+00, & 0.0112438862733455D+00, & 0.0051992199779198D+00, & 0.0051992199779198D+00, & 0.0051992199779198D+00, & 0.0051992199779198D+00, & 0.0051992199779198D+00, & 0.0051992199779198D+00, & 0.0102789491602273D+00, & 0.0102789491602273D+00, & 0.0102789491602273D+00, & 0.0102789491602273D+00, & 0.0102789491602273D+00, & 0.0102789491602273D+00, & 0.0043461072505006D+00, & 0.0043461072505006D+00, & 0.0043461072505006D+00, & 0.0043461072505006D+00, & 0.0043461072505006D+00, & 0.0043461072505006D+00, & 0.0022921742008679D+00, & 0.0022921742008679D+00, & 0.0022921742008679D+00, & 0.0022921742008679D+00, & 0.0022921742008679D+00, & 0.0022921742008679D+00, & 0.0130858129676685D+00, & 0.0130858129676685D+00, & 0.0130858129676685D+00, & 0.0130858129676685D+00, & 0.0130858129676685D+00, & 0.0130858129676685D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule18 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule18() returns the rule of precision 18. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 67 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00, & 0.3999556280675762D+00, & 0.2000887438648475D+00, & 0.3999556280675762D+00, & 0.4875803015748695D+00, & 0.0248393968502609D+00, & 0.4875803015748695D+00, & 0.4618095064064492D+00, & 0.0763809871871015D+00, & 0.4618095064064492D+00, & 0.2422647025142720D+00, & 0.5154705949714561D+00, & 0.2422647025142720D+00, & 0.0388302560886856D+00, & 0.9223394878226288D+00, & 0.0388302560886856D+00, & 0.0919477421216432D+00, & 0.8161045157567136D+00, & 0.0919477421216432D+00, & 0.0458049158598608D+00, & 0.7703723762146752D+00, & 0.1838227079254640D+00, & 0.7703723762146752D+00, & 0.1838227079254640D+00, & 0.0458049158598608D+00, & 0.2063492574338379D+00, & 0.6709539851942345D+00, & 0.1226967573719275D+00, & 0.6709539851942345D+00, & 0.1226967573719275D+00, & 0.2063492574338379D+00, & 0.0038976110334734D+00, & 0.6004189546342569D+00, & 0.3956834343322697D+00, & 0.6004189546342569D+00, & 0.3956834343322697D+00, & 0.0038976110334734D+00, & 0.0134620167414450D+00, & 0.8783421894675217D+00, & 0.1081957937910333D+00, & 0.8783421894675217D+00, & 0.1081957937910333D+00, & 0.0134620167414450D+00, & 0.0402602834699081D+00, & 0.6399880920047146D+00, & 0.3197516245253773D+00, & 0.6399880920047146D+00, & 0.3197516245253773D+00, & 0.0402602834699081D+00, & 0.0052983351866098D+00, & 0.7589294798551984D+00, & 0.2357721849581917D+00, & 0.7589294798551984D+00, & 0.2357721849581917D+00, & 0.0052983351866098D+00, & 0.0005483600420423D+00, & 0.9723607289627957D+00, & 0.0270909109951620D+00, & 0.9723607289627957D+00, & 0.0270909109951620D+00, & 0.0005483600420423D+00, & 0.1205876951639246D+00, & 0.5459187753861946D+00, & 0.3334935294498808D+00, & 0.5459187753861946D+00, & 0.3334935294498808D+00, & 0.1205876951639246D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00, & 0.2000887438648475D+00, & 0.3999556280675762D+00, & 0.3999556280675762D+00, & 0.0248393968502609D+00, & 0.4875803015748695D+00, & 0.4875803015748695D+00, & 0.0763809871871015D+00, & 0.4618095064064492D+00, & 0.4618095064064492D+00, & 0.5154705949714561D+00, & 0.2422647025142720D+00, & 0.2422647025142720D+00, & 0.9223394878226288D+00, & 0.0388302560886856D+00, & 0.0388302560886856D+00, & 0.8161045157567136D+00, & 0.0919477421216432D+00, & 0.0919477421216432D+00, & 0.7703723762146752D+00, & 0.0458049158598608D+00, & 0.7703723762146752D+00, & 0.1838227079254640D+00, & 0.0458049158598608D+00, & 0.1838227079254640D+00, & 0.6709539851942345D+00, & 0.2063492574338379D+00, & 0.6709539851942345D+00, & 0.1226967573719275D+00, & 0.2063492574338379D+00, & 0.1226967573719275D+00, & 0.6004189546342569D+00, & 0.0038976110334734D+00, & 0.6004189546342569D+00, & 0.3956834343322697D+00, & 0.0038976110334734D+00, & 0.3956834343322697D+00, & 0.8783421894675217D+00, & 0.0134620167414450D+00, & 0.8783421894675217D+00, & 0.1081957937910333D+00, & 0.0134620167414450D+00, & 0.1081957937910333D+00, & 0.6399880920047146D+00, & 0.0402602834699081D+00, & 0.6399880920047146D+00, & 0.3197516245253773D+00, & 0.0402602834699081D+00, & 0.3197516245253773D+00, & 0.7589294798551984D+00, & 0.0052983351866098D+00, & 0.7589294798551984D+00, & 0.2357721849581917D+00, & 0.0052983351866098D+00, & 0.2357721849581917D+00, & 0.9723607289627957D+00, & 0.0005483600420423D+00, & 0.9723607289627957D+00, & 0.0270909109951620D+00, & 0.0005483600420423D+00, & 0.0270909109951620D+00, & 0.5459187753861946D+00, & 0.1205876951639246D+00, & 0.5459187753861946D+00, & 0.3334935294498808D+00, & 0.1205876951639246D+00, & 0.3334935294498808D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00, & 0.3999556280675762D+00, & 0.3999556280675762D+00, & 0.2000887438648475D+00, & 0.4875803015748696D+00, & 0.4875803015748695D+00, & 0.0248393968502609D+00, & 0.4618095064064492D+00, & 0.4618095064064492D+00, & 0.0763809871871015D+00, & 0.2422647025142720D+00, & 0.2422647025142720D+00, & 0.5154705949714562D+00, & 0.0388302560886856D+00, & 0.0388302560886856D+00, & 0.9223394878226288D+00, & 0.0919477421216431D+00, & 0.0919477421216432D+00, & 0.8161045157567135D+00, & 0.1838227079254640D+00, & 0.1838227079254640D+00, & 0.0458049158598608D+00, & 0.0458049158598607D+00, & 0.7703723762146752D+00, & 0.7703723762146752D+00, & 0.1226967573719275D+00, & 0.1226967573719275D+00, & 0.2063492574338379D+00, & 0.2063492574338379D+00, & 0.6709539851942345D+00, & 0.6709539851942345D+00, & 0.3956834343322697D+00, & 0.3956834343322698D+00, & 0.0038976110334734D+00, & 0.0038976110334734D+00, & 0.6004189546342569D+00, & 0.6004189546342569D+00, & 0.1081957937910333D+00, & 0.1081957937910333D+00, & 0.0134620167414450D+00, & 0.0134620167414450D+00, & 0.8783421894675217D+00, & 0.8783421894675217D+00, & 0.3197516245253773D+00, & 0.3197516245253773D+00, & 0.0402602834699081D+00, & 0.0402602834699081D+00, & 0.6399880920047146D+00, & 0.6399880920047146D+00, & 0.2357721849581917D+00, & 0.2357721849581918D+00, & 0.0052983351866098D+00, & 0.0052983351866098D+00, & 0.7589294798551984D+00, & 0.7589294798551984D+00, & 0.0270909109951619D+00, & 0.0270909109951620D+00, & 0.0005483600420423D+00, & 0.0005483600420423D+00, & 0.9723607289627956D+00, & 0.9723607289627956D+00, & 0.3334935294498808D+00, & 0.3334935294498808D+00, & 0.1205876951639246D+00, & 0.1205876951639246D+00, & 0.5459187753861946D+00, & 0.5459187753861946D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0181778676507133D+00, & 0.0166522350166951D+00, & 0.0166522350166951D+00, & 0.0166522350166951D+00, & 0.0060233238169999D+00, & 0.0060233238169999D+00, & 0.0060233238169999D+00, & 0.0094745857533894D+00, & 0.0094745857533894D+00, & 0.0094745857533894D+00, & 0.0182375447044718D+00, & 0.0182375447044718D+00, & 0.0182375447044718D+00, & 0.0035646630098595D+00, & 0.0035646630098595D+00, & 0.0035646630098595D+00, & 0.0082795799760016D+00, & 0.0082795799760016D+00, & 0.0082795799760016D+00, & 0.0068798081174711D+00, & 0.0068798081174711D+00, & 0.0068798081174711D+00, & 0.0068798081174711D+00, & 0.0068798081174711D+00, & 0.0068798081174711D+00, & 0.0118909554500764D+00, & 0.0118909554500764D+00, & 0.0118909554500764D+00, & 0.0118909554500764D+00, & 0.0118909554500764D+00, & 0.0118909554500764D+00, & 0.0022652672511285D+00, & 0.0022652672511285D+00, & 0.0022652672511285D+00, & 0.0022652672511285D+00, & 0.0022652672511285D+00, & 0.0022652672511285D+00, & 0.0034200550598036D+00, & 0.0034200550598036D+00, & 0.0034200550598036D+00, & 0.0034200550598036D+00, & 0.0034200550598036D+00, & 0.0034200550598036D+00, & 0.0088737445510102D+00, & 0.0088737445510102D+00, & 0.0088737445510102D+00, & 0.0088737445510102D+00, & 0.0088737445510102D+00, & 0.0088737445510102D+00, & 0.0025053304372899D+00, & 0.0025053304372899D+00, & 0.0025053304372899D+00, & 0.0025053304372899D+00, & 0.0025053304372899D+00, & 0.0025053304372899D+00, & 0.0006114740634805D+00, & 0.0006114740634805D+00, & 0.0006114740634805D+00, & 0.0006114740634805D+00, & 0.0006114740634805D+00, & 0.0006114740634805D+00, & 0.0127410876559122D+00, & 0.0127410876559122D+00, & 0.0127410876559122D+00, & 0.0127410876559122D+00, & 0.0127410876559122D+00, & 0.0127410876559122D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule19 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule19() returns the rule of precision 19. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 73 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00, & 0.0525238903512090D+00, & 0.8949522192975821D+00, & 0.0525238903512090D+00, & 0.4925126750413369D+00, & 0.0149746499173263D+00, & 0.4925126750413369D+00, & 0.1114488733230214D+00, & 0.7771022533539573D+00, & 0.1114488733230214D+00, & 0.4591942010395437D+00, & 0.0816115979209127D+00, & 0.4591942010395437D+00, & 0.4039697225519012D+00, & 0.1920605548961976D+00, & 0.4039697225519012D+00, & 0.1781701047817643D+00, & 0.6436597904364714D+00, & 0.1781701047817643D+00, & 0.0116394611837894D+00, & 0.9767210776324211D+00, & 0.0116394611837894D+00, & 0.2551616329136077D+00, & 0.4896767341727846D+00, & 0.2551616329136077D+00, & 0.1306976762680324D+00, & 0.8301564644002754D+00, & 0.0391458593316922D+00, & 0.8301564644002754D+00, & 0.0391458593316922D+00, & 0.1306976762680324D+00, & 0.3113176298095413D+00, & 0.5593698057203009D+00, & 0.1293125644701578D+00, & 0.5593698057203009D+00, & 0.1293125644701578D+00, & 0.3113176298095413D+00, & 0.0020689258966048D+00, & 0.6333132931287841D+00, & 0.3646177809746111D+00, & 0.6333132931287841D+00, & 0.3646177809746111D+00, & 0.0020689258966048D+00, & 0.0745602946016267D+00, & 0.7040048199660421D+00, & 0.2214348854323312D+00, & 0.7040048199660421D+00, & 0.2214348854323312D+00, & 0.0745602946016267D+00, & 0.0050072882573545D+00, & 0.8525669543768892D+00, & 0.1424257573657563D+00, & 0.8525669543768892D+00, & 0.1424257573657563D+00, & 0.0050072882573545D+00, & 0.0408880111960169D+00, & 0.6050839790687079D+00, & 0.3540280097352752D+00, & 0.6050839790687079D+00, & 0.3540280097352752D+00, & 0.0408880111960169D+00, & 0.2418945789605796D+00, & 0.7431813689574364D+00, & 0.0149240520819841D+00, & 0.7431813689574364D+00, & 0.0149240520819841D+00, & 0.2418945789605796D+00, & 0.0600862753223067D+00, & 0.9301376988768051D+00, & 0.0097760258008882D+00, & 0.9301376988768051D+00, & 0.0097760258008882D+00, & 0.0600862753223067D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00, & 0.8949522192975821D+00, & 0.0525238903512090D+00, & 0.0525238903512090D+00, & 0.0149746499173263D+00, & 0.4925126750413369D+00, & 0.4925126750413369D+00, & 0.7771022533539573D+00, & 0.1114488733230214D+00, & 0.1114488733230214D+00, & 0.0816115979209127D+00, & 0.4591942010395437D+00, & 0.4591942010395437D+00, & 0.1920605548961976D+00, & 0.4039697225519012D+00, & 0.4039697225519012D+00, & 0.6436597904364714D+00, & 0.1781701047817643D+00, & 0.1781701047817643D+00, & 0.9767210776324211D+00, & 0.0116394611837894D+00, & 0.0116394611837894D+00, & 0.4896767341727846D+00, & 0.2551616329136077D+00, & 0.2551616329136077D+00, & 0.8301564644002754D+00, & 0.1306976762680324D+00, & 0.8301564644002754D+00, & 0.0391458593316922D+00, & 0.1306976762680324D+00, & 0.0391458593316922D+00, & 0.5593698057203009D+00, & 0.3113176298095413D+00, & 0.5593698057203009D+00, & 0.1293125644701578D+00, & 0.3113176298095413D+00, & 0.1293125644701578D+00, & 0.6333132931287841D+00, & 0.0020689258966048D+00, & 0.6333132931287841D+00, & 0.3646177809746111D+00, & 0.0020689258966048D+00, & 0.3646177809746111D+00, & 0.7040048199660421D+00, & 0.0745602946016267D+00, & 0.7040048199660421D+00, & 0.2214348854323312D+00, & 0.0745602946016267D+00, & 0.2214348854323312D+00, & 0.8525669543768892D+00, & 0.0050072882573545D+00, & 0.8525669543768892D+00, & 0.1424257573657563D+00, & 0.0050072882573545D+00, & 0.1424257573657563D+00, & 0.6050839790687079D+00, & 0.0408880111960169D+00, & 0.6050839790687079D+00, & 0.3540280097352752D+00, & 0.0408880111960169D+00, & 0.3540280097352752D+00, & 0.7431813689574364D+00, & 0.2418945789605796D+00, & 0.7431813689574364D+00, & 0.0149240520819841D+00, & 0.2418945789605796D+00, & 0.0149240520819841D+00, & 0.9301376988768051D+00, & 0.0600862753223067D+00, & 0.9301376988768051D+00, & 0.0097760258008882D+00, & 0.0600862753223067D+00, & 0.0097760258008882D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00, & 0.0525238903512090D+00, & 0.0525238903512090D+00, & 0.8949522192975821D+00, & 0.4925126750413368D+00, & 0.4925126750413369D+00, & 0.0149746499173262D+00, & 0.1114488733230212D+00, & 0.1114488733230213D+00, & 0.7771022533539571D+00, & 0.4591942010395437D+00, & 0.4591942010395437D+00, & 0.0816115979209127D+00, & 0.4039697225519012D+00, & 0.4039697225519012D+00, & 0.1920605548961976D+00, & 0.1781701047817643D+00, & 0.1781701047817643D+00, & 0.6436597904364714D+00, & 0.0116394611837894D+00, & 0.0116394611837894D+00, & 0.9767210776324211D+00, & 0.2551616329136077D+00, & 0.2551616329136077D+00, & 0.4896767341727846D+00, & 0.0391458593316922D+00, & 0.0391458593316922D+00, & 0.1306976762680324D+00, & 0.1306976762680324D+00, & 0.8301564644002754D+00, & 0.8301564644002754D+00, & 0.1293125644701578D+00, & 0.1293125644701578D+00, & 0.3113176298095413D+00, & 0.3113176298095413D+00, & 0.5593698057203009D+00, & 0.5593698057203009D+00, & 0.3646177809746111D+00, & 0.3646177809746111D+00, & 0.0020689258966048D+00, & 0.0020689258966048D+00, & 0.6333132931287841D+00, & 0.6333132931287841D+00, & 0.2214348854323313D+00, & 0.2214348854323312D+00, & 0.0745602946016267D+00, & 0.0745602946016267D+00, & 0.7040048199660423D+00, & 0.7040048199660423D+00, & 0.1424257573657564D+00, & 0.1424257573657563D+00, & 0.0050072882573544D+00, & 0.0050072882573545D+00, & 0.8525669543768892D+00, & 0.8525669543768892D+00, & 0.3540280097352753D+00, & 0.3540280097352752D+00, & 0.0408880111960169D+00, & 0.0408880111960169D+00, & 0.6050839790687079D+00, & 0.6050839790687080D+00, & 0.0149240520819841D+00, & 0.0149240520819841D+00, & 0.2418945789605796D+00, & 0.2418945789605796D+00, & 0.7431813689574365D+00, & 0.7431813689574365D+00, & 0.0097760258008882D+00, & 0.0097760258008882D+00, & 0.0600862753223068D+00, & 0.0600862753223068D+00, & 0.9301376988768052D+00, & 0.9301376988768051D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0172346988520062D+00, & 0.0035546282988991D+00, & 0.0035546282988991D+00, & 0.0035546282988991D+00, & 0.0051608775714721D+00, & 0.0051608775714721D+00, & 0.0051608775714721D+00, & 0.0076171755465091D+00, & 0.0076171755465091D+00, & 0.0076171755465091D+00, & 0.0114917950133708D+00, & 0.0114917950133708D+00, & 0.0114917950133708D+00, & 0.0157687674465775D+00, & 0.0157687674465775D+00, & 0.0157687674465775D+00, & 0.0123259574240954D+00, & 0.0123259574240954D+00, & 0.0123259574240954D+00, & 0.0008826613882214D+00, & 0.0008826613882214D+00, & 0.0008826613882214D+00, & 0.0158765096830015D+00, & 0.0158765096830015D+00, & 0.0158765096830015D+00, & 0.0048477422434275D+00, & 0.0048477422434275D+00, & 0.0048477422434275D+00, & 0.0048477422434275D+00, & 0.0048477422434275D+00, & 0.0048477422434275D+00, & 0.0131731609886954D+00, & 0.0131731609886954D+00, & 0.0131731609886954D+00, & 0.0131731609886954D+00, & 0.0131731609886954D+00, & 0.0131731609886954D+00, & 0.0016410382759179D+00, & 0.0016410382759179D+00, & 0.0016410382759179D+00, & 0.0016410382759179D+00, & 0.0016410382759179D+00, & 0.0016410382759179D+00, & 0.0090539724656062D+00, & 0.0090539724656062D+00, & 0.0090539724656062D+00, & 0.0090539724656062D+00, & 0.0090539724656062D+00, & 0.0090539724656062D+00, & 0.0014631575517351D+00, & 0.0014631575517351D+00, & 0.0014631575517351D+00, & 0.0014631575517351D+00, & 0.0014631575517351D+00, & 0.0014631575517351D+00, & 0.0080510813820121D+00, & 0.0080510813820121D+00, & 0.0080510813820121D+00, & 0.0080510813820121D+00, & 0.0080510813820121D+00, & 0.0080510813820121D+00, & 0.0042279437497682D+00, & 0.0042279437497682D+00, & 0.0042279437497682D+00, & 0.0042279437497682D+00, & 0.0042279437497682D+00, & 0.0042279437497682D+00, & 0.0016636006814297D+00, & 0.0016636006814297D+00, & 0.0016636006814297D+00, & 0.0016636006814297D+00, & 0.0016636006814297D+00, & 0.0016636006814297D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end subroutine rule20 ( n, a, b, c, w ) !*****************************************************************************80 ! !! rule20() returns the rule of precision 20. ! ! Discussion: ! ! The data is given for the following triangle: ! ! (0,1) ! | \ ! | \ ! | \ ! | \ ! (0,0)----(1,0) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 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 points and weights. ! ! Output: ! ! real ( kind = rk ) a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real ( kind = rk ) w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) integer n integer, parameter :: n_save = 79 real ( kind = rk ) a(n) real ( kind = rk ), save, dimension ( n_save ) :: a_save = (/ & 0.3333333333333334D+00, & 0.2545792676733391D+00, & 0.4908414646533218D+00, & 0.2545792676733391D+00, & 0.0109761410283978D+00, & 0.9780477179432046D+00, & 0.0109761410283978D+00, & 0.1093835967117146D+00, & 0.7812328065765708D+00, & 0.1093835967117146D+00, & 0.1862949977445409D+00, & 0.6274100045109181D+00, & 0.1862949977445409D+00, & 0.4455510569559248D+00, & 0.1088978860881504D+00, & 0.4455510569559248D+00, & 0.0373108805988847D+00, & 0.9253782388022307D+00, & 0.0373108805988847D+00, & 0.3934253478170999D+00, & 0.2131493043658003D+00, & 0.3934253478170999D+00, & 0.4762456115404990D+00, & 0.0475087769190020D+00, & 0.4762456115404990D+00, & 0.0075707805046965D+00, & 0.8332955118382362D+00, & 0.1591337076570672D+00, & 0.8332955118382362D+00, & 0.1591337076570672D+00, & 0.0075707805046965D+00, & 0.0465603649076643D+00, & 0.7549215028635474D+00, & 0.1985181322287882D+00, & 0.7549215028635474D+00, & 0.1985181322287882D+00, & 0.0465603649076643D+00, & 0.0640905856084341D+00, & 0.9310544767839422D+00, & 0.0048549376076237D+00, & 0.9310544767839422D+00, & 0.0048549376076237D+00, & 0.0640905856084341D+00, & 0.0549874791429868D+00, & 0.6118777035474257D+00, & 0.3331348173095875D+00, & 0.6118777035474257D+00, & 0.3331348173095875D+00, & 0.0549874791429868D+00, & 0.0999522962881387D+00, & 0.8616840189364867D+00, & 0.0383636847753746D+00, & 0.8616840189364867D+00, & 0.0383636847753746D+00, & 0.0999522962881387D+00, & 0.1062272047202700D+00, & 0.6781657378896355D+00, & 0.2156070573900944D+00, & 0.6781657378896355D+00, & 0.2156070573900944D+00, & 0.1062272047202700D+00, & 0.4200237588162241D+00, & 0.5701446928909734D+00, & 0.0098315482928026D+00, & 0.5701446928909734D+00, & 0.0098315482928026D+00, & 0.4200237588162241D+00, & 0.3178601238357720D+00, & 0.5423318041724281D+00, & 0.1398080719917999D+00, & 0.5423318041724281D+00, & 0.1398080719917999D+00, & 0.3178601238357720D+00, & 0.0107372128560111D+00, & 0.7086813757203236D+00, & 0.2805814114236652D+00, & 0.7086813757203236D+00, & 0.2805814114236652D+00, & 0.0107372128560111D+00 /) real ( kind = rk ) b(n) real ( kind = rk ), save, dimension ( n_save ) :: b_save = (/ & 0.3333333333333334D+00, & 0.4908414646533218D+00, & 0.2545792676733391D+00, & 0.2545792676733391D+00, & 0.9780477179432046D+00, & 0.0109761410283978D+00, & 0.0109761410283978D+00, & 0.7812328065765708D+00, & 0.1093835967117146D+00, & 0.1093835967117146D+00, & 0.6274100045109181D+00, & 0.1862949977445409D+00, & 0.1862949977445409D+00, & 0.1088978860881504D+00, & 0.4455510569559248D+00, & 0.4455510569559248D+00, & 0.9253782388022307D+00, & 0.0373108805988847D+00, & 0.0373108805988847D+00, & 0.2131493043658003D+00, & 0.3934253478170999D+00, & 0.3934253478170999D+00, & 0.0475087769190020D+00, & 0.4762456115404990D+00, & 0.4762456115404990D+00, & 0.8332955118382362D+00, & 0.0075707805046965D+00, & 0.8332955118382362D+00, & 0.1591337076570672D+00, & 0.0075707805046965D+00, & 0.1591337076570672D+00, & 0.7549215028635474D+00, & 0.0465603649076643D+00, & 0.7549215028635474D+00, & 0.1985181322287882D+00, & 0.0465603649076643D+00, & 0.1985181322287882D+00, & 0.9310544767839422D+00, & 0.0640905856084341D+00, & 0.9310544767839422D+00, & 0.0048549376076237D+00, & 0.0640905856084341D+00, & 0.0048549376076237D+00, & 0.6118777035474257D+00, & 0.0549874791429868D+00, & 0.6118777035474257D+00, & 0.3331348173095875D+00, & 0.0549874791429868D+00, & 0.3331348173095875D+00, & 0.8616840189364867D+00, & 0.0999522962881387D+00, & 0.8616840189364867D+00, & 0.0383636847753746D+00, & 0.0999522962881387D+00, & 0.0383636847753746D+00, & 0.6781657378896355D+00, & 0.1062272047202700D+00, & 0.6781657378896355D+00, & 0.2156070573900944D+00, & 0.1062272047202700D+00, & 0.2156070573900944D+00, & 0.5701446928909734D+00, & 0.4200237588162241D+00, & 0.5701446928909734D+00, & 0.0098315482928026D+00, & 0.4200237588162241D+00, & 0.0098315482928026D+00, & 0.5423318041724281D+00, & 0.3178601238357720D+00, & 0.5423318041724281D+00, & 0.1398080719917999D+00, & 0.3178601238357720D+00, & 0.1398080719917999D+00, & 0.7086813757203236D+00, & 0.0107372128560111D+00, & 0.7086813757203236D+00, & 0.2805814114236652D+00, & 0.0107372128560111D+00, & 0.2805814114236652D+00 /) real ( kind = rk ) c(n) real ( kind = rk ), save, dimension ( n_save ) :: c_save = (/ & 0.3333333333333333D+00, & 0.2545792676733390D+00, & 0.2545792676733391D+00, & 0.4908414646533217D+00, & 0.0109761410283977D+00, & 0.0109761410283977D+00, & 0.9780477179432046D+00, & 0.1093835967117146D+00, & 0.1093835967117146D+00, & 0.7812328065765708D+00, & 0.1862949977445409D+00, & 0.1862949977445409D+00, & 0.6274100045109181D+00, & 0.4455510569559248D+00, & 0.4455510569559248D+00, & 0.1088978860881504D+00, & 0.0373108805988847D+00, & 0.0373108805988846D+00, & 0.9253782388022307D+00, & 0.3934253478170998D+00, & 0.3934253478170999D+00, & 0.2131493043658002D+00, & 0.4762456115404991D+00, & 0.4762456115404991D+00, & 0.0475087769190021D+00, & 0.1591337076570672D+00, & 0.1591337076570672D+00, & 0.0075707805046965D+00, & 0.0075707805046965D+00, & 0.8332955118382362D+00, & 0.8332955118382362D+00, & 0.1985181322287883D+00, & 0.1985181322287883D+00, & 0.0465603649076645D+00, & 0.0465603649076644D+00, & 0.7549215028635476D+00, & 0.7549215028635476D+00, & 0.0048549376076238D+00, & 0.0048549376076238D+00, & 0.0640905856084342D+00, & 0.0640905856084341D+00, & 0.9310544767839423D+00, & 0.9310544767839422D+00, & 0.3331348173095876D+00, & 0.3331348173095875D+00, & 0.0549874791429869D+00, & 0.0549874791429869D+00, & 0.6118777035474257D+00, & 0.6118777035474258D+00, & 0.0383636847753746D+00, & 0.0383636847753746D+00, & 0.0999522962881387D+00, & 0.0999522962881387D+00, & 0.8616840189364867D+00, & 0.8616840189364867D+00, & 0.2156070573900944D+00, & 0.2156070573900944D+00, & 0.1062272047202701D+00, & 0.1062272047202701D+00, & 0.6781657378896356D+00, & 0.6781657378896355D+00, & 0.0098315482928025D+00, & 0.0098315482928025D+00, & 0.4200237588162240D+00, & 0.4200237588162241D+00, & 0.5701446928909732D+00, & 0.5701446928909732D+00, & 0.1398080719917999D+00, & 0.1398080719917999D+00, & 0.3178601238357720D+00, & 0.3178601238357720D+00, & 0.5423318041724281D+00, & 0.5423318041724281D+00, & 0.2805814114236653D+00, & 0.2805814114236653D+00, & 0.0107372128560111D+00, & 0.0107372128560111D+00, & 0.7086813757203236D+00, & 0.7086813757203237D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0139101107014531D+00, & 0.0140832013075202D+00, & 0.0140832013075202D+00, & 0.0140832013075202D+00, & 0.0007988407910666D+00, & 0.0007988407910666D+00, & 0.0007988407910666D+00, & 0.0078302307760745D+00, & 0.0078302307760745D+00, & 0.0078302307760745D+00, & 0.0091734629742529D+00, & 0.0091734629742529D+00, & 0.0091734629742529D+00, & 0.0094523999332324D+00, & 0.0094523999332324D+00, & 0.0094523999332324D+00, & 0.0021612754106656D+00, & 0.0021612754106656D+00, & 0.0021612754106656D+00, & 0.0137880506290705D+00, & 0.0137880506290705D+00, & 0.0137880506290705D+00, & 0.0071018253034084D+00, & 0.0071018253034084D+00, & 0.0071018253034084D+00, & 0.0022028974185585D+00, & 0.0022028974185585D+00, & 0.0022028974185585D+00, & 0.0022028974185585D+00, & 0.0022028974185585D+00, & 0.0022028974185585D+00, & 0.0059863985789547D+00, & 0.0059863985789547D+00, & 0.0059863985789547D+00, & 0.0059863985789547D+00, & 0.0059863985789547D+00, & 0.0059863985789547D+00, & 0.0011298696021259D+00, & 0.0011298696021259D+00, & 0.0011298696021259D+00, & 0.0011298696021259D+00, & 0.0011298696021259D+00, & 0.0011298696021259D+00, & 0.0086672255672193D+00, & 0.0086672255672193D+00, & 0.0086672255672193D+00, & 0.0086672255672193D+00, & 0.0086672255672193D+00, & 0.0086672255672193D+00, & 0.0041457115276139D+00, & 0.0041457115276139D+00, & 0.0041457115276139D+00, & 0.0041457115276139D+00, & 0.0041457115276139D+00, & 0.0041457115276139D+00, & 0.0077226078220992D+00, & 0.0077226078220992D+00, & 0.0077226078220992D+00, & 0.0077226078220992D+00, & 0.0077226078220992D+00, & 0.0077226078220992D+00, & 0.0036956815002553D+00, & 0.0036956815002553D+00, & 0.0036956815002553D+00, & 0.0036956815002553D+00, & 0.0036956815002553D+00, & 0.0036956815002553D+00, & 0.0116917457318277D+00, & 0.0116917457318277D+00, & 0.0116917457318277D+00, & 0.0116917457318277D+00, & 0.0116917457318277D+00, & 0.0116917457318277D+00, & 0.0035782002384577D+00, & 0.0035782002384577D+00, & 0.0035782002384577D+00, & 0.0035782002384577D+00, & 0.0035782002384577D+00, & 0.0035782002384577D+00 /) a = a_save b = b_save c = c_save w = w_save ! ! The weights should sum to 1, not 1/2! ! w = w * 2.0D+00 return end function triangle_area ( vert1, vert2, vert3 ) !*****************************************************************************80 ! !! triangle_area() returns the area of a triangle. ! ! Licensing: ! ! This code is distributed under the GNU GPL license. ! ! Modified: ! ! 22 April 2023 ! ! Author: ! ! John Burkardt. ! ! Input: ! ! real ( kind = rk ) vert1(2), vert2(2), vert3(2): the vertex coordinates. ! ! Output: ! ! real ( kind = rk ) triangle_area: the area of the triangle. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) real ( kind = rk ) triangle_area real ( kind = rk ) vert1(2) real ( kind = rk ) vert2(2) real ( kind = rk ) vert3(2) triangle_area = 0.5D+00 * & ( & ( vert2(1) - vert1(1) ) * ( vert3(2) - vert1(2) ) & - ( vert3(1) - vert1(1) ) * ( vert2(2) - vert1(2) ) & ) return end function triangle_unit_area ( ) !*****************************************************************************80 ! !! triangle_unit_area() returns the area of a unit triangle. ! ! Licensing: ! ! This code is distributed under the GNU GPL license. ! ! Modified: ! ! 24 May 2023 ! ! Author: ! ! John Burkardt. ! ! Output: ! ! real ( kind = rk ) triangle_unit_area: the area of the triangle. ! implicit none integer, parameter :: rk = kind ( 1.0D+00) real ( kind = rk ) triangle_unit_area triangle_unit_area = 0.5D+00 return end function triangle_unit_monomial_integral ( expon ) !*****************************************************************************80 ! !! triangle_unit_monomial_integral(): integral of a monomial over the unit triangle. ! ! Discussion: ! ! This routine integrates a monomial of the form ! ! product ( 1 <= dim <= dim_num ) x(dim)^expon(dim) ! ! where the exponents are nonnegative integers. Note that ! if the combination 0^0 is encountered, it should be treated ! as 1. ! ! Integral ( over unit triangle ) x^m y^n dx dy = m! * n! / ( m + n + 2 )! ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 July 2007 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer EXPON(3), the exponents. ! ! Output: ! ! real ( kind = rk ) triangle_unit_monomial_integral: the value of the integral. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer expon(3) integer i integer k real ( kind = rk ) triangle_unit_monomial_integral real ( kind = rk ) value ! ! The first computation ends with VALUE = 1.0; ! value = 1.0D+00 k = 0 do i = 1, expon(1) k = k + 1 ! value = value * real ( i, kind = rk ) / real ( k, kind = rk ) end do do i = 1, expon(2) k = k + 1 value = value * real ( i, kind = rk ) / real ( k, kind = rk ) end do k = k + 1 value = value / real ( k, kind = rk ) k = k + 1 value = value / real ( k, kind = rk ) triangle_unit_monomial_integral = value return end subroutine triangle_witherden_rule ( p, n, a, b, c, w ) !*****************************************************************************80 ! !! triangle_witherden_rule() returns a triangle quadrature rule of given precision. ! ! 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: ! ! integer p: the precision, 0 <= p <= 20. ! ! integer n: the order of the rule. ! ! Output: ! ! real a(n), b(n), c(n): the barycentric coordinates of quadrature points. ! ! real w(n): the weights of quadrature points. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a(n) real ( kind = rk ) b(n) real ( kind = rk ) c(n) integer p real ( kind = rk ) w(n) if ( p < 0 ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'triangle_witherden_rule(): Fatal error!' write ( *, '(a)' ) ' Input p < 0.' call exit ( 1 ) end if if ( 20 < p ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'triangle_witherden_rule(): Fatal error!' write ( *, '(a)' ) ' Input 20 < p.' call exit ( 1 ) end if if ( p == 0 ) then call rule00 ( n, a, b, c, w ) else if ( p == 1 ) then call rule01 ( n, a, b, c, w ) else if ( p == 2 ) then call rule02 ( n, a, b, c, w ) else if ( p == 3 ) then call rule03 ( n, a, b, c, w ) else if ( p == 4 ) then call rule04 ( n, a, b, c, w ) else if ( p == 5 ) then call rule05 ( n, a, b, c, w ) else if ( p == 6 ) then call rule06 ( n, a, b, c, w ) else if ( p == 7 ) then call rule07 ( n, a, b, c, w ) else if ( p == 8 ) then call rule08 ( n, a, b, c, w ) else if ( p == 9 ) then call rule09 ( n, a, b, c, w ) else if ( p == 10 ) then call rule10 ( n, a, b, c, w ) else if ( p == 11 ) then call rule11 ( n, a, b, c, w ) else if ( p == 12 ) then call rule12 ( n, a, b, c, w ) else if ( p == 13 ) then call rule13 ( n, a, b, c, w ) else if ( p == 14 ) then call rule14 ( n, a, b, c, w ) else if ( p == 15 ) then call rule15 ( n, a, b, c, w ) else if ( p == 16 ) then call rule16 ( n, a, b, c, w ) else if ( p == 17 ) then call rule17 ( n, a, b, c, w ) else if ( p == 18 ) then call rule18 ( n, a, b, c, w ) else if ( p == 19 ) then call rule19 ( n, a, b, c, w ) else if ( p == 20 ) then call rule20 ( n, a, b, c, w ) end if return end