30318764346

Committed 28 Jul 2026 12:55AM UTC coverage: 76.972% (+1.1%) from 75.909%

Build # 30318764346

Build Type

Pull #166

github

Committed by

web-flow

Commit Message

Merge 3346edfd3 into 955f92dca

Pull Request Pull Request #166: Add L2/H1 boundary optimization

Coverage Stats

1631 of 1964 new or added lines in 40 files covered. (83.04%)

48 existing lines in 4 files now uncovered.

9563 of 12424 relevant lines covered (76.97%)

634449.32 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

94.86

/Source/BoundaryOptimization/LegendreAlgorithms.f90

! MIT License
!
! Copyright (c) 2010-present David A. Kopriva and other contributors: AUTHORS.md
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!
! FTObjectLibrary contains code that, to the best of our knowledge, has been released as
! public domain software:
! * `b3hs_hash_key_jenkins`: originally by Rich Townsend,
! https://groups.google.com/forum/#!topic/comp.lang.fortran/RWoHZFt39ng, 2005
!
! --- End License
!
!////////////////////////////////////////////////////////////////////////
!
!      LegendreAlgorithms.f90
!      Created: June 3, 2025 at 1:57 PM
!      By: David Kopriva
!
!////////////////////////////////////////////////////////////////////////
!
   Module LegendreAlgorithms
      USE ProgramGlobals
      IMPLICIT NONE

      INTEGER, PARAMETER     :: LA_EVALUATE_FUNCTION = 0, LA_EVALUATE_DERIVATIVE = 1

      REAL(KIND=RP), PRIVATE :: testCoefs(0:4) = [2.2_RP, 1.3_RP, 4.7_RP, 3.14_RP, 0.7_RP]
!
!     ========
      CONTAINS
!     ========
!
!      Main entry points:
!
!      SUBROUTINE LegendreLobattoNodesAndWeights( N, x, w )
!         Computes the Lobatto points,x, and weights, w, for polynomial order N
!
!      SUBROUTINE LegendrePolyAndDerivative( N, x, L_N, LPrime_N )
!         Compute the Legendre polynomial of degree N and its derivative at a point x
!
!      SUBROUTINE LegendrePolysAndDerivatives( N, x, L, LPrime )
!         Same as the above, but returns the arrays of L_k(x) and L'_k(x), k = 0, ..., N
!
!      REAL(KIND=RP) FUNCTION LegendreSeries( x, N, coefs, which )
!         Given the coefficients computes either \sum_{k=0}^N coefs(k)L_k(x) for which = LA_EVALUATE_FUNCTION
!         or the derivative \sum_{k=0}^N coefs(k)L'_k(x) for which = LA_EVALUATE_DERIVATIVE
!
!      FUNCTION gaussInnerProduct(a, b, weights, N)
!         Compute the Gauss-Lobatto innerproduct <a,b>_N for the weights compute from
!         LegendreLobattoNodesAndWeights
!
!      SUBROUTINE LegendreCoefficients( N, vals, coefs)
!         Given an array of values, compute the Lagrange modal coefficients
!
!
!      Testing Routines:
!
!      LOGICAL FUNCTION legendreQuadratureIsOK()
!      LOGICAL FUNCTION gaussInnerProductIsOK()
!      LOGICAL FUNCTION LegendrePolyAndDerivIsOK()
!      LOGICAL FUNCTION legendreCoefsForTestAreOK()
!      LOGICAL FUNCTION legendreSeriesIsOK()
!
!////////////////////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE LegendrePolysAndDerivatives( N, x, L, LPrime )
!
!     ---------------------------------------------------------------------
!     Compute the Legendre Polynomials up to degree N and their derivatives
!     at the point x.
!     ---------------------------------------------------------------------
!
!     -----------------
!     Input parameters:
!     -----------------
!
      INTEGER      , INTENT(IN) :: N
      REAL(KIND=RP), INTENT(IN) :: x
!
!     ------------------
!     Output parameters:
!     ------------------
!
      REAL(KIND=RP), INTENT(OUT) :: L(0:N), LPrime(0:N)
!
!     ----------------
!     Local Variables:
!     ----------------
!
      INTEGER       :: k

      L(0)      = 1.0_rp
      LPrime(0) = 0.0_rp
      IF(N==0)    RETURN

      L(1)      = x
      LPrime(1) = 1.0_rp
      IF( N == 1)  RETURN

      DO k = 2, N
         L(k)        = ((2*k-1)*x*L(k-1) - (k-1)*L(k-2))/k
         LPrime(k)   = LPrime(k-2) + (2*k-1)*L(k-1)
      END DO

      END SUBROUTINE LegendrePolysAndDerivatives
!
!////////////////////////////////////////////////////////////////////////
!
      REAL(KIND=RP) FUNCTION LegendreSeries( x, N, coefs, which )
!
!     --------------------------------------------------------------------
!     Compute the value of the Legendre series \sum_{k=0}^N coefs_k \phi_k
!     for \phi_k = L_k if which = LA_EVALUATE_FUNCTION and \phi_k = L'_k
!     if it is LA_EVALUATE_DERIVATIVE
!     --------------------------------------------------------------------
!
         IMPLICIT NONE
         INTEGER      , INTENT(IN) :: N
         REAL(KIND=RP), INTENT(IN) :: x
         REAL(KIND=RP), INTENT(IN) :: coefs(0:N)
         INTEGER                   :: which

         REAL(KIND=RP) :: L(0:N), LPrime(0:N)

         CALL LegendrePolysAndDerivatives( N, x, L, LPrime )

         IF ( which == LA_EVALUATE_FUNCTION  )     THEN
            LegendreSeries = SumSeries( N, coefs, L)
         ELSE
            LegendreSeries = SumSeries( N, coefs, LPrime)
         END IF

      END FUNCTION LegendreSeries
!
!////////////////////////////////////////////////////////////////////////
!
      REAL(KIND=RP) FUNCTION SumSeries( N, coefs, basis )
         IMPLICIT NONE
         INTEGER      , INTENT(IN) :: N
         REAL(KIND=RP), INTENT(IN) :: coefs(0:N), basis(0:N)

         REAL(KIND=RP) :: s
         INTEGER       :: k

         s = 0.0_RP
         DO k = 0, N
            s = s +  coefs(k)*basis(k)
         END DO
         SumSeries = s

      END FUNCTION SumSeries
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION gaussInnerProduct(a, b, weights, N)
!
!     -----------------------------------------
!     Compute the discrete inner product
!     <a,b>_N = \sum_{k=0}^N a(k)b(k)weights(k)
!     -----------------------------------------
!
         IMPLICIT NONE
         INTEGER :: N
         REAL(KIND=RP) :: a(0:N), b(0:N), weights(0:N)
         REAL(KIND=RP) :: gaussInnerProduct

         INTEGER       :: j
         REAL(KIND=RP) :: s

         s = 0.0_RP
         DO j = 0,N
            s = s + a(j)*b(j)*weights(j)
         END DO

         gaussInnerProduct = s

      END FUNCTION gaussInnerProduct
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION gaussQuadrature(a, weights, N)
!
!     -----------------------------------------
!     Compute the discrete inner product
!     <a,b>_N = \sum_{k=0}^N a(k)b(k)weights(k)
!     -----------------------------------------
!
         IMPLICIT NONE
         INTEGER :: N
         REAL(KIND=RP) :: a(0:N), weights(0:N)
         REAL(KIND=RP) :: gaussQuadrature

         INTEGER       :: j
         REAL(KIND=RP) :: s

         s = 0.0_RP
         DO j = 0,N
            s = s + a(j)*weights(j)
         END DO

         gaussQuadrature = s

      END FUNCTION gaussQuadrature
!
!////////////////////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE GaussLegendreNodesAndWeights( N, x, w )
!
!     Compute the Gauss-legendre quadrature nodes and
!     weights
!
!     -----------------
!     Input parameters:
!     -----------------
!
      INTEGER, INTENT(IN) :: N
!
!     ------------------
!     Output parameters:
!     ------------------
!
      REAL(KIND=RP), DIMENSION(0:N), INTENT(OUT) :: x, w
!
!     ----------------
!     Local Variables:
!     ----------------
!
      REAL(KIND=RP) :: xj, L_NP1, LPrime_NP1, delta, tolerance
      INTEGER       :: j, k, NDiv2
!
!     ----------
!     Constants:
!     ----------
!
      INTEGER, PARAMETER       :: numNewtonIterations = 10
      REAL(KIND=RP), PARAMETER :: toleranceFactor     = 4.0_RP

      tolerance = toleranceFactor*EPSILON(L_NP1)
      IF( N == 0 )     THEN
         x(0) = 0.0_RP
         w(0) = 2.0_RP
         RETURN
      ELSE IF( N == 1 )     THEN
         x(0) = -SQRT(1.0_RP/3.0_RP)
         w(0) =  1.0_RP
         x(1) = -x(0)
         w(1) =  w(0)
      ELSE
!
!        ----------------------------------
!        Iterate on half the interior nodes
!        ----------------------------------
!
         NDiv2 = (N+1)/2
         DO j = 0, NDiv2-1
            xj = -COS( (2*j+1)*PI/(2*N+2) )
            DO k = 0, numNewtonIterations
               CALL LegendrePolyAndDerivative( N+1, xj, L_NP1, LPrime_NP1 )
               delta = -L_NP1/LPrime_NP1
               xj = xj + delta
               IF( ABS(delta) <=  tolerance*ABS(xj) )     EXIT
            END DO
            CALL LegendrePolyAndDerivative( N+1, xj, L_NP1, LPrime_NP1 )
            x(j)   = xj
            w(j)   = 2.0_RP/( (1.0_RP - xj**2)*LPrime_NP1**2 )
            x(N-j) = -xj
            w(N-j) = w(j)
         END DO
      END IF
!
!     ---------------------------
!     Fill in middle if necessary
!     ---------------------------
!
      IF( MOD(N,2) == 0 )     THEN
         CALL LegendrePolyAndDerivative( N+1, 0.0_RP, L_NP1, LPrime_NP1 )
         x(N/2) = 0.0_RP
         w(N/2) = 2.0_RP/LPrime_NP1**2
      END IF

      END SUBROUTINE GaussLegendreNodesAndWeights
!
!////////////////////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE LegendrePolyAndDerivative( N, x, L_N, LPrime_N )
!
!     --------------------------------------------------------------
!     Compute the Legendre Polynomial of degree k and its derivative
!     --------------------------------------------------------------
!
!     -----------------
!     Input parameters:
!     -----------------
!
      INTEGER      , INTENT(IN) :: N
      REAL(KIND=RP), INTENT(IN) :: x
!
!     ------------------
!     Output parameters:
!     ------------------
!
      REAL(KIND=RP), INTENT(OUT) :: L_N, LPrime_N
!
!     ----------------
!     Local Variables:
!     ----------------
!
      INTEGER       :: k
      REAL(KIND=RP) :: L_NM1, L_NM2, LPrime_NM2, LPrime_NM1

      IF( N == 0 )     THEN
         L_N      = 1.0_rp
         LPrime_N = 0.0_rp
      ELSE IF ( N == 1 )     THEN
         L_N      = x
         LPrime_N = 1.0_rp
      ELSE
         L_NM2 = 1.0_rp
         LPrime_NM2 = 0.0_rp
         L_NM1 = x
         LPrime_NM1 = 1.0_rp
         DO k = 2, N
            L_N        = ((2*k-1)*x*L_NM1 - (k-1)*L_NM2)/k
            LPrime_N   = LPrime_NM2 + (2*k-1)*L_NM1
            L_NM2      = L_NM1
            L_NM1      = L_N
            LPrime_NM2 = LPrime_NM1
            LPrime_NM1 = LPrime_N
         END DO
      END IF

      END SUBROUTINE LegendrePolyAndDerivative
!
!     -------
!     Testing
!     -------
!
!
!////////////////////////////////////////////////////////////////////////
!
      LOGICAL FUNCTION legendreQuadratureIsOK()
!
!     ---------------------------------------------------------------
!     The Legendre quadrature will be exact for polynomials of degree
!     2N-1. Checks up to N = 20, but we know it should be good at
!     least to N = 2000 from other tests.
!     ---------------------------------------------------------------
!
         IMPLICIT NONE
         REAL(KIND=RP), DIMENSION(:), ALLOCATABLE :: x, w, f, coefs
         INTEGER                                  :: N, j, m, nMax
         REAL(KIND=RP)                            :: sm, fInt, exact
         REAL(KIND=RP)                            :: maxError, maxSumError
         REAL(KIND=RP)                            :: tol = 1.0d-9

         nMax = 20
         ALLOCATE(x(0:nMax),w(0:nMax), f(0:nMax), coefs(0:2*nMax-1))

         maxError    = -1.0_RP
         maxSumError = -1.0_RP

         legendreQuadratureIsOK = .TRUE.

         DO N = 2, nMax

            CALL GaussLegendreNodesAndWeights( N, x, w )
            CALL RANDOM_NUMBER(HARVEST = coefs)
