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 hexahedron_unit_monomial_integral ( expon, value ) !*****************************************************************************80 ! !! hexahedron_unit_monomial_integral(): monomial integral in a unit hexahedron. ! ! Discussion: ! ! This routine returns the integral of ! ! product ( 1 <= I <= 3 ) X(I)^EXPON(I) ! ! over the unit hexahedron. ! ! The unit hexahedron H is bounded by the vertices ! (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 May 2023 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer EXPON(3): the exponents. ! ! Output: ! ! real VALUE: the integral of the monomial. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer expon(3) real ( kind = rk ) value value = 1.0D+00 / ( expon(1) + 1 ) & * 1.0D+00 / ( expon(2) + 1 ) & * 1.0D+00 / ( expon(3) + 1 ) return end function hexahedron_unit_volume ( ) !*****************************************************************************80 ! !! hexahedron_unit_volume() returns the volume of a unit hexahedron. ! ! Discussion: ! ! The unit hexahedron has vertices ! (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 May 2023 ! ! Author: ! ! John Burkardt ! ! Output: ! ! real hexahedron_unit_volume: the volume. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) hexahedron_unit_volume hexahedron_unit_volume = 1.0 return end subroutine hexahedron_witherden_rule ( p, n, x, y, z, w ) !*****************************************************************************80 ! !! hexahedron_witherden_rule() returns a hexahedron quadrature rule of given precision. ! ! Discussion: ! ! The unit hexahedron H is bounded by the vertices ! (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! int p: the precision, 0 <= p <= 11. ! ! int n: the order of the rule. ! ! Output: ! ! double x(n), y(n), z(n): the coordinates of quadrature points. ! ! double w(n): the quadrature weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer p real ( kind = rk ) w(n) real ( kind = rk ) x(n) real ( kind = rk ) y(n) real ( kind = rk ) z(n) if ( p < 0 ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'hexahedron_witherden_rule(): Fatal error!' write ( *, '(a)' ) ' Input p < 0.' stop ( 1 ) end if if ( 11 < p ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'hexahedron_witherden_rule(): Fatal error!' write ( *, '(a)' ) ' Input 11 < p.' stop ( 1 ) end if if ( p <= 1 ) then call rule01 ( n, x, y, z, w ) else if ( p <= 3 ) then call rule03 ( n, x, y, z, w ) else if ( p <= 5 ) then call rule05 ( n, x, y, z, w ) else if ( p <= 7 ) then call rule07 ( n, x, y, z, w ) else if ( p <= 9 ) then call rule09 ( n, x, y, z, w ) else if ( p <= 11 ) then call rule11 ( n, x, y, z, w ) end if 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 hexahedron quadrature rule of given precision. ! ! Discussion: ! ! The "order" of a quadrature rule is the number of points involved. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer p: the precision, 0 <= p <= 11. ! ! Output: ! ! integer order: the order of the rule. ! implicit none integer order integer, dimension ( 0 : 11 ) :: order_vec = (/ & 1, & 1, 6, 6, 14, 14, 34, 34, 58, 58, 90, & 90 /) integer p if ( p < 0 ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'rule_order(): Fatal error.' write ( *, '(a)' ) ' Input p < 0.' end if if ( 11 < p ) then write ( *, '(a)' ) '' write ( *, '(a)' ) 'rule_order(): Fatal error.' write ( *, '(a)' ) ' Input 11 < p.' end if order = order_vec(p) return end subroutine rule01 ( n, x, y, z, w ) !*****************************************************************************80 ! !! rule01() returns the hexahedron rule of precision 1. ! ! Discussion: ! ! We suppose we are given a hexahedron H with vertices ! (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1), ! We call a rule with n points, returning coordinates ! x, y, z, and weights w. Then the integral I of f(x,y,z) over H is ! approximated by Q as follows: ! ! Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n), z(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 1 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.5000000000000000D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00 /) real ( kind = rk ) z(n) real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ & 0.5000000000000000D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 1.0000000000000000D+00 /) x = x_save y = y_save z = z_save w = w_save return end subroutine rule03 ( n, x, y, z, w ) !*****************************************************************************80 ! !! rule03() returns the hexahedron rule of precision 3. ! ! Discussion: ! ! We suppose we are given a hexahedron H with vertices ! (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1), ! We call a rule with n points, returning coordinates ! x, y, z, and weights w. Then the integral I of f(x,y,z) over H is ! approximated by Q as follows: ! ! Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n), z(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 6 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.0000000000000000D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 1.0000000000000000D+00, & 0.5000000000000000D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 1.0000000000000000D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.0000000000000000D+00 /) real ( kind = rk ) z(n) real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ & 0.5000000000000000D+00, & 1.0000000000000000D+00, & 0.5000000000000000D+00, & 0.0000000000000000D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.1666666666666667D+00, & 0.1666666666666667D+00, & 0.1666666666666667D+00, & 0.1666666666666667D+00, & 0.1666666666666667D+00, & 0.1666666666666667D+00 /) x = x_save y = y_save z = z_save w = w_save return end subroutine rule05 ( n, x, y, z, w ) !*****************************************************************************80 ! !! rule05() returns the hexahedron rule of precision 5. ! ! Discussion: ! ! We suppose we are given a hexahedron H with vertices ! (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1), ! We call a rule with n points, returning coordinates ! x, y, z, and weights w. Then the integral I of f(x,y,z) over H is ! approximated by Q as follows: ! ! Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n), z(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 14 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.1020887871228893D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.8979112128771107D+00, & 0.5000000000000000D+00, & 0.8793934553196641D+00, & 0.1206065446803359D+00, & 0.1206065446803359D+00, & 0.1206065446803359D+00, & 0.1206065446803359D+00, & 0.8793934553196641D+00, & 0.8793934553196641D+00, & 0.8793934553196641D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.8979112128771107D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.1020887871228893D+00, & 0.1206065446803359D+00, & 0.8793934553196641D+00, & 0.8793934553196641D+00, & 0.1206065446803359D+00, & 0.1206065446803359D+00, & 0.8793934553196641D+00, & 0.8793934553196641D+00, & 0.1206065446803359D+00 /) real ( kind = rk ) z(n) real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ & 0.5000000000000000D+00, & 0.8979112128771107D+00, & 0.5000000000000000D+00, & 0.1020887871228893D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.1206065446803359D+00, & 0.8793934553196641D+00, & 0.1206065446803359D+00, & 0.1206065446803359D+00, & 0.8793934553196641D+00, & 0.1206065446803359D+00, & 0.8793934553196641D+00, & 0.8793934553196641D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.1108033240997230D+00, & 0.1108033240997230D+00, & 0.1108033240997230D+00, & 0.1108033240997230D+00, & 0.1108033240997230D+00, & 0.1108033240997230D+00, & 0.0418975069252078D+00, & 0.0418975069252078D+00, & 0.0418975069252078D+00, & 0.0418975069252078D+00, & 0.0418975069252078D+00, & 0.0418975069252078D+00, & 0.0418975069252078D+00, & 0.0418975069252078D+00 /) x = x_save y = y_save z = z_save w = w_save return end subroutine rule07 ( n, x, y, z, w ) !*****************************************************************************80 ! !! rule07() returns the hexahedron rule of precision 7. ! ! Discussion: ! ! We suppose we are given a hexahedron H with vertices ! (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1), ! We call a rule with n points, returning coordinates ! x, y, z, and weights w. Then the integral I of f(x,y,z) over H is ! approximated by Q as follows: ! ! Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n), z(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 34 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.0059569194056617D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.9940430805943383D+00, & 0.5000000000000000D+00, & 0.7039775836791597D+00, & 0.2960224163208403D+00, & 0.2960224163208403D+00, & 0.2960224163208403D+00, & 0.2960224163208403D+00, & 0.7039775836791597D+00, & 0.7039775836791597D+00, & 0.7039775836791597D+00, & 0.8905514105020593D+00, & 0.1094485894979407D+00, & 0.1094485894979407D+00, & 0.1094485894979407D+00, & 0.1094485894979407D+00, & 0.8905514105020593D+00, & 0.8905514105020593D+00, & 0.8905514105020593D+00, & 0.0759738621579806D+00, & 0.9240261378420194D+00, & 0.5000000000000000D+00, & 0.9240261378420194D+00, & 0.9240261378420194D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.0759738621579806D+00, & 0.0759738621579806D+00, & 0.9240261378420194D+00, & 0.5000000000000000D+00, & 0.0759738621579806D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.9940430805943383D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.0059569194056617D+00, & 0.2960224163208403D+00, & 0.7039775836791597D+00, & 0.7039775836791597D+00, & 0.2960224163208403D+00, & 0.2960224163208403D+00, & 0.7039775836791597D+00, & 0.7039775836791597D+00, & 0.2960224163208403D+00, & 0.1094485894979407D+00, & 0.8905514105020593D+00, & 0.8905514105020593D+00, & 0.1094485894979407D+00, & 0.1094485894979407D+00, & 0.8905514105020593D+00, & 0.8905514105020593D+00, & 0.1094485894979407D+00, & 0.0759738621579806D+00, & 0.5000000000000000D+00, & 0.9240261378420194D+00, & 0.9240261378420194D+00, & 0.5000000000000000D+00, & 0.0759738621579806D+00, & 0.0759738621579806D+00, & 0.5000000000000000D+00, & 0.9240261378420194D+00, & 0.0759738621579806D+00, & 0.9240261378420194D+00, & 0.5000000000000000D+00 /) real ( kind = rk ) z(n) real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ & 0.5000000000000000D+00, & 0.9940430805943383D+00, & 0.5000000000000000D+00, & 0.0059569194056617D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.2960224163208403D+00, & 0.7039775836791597D+00, & 0.2960224163208403D+00, & 0.2960224163208403D+00, & 0.7039775836791597D+00, & 0.2960224163208403D+00, & 0.7039775836791597D+00, & 0.7039775836791597D+00, & 0.1094485894979407D+00, & 0.8905514105020593D+00, & 0.1094485894979407D+00, & 0.1094485894979407D+00, & 0.8905514105020593D+00, & 0.1094485894979407D+00, & 0.8905514105020593D+00, & 0.8905514105020593D+00, & 0.5000000000000000D+00, & 0.0759738621579806D+00, & 0.0759738621579806D+00, & 0.5000000000000000D+00, & 0.9240261378420194D+00, & 0.9240261378420194D+00, & 0.0759738621579806D+00, & 0.9240261378420194D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.9240261378420194D+00, & 0.0759738621579806D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0250162363729459D+00, & 0.0250162363729459D+00, & 0.0250162363729459D+00, & 0.0250162363729459D+00, & 0.0250162363729459D+00, & 0.0250162363729459D+00, & 0.0571442320087316D+00, & 0.0571442320087316D+00, & 0.0571442320087316D+00, & 0.0571442320087316D+00, & 0.0571442320087316D+00, & 0.0571442320087316D+00, & 0.0571442320087316D+00, & 0.0571442320087316D+00, & 0.0192245175082448D+00, & 0.0192245175082448D+00, & 0.0192245175082448D+00, & 0.0192245175082448D+00, & 0.0192245175082448D+00, & 0.0192245175082448D+00, & 0.0192245175082448D+00, & 0.0192245175082448D+00, & 0.0199127154688761D+00, & 0.0199127154688761D+00, & 0.0199127154688761D+00, & 0.0199127154688761D+00, & 0.0199127154688761D+00, & 0.0199127154688761D+00, & 0.0199127154688761D+00, & 0.0199127154688761D+00, & 0.0199127154688761D+00, & 0.0199127154688761D+00, & 0.0199127154688761D+00, & 0.0199127154688761D+00 /) x = x_save y = y_save z = z_save w = w_save return end subroutine rule09 ( n, x, y, z, w ) !*****************************************************************************80 ! !! rule09() returns the hexahedron rule of precision 9. ! ! Discussion: ! ! We suppose we are given a hexahedron H with vertices ! (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1), ! We call a rule with n points, returning coordinates ! x, y, z, and weights w. Then the integral I of f(x,y,z) over H is ! approximated by Q as follows: ! ! Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n), z(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 58 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.1931592652041455D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.8068407347958545D+00, & 0.5000000000000000D+00, & 0.9350498923309880D+00, & 0.0649501076690120D+00, & 0.0649501076690120D+00, & 0.0649501076690120D+00, & 0.0649501076690120D+00, & 0.9350498923309880D+00, & 0.9350498923309880D+00, & 0.9350498923309880D+00, & 0.7820554035100150D+00, & 0.2179445964899850D+00, & 0.2179445964899850D+00, & 0.2179445964899850D+00, & 0.2179445964899850D+00, & 0.7820554035100150D+00, & 0.7820554035100150D+00, & 0.7820554035100150D+00, & 0.0611564383711609D+00, & 0.9388435616288391D+00, & 0.5000000000000000D+00, & 0.9388435616288391D+00, & 0.9388435616288391D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.0611564383711609D+00, & 0.0611564383711609D+00, & 0.9388435616288391D+00, & 0.5000000000000000D+00, & 0.0611564383711609D+00, & 0.7161339513154311D+00, & 0.2838660486845689D+00, & 0.2838660486845689D+00, & 0.0307347890676641D+00, & 0.2838660486845689D+00, & 0.9692652109323359D+00, & 0.7161339513154311D+00, & 0.7161339513154311D+00, & 0.2838660486845689D+00, & 0.7161339513154311D+00, & 0.0307347890676641D+00, & 0.2838660486845689D+00, & 0.2838660486845689D+00, & 0.2838660486845689D+00, & 0.7161339513154311D+00, & 0.9692652109323359D+00, & 0.7161339513154311D+00, & 0.9692652109323359D+00, & 0.9692652109323359D+00, & 0.0307347890676641D+00, & 0.7161339513154311D+00, & 0.2838660486845689D+00, & 0.7161339513154311D+00, & 0.0307347890676641D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.8068407347958545D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.1931592652041455D+00, & 0.0649501076690120D+00, & 0.9350498923309880D+00, & 0.9350498923309880D+00, & 0.0649501076690120D+00, & 0.0649501076690120D+00, & 0.9350498923309880D+00, & 0.9350498923309880D+00, & 0.0649501076690120D+00, & 0.2179445964899850D+00, & 0.7820554035100150D+00, & 0.7820554035100150D+00, & 0.2179445964899850D+00, & 0.2179445964899850D+00, & 0.7820554035100150D+00, & 0.7820554035100150D+00, & 0.2179445964899850D+00, & 0.0611564383711609D+00, & 0.5000000000000000D+00, & 0.9388435616288391D+00, & 0.9388435616288391D+00, & 0.5000000000000000D+00, & 0.0611564383711609D+00, & 0.0611564383711609D+00, & 0.5000000000000000D+00, & 0.9388435616288391D+00, & 0.0611564383711609D+00, & 0.9388435616288391D+00, & 0.5000000000000000D+00, & 0.9692652109323359D+00, & 0.0307347890676641D+00, & 0.9692652109323359D+00, & 0.2838660486845689D+00, & 0.2838660486845689D+00, & 0.7161339513154311D+00, & 0.0307347890676641D+00, & 0.0307347890676641D+00, & 0.0307347890676641D+00, & 0.7161339513154311D+00, & 0.7161339513154311D+00, & 0.7161339513154311D+00, & 0.7161339513154311D+00, & 0.2838660486845689D+00, & 0.2838660486845689D+00, & 0.2838660486845689D+00, & 0.9692652109323359D+00, & 0.7161339513154311D+00, & 0.2838660486845689D+00, & 0.7161339513154311D+00, & 0.7161339513154311D+00, & 0.9692652109323359D+00, & 0.2838660486845689D+00, & 0.2838660486845689D+00 /) real ( kind = rk ) z(n) real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ & 0.5000000000000000D+00, & 0.8068407347958545D+00, & 0.5000000000000000D+00, & 0.1931592652041455D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.0649501076690120D+00, & 0.9350498923309880D+00, & 0.0649501076690120D+00, & 0.0649501076690120D+00, & 0.9350498923309880D+00, & 0.0649501076690120D+00, & 0.9350498923309880D+00, & 0.9350498923309880D+00, & 0.2179445964899850D+00, & 0.7820554035100150D+00, & 0.2179445964899850D+00, & 0.2179445964899850D+00, & 0.7820554035100150D+00, & 0.2179445964899850D+00, & 0.7820554035100150D+00, & 0.7820554035100150D+00, & 0.5000000000000000D+00, & 0.0611564383711609D+00, & 0.0611564383711609D+00, & 0.5000000000000000D+00, & 0.9388435616288391D+00, & 0.9388435616288391D+00, & 0.0611564383711609D+00, & 0.9388435616288391D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.9388435616288391D+00, & 0.0611564383711609D+00, & 0.2838660486845689D+00, & 0.7161339513154311D+00, & 0.2838660486845689D+00, & 0.7161339513154311D+00, & 0.9692652109323359D+00, & 0.2838660486845689D+00, & 0.2838660486845689D+00, & 0.7161339513154311D+00, & 0.2838660486845689D+00, & 0.0307347890676641D+00, & 0.2838660486845689D+00, & 0.0307347890676641D+00, & 0.9692652109323359D+00, & 0.0307347890676641D+00, & 0.9692652109323359D+00, & 0.2838660486845689D+00, & 0.7161339513154311D+00, & 0.7161339513154311D+00, & 0.7161339513154311D+00, & 0.7161339513154311D+00, & 0.9692652109323359D+00, & 0.7161339513154311D+00, & 0.0307347890676641D+00, & 0.2838660486845689D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0541593744687068D+00, & 0.0541593744687068D+00, & 0.0541593744687068D+00, & 0.0541593744687068D+00, & 0.0541593744687068D+00, & 0.0541593744687068D+00, & 0.0062685994124186D+00, & 0.0062685994124186D+00, & 0.0062685994124186D+00, & 0.0062685994124186D+00, & 0.0062685994124186D+00, & 0.0062685994124186D+00, & 0.0062685994124186D+00, & 0.0062685994124186D+00, & 0.0248574797680029D+00, & 0.0248574797680029D+00, & 0.0248574797680029D+00, & 0.0248574797680029D+00, & 0.0248574797680029D+00, & 0.0248574797680029D+00, & 0.0248574797680029D+00, & 0.0248574797680029D+00, & 0.0114737257670222D+00, & 0.0114737257670222D+00, & 0.0114737257670222D+00, & 0.0114737257670222D+00, & 0.0114737257670222D+00, & 0.0114737257670222D+00, & 0.0114737257670222D+00, & 0.0114737257670222D+00, & 0.0114737257670222D+00, & 0.0114737257670222D+00, & 0.0114737257670222D+00, & 0.0114737257670222D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00, & 0.0120146004391717D+00 /) x = x_save y = y_save z = z_save w = w_save return end subroutine rule11 ( n, x, y, z, w ) !*****************************************************************************80 ! !! rule11() returns the hexahedron rule of precision 11. ! ! Discussion: ! ! We suppose we are given a hexahedron H with vertices ! (0,0,0), (1,0,0), (1,1,0), (0,1,0), (0,0,1), (1,0,1), (1,1,1), (0,1,1), ! We call a rule with n points, returning coordinates ! x, y, z, and weights w. Then the integral I of f(x,y,z) over H is ! approximated by Q as follows: ! ! Q = volume(H) * sum ( 1 <= i <= n ) w(i) * f(x(i),y(i),z(i)) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 May 2023 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Freddie Witherden, Peter Vincent, ! On the identification of symmetric quadrature rules for finite element methods, ! Computers and Mathematics with Applications, ! Volume 69, pages 1232-1241, 2015. ! ! Input: ! ! integer n: the number of quadrature points. ! ! Output: ! ! real x(n), y(n), z(n): the coordinates of quadrature points. ! ! real w(n): the quadrature weights, which sum to 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer, parameter :: n_save = 90 real ( kind = rk ) x(n) real ( kind = rk ), save, dimension ( n_save ) :: x_save = (/ & 0.0936928329501868D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.9063071670498133D+00, & 0.5000000000000000D+00, & 0.8008376320991313D+00, & 0.1991623679008687D+00, & 0.1991623679008687D+00, & 0.1991623679008687D+00, & 0.1991623679008687D+00, & 0.8008376320991313D+00, & 0.8008376320991313D+00, & 0.8008376320991313D+00, & 0.9277278805088800D+00, & 0.0722721194911200D+00, & 0.0722721194911200D+00, & 0.0722721194911200D+00, & 0.0722721194911200D+00, & 0.9277278805088800D+00, & 0.9277278805088800D+00, & 0.9277278805088800D+00, & 0.6566967022580273D+00, & 0.3433032977419727D+00, & 0.3433032977419727D+00, & 0.3433032977419727D+00, & 0.3433032977419727D+00, & 0.6566967022580273D+00, & 0.6566967022580273D+00, & 0.6566967022580273D+00, & 0.1326658565014960D+00, & 0.8673341434985040D+00, & 0.5000000000000000D+00, & 0.8673341434985040D+00, & 0.8673341434985040D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.1326658565014960D+00, & 0.1326658565014960D+00, & 0.8673341434985040D+00, & 0.5000000000000000D+00, & 0.1326658565014960D+00, & 0.7253999675572547D+00, & 0.2746000324427453D+00, & 0.2746000324427453D+00, & 0.0174501672436448D+00, & 0.2746000324427453D+00, & 0.9825498327563551D+00, & 0.7253999675572547D+00, & 0.7253999675572547D+00, & 0.2746000324427453D+00, & 0.7253999675572547D+00, & 0.0174501672436448D+00, & 0.2746000324427453D+00, & 0.2746000324427453D+00, & 0.2746000324427453D+00, & 0.7253999675572547D+00, & 0.9825498327563551D+00, & 0.7253999675572547D+00, & 0.9825498327563551D+00, & 0.9825498327563551D+00, & 0.0174501672436448D+00, & 0.7253999675572547D+00, & 0.2746000324427453D+00, & 0.7253999675572547D+00, & 0.0174501672436448D+00, & 0.9706224286053016D+00, & 0.0293775713946984D+00, & 0.0293775713946984D+00, & 0.3230485927016850D+00, & 0.0293775713946984D+00, & 0.6769514072983150D+00, & 0.9706224286053016D+00, & 0.9706224286053016D+00, & 0.0293775713946984D+00, & 0.9706224286053016D+00, & 0.3230485927016850D+00, & 0.0293775713946984D+00, & 0.0293775713946984D+00, & 0.0293775713946984D+00, & 0.9706224286053016D+00, & 0.6769514072983150D+00, & 0.9706224286053016D+00, & 0.6769514072983150D+00, & 0.6769514072983150D+00, & 0.3230485927016850D+00, & 0.9706224286053016D+00, & 0.0293775713946984D+00, & 0.9706224286053016D+00, & 0.3230485927016850D+00 /) real ( kind = rk ) y(n) real ( kind = rk ), save, dimension ( n_save ) :: y_save = (/ & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.9063071670498133D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.0936928329501868D+00, & 0.1991623679008687D+00, & 0.8008376320991313D+00, & 0.8008376320991313D+00, & 0.1991623679008687D+00, & 0.1991623679008687D+00, & 0.8008376320991313D+00, & 0.8008376320991313D+00, & 0.1991623679008687D+00, & 0.0722721194911200D+00, & 0.9277278805088800D+00, & 0.9277278805088800D+00, & 0.0722721194911200D+00, & 0.0722721194911200D+00, & 0.9277278805088800D+00, & 0.9277278805088800D+00, & 0.0722721194911200D+00, & 0.3433032977419727D+00, & 0.6566967022580273D+00, & 0.6566967022580273D+00, & 0.3433032977419727D+00, & 0.3433032977419727D+00, & 0.6566967022580273D+00, & 0.6566967022580273D+00, & 0.3433032977419727D+00, & 0.1326658565014960D+00, & 0.5000000000000000D+00, & 0.8673341434985040D+00, & 0.8673341434985040D+00, & 0.5000000000000000D+00, & 0.1326658565014960D+00, & 0.1326658565014960D+00, & 0.5000000000000000D+00, & 0.8673341434985040D+00, & 0.1326658565014960D+00, & 0.8673341434985040D+00, & 0.5000000000000000D+00, & 0.9825498327563551D+00, & 0.0174501672436448D+00, & 0.9825498327563551D+00, & 0.2746000324427453D+00, & 0.2746000324427453D+00, & 0.7253999675572547D+00, & 0.0174501672436448D+00, & 0.0174501672436448D+00, & 0.0174501672436448D+00, & 0.7253999675572547D+00, & 0.7253999675572547D+00, & 0.7253999675572547D+00, & 0.7253999675572547D+00, & 0.2746000324427453D+00, & 0.2746000324427453D+00, & 0.2746000324427453D+00, & 0.9825498327563551D+00, & 0.7253999675572547D+00, & 0.2746000324427453D+00, & 0.7253999675572547D+00, & 0.7253999675572547D+00, & 0.9825498327563551D+00, & 0.2746000324427453D+00, & 0.2746000324427453D+00, & 0.6769514072983150D+00, & 0.3230485927016850D+00, & 0.6769514072983150D+00, & 0.0293775713946984D+00, & 0.0293775713946984D+00, & 0.9706224286053016D+00, & 0.3230485927016850D+00, & 0.3230485927016850D+00, & 0.3230485927016850D+00, & 0.9706224286053016D+00, & 0.9706224286053016D+00, & 0.9706224286053016D+00, & 0.9706224286053016D+00, & 0.0293775713946984D+00, & 0.0293775713946984D+00, & 0.0293775713946984D+00, & 0.6769514072983150D+00, & 0.9706224286053016D+00, & 0.0293775713946984D+00, & 0.9706224286053016D+00, & 0.9706224286053016D+00, & 0.6769514072983150D+00, & 0.0293775713946984D+00, & 0.0293775713946984D+00 /) real ( kind = rk ) z(n) real ( kind = rk ), save, dimension ( n_save ) :: z_save = (/ & 0.5000000000000000D+00, & 0.9063071670498133D+00, & 0.5000000000000000D+00, & 0.0936928329501868D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.1991623679008687D+00, & 0.8008376320991313D+00, & 0.1991623679008687D+00, & 0.1991623679008687D+00, & 0.8008376320991313D+00, & 0.1991623679008687D+00, & 0.8008376320991313D+00, & 0.8008376320991313D+00, & 0.0722721194911200D+00, & 0.9277278805088800D+00, & 0.0722721194911200D+00, & 0.0722721194911200D+00, & 0.9277278805088800D+00, & 0.0722721194911200D+00, & 0.9277278805088800D+00, & 0.9277278805088800D+00, & 0.3433032977419727D+00, & 0.6566967022580273D+00, & 0.3433032977419727D+00, & 0.3433032977419727D+00, & 0.6566967022580273D+00, & 0.3433032977419727D+00, & 0.6566967022580273D+00, & 0.6566967022580273D+00, & 0.5000000000000000D+00, & 0.1326658565014960D+00, & 0.1326658565014960D+00, & 0.5000000000000000D+00, & 0.8673341434985040D+00, & 0.8673341434985040D+00, & 0.1326658565014960D+00, & 0.8673341434985040D+00, & 0.5000000000000000D+00, & 0.5000000000000000D+00, & 0.8673341434985040D+00, & 0.1326658565014960D+00, & 0.2746000324427453D+00, & 0.7253999675572547D+00, & 0.2746000324427453D+00, & 0.7253999675572547D+00, & 0.9825498327563551D+00, & 0.2746000324427453D+00, & 0.2746000324427453D+00, & 0.7253999675572547D+00, & 0.2746000324427453D+00, & 0.0174501672436448D+00, & 0.2746000324427453D+00, & 0.0174501672436448D+00, & 0.9825498327563551D+00, & 0.0174501672436448D+00, & 0.9825498327563551D+00, & 0.2746000324427453D+00, & 0.7253999675572547D+00, & 0.7253999675572547D+00, & 0.7253999675572547D+00, & 0.7253999675572547D+00, & 0.9825498327563551D+00, & 0.7253999675572547D+00, & 0.0174501672436448D+00, & 0.2746000324427453D+00, & 0.0293775713946984D+00, & 0.9706224286053016D+00, & 0.0293775713946984D+00, & 0.9706224286053016D+00, & 0.6769514072983150D+00, & 0.0293775713946984D+00, & 0.0293775713946984D+00, & 0.9706224286053016D+00, & 0.0293775713946984D+00, & 0.3230485927016850D+00, & 0.0293775713946984D+00, & 0.3230485927016850D+00, & 0.6769514072983150D+00, & 0.3230485927016850D+00, & 0.6769514072983150D+00, & 0.0293775713946984D+00, & 0.9706224286053016D+00, & 0.9706224286053016D+00, & 0.9706224286053016D+00, & 0.9706224286053016D+00, & 0.6769514072983150D+00, & 0.9706224286053016D+00, & 0.3230485927016850D+00, & 0.0293775713946984D+00 /) real ( kind = rk ) w(n) real ( kind = rk ), save, dimension ( n_save ) :: w_save = (/ & 0.0253096342016000D+00, & 0.0253096342016000D+00, & 0.0253096342016000D+00, & 0.0253096342016000D+00, & 0.0253096342016000D+00, & 0.0253096342016000D+00, & 0.0146922934945570D+00, & 0.0146922934945570D+00, & 0.0146922934945570D+00, & 0.0146922934945570D+00, & 0.0146922934945570D+00, & 0.0146922934945570D+00, & 0.0146922934945570D+00, & 0.0146922934945570D+00, & 0.0055804890098537D+00, & 0.0055804890098537D+00, & 0.0055804890098537D+00, & 0.0055804890098537D+00, & 0.0055804890098537D+00, & 0.0055804890098537D+00, & 0.0055804890098537D+00, & 0.0055804890098537D+00, & 0.0269990056568711D+00, & 0.0269990056568711D+00, & 0.0269990056568711D+00, & 0.0269990056568711D+00, & 0.0269990056568711D+00, & 0.0269990056568711D+00, & 0.0269990056568711D+00, & 0.0269990056568711D+00, & 0.0181499182325145D+00, & 0.0181499182325145D+00, & 0.0181499182325145D+00, & 0.0181499182325145D+00, & 0.0181499182325145D+00, & 0.0181499182325145D+00, & 0.0181499182325145D+00, & 0.0181499182325145D+00, & 0.0181499182325145D+00, & 0.0181499182325145D+00, & 0.0181499182325145D+00, & 0.0181499182325145D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0076802492622294D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00, & 0.0028267870173527D+00 /) x = x_save y = y_save z = z_save w = w_save return end