!= テストの物理場を作成し, またその場の微分についての情報を解析的に与えるモジュール.  
!
! Authors::   Yuta Kawai
! Version::   $Id: test_diff_field.f90,v 1.3 2011-06-22 18:17:02 ykawai Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
module test_diff_field
  != テストの物理場を作成し, またその場の微分についての情報を解析的に与えるモジュール.  
  !
  !
  use dc_types, only: DP

  use igmcore_math
  !
  !
  use igmcore_coordinate_conversion
  
  !
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & get_IdMin, get_IdMax, get_EffSize_Min, get_EffSize_Max, &
    & get_IcRadius, &
    & get_LocalProcRectInfo, RC_REGIONS_NUM

  !
  !
  !
  use LocalProcRect_Manager, only: &
    & LocalProcRectInfo, &
    & get_InteriorRegion_IndexInfo
  !
  !
  use Field_IcGrid2D_Manager, only: &
    & Field_IcGrid2D, &
    & Field_IcGrid2D_Init, Field_IcGrid2D_Final, &
    & get_icgrid

  ! 宣言文 ; Declaration statements
  !	
  implicit none
  private
  
  type(Field_IcGrid2D), save, target, public :: field_scalar, field_vector
  type(Field_IcGrid2D), save, target, public :: dfield_scalar, dfield_vector
  type(Field_IcGrid2D), save, target, public :: ana_dfield_scalar, ana_dfield_vector
  
  integer, parameter, public :: DIV_OPERATOR = 1
  integer, parameter, public :: ROT_OPERATOR = 2
  integer, parameter, public :: GRADIENT_OPERATOR = 3
  integer, parameter, public :: LAPLACIAN_OPERATOR = 4

  real(DP) :: ic_radius

  public init_test_diff_field,  final_test_diff_field, &
    &    create_field, calc_ana_dfield
  
contains

subroutine init_test_diff_field(ic_grid)
  !
  !
  !
  type(IcGrid2D_FVM), intent(inout) :: ic_grid
  
  ic_radius = get_IcRadius(ic_grid)

  call Field_IcGrid2D_Init(self=field_scalar, icgrid=ic_grid, name='field', rank=1)
  call Field_IcGrid2D_Init(self=dfield_scalar, icgrid=ic_grid, name='dfield', rank=1)
  call Field_IcGrid2D_Init(self=ana_dfield_scalar, icgrid=ic_grid, name='ana_dfield', rank=1)

  call Field_IcGrid2D_Init(self=field_vector, icgrid=ic_grid, name='field', rank=3)
  call Field_IcGrid2D_Init(self=dfield_vector, icgrid=ic_grid, name='dfield', rank=3)
  call Field_IcGrid2D_Init(self=ana_dfield_vector, icgrid=ic_grid, name='ana_dfield', rank=3)

end subroutine init_test_diff_field

subroutine final_test_diff_field()
  
  call Field_IcGrid2D_Final(field_scalar)
  call Field_IcGrid2D_Final(dfield_scalar)
  call Field_IcGrid2D_Final(ana_dfield_scalar)
  call Field_IcGrid2D_Final(field_vector)
  call Field_IcGrid2D_Final(dfield_vector)
  call Field_IcGrid2D_Final(ana_dfield_vector)
  
end subroutine final_test_diff_field

subroutine create_field(field, diff_optr_type, m, n)
  !
  !
  !
  type(Field_IcGrid2D), intent(inout) :: field
  integer, intent(in) :: diff_Optr_type
  integer, intent(in) :: m, n

  !
  !
  !
  type(IcGrid2D_FVM), pointer :: ic_grid

  integer rcID, i,j
  integer :: glEMin, glEMax
  type(LocalProcRectInfo), pointer :: lcRCInfo
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax
  
  real(DP) :: s_pos(3), s_vec(3)

  !
  !
  ic_grid => get_icgrid(field)
  glEMin = get_EffSize_Min(ic_grid)
  glEMax = get_EffSize_Max(ic_grid)

  lcRCInfo => get_LocalProcRectInfo(ic_grid)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out) 
    
  do rcID=rcIDMin, rcIDMax
    do j=jEMin,jEMax
      do i=iEMin,iEMax
        if( ( rcID <= 5 .and. i==glEMin .and. j==glEMin ) .or. &
          & ( rcID >  5 .and. i==glEMax .and. j==glEMax ) ) then
           !
          field%val(:, i,j,rcID) = 0.0d0
        else
           !
          s_pos(:) = orth_to_geo_pos(ic_grid%rcs_AGrid(:, i,j,rcID))
          field%val(:, i,j,rcID) = calc_field_on_GP( diff_Optr_type, s_pos(1), s_pos(2), m, n)
          
        end if
      end do
    end do
  end do
   
