!> \~japanese
!! \brief 正二十面体格子の一つの格子点に対して, 物理場の微分(勾配, 発散, 回転)を評価するためのサブルーチンを提供するモジュール.
!! \~english
!! \brief 正二十面体格子の一つの格子点に対して, 物理場の微分(勾配, 発散, 回転)を評価するためのサブルーチンを提供するモジュール.
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module FVM_HDiff_operator

  ! モジュール引用 ; Use statement
  !
   
  ! 種類型パラメータ
  ! Kind type parameter
  !
  use dc_types, only: DP  ! 倍精度実数型. Double precision. 
  
  ! 線形代数
  ! Linear algebra
  !
  use igmcore_linear_algebra, only: &
    & vec_length, vec_normalize, &
    & cross, dot

  ! 球面三角法
  ! Spherical trigonometry
  !
  use igmcore_spherical_trigonometry, only: &
    & geodesic_arc_length

  ! 宣言文 ; Declaration statement
  !
  implicit none
  private

  ! 公開手続き 
  ! Public procedure
  !
  public :: eval_div_operator, eval_rot_operator, eval_grad_operator

  type, public :: FVM_HDerivate_CVInfo
    real(DP) :: geoArcLength_CV(6)
    real(DP) :: normalUniVec1_CV(3, 6)
    real(DP) :: normalUniVec2_CV(3, 6)
    real(DP) :: tangentUniVec_CV(3, 6)

    real(DP) :: geoArcLength_Tri(3, 6)
    real(DP) :: CVTriArea(6)
    real(DP) :: normalUniVec2_Tri(3, 3, 6)
    
    integer :: CV_num
  end type FVM_HDerivate_CVInfo

contains


!
!>
!!
!!
function eval_div_operator( &
  & val_GP0, val_CV, GP0, CVInfo,                   &  ! (in)
  & CVArea, val_dim, ret_dim &  ! (in)
  & ) result(val)

  ! 宣言文 ; Declaration statements
  !

  integer, intent(in) :: val_dim
  integer, intent(in) :: ret_dim
  type(FVM_HDerivate_CVInfo), intent(in) :: CVInfo
  real(DP), intent(in) :: val_GP0(val_dim)
  real(DP), intent(in) :: GP0(3)
  real(DP), intent(in) :: val_CV(3, CVInfo%CV_num)
  real(DP), intent(in) :: CVArea
  real(DP) val(ret_dim)

  ! 実行文 ; Executable statements
  !

  val(1) =  div_operator( &
             & vec_CV= val_CV, CVInfo=CVInfo, CVArea=CVArea )  ! (in)

end function eval_div_operator

!
!>
!!
!!
function eval_rot_operator( &
  & val_GP0, val_CV, GP0, CVInfo,               &  ! (in)
  & CVArea, val_dim, ret_dim &  ! (in)
  & ) result(val)

  ! 宣言文 ; Declaration statements
  !
  integer, intent(in) :: val_dim
  integer, intent(in) :: ret_dim
  type(FVM_HDerivate_CVInfo), intent(in) :: CVInfo
  real(DP), intent(in) :: val_GP0(val_dim)
  real(DP), intent(in) :: GP0(3)
  real(DP), intent(in) :: val_CV(3, CVInfo%CV_num)
  real(DP), intent(in) :: CVArea
  real(DP) val(ret_dim)

  ! 実行文 ; Executable statements
  !

  val(1) =  rot_operator( &
             & vec_CV= val_CV, CVInfo=CVInfo, CVArea=CVArea &  ! (in)
             & )
end function eval_rot_operator

!
!>
!!
!!
function eval_grad_operator( &
  & val_GP0, val_CV, GP0,             &  ! (in)
  & CVInfo, CVArea, val_dim, ret_dim  &  ! (in)
  & ) result(val)

  ! 宣言文 ; Declaration statements
  !
  integer, intent(in) :: val_dim
  integer, intent(in) :: ret_dim
  type(FVM_HDerivate_CVInfo), intent(in) :: CVInfo
  real(DP), intent(in) :: val_GP0(val_dim)
  real(DP), intent(in) :: GP0(3)
  real(DP), intent(in) :: val_CV(1, CVInfo%CV_num)
  real(DP), intent(in) :: CVArea
  real(DP) val(ret_dim)

  ! 実行文 ; Executable statements
  !

!  val(:) =  grad_operator( &
!             & q_GP0=val_GP0(1), q_CV= val_CV(:,1), pt_CV=CV, CVArea=CVArea, pt_size=CV_num, ic_radius=ic_radius )


  val(:) =  grad_operator2( &
             & q_GP0=val_GP0(1), q_CV= val_CV(1,:), GP=GP0, CVInfo=CVInfo, &  ! (in)
             & CVArea=CVArea )      ! (in)

end function eval_grad_operator


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 非公開手続き
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!> ベクトル場の水平発散を計算する. 
!!
!! @param[in] vec_CV    コントロールボリュームの頂点上のベクトル場の値を保持する配列. 
!! @param[in] CVInfo
!! @param[in] CVArea    コントロールボリュームの面積. 
!! @return    val       ベクトル場の水平発散. 
!!
function div_operator( vec_CV, CVInfo, CVArea ) result(val)

  ! 宣言文 ; Declaration statement
  !
  type(FVM_HDerivate_CVInfo), intent(in) :: CVInfo
  real(DP), intent(in) :: vec_CV(3, CVInfo%CV_num)
  real(DP), intent(in) :: CVArea
  real(DP) val

  ! 作業変数 
  ! Work variables
  !
  integer :: cvID, cvIDup, cv_num
  real(DP) :: n(3)

  ! 実行文 ; Executable statement
  !

  val = 0.0d0
  cv_num = CVInfo%CV_num

  do cvID=1, cv_num
    cvIDup = 1+mod(cvID, cv_num)

    val = val + CVInfo%geoArcLength_CV(cvID) * &
      &   dot( 0.5d0 * ( vec_CV(1:3, cvIDup) + vec_CV(1:3, cvID) ), CVInfo%normalUniVec1_CV(1:3, cvID) )
  end do

  val = val / CVArea

end function div_operator

!
!> 水平ベクトル場の回転の鉛直成分を計算する. 
!!
!! @param[in] vec_CV
!! @param[in] pt_CV
!! @param[in] CVArea
!! @param[in] pt_size
!! @param[in] ic_radius
!! @return    val
!!
function rot_operator( vec_CV, CVInfo, CVArea) result(val)
  
  ! 宣言文 ; Declaration statement
  !
  type(FVM_HDerivate_CVInfo), intent(in) :: CVInfo
  real(DP), intent(in) :: vec_CV(3, CVInfo%CV_num)
  real(DP), intent(in) :: CVArea
  real(DP) val

  ! 作業変数
  ! Work variable
  !
  integer :: cvID, CV_num, cvIDup

  ! 実行文 ; Executable statement
  !
  val = 0.0d0
  CV_num = CVInfo%CV_num

  do cvID=1, CV_num
    cvIDup = 1 + mod(cvID, CV_num)

    val = val + cvInfo%geoArcLength_CV(cvID) * &
      &   dot( 0.5d0 * ( vec_CV(1:3, cvIDup) + vec_CV(1:3, cvID) ), CVInfo%tangentUniVec_CV(1:3, cvID) )
  end do

  val = val / CVArea

end function rot_operator

!
!> 2 次元のスカラー場の水平勾配を計算する. 
!!
!! @param[in] q_GP0
!! @param[in] q_CV
!! @param[in] pt_CV
!! @param[in] CVArea
!! @param[in] pt_size
!! @param[in] ic_radius
!! @return    val
!!
function grad_operator( q_GP0, q_CV, pt_CV, CVArea, pt_size, ic_radius ) result(val)

  ! 宣言文 ; Declaration statement
  !
  integer, intent(in) :: pt_size
  real(DP), intent(in) :: q_GP0
  real(DP), intent(in) :: q_CV(pt_size)
  real(DP), intent(in) :: pt_CV(pt_size, 3)
  real(DP), intent(in) :: CVArea
  real(DP), intent(in) :: ic_radius
  real(DP) val(3)

  ! 作業変数　
  ! Work variable
  !
  integer :: i, iup
  real(DP) :: n(3)

  ! 実行文 ; Executable statement
  !

!  val(:) = 0.0d0
!
!  do i=1, pt_size
!    n(:) = cross( pt_CV(1+mod(i, pt_size), :), pt_CV(i,:) )
!    call vec_normalize(n)

!    val(:) = val(:) + geodesic_arc_length(ic_radius, pt_CV(1+mod(i, pt_size), :), pt_CV(i,:)) * &
!      &   ( 0.5d0 * ( q_CV(1+mod(i, pt_size)) + q_CV(i) ) - q_GP0 ) * n(:)
!  end do

!  val(:) = val(:) / CVArea

end function grad_operator

!
!>
!!
!! @param[in] q_GP0
!! @param[in] q_CV
!! @param[in] pt_CV
!! @param[in] CVArea
!! @param[in] pt_size
!! @param[in] ic_radius
!! @return    val
!!
function grad_operator2( q_GP0, q_CV, GP, CVInfo, CVArea) result(val)
  
  ! 宣言文  ; Declaration statement
  !
  type(FVM_HDerivate_CVInfo), intent(in) :: CVInfo
  real(DP), intent(in) :: q_GP0
  real(DP), intent(in) :: q_CV(CVInfo%CV_num)
  real(DP), intent(in) :: GP(3)
  real(DP), intent(in) :: CVArea
  real(DP) val(3)

  ! 作業変数 
  ! Work variables
  !
  integer :: cvID, cvIDup, CV_num

  ! 実行文 ; Executable statement
  !

  val(:) = 0.0d0
  CV_num = CVInfo%CV_num

  do cvID=1, CV_num
    cvIDup = 1+mod(cvID, CV_num)

    val(:) = val(:) + CVInfo%geoArcLength_CV(cvID) * &
      &   ( 0.5d0 * ( q_CV(cvIDup) + q_CV(cvID) ) - q_GP0 ) * CVInfo%normalUniVec2_CV(1:3, cvID)
  end do

  val(:) = val(:) / CVArea

end function grad_operator2


end module FVM_HDiff_operator
