MthPolynomialDerivativeMatrix Subroutine

public subroutine MthPolynomialDerivativeMatrix(N_in, xGP, deriv, D)

Computes mth polynomial differentiation matrix for interpolation polynomial given by set of nodes. (Algorithm 38, Kopriva book)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: N_in

polynomial degree
real(kind=wp), intent(in) :: xGP(0:N_in)

Gauss point positions for the reference interval [-1,1]
integer, intent(in) :: deriv

derivative (starting at 1)
real(kind=wp), intent(out) :: D(0:N_in,0:N_in)

differentiation Matrix

Calls

proc~~mthpolynomialderivativematrix~~CallsGraph proc~mthpolynomialderivativematrix MthPolynomialDerivativeMatrix interface~barycentricweights BarycentricWeights proc~mthpolynomialderivativematrix->interface~barycentricweights interface~polynomialderivativematrix PolynomialDerivativeMatrix proc~mthpolynomialderivativematrix->interface~polynomialderivativematrix interface~barycentricweights->interface~barycentricweights interface~polynomialderivativematrix->interface~polynomialderivativematrix

Source Code

SUBROUTINE MthPolynomialDerivativeMatrix(N_in,xGP,deriv,D)
IMPLICIT NONE
!----------------------------------------------------------------------------------------------------------------------------------
! INPUT/OUTPUT VARIABLES
INTEGER,INTENT(IN)     :: N_in              !! polynomial degree
REAL(wp),INTENT(IN)    :: xGP(0:N_in)       !! Gauss point positions for the reference interval [-1,1]
INTEGER,INTENT(IN)     :: deriv             !! derivative (starting at 1)
REAL(wp),INTENT(OUT)   :: D(0:N_in,0:N_in)  !! differentiation Matrix
!----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
INTEGER            :: ideriv,iGP,iLagrange
REAL(wp)           :: wBary(0:N_in),Dii_last
!==================================================================================================================================
D=0.0_wp
IF(deriv.LT.1) CALL abort(__STAMP__, &
                          'deriv in MthPolyDerivativeMatrix must be >=1!')
IF(deriv.GT.N_in) RETURN
CALL BarycentricWeights(N_in,xGP,wBary)
CALL PolynomialDerivativeMatrix(N_in,xGP,D)
IF(deriv.EQ.1) RETURN
DO ideriv=2,deriv
  DO iGP=0,N_in
    Dii_last = D(iGP,iGP)
    D(iGP,iGP)=0.0_wp
    DO iLagrange=0,N_in
      IF(iLagrange.NE.iGP)THEN
        D(iGP,iLagrange)=REAL(ideriv,wp)/(xGP(iGP)-xGP(iLagrange))*(wBary(iLagrange)/wBary(iGP)*Dii_last-D(iGP,iLagrange))
        D(iGP,iGP)=D(iGP,iGP)-D(iGP,iLagrange)
      END IF ! (iLagrange.NE.iGP)
    END DO ! iLagrange
  END DO ! iGP
END DO !iDeriv
END SUBROUTINE MthPolynomialDerivativeMatrix