scale_polynominal.F90

Go to the documentation of this file.

 
#include "scaleFElib.h"
module scale_polynominal
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_const, only: &
    pi => const_pi
  use scale_io
  
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !  
 
  public :: polynominal_genlegendrepoly
  public :: polynominal_genlegendrepoly_sub
  public :: polynominal_gendlegendrepoly
 
  public :: polynominal_gengausslobattopt
  public :: polynominal_gengausslobattoptintweight
 
  public :: polynominal_gengausslegendrept
  public :: polynominal_gengausslegendreptintweight
 
  public :: polynominal_genlagrangepoly
  public :: polynominal_gendlagrangepoly_lglpt
  
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  private :: gen_jacobigaussquadraturepts
 
  !-----------------------------------------------------------------------------
 
contains
!OCL SERIAL
  function polynominal_genlagrangepoly(Nord, x_lgl, x) result(l)
    implicit none
    
    integer, intent(in) :: nord
    real(rp), intent(in) :: x_lgl(nord+1)
    real(rp), intent(in) :: x(:)
    real(rp) :: l(size(x), nord+1)
 
    integer :: n, i
    real(rp) :: p_lgl(nord+1,nord+1)    
    real(rp) :: p(size(x),nord+1)
    real(rp) :: pr(size(x),nord+1)
    !---------------------------------------------------------------------------
 
    p_lgl(:,:)  = polynominal_genlegendrepoly(nord, x_lgl)
    p(:,:)  = polynominal_genlegendrepoly(nord, x)
    pr(:,:) = polynominal_gendlegendrepoly(nord, x, p)
 
    do n=1, nord+1
      do i=1, size(x)
        if ( abs(x(i)-x_lgl(n)) < 1e-15_rp ) then
          l(i,n) = 1.0_rp
        else
          l(i,n) = &
            (x(i) - 1.0_rp)*(x(i) + 1.0_rp)*pr(i,nord+1)              &
            / (dble(nord*(nord+1))*p_lgl(n,nord+1)*(x(i) - x_lgl(n)))
        end if
      end do
    end do
 
    return
  end function polynominal_genlagrangepoly
 
!OCL SERIAL
  function polynominal_gendlagrangepoly_lglpt(Nord, x_lgl) result(lr)
    implicit none
 
    integer, intent(in) :: nord
    real(rp), intent(in) :: x_lgl(nord+1)
    real(rp) :: lr(nord+1, nord+1)
 
    integer :: n, k
    real(rp) :: p(nord+1,nord+1)
    real(rp) :: lr_nn
    !---------------------------------------------------------------------------
 
    p(:,:)  = polynominal_genlegendrepoly(nord, x_lgl)
 
    do n=1, nord+1
      lr_nn = 0.0_rp
      do k=1, nord+1
        if (k==1 .and. n==1) then
          lr(k,n) = - 0.25_rp*dble(nord*(nord+1))
        else if (k==nord+1 .and. n==nord+1) then
          lr(k,n) = + 0.25_rp*dble(nord*(nord+1))
        else if (k==n) then
          lr(k,n) = 0.0_rp
        else
          lr(k,n) = p(n,nord+1)/(p(k,nord+1)*(x_lgl(n) - x_lgl(k)))
        end if
 
        if ( k /= n ) then
          lr_nn = lr_nn + lr(k,n)
        end if        
      end do
      lr(n,n) = - lr_nn
    end do
    
    return
  end function polynominal_gendlagrangepoly_lglpt
 
!OCL SERIAL
  subroutine polynominal_genlegendrepoly_sub(Nord, x, P)
    implicit none
 
    integer, intent(in) :: nord
    real(rp), intent(in) :: x(:)
    real(rp), intent(out) :: p(size(x), nord+1)
 
    integer :: n
    !---------------------------------------------------------------------------
 
    if (nord < 0) then
       log_error("Polynominal_GenLegendrePoly",*)   "Nord must be larger than 0."
    end if
 
    p(:,1) = 1.0_rp
    if (nord==0) return
 
    p(:,2) = x(:)
    do n=2, nord
      p(:,n+1) = ( dble(2*n-1)*x(:)*p(:,n) - dble(n-1)*p(:,n-1) )/dble(n)
    end do
    
    return    
  end subroutine polynominal_genlegendrepoly_sub
 
!OCL SERIAL
  function polynominal_genlegendrepoly(Nord, x) result(P)
    implicit none
 
    integer, intent(in) :: nord
    real(rp), intent(in) :: x(:)
    real(rp) :: p(size(x), nord+1)
    !---------------------------------------------------------------------------
 
    call polynominal_genlegendrepoly_sub( nord, x(:), & ! (in)
       p(:,:)                                         ) ! (out)
    return    
  end function polynominal_genlegendrepoly
 
