program fempack_prb
!
!*******************************************************************************
!
!! FEMPACK_PRB calls the various FEMPACK tests.
!
  implicit none
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FEMPACK_PRB'
  write ( *, '(a)' ) '  Begin FEMPACK tests.'
 
  call test01
  call test02
  call test03
  call test04
  call test05

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'FEMPACK_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine test01
!
!*******************************************************************************
!
!! TEST01 tests the shape routines.
!
  implicit none
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST01'
  write ( *, '(a)' ) '  SHAPE_TEST tests the shape routines.'

  call shape_test ( 'Q4' )

  call shape_test ( 'Q8' )

  call shape_test ( 'Q9' )

  call shape_test ( 'Q12' )

  call shape_test ( 'Q16' )

  call shape_test ( 'QL' )

  call shape_test ( 'T3' )

  call shape_test ( 'T6' )

  call shape_test ( 'T10' )

  return
end
subroutine test02
!
!*******************************************************************************
!
!! TEST02 tests the grid routines.
!
  implicit none
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST02'
  write ( *, '(a)' ) '  Test the grid routines.'

  call grid_test ( 'Q4' )

  call grid_test ( 'Q8' )

  call grid_test ( 'Q9' )

  call grid_test ( 'Q12' )

  call grid_test ( 'Q16' )

  call grid_test ( 'QL' )

  call grid_test ( 'T3' )

  call grid_test ( 'T6' )

  call grid_test ( 'T10' )

  return
end
subroutine test03
!
!*******************************************************************************
!
!! TEST03 tests the map routines.
!
  implicit none
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST03'
  write ( *, '(a)' ) '  Test the map routines.'

  call map_test ( 'Q4' )

  call map_test ( 'Q8' )

  call map_test ( 'Q9' )

  call map_test ( 'Q12' )

  call map_test ( 'Q16' )

  call map_test ( 'QL' )

  call map_test ( 'T3' )

  call map_test ( 'T6' )

  call map_test ( 'T10' )

  return
end
subroutine test04
!
!*******************************************************************************
!
!! TEST04 tests DIV_Q4.
!
  implicit none
!
  integer, parameter :: m = 21
  integer, parameter :: n = 13
!
  real diff
  real div(m-1,n-1)
  real dudx
  real dudy
  real dvdx
  real dvdy
  integer i
  integer j
  real temp
  real u(m,n)
  real v(m,n)
  real vort(m-1,n-1)
  real x
  real xhi
  real xlo
  real xm
  real y
  real yhi
  real ylo
  real ym
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST04'
  write ( *, '(a)' ) '  DIV_Q4 estimates divergence and vorticity'
  write ( *, '(a)' ) '  using 4 node quadrilateral elements.'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6,a)' ) '  Original U, V data forms ', m, ' rows and '
  write ( *, '(a,i6,a)' ) '  ', n, ' columns.'
!
!  Set limits of the data points.
!
  xlo = 0.0E+00
  xhi = 1.0E+00
  ylo = 0.0E+00
  yhi = 2.0E+00
!
!  Put dummy data into U and V at the data nodes.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  I, J, X, Y, U(I,J), V(I,J)'
  write ( *, '(a)' ) ' '

  do i = 1, m
    y = ( real ( m - i ) * ylo + real ( i - 1 ) * yhi ) / real ( m - 1 )
    do j = 1, n
      x = ( real ( n - j ) * xlo + real ( j - 1 ) * xhi ) / real ( n - 1 )
      u(i,j) = x * y
      v(i,j) = sin ( x**2 + y**2 )
      write ( *, '(2i4,4g14.6)' ) i, j, x, y, u(i,j), v(i,j)
    end do
    write ( *, '(a)' ) ' '
  end do
!
!  Get DIV and VORT.
!
  call div_q4 ( div, m, n, u, v, vort, xhi, xlo, yhi, ylo )
!
!  Compare computed and known values at the centers of the elements.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'I, J, XM(I,J), YM(I,J)'
  write ( *, '(a)' ) ' '

  do i = 1, m-1

    ym = ( real(2*m-2*i-1) * ylo + real(2*i-1) * yhi ) / real ( 2 * m - 2 )

    do j = 1, n-1

      xm = ( real(2*n-2*j-1) * xlo + real(2*j-1) * xhi ) / real( 2 * n - 2 )

      write ( *, '(2i4,3g14.6)' ) i, j, xm, ym

    end do
    write ( *, '(a)' ) ' '
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  I, J, DIV(I,J), Exact Divergence, Difference'
  write ( *, '(a)' ) ' '

  do i = 1, m - 1

    ym = ( real(2*m-2*i-1) * ylo + real(2*i-1) * yhi ) / real ( 2 * m - 2 )

    do j = 1, n - 1

      xm = ( real(2*n-2*j-1) * xlo + real(2*j-1) * xhi ) / real ( 2 * n - 2 )

      dudx = ym
      dudy = xm
      dvdx = 2.0E+00 * xm * cos ( xm**2 + ym**2 )
      dvdy = 2.0E+00 * ym * cos ( xm**2 + ym**2 )

      write ( *, '(2i4,3g14.6)' ) &
        i, j, div(i,j), dudx + dvdy, div(i,j) - ( dudx + dvdy )

    end do
    write ( *, '(a)' ) ' '
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  I, J, VORT(I, j), Exact Vorticity, Difference'
  write ( *, '(a)' ) ' '

  do i = 1, m - 1

    y = ( real(2*m-2*i-1) * ylo + real(2*i-1) * yhi ) / real ( 2 * m - 2 )

    do j = 1, n-1

      x = ( real(2*n-2*j-1) * xlo + real(2*j-1) * xhi ) / real ( 2 * n - 2 )

      dudx = y
      dudy = x
      dvdx = 2.0E+00 * x * cos ( x**2 + y**2 )
      dvdy = 2.0E+00 * y * cos ( x**2 + y**2 )

      write ( *, '(2i4,3g14.6)' ) &
        i, j, vort(i,j), dvdx - dudy, vort(i,j) - ( dvdx - dudy )

    end do
    write ( *, '(a)' ) ' '
  end do

  return
end
subroutine test05
!
!*******************************************************************************
!
!! TEST05 demonstrates S_L2NORM.
!
  implicit none
!
  integer, parameter :: element_num = 40
  integer, parameter :: psi_num = element_num + 1
  integer, parameter :: quad_num = 5
!
  integer element
  real element_area(element_num)
  real exact
  integer i
  integer j
  integer k
  real l2norm
  real node_x(psi_num)
  integer psi
  real psi_quad(psi_num,element_num,quad_num)
  integer quad
  real quad_weight(quad_num)
  real quad_x(quad_num)
  real s_coef(psi_num)
  real x
  real, parameter :: x_max = 1.0E+00
  real, parameter :: x_min = 0.0E+00
  real xl
  real xr
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'TEST05'
  write ( *, '(a)' ) '  S_L2NORM computes the L2 norm of a scalar function'
  write ( *, '(a)' ) '  S(X) over a region (of any dimension), assuming:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  * that there is a set of finite element basis'
  write ( *, '(a)' ) '    functions PSI;'
  write ( *, '(a)' ) '  * that the region is decomposed into a number of'
  write ( *, '(a)' ) '    elements whose areas are known;'
  write ( *, '(a)' ) '  * that the integral is to be computed by a quadrature'
  write ( *, '(a)' ) '    rule applied in the same way to each element;'
  write ( *, '(a)' ) '  * that the value of the basis functions is given'
  write ( *, '(a)' ) '    at every quadrature node in every element;'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Our example will have one spatial dimension.'
  write ( *, '(a,g14.6)' ) '  XMIN = ', x_min
  write ( *, '(a,g14.6)' ) '  XMAX = ', x_max
  write ( *, '(a,i6)' ) '  The number of intervals will be    ', element_num
  write ( *, '(a,i6)' ) '  The number of basis functions is   ', psi_num
  write ( *, '(a,i6)' ) '  The number of quadrature points is ', quad_num
!
!  Set the nodes.
!
  do i = 1, psi_num
    node_x(i) = ( real ( psi_num - i )     * x_min &
                + real (           i - 1 ) * x_max ) / real ( psi_num - 1 )
  end do
!
!  Set the element "areas".
!
  do i = 1, element_num
    element_area(i) = node_x(i+1) - node_x(i)
  end do
!
!  Set the quadrature weights and abscissas.
!
  quad_weight(1) = 0.236926885056189087514264040720E+00
  quad_weight(2) = 0.478628670499366468041291514836E+00
  quad_weight(3) = 0.568888888888888888888888888889E+00
  quad_weight(4) = 0.478628670499366468041291514836E+00
  quad_weight(5) = 0.236926885056189087514264040720E+00

  quad_weight(1:5) = quad_weight(1:5) / 2.0E+00

  quad_x(1) = - 0.906179845938663992797626878299E+00
  quad_x(2) = - 0.538469310105683091036314420700E+00
  quad_x(3) =   0.0E+00
  quad_x(4) =   0.538469310105683091036314420700E+00
  quad_x(5) =   0.906179845938663992797626878299E+00
!
!  Set the finite element coefficients of S.  In our formulation,
!  these finite element coefficients are simply the function values
!  at the nodes.
!
  do i = 1, psi_num
    s_coef(i) = sin ( node_x(i) )
  end do
!
!  For each basis function I,
!    for each element J,
!      for each quadrature point K,
!
!  ...evaluate the basis function in the element at the quadrature point.
!
  psi_quad(1:psi_num,1:element_num,1:quad_num) = 0.0E+00

  do psi = 1, psi_num

    do element = 1, element_num

      if ( element < psi - 1 ) then

        cycle

      else if ( element == psi - 1 ) then

        xl = node_x(element)
        xr = node_x(element+1)

        do quad = 1, quad_num

          x = ( ( 1.0E+00 - quad_x(quad) ) * xl &
              + ( 1.0E+00 + quad_x(quad) ) * xr ) / 2.0E+00

          psi_quad(psi,element,quad) = ( x - xl ) / ( xr - xl )

        end do

      else if ( element == psi ) then

        xl = node_x(element)
        xr = node_x(element+1)

        do quad = 1, quad_num

          x = ( ( 1.0E+00 - quad_x(quad) ) * xl &
              + ( 1.0E+00 + quad_x(quad) ) * xr ) / 2.0E+00

          psi_quad(psi,element,quad) = ( xr - x ) / ( xr - xl )

        end do

      else if ( element > psi ) then
        cycle
      end if

    end do

  end do
!
!  Ask S_L2NORM to compute the integral of L2 norm of S.
!
  call s_l2norm ( psi_num, element_num, quad_num, element_area, &
    quad_weight, psi_quad, s_coef, l2norm )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  The computed L2 norm is ', l2norm
!
!  The integral of (sin(x))**2 is x/2 - sin(2x)/4
!
  exact = sqrt ( &
      ( 0.5E+00 * x_max - 0.25E+00 * sin ( 2.0E+00 * x_max ) ) &
    - ( 0.5E+00 * x_min - 0.25E+00 * sin ( 2.0E+00 * x_min ) ) )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  The exact value is      ', exact

  return
end