!
!           --------------------------------
!           Check sum of weights should be 2
!           --------------------------------
!
            sm = 0.0_RP
            DO j = 0, N
               sm = sm + w(j)
            END DO
            maxSumError = MAX(maxSumError,ABS(sm-2.0_RP))
!
!           ----------------------------------------------------------
!           Lobatto quadrature is exact for polynomials to degree 2N-1
!           ----------------------------------------------------------
!
            DO j = 0,N
               f(j) = 0.0_RP
               DO m = 0, 2*N-1
                  f(j) = f(j) + coefs(m)*x(j)**m
               END DO
            END DO
!
!           ---------------------------
!           Exact value of the integral
!           ---------------------------
!
            exact = 0.0_RP
            DO m = 0, 2*N-1
               exact = exact + coefs(m)*(1.0_RP - (-1.0)**(m+1))/(m+1)
            END DO
!
!           ----------------------
!           Compute the quadrature
!           ----------------------
!
            fInt = gaussQuadrature(f, w, N)

            maxError = MAX(maxError, ABS(fInt-exact))

            legendreQuadratureIsOK = legendreQuadratureIsOK .AND. maxError < tol .AND. maxSumError < tol

         END DO

      END FUNCTION legendreQuadratureIsOK
!
!////////////////////////////////////////////////////////////////////////
!
      LOGICAL FUNCTION gaussInnerProductIsOK()
        IMPLICIT NONE
        INTEGER, PARAMETER :: N = 4
        REAL(KIND=RP)      :: a(0:N), b(0:N)
        REAL(KIND=RP)      :: nodes(0:N), weights(0:N)
        REAL(KIND=RP)      :: nrm, tol = 1.0d-8
        INTEGER            :: j

        CALL GaussLegendreNodesAndWeights( N, nodes, weights )
!
!       --------------------------------------
!       Orthogonal polynomials shall give zero
!       --------------------------------------
!
        DO j = 0,N
           a(j) = P3(nodes(j))
           b(j) = P2(nodes(j))
        END DO
        nrm = gaussInnerProduct(a,b,weights,N)
        gaussInnerProductIsOK = nrm < tol
!
!       -------------------------------
!       Same polynomial shall give norm
!       -------------------------------
!
        DO j = 0,N
           a(j) = P3(nodes(j))
           b(j) = a(j)
        END DO
        nrm = gaussInnerProduct(a,b,weights,N)
        gaussInnerProductIsOK = ABS(nrm - 2.0_RP/(2.0_RP*3.0_RP+1.0_RP)) < tol &
                                .AND. gaussInnerProductIsOK

      END FUNCTION gaussInnerProductIsOK
!
!////////////////////////////////////////////////////////////////////////
!
      LOGICAL FUNCTION LegendrePolyAndDerivIsOK()
         IMPLICIT NONE
         INTEGER      , PARAMETER  :: N = 4
         REAL(KIND=RP)             :: x, dx, tol = 1.0d-8
         REAL(KIND=RP)             :: L(0:N)     , LPrime(0:N)
         REAL(KIND=RP)             :: LExact(0:N), LPrimeExact(0:N)
         INTEGER                   :: j

         dx = 2.0_RP/5.0_RP
         LegendrePolyAndDerivIsOK = .TRUE.
         DO j = 0, 5
            x = -1.0_RP + j*dx
            CALL LegendrePolysAndDerivatives(N,x,L,LPrime)

            LExact(0) = P0(x)
            LExact(1) = P1(x)
            LExact(2) = P2(x)
            LExact(3) = P3(x)
            LExact(4) = P4(x)

            LegendrePolyAndDerivIsOK = MAXVAL(ABS(L - LExact)) < tol &
                                       .AND. LegendrePolyAndDerivIsOK

            LPrimeExact(0) = P0Prime(x)
            LPrimeExact(1) = P1Prime(x)
            LPrimeExact(2) = P2Prime(x)
            LPrimeExact(3) = P3Prime(x)
            LPrimeExact(4) = P4Prime(x)

            LegendrePolyAndDerivIsOK = MAXVAL(ABS(LPrime - LPrimeExact)) < tol &
                                       .AND. LegendrePolyAndDerivIsOK
         END DO

      END FUNCTION LegendrePolyAndDerivIsOK
!
!////////////////////////////////////////////////////////////////////////
!
      LOGICAL FUNCTION legendreCoefsForTestAreOK()
!
!     -------------------------------------------------------------------
!     The Legendre series will be exact for any polynomial of degree <= N
!     -------------------------------------------------------------------
!
         IMPLICIT NONE
         INTEGER, PARAMETER :: N = 4
         REAL(KIND=RP)      :: tol = 1.0d-8
         REAL(KIND=RP)      :: coefs(0:N)

         CALL LegendreCoefsForTestPolynomial(coefs)
         legendreCoefsForTestAreOK = MAXVAL(ABS(coefs - testCoefs)) < tol

      END FUNCTION legendreCoefsForTestAreOK
!
!////////////////////////////////////////////////////////////////////////
!
      LOGICAL FUNCTION legendreSeriesIsOK()
         IMPLICIT NONE
         REAL(KIND=RP) :: tol = 1.0d-8
         REAL(KIND=RP) :: x, dx
         INTEGER       :: j

         legendreSeriesIsOK = .TRUE.
         dx = 2.0_RP/5.0_RP
         DO j = 0, 5
            x = -1.0_RP + j*dx
            legendreSeriesIsOK = ABS(LegendreSeries(x,4,testCoefs,LA_EVALUATE_FUNCTION) &
                                    - testPolynomial(x)) < tol &
                                    .AND. legendreSeriesIsOK
         END DO

      END FUNCTION legendreSeriesIsOK
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE LegendreCoefsForTestPolynomial(coefs)
         IMPLICIT NONE
         INTEGER, PARAMETER :: N = 6
         REAL(KIND=RP)      :: coefs(0:4)
         REAL(KIND=RP)      :: a(0:N), b(0:N)
         REAL(KIND=RP)      :: nodes(0:N), weights(0:N)
         INTEGER            :: j,k

        CALL GaussLegendreNodesAndWeights( N, nodes, weights )