end subroutine create_field

subroutine calc_ana_dfield(ana_dfield, diff_Optr_type, ret_dim, m, n)
  !
  !
  !
  type(Field_IcGrid2D), intent(inout) :: ana_dfield
  integer, intent(in) :: diff_Optr_type
  integer, intent(in) :: ret_dim
  integer, intent(in) :: m, n
  
  !
  !
  !  
  type(IcGrid2D_FVM), pointer :: ic_grid

  integer :: rcID, i, j
  real(DP) :: s_pos(3)
  integer :: glEMin, glEMax

  type(LocalProcRectInfo), pointer :: lcRCInfo
  integer :: iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax

  !
  !
  ic_grid => get_icgrid(ana_dfield)  
  glEMin = get_EffSize_Min(ic_grid)
  glEMax = get_EffSize_Max(ic_grid)

  lcRCInfo => get_LocalProcRectInfo(ic_grid)

  call get_InteriorRegion_IndexInfo( &
    & lcRCInfo,            &                         ! (in)
    & iEMin, iEMax, jEMin, jEMax, rcIDMin, rcIDMax ) ! (out) 
    
  do rcID=rcIDMin, rcIDMax
    do j=jEMin,jEMax
      do i=iEMin,iEMax
        if ( ( rcID <= 5 .and. i==glEMin .and. j==glEMin ) .or. &
          &  ( rcID >  5 .and. i==glEMax .and. j==glEMax ) ) then
          ana_dfield%val(:,i,j,rcID) = 0.0d0

        else
          s_pos(:) = orth_to_geo_pos(ic_grid%rcs_AGrid(:, i,j,rcID))
          ana_dfield%val(1:ret_dim, i,j,rcID) = &
            & eval_ana_diff_operator(diff_Optr_type, ret_dim, m, n, s_pos(1), s_pos(2))
        
        end if
      end do
    end do
  end do

end subroutine calc_ana_dfield

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! モジュール外非公開手続き
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

function calc_field_on_GP(diff_Optr_type, lambda, theta, m, n) result(val)
  !
  !
  !
  integer, intent(in) :: diff_Optr_type
  real(DP), intent(in) :: lambda
  real(DP), intent(in) :: theta
  integer, intent(in) :: m, n
  real(DP) val(3)
  

  select case(diff_Optr_type)
  case(GRADIENT_OPERATOR)
    val(1:1) = calc_betafield(lambda, theta, m, n)
    return
  case(DIV_OPERATOR)
    val(1:3) = geo_to_orth_vec( &
      & calc_vecfield_components(lambda, theta, m, n), (/ lambda, theta, ic_radius /) )
    return
  case(ROT_OPERATOR)
    val(1:3) = geo_to_orth_vec( &
      & calc_vecfield_components(lambda, theta, m, n), (/ lambda, theta, ic_radius /) )
    return
  end select

end function calc_field_on_GP

function calc_betafield(lambda, theta, m, n) result(val)
  !
  !
  !
  real(DP), intent(in) :: lambda
  real(DP), intent(in) :: theta
  integer, intent(in) :: m, n

  real(DP) val(1)

  val(1) = cos(m*lambda) * cos(n*theta)**4
  
end function calc_betafield

function calc_vecfield_components(lambda, theta, m, n) result(s_val)
  !
  !
  !
  real(DP), intent(in) :: lambda
  real(DP), intent(in) :: theta
  integer, intent(in) :: m, n

  real(DP) s_val(3)
  
  s_val(1) = - m * cos(n*theta)**4 * sin(lambda) * sin(m*lambda) / cos(theta)
  s_val(2) = - 4.0d0 * n * cos(n*theta)**3 * sin(n*theta) * sin(lambda) * cos(m*lambda)
  s_val(3) = 0.0d0

end function calc_vecfield_components

function eval_ana_diff_operator(diff_optr_type, ret_dim, m, n, lambda, theta) result(val)
  !
  !
  !
  integer, intent(in) :: diff_optr_type
  integer, intent(in) :: m, n
  real(8), intent(in) :: theta, lambda
  integer, intent(in) :: ret_dim

  real(8) val(ret_dim)

  select case (diff_optr_type)
    case (DIV_OPERATOR)
      val(1:ret_dim) = div_field(m, n, lambda, theta)
      return

    case (ROT_OPERATOR)
      val(1:ret_dim) = rot_field(m, n, lambda, theta)
      return

    case (GRADIENT_OPERATOR)
      val(1:ret_dim) = geo_to_orth_vec( &
        & grad_beta_field(m, n, lambda, theta), (/ lambda, theta, ic_radius /) )

      return

    case(LAPLACIAN_OPERATOR)
      val(1:ret_dim) = laplacian_betafield(lambda, theta, m, n)
    return  

  end select

end function eval_ana_diff_operator


function laplacian_betafield(lambda, theta_, m, n) result(val)
  !
  !
  !
  real(DP), intent(in) :: lambda
  real(DP), intent(in) :: theta_
  integer, intent(in) :: m, n

  real(DP) val(1)
  real(DP) :: theta

  theta = theta_
  val(1) = ( - m**2 * cos(n*theta)**4 / cos(theta)**2   &
    &        + 2d0 * n *tan(theta) * sin(2d0*n*theta) * cos(n*theta)**2                  &
    &        - 2d0 * n**2 * (2d0*cos(2d0*n*theta)*cos(n*theta)**2 - sin(2d0*n*theta)**2) &
    &       ) * cos(m*lambda)  
  
end function laplacian_betafield

function div_field(m, n, lambda, theta) result(val)
  !
  !
  !
  integer, intent(in) :: m, n
  real(8), intent(in) :: theta, lambda
  real(8) val(1)

  val(:) = &
    &    - m * cos(n*theta)**4 * ( cos(lambda) * sin(m*lambda) + m * sin(lambda) * cos(m*lambda) ) / cos(theta) **2 &
    &    - 4.0d0 * n * cos(n*theta)**2 * sin(lambda) * cos(m*lambda) * &
    &      ( &
    &        - 3.0d0 * n * sin(n*theta)**2 + n * cos(n*theta)**2 - sin(n*theta) * cos(n*theta) * tan(theta) &
    &      )

  val = val / ic_radius
end function div_field

function rot_field(m, n, lambda, theta) result(val)
  !
  !
  !
  integer, intent(in) :: m, n
  real(8), intent(in) :: lambda, theta
  real(8) val(1)

  val(:) = &
    &     - 4.0d0 * n * cos(n*theta)**3 * sin(n*theta) / cos(theta) * &
    &       ( &
    &         cos(lambda) * cos(m*lambda) - m * sin(lambda) * sin(m*lambda) &
    &       ) &
    &     - 4.0d0 * n * m * sin(lambda) * sin(m*lambda) * cos(n*theta)**3 * sin(n*theta) / cos(theta)

  val = val / ic_radius
end function rot_field

function grad_beta_field(m, n, lambda, theta) result(val)
  !
  !
  !
  integer, intent(in) :: m, n
  real(8), intent(in) :: theta, lambda
  real(8) val(3)

  val(1) = - m * cos(n * theta)**4 * sin(m * lambda) / cos(theta)
  val(2) = - 4.0d0 * n * cos(n * theta)**3 * sin(n * theta) * cos(m * lambda)
  val(3) = 0.0d0

  val = val / ic_radius
end function grad_beta_field

end module test_diff_field