!OCL SERIAL
  function polynominal_gendlegendrepoly(Nord, x, P) result(GradP)
    implicit none
 
    integer, intent(in) :: nord
    real(rp), intent(in) :: x(:)
    real(rp), intent(in) :: p(:,:)
    real(rp) :: gradp(size(x), nord+1)
 
    integer :: n
    !---------------------------------------------------------------------------
 
    if (nord < 0) then
      log_error("Polynominal_GenDLegendrePoly",*)   "Nord must be larger than 0."
    end if
 
    gradp(:,1) = 0.0_rp
    if (nord == 0) return
 
    gradp(:,2) = 1.0_rp
    do n=2, nord
      gradp(:,n+1) = 2.0_rp*x(:)*gradp(:,n) - gradp(:,n-1) + p(:,n)
    end do
 
    return    
  end function polynominal_gendlegendrepoly
 
!OCL SERIAL
  function polynominal_gengausslobattopt(Nord) result(pts)
    implicit none
 
    integer, intent(in) :: nord
    real(rp) :: pts(nord+1)
 
    integer :: n1
    !---------------------------------------------------------------------------
 
    pts(1) = -1.0_rp; pts(nord+1) = 1.0_rp
    if (nord==1) return
 
    call gen_jacobigaussquadraturepts( 1, 1, nord-2, pts(2:nord) )
    return   
  end function polynominal_gengausslobattopt
 
!OCL SERIAL
  function polynominal_gengausslobattoptintweight(Nord) result(int_weight_lgl)
    implicit none
 
    integer, intent(in) :: nord
    real(rp) :: int_weight_lgl(nord+1)
 
    real(rp) :: lglpts1d(nord+1)
    real(rp) :: p1d_ori(nord+1, nord+1)
    !---------------------------------------------------------------------------
 
    lglpts1d(:)      = polynominal_gengausslobattopt( nord )
    p1d_ori(:,:)     = polynominal_genlegendrepoly( nord, lglpts1d )
 
    int_weight_lgl(:) = 2.0_rp/(dble(nord*(nord+1))*p1d_ori(:,nord+1)**2)
 
    return
  end function polynominal_gengausslobattoptintweight
 
!OCL SERIAL
  function polynominal_gengausslegendrept(Nord) result(pts)
    implicit none
 
    integer, intent(in) :: nord
    real(rp) :: pts(nord)
    !---------------------------------------------------------------------------
 
    call gen_jacobigaussquadraturepts( 0, 0, nord-1, pts(:) )
    return   
  end function polynominal_gengausslegendrept
 
!OCL SERIAL
  function polynominal_gengausslegendreptintweight(Nord) result(int_weight_gl)
    implicit none
 
    integer, intent(in) :: nord
    real(rp) :: int_weight_gl(nord)
 
    real(rp) :: glpts1d(nord)
    real(rp) :: p1d_ori(nord, nord+1)
    real(rp) :: dp1d_ori(nord, nord+1)
    !---------------------------------------------------------------------------
 
    glpts1d(:)    = polynominal_gengausslegendrept( nord )
    p1d_ori(:,:)    = polynominal_genlegendrepoly( nord, glpts1d )   
    dp1d_ori(:,:) = polynominal_gendlegendrepoly( nord, glpts1d, p1d_ori )
 
    int_weight_gl(:) = 2.0_rp / ( (1.0_rp - glpts1d(:)**2) * dp1d_ori(:,nord+1)**2 )
 
    return
  end function polynominal_gengausslegendreptintweight
 
  !- private -------------------------------
 
!OCL SERIAL
  subroutine gen_jacobigaussquadraturepts( alpha, beta, N, &
      x )
 
    implicit none
    integer, intent(in) :: alpha
    integer, intent(in) :: beta
    integer, intent(in) :: n
    real(rp), intent(out) :: x(n+1)
 
    integer :: i
    real(dp) :: d(n+1), e(n)
    real(dp) :: work(2*(n+1)-2), z(n+1,n+1)
    real(dp) :: h1(n+1)
    integer :: info
    !--------------------------------------------------------------
 
    if (n==0) then
      x(1) = - dble(alpha - beta) / dble(alpha + beta + 2)
      return
    end if
 
    do i=0, n
      h1(i+1) = dble( 2*i + alpha + beta )
    end do
 
    do i=1, n+1
      d(i) = - dble(alpha**2 - beta**2) / (h1(i) * (h1(i) + 2d0))
    end do
    do i=1, n
      e(i) = 2d0 / (h1(i) + 2d0) &
           * sqrt( dble(i * (i + alpha + beta) * (i + alpha) * (i + beta)) &
                   / ((h1(i) + 1d0) * (h1(i) + 3d0))                       )
    end do
    if ( dble(alpha + beta) < 1d-16) d(1) = 0.0_rp
    
    call dstev( 'Vectors', n+1, d, e, z, n+1, work, info )
    x(:) = d(:) 
 
    return
  end subroutine gen_jacobigaussquadraturepts
 
end module scale_polynominal