!
!       -----------------------
!       Select the coefficients
!       -----------------------
!
        k = 0
        DO j = 0,N
           a(j) = testPolynomial(nodes(j))
           b(j) = P0(nodes(j))
        END DO
        coefs(k) = gaussInnerProduct(a,b,weights,N)*(2.0_RP*k + 1.0_RP)/2.0_RP
        k = 1
        DO j = 0,N
           b(j) = P1(nodes(j))
        END DO
        coefs(k) = gaussInnerProduct(a,b,weights,N)*(2.0_RP*k + 1.0_RP)/2.0_RP
        k = 2
        DO j = 0,N
           b(j) = P2(nodes(j))
        END DO
        coefs(k) = gaussInnerProduct(a,b,weights,N)*(2.0_RP*k + 1.0_RP)/2.0_RP
        k = 3
        DO j = 0,N
           b(j) = P3(nodes(j))
        END DO
        coefs(k) = gaussInnerProduct(a,b,weights,N)*(2.0_RP*k + 1.0_RP)/2.0_RP
        k = 4
        DO j = 0,N
           b(j) = P4(nodes(j))
        END DO
        coefs(k) = gaussInnerProduct(a,b,weights,N)*(2.0_RP*k + 1.0_RP)/2.0_RP

      END SUBROUTINE LegendreCoefsForTestPolynomial
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION testPolynomial(x) RESULT(s)
!
!     ---------------------------------
!     Evaluate a polynomial of degree 4
!     ---------------------------------
!
         IMPLICIT NONE
         REAL(KIND=RP) :: x
         REAL(KIND=RP) :: P(0:4), s
         INTEGER       :: k

         P(0) = P0(x)
         P(1) = P1(x)
         P(2) = P2(x)
         P(3) = P3(x)
         P(4) = P4(x)

         s = 0.0_RP
         DO k = 0, 4
            s = s + testCoefs(k)*P(k)
         END DO

      END FUNCTION testPolynomial
!
!     -------------------------
!     Some Legendre polynomials
!     -------------------------
!
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION P0(x)
         IMPLICIT NONE
         REAL(KIND=RP) :: x, P0
         P0 = 1.0_RP
      END FUNCTION P0
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION P0Prime(x)
         IMPLICIT NONE
         REAL(KIND=RP) :: x, P0Prime
         P0Prime = 0.0_RP
      END FUNCTION P0Prime
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION P1(x)
         IMPLICIT NONE
         REAL(KIND=RP) :: x, P1
         P1 = x
      END FUNCTION P1
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION P1Prime(x)
         IMPLICIT NONE
         REAL(KIND=RP) :: x, P1Prime
         P1Prime = 1.0_RP
      END FUNCTION P1Prime
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION P2(x)
         IMPLICIT NONE
         REAL(KIND=RP) :: x, P2
         P2 = 0.5_RP*(3.0_RP*x**2 - 1.0_RP)
      END FUNCTION P2
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION P2Prime(x)
         IMPLICIT NONE
         REAL(KIND=RP) :: x, P2Prime
         P2Prime = 0.5_RP*(2.0_RP*3.0_RP*x)
      END FUNCTION P2Prime
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION P3(x)
         IMPLICIT NONE
         REAL(KIND=RP) :: x, P3
         P3 = 0.5_RP*(5.0_RP*x**3 - 3.0_RP*x)
      END FUNCTION P3
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION P3Prime(x)
         IMPLICIT NONE
         REAL(KIND=RP) :: x, P3Prime
         P3Prime = 0.5_RP*(3.0_RP*5.0_RP*x**2 - 3.0_RP)
      END FUNCTION P3Prime
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION P4(x)
         IMPLICIT NONE
         REAL(KIND=RP) :: x, P4
         P4 = 0.125_RP*(35.0_RP*x**4 - 30.0_RP*x**2 + 3.0_RP)
      END FUNCTION P4
!
!////////////////////////////////////////////////////////////////////////
!
      FUNCTION P4Prime(x)
         IMPLICIT NONE
         REAL(KIND=RP) :: x, P4Prime
         P4Prime = 0.125_RP*(4.0_RP*35.0_RP*x**3 - 2.0_RP*30.0_RP*x)
      END FUNCTION P4Prime
   END Module LegendreAlgorithms

