MAKE_GL_NW.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	make_gl_nw (ngp, xabsc, weig)
	Makes Gauss Legendre nodes and weights.

Function/Subroutine Documentation

make_gl_nw()

subroutine make_gl_nw	(	integer*4, intent(in)	ngp,
		real(dbp), dimension(ngp), intent(out)	xabsc,
		real(dbp), dimension(ngp), intent(out)	weig
	)

Makes Gauss Legendre nodes and weights.

Author

Ilario Mazzieri

Date

September, 2013

Version

1.0

Parameters

[in]	ngp	number of Gauss points
[out]	xabsc	Gauss points
[out]	weig	weights for the Gauss quadrature rule

Definition at line 28 of file MAKE_GL_NW.f90.

 
      implicit none
 
      INTEGER, parameter :: dbp = selected_real_kind (15,307)
      REAL(dbp) :: newv
      REAL(dbp) :: EPS, M_PI
      parameter(eps=3.0d-15)               !EPS is the relative precision
      parameter(m_pi=3.141592654d0) 
      
      INTEGER*4  i, j, m
      REAL(dbp)  p1, p2, p3, pp, z, z1
      INTEGER*4, INTENT(IN) :: ngp            ! # of Gauss Points
      REAL(dbp), INTENT(OUT) :: xabsc(ngp), weig(ngp)
 
 
           m = (ngp + 1) / 2
!* Roots are symmetric in the interval - so only need to find half of them  */
 
           do i = 1, m                                ! Loop over the desired roots */
 
                     z = cos( m_pi * (i-0.25d0) / (ngp+0.5d0) )
!*   Starting with the above approximation to the ith root,
!*          we enter the main loop of refinement by NEWTON'S method   */
 
100             p1 = 1.0d0
                p2 = 0.0d0
!*  Loop up the recurrence relation to get the Legendre
!*  polynomial evaluated at z                 */
 
                do j = 1, ngp
                   p3 = p2
                   p2 = p1
                   p1 = ((2.0d0*j-1.0d0) * z * p2 - (j-1.0d0)*p3) / j
                enddo
 
!* p1 is now the desired Legendre polynomial. We next compute pp,
!* its derivative, by a standard relation involving also p2, the
!* polynomial of one lower order.      */
                pp = ngp*(z*p1-p2)/(z*z-1.0d0)
                z1 = z
                z = z1 - p1/pp             ! Newton's Method  */
 
                if (dabs(z-z1) .gt. eps) GOTO  100
 
              xabsc(i) =  - z                            ! Roots will be bewteen -1.0 & 1.0 */
              xabsc(ngp+1-i) =  + z                        ! and symmetric about the origin   */
              weig(i) = 2.0d0/((1.0d0-z*z)*pp*pp)     ! Compute the weight and its       */
              weig(ngp+1-i) = weig(i)                 ! symmetric counterpart            */
 
      end do     ! i loop
 

Referenced by setup_dg_elem().

Here is the caller graph for this function: