test_labelize_point Subroutine

   subroutine test_labelize_point()
   !================================================================================================
   !! Function to test the function "labelize_point" on a QU9 domain with a 2nd order polynomial
   !!        along x and y
   implicit none

      real(kind=R8) :: x, y, x0, y0
      real(kind=R8) :: aax, aay

      character(len=1) :: point_kind, label

      real(kind=R8), parameter :: xx0 = -0.1_R8, fxx0 = +2._R8
      real(kind=R8), parameter :: yy0 = +0.2_R8, fyy0 = +1._R8

      real(kind=R8), dimension(1:9) :: Ni, dNidx, dNidy, d2Nidx2, d2Nidxdy, d2Nidy2, tab

      real(kind=R8), dimension(1:9, 1:6) :: derivatives
      real(kind=R8), dimension(1:3, 1:3) :: height

      ! select what kind of point to check
      point_kind = 'S'

      selectcase( point_kind )
         case( 'V' )  ; aax = +1._R8 ; aay = +1._R8
         case( 'P' )  ; aax = -1._R8 ; aay = -1._R8
         case( 'S' )  ; aax = +1._R8 ; aay = -1._R8
         case default ; stop
      endselect

      ! "tab" is given row wise:
      !     Nodes order:
      !
      !     4---7---3
      !     |   |   |
      !     8---9-- 6
      !     |   |   |
      !     1---5---2
      !
      tab(1:9) = [ f(x = -1._R8, x0 = xx0, a = aax, fx0 = fxx0) * f(x = -1._R8, x0 = yy0, a = aay, fx0 = fyy0),   &  ! 1
                   f(x = +1._R8, x0 = xx0, a = aax, fx0 = fxx0) * f(x = -1._R8, x0 = yy0, a = aay, fx0 = fyy0),   &  ! 5
                   f(x = +1._R8, x0 = xx0, a = aax, fx0 = fxx0) * f(x = +1._R8, x0 = yy0, a = aay, fx0 = fyy0),   &  ! 2
                   f(x = -1._R8, x0 = xx0, a = aax, fx0 = fxx0) * f(x = +1._R8, x0 = yy0, a = aay, fx0 = fyy0),   &  ! 8
                   f(x =  0._R8, x0 = xx0, a = aax, fx0 = fxx0) * f(x = -1._R8, x0 = yy0, a = aay, fx0 = fyy0),   &  ! 9
                   f(x = +1._R8, x0 = xx0, a = aax, fx0 = fxx0) * f(x =  0._R8, x0 = yy0, a = aay, fx0 = fyy0),   &  ! 6
                   f(x =  0._R8, x0 = xx0, a = aax, fx0 = fxx0) * f(x = +1._R8, x0 = yy0, a = aay, fx0 = fyy0),   &  ! 4
                   f(x = -1._R8, x0 = xx0, a = aax, fx0 = fxx0) * f(x =  0._R8, x0 = yy0, a = aay, fx0 = fyy0),   &  ! 7
                   f(x =  0._R8, x0 = xx0, a = aax, fx0 = fxx0) * f(x =  0._R8, x0 = yy0, a = aay, fx0 = fyy0) ]     ! 3

      ! stack is done column wise for reshape
      height(1:3, 1:3) = reshape( [ tab(1),   &  !
                                    tab(8),   &  !
                                    tab(4),   &  !
                                    tab(5),   &  !
                                    tab(9),   &  !
                                    tab(7),   &  !
                                    tab(2),   &  !
                                    tab(6),   &  !
                                    tab(3) ], [3, 3] )

      ! random point to compare the analytic function to the quadratic approx
      x = -0.35_R8
      y = +0.52_R8

      call deriv_N(x = x, y = y, mat_d = derivatives(1:9, 1:6))

      Ni(1:9)       = derivatives(1:9, 1)
      dNidx(1:9)    = derivatives(1:9, 2)
      dNidy(1:9)    = derivatives(1:9, 3)
      d2Nidx2(1:9)  = derivatives(1:9, 4)
      d2Nidxdy(1:9) = derivatives(1:9, 5)
      d2Nidy2(1:9)  = derivatives(1:9, 6)

      write(*,*) 'numerical h      : ', sum(tab * Ni)        ,  ' ; theoretical h      : ',   f(x = x, x0 = xx0, a = aax, fx0 = fxx0) *   f(x = y, x0 = yy0, a = aay, fx0 = fyy0)
      write(*,*) 'numerical dhdx   : ', sum(tab * dNidx)     ,  ' ; theoretical dhdx   : ',  df(x = x, x0 = xx0, a = aax            ) *   f(x = y, x0 = yy0, a = aay, fx0 = fyy0)
      write(*,*) 'numerical dhdy   : ', sum(tab * dNidy)     ,  ' ; theoretical dhdy   : ',   f(x = x, x0 = xx0, a = aax, fx0 = fxx0) *  df(x = y, x0 = yy0, a = aay            )
      write(*,*) 'numerical d2hdx2 : ', sum(tab * d2Nidx2)   ,  ' ; theoretical d2hdx2 : ', d2f(                 a = aax            ) *   f(x = y, x0 = yy0, a = aay, fx0 = fyy0)
      write(*,*) 'numerical d2hdy2 : ', sum(tab * d2Nidy2)   ,  ' ; theoretical d2hdy2 : ',   f(x = x, x0 = xx0, a = aax, fx0 = fxx0) * d2f(                 a = aay            )
      write(*,*) 'numerical d2hdxdy: ', sum(tab * d2Nidxdy)  ,  ' ; theoretical d2hdxdy: ',  df(x = x, x0 = xx0, a = aax            ) *  df(x = y, x0 = yy0, a = aay            )

      !===========================================

      call labelize_point(height = height(1:3, 1:3), label = label, x = x0, y = y0)

      write(*,*) 'theoretical xx0: ', xx0, 'numerical xx0: ', x0
      write(*,*) 'theoretical yy0: ', yy0, 'numerical yy0: ', y0

      write(*,*) 'theoretical point: ', point_kind, ' numerical point: ', label

   contains

      real(kind=R8) function f(x, x0, a, fx0)
      implicit none
      real(kind=R8), intent(in) :: x
      real(kind=R8), intent(in) :: x0
      real(kind=R8), intent(in) :: a
      real(kind=R8), intent(in) :: fx0

         real(kind=R8) :: b, c

         ! 2nd order polynomial defined by:
         ! ->   a : +1 or -1 (curvature sign)
         ! ->  x0 : the extremum abscissa
         ! -> fx0 : the extremum value
         b = -2 * a * x0
         c = fx0 + a * x0**2

         f = a * x**2 + b * x + c

      return
      endfunction f

      real(kind=R8) function df(x, x0, a)
      implicit none
      real(kind=R8), intent(in) :: x
      real(kind=R8), intent(in) :: x0
      real(kind=R8), intent(in) :: a

         real(kind=R8) :: b

         b = -2 * a * x0

         df = 2 * a * x + b

      return
      endfunction df

      real(kind=R8) function d2f(a)
      implicit none
      real(kind=R8), intent(in) :: a

         d2f = 2 * a

      return
      endfunction d2f

   endsubroutine test_labelize_point