1	! MIT License
2	!
3	! Copyright (c) 2010-present David A. Kopriva and other contributors: AUTHORS.md
4	!
5	! Permission is hereby granted, free of charge, to any person obtaining a copy
6	! of this software and associated documentation files (the "Software"), to deal
7	! in the Software without restriction, including without limitation the rights
8	! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9	! copies of the Software, and to permit persons to whom the Software is
10	! furnished to do so, subject to the following conditions:
11	!
12	! The above copyright notice and this permission notice shall be included in all
13	! copies or substantial portions of the Software.
14	!
15	! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16	! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17	! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18	! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19	! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20	! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
21	! SOFTWARE.
22	!
23	! FTObjectLibrary contains code that, to the best of our knowledge, has been released as
24	! public domain software:
25	! * `b3hs_hash_key_jenkins`: originally by Rich Townsend,
26	! https://groups.google.com/forum/#!topic/comp.lang.fortran/RWoHZFt39ng, 2005
27	!
28	! --- End License
29	!
30	!////////////////////////////////////////////////////////////////////////
31	!
32	! LegendreAlgorithms.f90
33	! Created: June 3, 2025 at 1:57 PM
34	! By: David Kopriva
35	!
36	!////////////////////////////////////////////////////////////////////////
37	!
38	Module LegendreAlgorithms
39	USE ProgramGlobals
40	IMPLICIT NONE
41
42	INTEGER, PARAMETER :: LA_EVALUATE_FUNCTION = 0, LA_EVALUATE_DERIVATIVE = 1
43
44	REAL(KIND=RP), PRIVATE :: testCoefs(0:4) = [2.2_RP, 1.3_RP, 4.7_RP, 3.14_RP, 0.7_RP]
45	!
46	! ========
47	CONTAINS
48	! ========
49	!
50	! Main entry points:
51	!
52	! SUBROUTINE LegendreLobattoNodesAndWeights( N, x, w )
53	! Computes the Lobatto points,x, and weights, w, for polynomial order N
54	!
55	! SUBROUTINE LegendrePolyAndDerivative( N, x, L_N, LPrime_N )
56	! Compute the Legendre polynomial of degree N and its derivative at a point x
57	!
58	! SUBROUTINE LegendrePolysAndDerivatives( N, x, L, LPrime )
59	! Same as the above, but returns the arrays of L_k(x) and L'_k(x), k = 0, ..., N
60	!
61	! REAL(KIND=RP) FUNCTION LegendreSeries( x, N, coefs, which )
62	! Given the coefficients computes either \sum_{k=0}^N coefs(k)L_k(x) for which = LA_EVALUATE_FUNCTION
63	! or the derivative \sum_{k=0}^N coefs(k)L'_k(x) for which = LA_EVALUATE_DERIVATIVE
64	!
65	! FUNCTION gaussInnerProduct(a, b, weights, N)
66	! Compute the Gauss-Lobatto innerproduct <a,b>_N for the weights compute from
67	! LegendreLobattoNodesAndWeights
68	!
69	! SUBROUTINE LegendreCoefficients( N, vals, coefs)
70	! Given an array of values, compute the Lagrange modal coefficients
71	!
72	!
73	! Testing Routines:
74	!
75	! LOGICAL FUNCTION legendreQuadratureIsOK()
76	! LOGICAL FUNCTION gaussInnerProductIsOK()
77	! LOGICAL FUNCTION LegendrePolyAndDerivIsOK()
78	! LOGICAL FUNCTION legendreCoefsForTestAreOK()
79	! LOGICAL FUNCTION legendreSeriesIsOK()
80	!
81	!////////////////////////////////////////////////////////////////////////////////////////
82	!
83	SUBROUTINE LegendrePolysAndDerivatives( N, x, L, LPrime )	63,124✔
84	!
85	! ---------------------------------------------------------------------
86	! Compute the Legendre Polynomials up to degree N and their derivatives
87	! at the point x.
88	! ---------------------------------------------------------------------
89	!
90	! -----------------
91	! Input parameters:
92	! -----------------
93	!
94	INTEGER , INTENT(IN) :: N
95	REAL(KIND=RP), INTENT(IN) :: x
96	!
97	! ------------------
98	! Output parameters:
99	! ------------------
100	!
101	REAL(KIND=RP), INTENT(OUT) :: L(0:N), LPrime(0:N)
102	!
103	! ----------------
104	! Local Variables:
105	! ----------------
106	!
107	INTEGER :: k
108
109	L(0) = 1.0_rp	63,124✔
110	LPrime(0) = 0.0_rp	63,124✔
111	IF(N==0) RETURN	63,124✔
112
113	L(1) = x	63,124✔
114	LPrime(1) = 1.0_rp	63,124✔
115	IF( N == 1) RETURN	63,124✔
116
117	DO k = 2, N	369,293✔
118	L(k) = ((2k-1)xL(k-1) - (k-1)L(k-2))/k	306,169✔
119	LPrime(k) = LPrime(k-2) + (2k-1)L(k-1)	306,169✔
120	END DO
121
122	END SUBROUTINE LegendrePolysAndDerivatives
123	!
124	!////////////////////////////////////////////////////////////////////////
125	!
126	REAL(KIND=RP) FUNCTION LegendreSeries( x, N, coefs, which )	56,722✔
127	!
128	! --------------------------------------------------------------------
129	! Compute the value of the Legendre series \sum_{k=0}^N coefs_k \phi_k
130	! for \phi_k = L_k if which = LA_EVALUATE_FUNCTION and \phi_k = L'_k
131	! if it is LA_EVALUATE_DERIVATIVE
132	! --------------------------------------------------------------------
133	!
134	IMPLICIT NONE
135	INTEGER , INTENT(IN) :: N
136	REAL(KIND=RP), INTENT(IN) :: x
137	REAL(KIND=RP), INTENT(IN) :: coefs(0:N)
138	INTEGER :: which
139
140	REAL(KIND=RP) :: L(0:N), LPrime(0:N)	113,444✔
141
142	CALL LegendrePolysAndDerivatives( N, x, L, LPrime )	56,722✔
143
144	IF ( which == LA_EVALUATE_FUNCTION ) THEN	56,722✔
145	LegendreSeries = SumSeries( N, coefs, L)	29,882✔
146	ELSE
147	LegendreSeries = SumSeries( N, coefs, LPrime)	26,840✔
148	END IF
149
150	END FUNCTION LegendreSeries	56,722✔
151	!
152	!////////////////////////////////////////////////////////////////////////
153	!
154	REAL(KIND=RP) FUNCTION SumSeries( N, coefs, basis )	56,722✔
155	IMPLICIT NONE
156	INTEGER , INTENT(IN) :: N
157	REAL(KIND=RP), INTENT(IN) :: coefs(0:N), basis(0:N)
158
159	REAL(KIND=RP) :: s
160	INTEGER :: k
161
162	s = 0.0_RP	56,722✔
163	DO k = 0, N	445,342✔
164	s = s + coefs(k)*basis(k)	445,342✔
165	END DO
166	SumSeries = s	56,722✔
167
168	END FUNCTION SumSeries	56,722✔
169	!
170	!////////////////////////////////////////////////////////////////////////
171	!
172	FUNCTION gaussInnerProduct(a, b, weights, N)	65,725✔
173	!
174	! -----------------------------------------
175	! Compute the discrete inner product
176	! <a,b>_N = \sum_{k=0}^N a(k)b(k)weights(k)
177	! -----------------------------------------
178	!
179	IMPLICIT NONE
180	INTEGER :: N
181	REAL(KIND=RP) :: a(0:N), b(0:N), weights(0:N)
182	REAL(KIND=RP) :: gaussInnerProduct
183
184	INTEGER :: j
185	REAL(KIND=RP) :: s
186
187	s = 0.0_RP	65,725✔
188	DO j = 0,N	905,728✔
189	s = s + a(j)b(j)weights(j)	905,728✔
190	END DO
191
192	gaussInnerProduct = s	65,725✔
193
194	END FUNCTION gaussInnerProduct	65,725✔
195	!
196	!////////////////////////////////////////////////////////////////////////
197	!
198	FUNCTION gaussQuadrature(a, weights, N)	19✔
199	!
200	! -----------------------------------------
201	! Compute the discrete inner product
202	! <a,b>_N = \sum_{k=0}^N a(k)b(k)weights(k)
203	! -----------------------------------------
204	!
205	IMPLICIT NONE
206	INTEGER :: N
207	REAL(KIND=RP) :: a(0:N), weights(0:N)
208	REAL(KIND=RP) :: gaussQuadrature
209
210	INTEGER :: j
211	REAL(KIND=RP) :: s
212
213	s = 0.0_RP	19✔
214	DO j = 0,N	247✔
215	s = s + a(j)*weights(j)	247✔
216	END DO
217
218	gaussQuadrature = s	19✔
219
220	END FUNCTION gaussQuadrature	19✔
221	!
222	!////////////////////////////////////////////////////////////////////////////////////////
223	!
224	SUBROUTINE GaussLegendreNodesAndWeights( N, x, w )	1,081✔
225	!
226	! Compute the Gauss-legendre quadrature nodes and
227	! weights
228	!
229	! -----------------
230	! Input parameters:
231	! -----------------
232	!
233	INTEGER, INTENT(IN) :: N
234	!
235	! ------------------
236	! Output parameters:
237	! ------------------
238	!
239	REAL(KIND=RP), DIMENSION(0:N), INTENT(OUT) :: x, w
240	!
241	! ----------------
242	! Local Variables:
243	! ----------------
244	!
245	REAL(KIND=RP) :: xj, L_NP1, LPrime_NP1, delta, tolerance
246	INTEGER :: j, k, NDiv2
247	!
248	! ----------
249	! Constants:
250	! ----------
251	!
252	INTEGER, PARAMETER :: numNewtonIterations = 10
253	REAL(KIND=RP), PARAMETER :: toleranceFactor = 4.0_RP
254
255	tolerance = toleranceFactor*EPSILON(L_NP1)	1,081✔
256	IF( N == 0 ) THEN	1,081✔
NEW 257	x(0) = 0.0_RP	×
NEW 258	w(0) = 2.0_RP	×
NEW 259	RETURN	×
260	ELSE IF( N == 1 ) THEN	1,081✔
NEW 261	x(0) = -SQRT(1.0_RP/3.0_RP)	×
NEW 262	w(0) = 1.0_RP	×
NEW 263	x(1) = -x(0)	×
NEW 264	w(1) = w(0)	×
265	ELSE
266	!
267	! ----------------------------------
268	! Iterate on half the interior nodes
269	! ----------------------------------
270	!
271	NDiv2 = (N+1)/2	1,081✔
272	DO j = 0, NDiv2-1	7,211✔
273	xj = -COS( (2j+1)PI/(2*N+2) )	6,130✔
274	DO k = 0, numNewtonIterations	29,592✔
275	CALL LegendrePolyAndDerivative( N+1, xj, L_NP1, LPrime_NP1 )	29,592✔
276	delta = -L_NP1/LPrime_NP1	29,592✔
277	xj = xj + delta	29,592✔
278	IF( ABS(delta) <= tolerance*ABS(xj) ) EXIT	29,592✔
279	END DO
280	CALL LegendrePolyAndDerivative( N+1, xj, L_NP1, LPrime_NP1 )	6,130✔
281	x(j) = xj	6,130✔
282	w(j) = 2.0_RP/( (1.0_RP - xj*2)LPrime_NP1**2 )	6,130✔
283	x(N-j) = -xj	6,130✔
284	w(N-j) = w(j)	7,211✔
285	END DO
286	END IF
287	!
288	! ---------------------------
289	! Fill in middle if necessary
290	! ---------------------------
291	!
292	IF( MOD(N,2) == 0 ) THEN	1,081✔
293	CALL LegendrePolyAndDerivative( N+1, 0.0_RP, L_NP1, LPrime_NP1 )	1,063✔
294	x(N/2) = 0.0_RP	1,063✔
295	w(N/2) = 2.0_RP/LPrime_NP1**2	1,063✔
296	END IF
297
298	END SUBROUTINE GaussLegendreNodesAndWeights
299	!
300	!////////////////////////////////////////////////////////////////////////////////////////
301	!
302	SUBROUTINE LegendrePolyAndDerivative( N, x, L_N, LPrime_N )	36,785✔
303	!
304	! --------------------------------------------------------------
305	! Compute the Legendre Polynomial of degree k and its derivative
306	! --------------------------------------------------------------
307	!
308	! -----------------
309	! Input parameters:
310	! -----------------
311	!
312	INTEGER , INTENT(IN) :: N
313	REAL(KIND=RP), INTENT(IN) :: x
314	!
315	! ------------------
316	! Output parameters:
317	! ------------------
318	!
319	REAL(KIND=RP), INTENT(OUT) :: L_N, LPrime_N
320	!
321	! ----------------
322	! Local Variables:
323	! ----------------
324	!
325	INTEGER :: k
326	REAL(KIND=RP) :: L_NM1, L_NM2, LPrime_NM2, LPrime_NM1
327
328	IF( N == 0 ) THEN	36,785✔
NEW 329	L_N = 1.0_rp	×
NEW 330	LPrime_N = 0.0_rp	×
331	ELSE IF ( N == 1 ) THEN	36,785✔
NEW 332	L_N = x	×
NEW 333	LPrime_N = 1.0_rp	×
334	ELSE
335	L_NM2 = 1.0_rp	36,785✔
336	LPrime_NM2 = 0.0_rp	36,785✔
337	L_NM1 = x	36,785✔
338	LPrime_NM1 = 1.0_rp	36,785✔
339	DO k = 2, N	470,033✔
340	L_N = ((2k-1)xL_NM1 - (k-1)L_NM2)/k	433,248✔
341	LPrime_N = LPrime_NM2 + (2k-1)L_NM1	433,248✔
342	L_NM2 = L_NM1	433,248✔
343	L_NM1 = L_N	433,248✔
344	LPrime_NM2 = LPrime_NM1	433,248✔
345	LPrime_NM1 = LPrime_N	470,033✔
346	END DO
347	END IF
348
349	END SUBROUTINE LegendrePolyAndDerivative	36,785✔
350	!
351	! -------
352	! Testing
353	! -------
354	!
355	!
356	!////////////////////////////////////////////////////////////////////////
357	!
358	LOGICAL FUNCTION legendreQuadratureIsOK()	1✔
359	!
360	! ---------------------------------------------------------------
361	! The Legendre quadrature will be exact for polynomials of degree
362	! 2N-1. Checks up to N = 20, but we know it should be good at
363	! least to N = 2000 from other tests.
364	! ---------------------------------------------------------------
365	!
366	IMPLICIT NONE
367	REAL(KIND=RP), DIMENSION(:), ALLOCATABLE :: x, w, f, coefs	1✔
368	INTEGER :: N, j, m, nMax
369	REAL(KIND=RP) :: sm, fInt, exact
370	REAL(KIND=RP) :: maxError, maxSumError
371	REAL(KIND=RP) :: tol = 1.0d-9
372
373	nMax = 20	1✔
374	ALLOCATE(x(0:nMax),w(0:nMax), f(0:nMax), coefs(0:2*nMax-1))	1✔
375
376	maxError = -1.0_RP	1✔
377	maxSumError = -1.0_RP	1✔
378
379	legendreQuadratureIsOK = .TRUE.	1✔
380
381	DO N = 2, nMax	20✔
382
383	CALL GaussLegendreNodesAndWeights( N, x, w )	19✔
384	CALL RANDOM_NUMBER(HARVEST = coefs)	19✔
385	!
386	! --------------------------------
387	! Check sum of weights should be 2
388	! --------------------------------
389	!
390	sm = 0.0_RP	19✔
391	DO j = 0, N	247✔
392	sm = sm + w(j)	247✔
393	END DO
394	maxSumError = MAX(maxSumError,ABS(sm-2.0_RP))	19✔
395	!
396	! ----------------------------------------------------------
397	! Lobatto quadrature is exact for polynomials to degree 2N-1
398	! ----------------------------------------------------------
399	!
400	DO j = 0,N	247✔
401	f(j) = 0.0_RP	228✔
402	DO m = 0, 2*N-1	6,403✔
403	f(j) = f(j) + coefs(m)x(j)*m	6,384✔
404	END DO
405	END DO
406	!
407	! ---------------------------
408	! Exact value of the integral
409	! ---------------------------
410	!
411	exact = 0.0_RP	19✔
412	DO m = 0, 2*N-1	437✔
413	exact = exact + coefs(m)(1.0_RP - (-1.0)*(m+1))/(m+1)	437✔
414	END DO
415	!
416	! ----------------------
417	! Compute the quadrature
418	! ----------------------
419	!
420	fInt = gaussQuadrature(f, w, N)	19✔
421
422	maxError = MAX(maxError, ABS(fInt-exact))	19✔
423
424	legendreQuadratureIsOK = legendreQuadratureIsOK .AND. maxError < tol .AND. maxSumError < tol	20✔
425
426	END DO
427
428	END FUNCTION legendreQuadratureIsOK	1✔
429	!
430	!////////////////////////////////////////////////////////////////////////
431	!
432	LOGICAL FUNCTION gaussInnerProductIsOK()	1✔
433	IMPLICIT NONE
434	INTEGER, PARAMETER :: N = 4
435	REAL(KIND=RP) :: a(0:N), b(0:N)
436	REAL(KIND=RP) :: nodes(0:N), weights(0:N)
437	REAL(KIND=RP) :: nrm, tol = 1.0d-8
438	INTEGER :: j
439
440	CALL GaussLegendreNodesAndWeights( N, nodes, weights )	1✔
441	!
442	! --------------------------------------
443	! Orthogonal polynomials shall give zero
444	! --------------------------------------
445	!
446	DO j = 0,N	6✔
447	a(j) = P3(nodes(j))	5✔
448	b(j) = P2(nodes(j))	6✔
449	END DO
450	nrm = gaussInnerProduct(a,b,weights,N)	1✔
451	gaussInnerProductIsOK = nrm < tol	1✔
452	!
453	! -------------------------------
454	! Same polynomial shall give norm
455	! -------------------------------
456	!
457	DO j = 0,N	6✔
458	a(j) = P3(nodes(j))	5✔
459	b(j) = a(j)	6✔
460	END DO
461	nrm = gaussInnerProduct(a,b,weights,N)	1✔
462	gaussInnerProductIsOK = ABS(nrm - 2.0_RP/(2.0_RP*3.0_RP+1.0_RP)) < tol &
463	.AND. gaussInnerProductIsOK	1✔
464
465	END FUNCTION gaussInnerProductIsOK	1✔
466	!
467	!////////////////////////////////////////////////////////////////////////
468	!
469	LOGICAL FUNCTION LegendrePolyAndDerivIsOK()	1✔
470	IMPLICIT NONE
471	INTEGER , PARAMETER :: N = 4
472	REAL(KIND=RP) :: x, dx, tol = 1.0d-8
473	REAL(KIND=RP) :: L(0:N) , LPrime(0:N)
474	REAL(KIND=RP) :: LExact(0:N), LPrimeExact(0:N)
475	INTEGER :: j
476
477	dx = 2.0_RP/5.0_RP	1✔
478	LegendrePolyAndDerivIsOK = .TRUE.	1✔
479	DO j = 0, 5	7✔
480	x = -1.0_RP + j*dx	6✔
481	CALL LegendrePolysAndDerivatives(N,x,L,LPrime)	6✔
482
483	LExact(0) = P0(x)	6✔
484	LExact(1) = P1(x)	6✔
485	LExact(2) = P2(x)	6✔
486	LExact(3) = P3(x)	6✔
487	LExact(4) = P4(x)	6✔
488
489	LegendrePolyAndDerivIsOK = MAXVAL(ABS(L - LExact)) < tol &
490	.AND. LegendrePolyAndDerivIsOK	36✔
491
492	LPrimeExact(0) = P0Prime(x)	6✔
493	LPrimeExact(1) = P1Prime(x)	6✔
494	LPrimeExact(2) = P2Prime(x)	6✔
495	LPrimeExact(3) = P3Prime(x)	6✔
496	LPrimeExact(4) = P4Prime(x)	6✔
497
498	LegendrePolyAndDerivIsOK = MAXVAL(ABS(LPrime - LPrimeExact)) < tol &
499	.AND. LegendrePolyAndDerivIsOK	37✔
500	END DO
501
502	END FUNCTION LegendrePolyAndDerivIsOK	1✔
503	!
504	!////////////////////////////////////////////////////////////////////////
505	!
506	LOGICAL FUNCTION legendreCoefsForTestAreOK()	1✔
507	!
508	! -------------------------------------------------------------------
509	! The Legendre series will be exact for any polynomial of degree <= N
510	! -------------------------------------------------------------------
511	!
512	IMPLICIT NONE
513	INTEGER, PARAMETER :: N = 4
514	REAL(KIND=RP) :: tol = 1.0d-8
515	REAL(KIND=RP) :: coefs(0:N)
516
517	CALL LegendreCoefsForTestPolynomial(coefs)	1✔
518	legendreCoefsForTestAreOK = MAXVAL(ABS(coefs - testCoefs)) < tol	6✔
519
520	END FUNCTION legendreCoefsForTestAreOK	1✔
521	!
522	!////////////////////////////////////////////////////////////////////////
523	!
524	LOGICAL FUNCTION legendreSeriesIsOK()	1✔
525	IMPLICIT NONE
526	REAL(KIND=RP) :: tol = 1.0d-8
527	REAL(KIND=RP) :: x, dx
528	INTEGER :: j
529
530	legendreSeriesIsOK = .TRUE.	1✔
531	dx = 2.0_RP/5.0_RP	1✔
532	DO j = 0, 5	7✔
533	x = -1.0_RP + j*dx	6✔
534	legendreSeriesIsOK = ABS(LegendreSeries(x,4,testCoefs,LA_EVALUATE_FUNCTION) &
535	- testPolynomial(x)) < tol &
536	.AND. legendreSeriesIsOK	7✔
537	END DO
538
539	END FUNCTION legendreSeriesIsOK	1✔
540	!
541	!////////////////////////////////////////////////////////////////////////
542	!
543	SUBROUTINE LegendreCoefsForTestPolynomial(coefs)	1✔
544	IMPLICIT NONE
545	INTEGER, PARAMETER :: N = 6
546	REAL(KIND=RP) :: coefs(0:4)
547	REAL(KIND=RP) :: a(0:N), b(0:N)
548	REAL(KIND=RP) :: nodes(0:N), weights(0:N)
549	INTEGER :: j,k
550
551	CALL GaussLegendreNodesAndWeights( N, nodes, weights )	1✔
552	!
553	! -----------------------
554	! Select the coefficients
555	! -----------------------
556	!
557	k = 0	1✔
558	DO j = 0,N	8✔
559	a(j) = testPolynomial(nodes(j))	7✔
560	b(j) = P0(nodes(j))	8✔
561	END DO
562	coefs(k) = gaussInnerProduct(a,b,weights,N)(2.0_RPk + 1.0_RP)/2.0_RP	1✔
563	k = 1	1✔
564	DO j = 0,N	8✔
565	b(j) = P1(nodes(j))	8✔
566	END DO
567	coefs(k) = gaussInnerProduct(a,b,weights,N)(2.0_RPk + 1.0_RP)/2.0_RP	1✔
568	k = 2	1✔
569	DO j = 0,N	8✔
570	b(j) = P2(nodes(j))	8✔
571	END DO
572	coefs(k) = gaussInnerProduct(a,b,weights,N)(2.0_RPk + 1.0_RP)/2.0_RP	1✔
573	k = 3	1✔
574	DO j = 0,N	8✔
575	b(j) = P3(nodes(j))	8✔
576	END DO
577	coefs(k) = gaussInnerProduct(a,b,weights,N)(2.0_RPk + 1.0_RP)/2.0_RP	1✔
578	k = 4	1✔
579	DO j = 0,N	8✔
580	b(j) = P4(nodes(j))	8✔
581	END DO
582	coefs(k) = gaussInnerProduct(a,b,weights,N)(2.0_RPk + 1.0_RP)/2.0_RP	1✔
583
584	END SUBROUTINE LegendreCoefsForTestPolynomial	1✔
585	!
586	!////////////////////////////////////////////////////////////////////////
587	!
588	FUNCTION testPolynomial(x) RESULT(s)	13✔
589	!
590	! ---------------------------------
591	! Evaluate a polynomial of degree 4
592	! ---------------------------------
593	!
594	IMPLICIT NONE
595	REAL(KIND=RP) :: x
596	REAL(KIND=RP) :: P(0:4), s
597	INTEGER :: k
598
599	P(0) = P0(x)	13✔
600	P(1) = P1(x)	13✔
601	P(2) = P2(x)	13✔
602	P(3) = P3(x)	13✔
603	P(4) = P4(x)	13✔
604
605	s = 0.0_RP	13✔
606	DO k = 0, 4	78✔
607	s = s + testCoefs(k)*P(k)	78✔
608	END DO
609
610	END FUNCTION testPolynomial	13✔
611	!
612	! -------------------------
613	! Some Legendre polynomials
614	! -------------------------
615	!
616	!
617	!////////////////////////////////////////////////////////////////////////
618	!
619	FUNCTION P0(x)	26✔
620	IMPLICIT NONE
621	REAL(KIND=RP) :: x, P0
622	P0 = 1.0_RP	26✔
623	END FUNCTION P0	26✔
624	!
625	!////////////////////////////////////////////////////////////////////////
626	!
627	FUNCTION P0Prime(x)	6✔
628	IMPLICIT NONE
629	REAL(KIND=RP) :: x, P0Prime
630	P0Prime = 0.0_RP	6✔
631	END FUNCTION P0Prime	6✔
632	!
633	!////////////////////////////////////////////////////////////////////////
634	!
635	FUNCTION P1(x)	26✔
636	IMPLICIT NONE
637	REAL(KIND=RP) :: x, P1
638	P1 = x	26✔
639	END FUNCTION P1	26✔
640	!
641	!////////////////////////////////////////////////////////////////////////
642	!
643	FUNCTION P1Prime(x)	6✔
644	IMPLICIT NONE
645	REAL(KIND=RP) :: x, P1Prime
646	P1Prime = 1.0_RP	6✔
647	END FUNCTION P1Prime	6✔
648	!
649	!////////////////////////////////////////////////////////////////////////
650	!
651	FUNCTION P2(x)	31✔
652	IMPLICIT NONE
653	REAL(KIND=RP) :: x, P2
654	P2 = 0.5_RP(3.0_RPx**2 - 1.0_RP)	31✔
655	END FUNCTION P2	31✔
656	!
657	!////////////////////////////////////////////////////////////////////////
658	!
659	FUNCTION P2Prime(x)	6✔
660	IMPLICIT NONE
661	REAL(KIND=RP) :: x, P2Prime
662	P2Prime = 0.5_RP(2.0_RP3.0_RP*x)	6✔
663	END FUNCTION P2Prime	6✔
664	!
665	!////////////////////////////////////////////////////////////////////////
666	!
667	FUNCTION P3(x)	36✔
668	IMPLICIT NONE
669	REAL(KIND=RP) :: x, P3
670	P3 = 0.5_RP(5.0_RPx*3 - 3.0_RPx)	36✔
671	END FUNCTION P3	36✔
672	!
673	!////////////////////////////////////////////////////////////////////////
674	!
675	FUNCTION P3Prime(x)	6✔
676	IMPLICIT NONE
677	REAL(KIND=RP) :: x, P3Prime
678	P3Prime = 0.5_RP(3.0_RP5.0_RPx*2 - 3.0_RP)	6✔
679	END FUNCTION P3Prime	6✔
680	!
681	!////////////////////////////////////////////////////////////////////////
682	!
683	FUNCTION P4(x)	26✔
684	IMPLICIT NONE
685	REAL(KIND=RP) :: x, P4
686	P4 = 0.125_RP(35.0_RPx*4 - 30.0_RPx**2 + 3.0_RP)	26✔
687	END FUNCTION P4	26✔
688	!
689	!////////////////////////////////////////////////////////////////////////
690	!
691	FUNCTION P4Prime(x)	6✔
692	IMPLICIT NONE
693	REAL(KIND=RP) :: x, P4Prime
694	P4Prime = 0.125_RP(4.0_RP35.0_RPx3 - 2.0_RP30.0_RP*x)	6✔
695	END FUNCTION P4Prime	6✔
696	END Module LegendreAlgorithms

trixi-framework / HOHQMesh / 30318764346

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous