!> \~japanese
!! \brief STD-Grid をバネ力学を用いて改良した SPR-GC-grid を作成する手続きを提供するモジュール.
!!
!! - SPRGCGrid_Builder モジュールについて
!!  - 本モジュールで提供される手続き construct_SPRGCGrid によって生成される SPR-GC-grid は, STD-grid の格子点の位置を
!!    バネ力学を用いて準一様に再配置させたのち, さらに各格子点位置をコントールボリュームの重心に移動させることによって
!!    できる正二十面格子である.
!!
!! - 正二十面体格子 SPR-GC-grid について
!!   - このバネ力学を用いた修正により, コントロールボリュームの面積・歪みが, 全球領域で単調に変化するようになる.
!!   - STD-grid の場合, コントロールボリュームの面積・歪みがフラクタル的な分布をしているために,
!!     数値的に微分値を評価する際, 組織的な数値誤差が発生してしまうが, この修正により数値誤差は改善される.
!!
!! \~english
!! \brief This module provides some subroutines to construct the SPR-GC-grid which is obtained by modifying the STD-Grid using the spring dynamics.
!!
!! \par About SPRGCGrid_Builder module
!!  - 本モジュールで提供される手続き construct_SPRGCGrid によって生成される SPR-GC-grid は, STD-grid の格子点の位置を
!!    バネ力学を用いて準一様に再配置させたのち, さらに各格子点位置をコントールボリュームの重心に移動させることによって
!!    できる正二十面格子である.
!!
!! \par About icosahedral grid: 'SPR-GC-grid'
!!   - このバネ力学を用いた修正により, コントロールボリュームの面積・歪みが, 全球領域で単調に変化するようになる.
!!   - STD-grid の場合, コントロールボリュームの面積・歪みがフラクタル的な分布をしているために,
!!     数値的に微分値を評価する際, 組織的な数値誤差が発生してしまうが, この修正により数値誤差は改善される.
!!
!! \~
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module SPRGCGrid_Builder

  ! モジュール引用 ; Use statements
  !

  ! 種類型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP         ! 倍精度実数型. Double precision.

  ! 数学定数
  ! Mathematic
  !
  use igmcore_math, only: PI   ! 円周率

  ! 線形代数
  ! Linear algebra
  !
  use igmcore_linear_algebra, only: &
    & vec_length, vec_normalize, &
    & cross

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

  ! 正二十面体格子データの管理
  ! Manager of Icosahedral grid data
  !
  use IcGrid2D_FVM_Manager, only: &
    & IcGrid2D_FVM, &
    & IcGrid2D_FVM_Init, IcGrid2D_FVM_Final, &
    & malloc_GP_array, malloc_CV_array, &
    & paste_margin_width, generate_CV5_GPindex, &
    & get_idMin, get_idMax, get_EffSize_Min, get_EffSize_Max, &
    & get_IcRadius, get_glevel, get_IcGridMeanInterval, &
    & check_pole, NOT_POLE_FLAG, &
    & RC_REGIONS_NUM


  ! STD-GC-grid の構築
  ! STD-GC-grid construction
  !
  use STDGCGrid_Builder, only: &
    & calc_STDGrid_GP_vertex, calc_CV_vertex, calc_CV_Area, &
    & move_GP_into_CVGravPt, split_rectangle_region, project_on_sphere

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedures
  !
  public :: construct_SPRGCGrid, set_iteration_limiter

  ! 非公開変数
  ! Private variables
  !

  !> \~japanese 力が課されていない時のバネの長さ.
  !! \~english The length of spring without imposed force.
  real(DP) :: eq_d

  !> \~japanese 正二十面体を内包する球の半径.
  !! \~english The radius of a sphere which an icosahedron is embedded into.
  real(DP) :: radius

  !> \~japanese チューニングパラメータ.
  !! \~english Tuning parameter.
  real(DP) :: beta = 1.15d0

  real(DP) :: lambda 
  
  ! ばねの振動の減衰のふるまいがおよそ臨界減衰となるように,
  ! 摩擦定数とバネ定数を選ぶ.

  !> \~japanese 摩擦定数.
  !! \~english Frictional constant.
  real(DP) :: alpha = 2000.0d0

  !> \~japanese バネ定数.
  !! \~english Spring constant.
  real(DP) :: k = 1.0d06

  !> \~japanese 時間刻み幅.
  !! \~english Time step.
  real(DP) :: dt = 0.001d0

  !> \~japanese 時間積分のステップ数.
  !! \~english The number of time integration steps
  integer :: end_tstep = 20000

  !> \~japanese 時間積分に伴う逐次計算ループの反復回数を制限するか否かのフラグ.
  !! デフォルトは, ループ回数を制限しない.
  !! \~english 時間積分に伴う逐次計算ループの反復回数を制限するか否かのフラグ.
  !! デフォルトは, ループ回数を制限しない.
  logical :: tstep_limiter = .false.

  !>
  real(DP) :: EPS = 1.0d-06

  !> \~japanese 収束判定を行うステップ間隔
  !! \~english The step interval for which cheking the convergence is performed.
  integer :: conv_check_step = 100

  !> \~japanese ある格子点を取り囲む格子点(6 個 or 5 個)のインデックスを保持する配列.
  !! \~english ある格子点を取り囲む格子点(6 個 or 5 個)のインデックスを保持する配列.
  integer, allocatable :: GP_IJ(:,:,:,:,:)

  integer, parameter :: INITIAL_GLEVEL = 4

contains

!
!> \~japanese
!! \brief SPR-GC-grid を構築する.
!!
!! @param[in,out] icgrid 構造型 IcGrid2D_FVM の変数.
!!
!! \~english
!! \brief Constructs the SPR-GC-grid.
!!
!! @param[in,out] icgrid The variable of derived type IcGrid2D_FVM.
!!
subroutine construct_SPRGCGrid( &
  & icgrid,     &  ! (inout)
  & tstep_limit &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid
  integer, intent(in), optional :: tstep_limit

  ! 作業変数
  ! Work variables
  !
  integer :: glevel
  integer :: gleveli
  integer :: rcID, i, j
  type(IcGrid2D_FVM), pointer :: nextcreated_sprIcgrid
  integer :: nxEMin, nxEMax
  type(IcGrid2D_FVM), pointer :: created_sprIcgrid
  integer :: EMin, EMax
  

  ! 実行文 ; Executable statements
  !

  radius = get_IcRadius(icgrid)
  glevel = get_glevel(icgrid)

  !
  if ( present(tstep_limit) ) then
    tstep_limiter = .false.
    end_tstep = tstep_limit
  end if

  allocate(created_sprIcgrid)
  call IcGrid2D_FVM_Init(created_sprIcGrid, INITIAL_GLEVEL, radius)
  call malloc_GP_array(created_sprIcGrid)

  ! 最初に STD-grid の格子点座標を求める.
  call calc_STDGrid_GP_vertex(created_sprIcGrid)


  do gleveli=INITIAL_GLEVEL, glevel
     
     write(*,*) 'Create icosahedral grid of glevel', gleveli, &
          & '(for doing one whose glevel is ', glevel, ')...'
 
     ! つり合いの状態における, 近傍の格子点までの距離の推定値を定める.
     ! ( Tomita, etal(2002) の式 (3), (4) )
     
     lambda = 2d0 * PI * radius / ( 10d0 * 2**(gleveli-1) )
     eq_d = beta * lambda

     ! SPR-GC-grid を生成する.
     !

     ! バネ力学を用いて, STD-grid の格子の空間分布を準一様にする.
     call modify_STDGrid_with_spring_dyn(created_sprIcGrid)

     if ( gleveli /= glevel ) then
        
        allocate(nextcreated_sprIcgrid)
        call IcGrid2D_FVM_Init(nextcreated_sprIcGrid, gleveli+1, radius)
        call malloc_GP_array(nextcreated_sprIcGrid)
!        call malloc_CV_array(nextcreated_sprIcGrid)

        nxEMin = get_EffSize_Min(nextcreated_sprIcGrid)
        nxEMax = get_EffSize_Max(nextcreated_sprIcGrid)
        
        EMin = get_EffSize_Min(created_sprIcGrid)
        EMax = get_EffSize_Max(created_sprIcGrid)

        do rcID=1, RC_REGIONS_NUM
          do j=EMin, EMax
            do i=EMin, EMax

              nextcreated_sprIcGrid%rcs_AGrid(1:3, 2*i-1, 2*j-1, rcID) &
               &   = created_sprIcGrid%rcs_AGrid(1:3, i,j,rcID)

            end do
          end do
        end do

        do rcID=1, RC_REGIONS_NUM
           do j=nxEMin,nxEMax-2,2
              do i=nxEMin,nxEMax-2,2
              
                 call split_rectangle_region( &
                      & rc_AGrid=nextcreated_sprIcGrid%rcs_AGrid(1:3,nxEMin:nxEMax,nxEMin:nxEMax,rcID), &
                      & ic_radius=radius, &
                      & row1=i, row2=i+2, col1=j, col2=j+2 &
                      & )
           
              end do
           end do
        end do

        call paste_margin_width(nextcreated_sprIcgrid)

        !
        call IcGrid2D_FVM_Final(created_sprIcGrid)
        deallocate(created_sprIcgrid)

        !
        created_sprIcGrid => nextcreated_sprIcGrid

   
     else

        icgrid%rcs_AGrid(:,:,:,:) = created_sprIcGrid%rcs_AGrid(:,:,:,:)

        call IcGrid2D_FVM_Final(created_sprIcGrid)
        deallocate(created_sprIcgrid)

     end if

  end do

  ! 改良した格子でコントールボリュームの格子位置, 面積を求める.
  call calc_CV_vertex(icgrid)
  call calc_CV_Area(icgrid)

  ! 格子位置を, コントールボリュームの重心位置に移動させる.
  call move_GP_into_CVGravPt(icgrid)

end subroutine construct_SPRGCGrid

!
!> \~japanese
!! \brief 時間積分に伴う逐次計算ループの反復回数を制限するか否かのフラグを設定する.
!!
!! デフォルトではフラグ tstep_limiter は真であり, 最大 2 万回で逐次計算ループを抜ける.
!! このフラグを偽に設定した場合には, check_convergence による収束判定が真を返すまで, 逐次計算ループを続行される.
!!
!! @param[in] limiter_flag 反復回数を制限するか否かのフラグ.
!!
!! \~english
!! \brief 時間積分に伴う逐次計算ループの反復回数を制限するか否かのフラグを設定する.
!!
!! デフォルトではフラグ tstep_limiter は真であり, 最大 2 万回で逐次計算ループを抜ける.
!! このフラグを偽に設定した場合には, check_convergence による収束判定が真を返すまで, 逐次計算ループを続行される.
!!
!! @param[in] limiter_flag 反復回数を制限するか否かのフラグ.
!!
subroutine set_iteration_limiter( &
  & limiter_flag &  ! (in)
  & )

  ! 宣言文 ; Declaration statement
  !
  logical, intent(in)  :: limiter_flag

  ! 実行文 ; Executable statement
  !
  tstep_limiter = limiter_flag

end subroutine set_iteration_limiter

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

!
!> \~japanese
!! \brief STD-grid を修正することによって得られる, SPR-GC-grid を生成する.
!!
!! 具体的には, バネ力学を用いて STD-grid の配置を修正する.
!! その後, STD-GC-grid と同様の方法で, コントロールボリュームを定義して, コントロールボリュームの重心に
!! 格子点を移動させる.
!!
!! このバネ力学を用いた修正により, コントロールボリュームの面積・歪みが, 全球領域で単調に変化するようになる.
!! STD-grid の場合, コントロールボリュームの面積・歪みがフラクタル的な分布をしているために,
!! 数値的に微分値を評価する際, 組織的な数値誤差が発生してしまうが, この修正により数値誤差は改善される.
!!
!! @param[in,out] icgrid 構造体 IcGrid2D_FVM の変数.
!!
!! \~english
!! \brief Generates the SPR-GC-grid obtained by modifying the STD-grid.
!!
!! Specifically, the STD-grid configuration is modified by employing the spring dynamics.
!! After that, control volumes are defined and grid points are moved to the gravitational centers of the control volumes
!! in the same way as STD-GC-grid.
!!
!! このバネ力学を用いた修正により, コントロールボリュームの面積・歪みが, 全球領域で単調に変化するようになる.
!! STD-grid の場合, コントロールボリュームの面積・歪みがフラクタル的な分布をしているために,
!! 数値的に微分値を評価する際, 組織的な数値誤差が発生してしまうが, この修正により数値誤差は改善される.
!!
!! @param[in,out] icgrid A variable of derived type IcGrid2D_FVM.
!!
subroutine modify_STDGrid_with_spring_dyn( &
  & icgrid )                                    ! (inout)

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid

  ! 作業変数
  ! Work variables
  !
  integer :: idMin,  idMax
  integer :: EMin,  EMax
  integer :: rcID, i, j
  integer :: pt_num

  ! 実行文 ; Executable statements
  !

  !write(*,*) 'modify STD-GC-grid'
  ! 配列の境界インデックスを取得
  idMin = get_idMin(icgrid)
  idMax = get_idMax(icgrid)
  EMin = get_EffSize_Min(icgrid)
  EMax = get_EffSize_Max(icgrid)

  ! バネ力学を利用して, STD-grid の格子位置を改良する.
  !

  allocate(GP_IJ(RC_REGIONS_NUM, EMin:EMax, EMin:EMax, 6, 2))

  do rcID=1, RC_REGIONS_NUM
    do j=EMin, EMax
      do i=EMin, EMax
           if( ( rcID > 5 .and. i==EMin .and. j==EMin ) .or. ( i==EMin .and. j==EMax ) .or. &
            & ( i==EMax .and. j==EMin ) .or. ( rcID <= 5 .and. i==EMax .and. j==EMax )     &
            &       ) then
            pt_num = 5
            GP_IJ(rcID,i,j, 1:pt_num,:) = generate_near_GP5index(icgrid, rcID,i,j)
          else
            pt_num = 6
            GP_IJ(rcID,i,j, 1:pt_num,:) = generate_near_GP6index(rcID,i,j, EMin, EMax)
          end if

      end do
    end do
  end do

  call move_GP_with_spring_dyn(icgrid, idMin, idMax, EMin, EMax)
  deallocate(GP_IJ)

end subroutine modify_STDGrid_with_spring_dyn

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

!
!> \~japanese
!! \brief 全球領域の各格子点に対してバネ力学を適用し, 仮想的なつり合い位置を求める.
!!
!! 各格子点を仮想的な質点とみなし, それぞれ近傍の質点 6 個とを仮想的なバネでつなぐ.
!! そして, 全質点がつり合いの状態にある位置を求める.
!!
!! @param[in]     idMin 正二十面格子の座標データを保持する配列の袖領域も含めたインデックスの最小値.
!! @param[in]     idMax 正二十面格子の座標データを保持する配列の袖領域も含めたインデックスの最大値.
!! @param[in]     EMin 正二十面格子の座標データを保持する配列の袖領域も含めないインデックスの最小値.
!! @param[in]     EMax 正二十面格子の座標データを保持する配列の袖領域も含めないインデックスの最大値.
!!
!! \~english
!! \brief 全球領域の各格子点に対してバネ力学を適用し, 仮想的なつり合い位置を求める.
!!
!! 各格子点を仮想的な質点とみなし, それぞれ近傍の質点 6 個とを仮想的なバネでつなぐ.
!! そして, 全質点がつり合いの状態にある位置を求める.
!!
!! @param[in]     idMin 正二十面格子の座標データを保持する配列の袖領域も含めたインデックスの最小値.
!! @param[in]     idMax 正二十面格子の座標データを保持する配列の袖領域も含めたインデックスの最大値.
!! @param[in]     EMin 正二十面格子の座標データを保持する配列の袖領域も含めないインデックスの最小値.
!! @param[in]     EMax 正二十面格子の座標データを保持する配列の袖領域も含めないインデックスの最大値.
!!
subroutine move_GP_with_spring_dyn( &
  & icgrid,                           & ! (inout)
  & idMin, idMax, EMin, EMax )        ! (in)

  ! 宣言文 ; Declaration statements
  !

  type(IcGrid2D_FVM), intent(inout) :: icgrid
  integer, intent(in) :: idMin,  idMax
  integer, intent(in) :: EMin,  EMax


  ! 作業変数
  ! Work variables
  !
  integer :: tstep
  integer rcID, i,j, k, convrecID
  real(DP) :: GP(6,3)
  integer :: pt_num
  real(DP) :: w0(RC_REGIONS_NUM, idMin:idMax, idMin:idMax, 1:3)

  integer, parameter :: DCONV_REC_NUM = 20
  integer, parameter :: DCONV_PLUS_ULIMIT = 10
  real(DP) :: conv_max, conv, conv_rec, conv_total
  real(DP) :: dconv_recs(DCONV_REC_NUM);
  integer :: dconv_plus_counter = 0
  logical :: conv_flag

  ! 実行文 ; Executable statements
  !

  !
  tstep = 1
  conv_flag = .false.

  ! 初期速度の設定
  w0(:,:,:,:) = 0.0d0

  do
    do rcID=1, RC_REGIONS_NUM
      !$omp parallel do private(i, k, GP, pt_num)
      do j=EMin,EMax
        do i=EMin,EMax
          if( ( rcID > 5 .and. i==EMin .and. j==EMin ) .or. ( i==EMin .and. j==EMax ) .or. &
            & ( i==EMax .and. j==EMin ) .or. ( rcID <= 5 .and. i==EMax .and. j==EMax )     &
            &       ) then
             ! 特異点
             ! 格子点(rcID,i,j)の近傍にある格子点 5 個の格子点を求める.
            pt_num = 5
          else
             ! 格子点(rcID,i,j)の近傍にある格子点 6 個の格子点を求める.
            pt_num = 6
          end if

          if( check_pole(icgrid, rcID,i,j) == NOT_POLE_FLAG ) then
            do k=1,pt_num
              GP(k,1:3) = icgrid%rcs_AGrid(1:3, GP_IJ(rcID,i,j,k,1), GP_IJ(rcID,i,j,k,2), rcID)
            end do

             ! 常微分方程式系を一タイムスッテプ分時間積分して, 一タイムスッテプ後の格子点位置を計算する.
            call integ_ode_1step_spring_dynamics( &
                 & r0=icgrid%rcs_AGrid(:,i,j,rcID), w0=w0(rcID,i,j,:), &
                 & GP=GP(1:pt_num,:), pt_num=pt_num &
                 & )


            call project_on_sphere( radius, icgrid%rcs_AGrid(:,i,j,rcID))

          end if

        end do ! for i
      end do ! for j
    end do ! for rcID

    ! 次のタイムステップに移る前に, のりしろを貼り合わせる.
    call paste_margin_width(icgrid)

    ! 収束判定
    if ( mod(tstep, 50) == 0 ) then

      !
      call check_global_convergence(icgrid, conv_total, conv_max, EMin, EMax)
      write(*,*) 'conv_max=', conv_max
      write(*,*) get_IcGridMeanInterval(icgrid), EPS

      if (  conv_max / lambda < EPS ) then
        conv_flag = .true.
      end if

!      convrecID = mod(int(tstep / conv_check_step), DCONV_REC_NUM)+1
!      dconv_recs(convrecID) = ( conv_total - conv_rec ) / abs(conv_max)
!      conv_rec = conv_total



!
!      if ( int(tstep / conv_check_step) > DCONV_REC_NUM ) then
!        if ( sum(dconv_recs) > 0.0d0 ) then
!           dconv_plus_counter = dconv_plus_counter + 1
!             if ( dconv_plus_counter > DCONV_PLUS_ULIMIT ) then
!               conv_flag = .true.
!             end if
!        end if
!      end if

    end if

    if ( mod(tstep, 50) == 0 ) then
      write(*,*) '-- tstep:', tstep, '/', end_tstep
      write(*,*) '* conv_max=', conv_max, ', dconv=', sum(dconv_recs) / DCONV_REC_NUM
      write(*,*) '* conv_total=', conv_total,  ', dconv_plus_counter=', dconv_plus_counter
     ! write(*,*) conv_total / ( RC_REGIONS_NUM * ( (EMax - EMin + 1)**2 -4 ) )
    end if

    if ( ( tstep_limiter .and. ( tstep > end_tstep ) ) .or. conv_flag ) then
      write(*,*) tstep, tstep_limiter, end_tstep
      exit
    end if

    tstep = tstep + 1
  end do ! 時間積分ループ

  ! 反復回数と収束の程度を表す値を出力する.
  !

  write(*,*) 'SPR-GC-grid generation: the number of interation=', tstep
  write(*,*) 'The value for checking convergence:', conv_max / get_IcGridMeanInterval(icgrid)

end subroutine move_GP_with_spring_dyn

!
!> \~japanese
!! \brief バネでつながれた格子点の振動が, 全球領域で収束したかを調べる.
!!
!! @param[in]     conv_toal_ave
!! @param[in]     conv_max
!! @param[in]     EMin 正二十面格子の座標データを保持する配列の袖領域も含めないインデックスの最小値.
!! @param[in]     EMax 正二十面格子の座標データを保持する配列の袖領域も含めないインデックスの最大値.
!!
!! \~english
!! \brief バネでつながれた格子点の振動が, 全球で収束したかを調べる.
!!
!! @param[in]     conv_toal_ave
!! @param[in]     conv_max
!! @param[in]     EMin 正二十面格子の座標データを保持する配列の袖領域も含めないインデックスの最小値.
!! @param[in]     EMax 正二十面格子の座標データを保持する配列の袖領域も含めないインデックスの最大値.
!!
subroutine check_global_convergence( &
  & icgrid,                   &  ! (in)
  & conv_total_ave, conv_max, &  ! (out)
  & EMin, EMax                &  ! (in)
  & )

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(in) :: icgrid
  real(DP), intent(out) :: conv_total_ave
  real(DP), intent(out) :: conv_max
  integer, intent(in) :: EMin
  integer, intent(in) :: EMax

  ! 作業変数
  ! Work variables
  !
  integer rcID, i, j, k
  real(DP) :: GP(6,3)
  real(DP) :: conv
  integer :: pt_num

  ! 実行文 ; Executable statements
  !

  conv_max = 0.0d0
  conv_total_ave = 0.0d0

  do rcID=1, RC_REGIONS_NUM
    !$omp parallel do private(i, k, GP, conv, pt_num) reduction(+:conv_total_ave)
    do j=EMin,EMax
      do i=EMin, EMax
        
        if( check_pole(icgrid, rcID,i,j) /= NOT_POLE_FLAG ) then
          cycle
        end if
        
        if(  ( ( i==EMin .and. j==EMin ) .or. ( i==EMin .and. j==EMax ) .or. &
          &    ( i==EMax .and. j==EMin ) .or. ( i==EMax .and. j==EMax ) )    &
          &   ) then
          pt_num = 5
        else
          pt_num = 6
        end if
        
          do k=1,pt_num
            GP(k,1:3) = icgrid%rcs_AGrid(1:3, GP_IJ(rcID,i,j,k,1), GP_IJ(rcID,i,j,k,2), rcID)
          end do

          conv = check_convergence(icgrid%rcs_AGrid(:,i,j,rcID), GP(1:pt_num,:), pt_num)
          conv_total_ave = conv_total_ave + conv

          if( conv  > conv_max ) then
            conv_max = conv
          end if
        
      end do
    end  do
  end do

  conv_total_ave = conv_total_ave / ( RC_REGIONS_NUM * ( (EMax - EMin + 1)**2 -4 ) )

end subroutine check_global_convergence

!
!> \~japanese
!! \brief 収束度合いを, Tomita etal, 2001 の式(22) を用いてチェックする.
!!
!! 戻りがゼロに近ければ近いほど, 質点はつり合いの位置に存在することになる.
!!
!! @param[in]     r0  点 P0 の位置ベクトル.
!! @param[in]     pts 点 P0 を取り囲む点の位置ベクトル.
!! @param[in]     pt_num 点 P0 を取り囲む近傍の格子点数.
!! @return Tomita etal, 2001 の式(22)の左辺の値.
!!
!! \~english
!! \brief 収束度合いを, Tomita etal, 2001 の式(22) を用いてチェックする.
!!
!! 戻りがゼロに近ければ近いほど, 質点はつり合いの位置に存在することになる.
!!
!! @param[in]     r0  The position vector of P0.
!! @param[in]     pts The position vectors of the points surrounding point P0.
!! @param[in]     pt_num The number of the grid points surrounding point P0.
!! @return The left-hand side value of the equation described as the equation (22) in Tomita, et al(2001).
!!
function check_convergence( &
  & r0, pts, pt_num &  ! (in)
  & ) result(val)

  ! 宣言文 ; Declaration statement
  !
  real(8), intent(in) :: r0(3)    !
  integer, intent(in) :: pt_num
  real(8), intent(in) :: pts(pt_num, 3) !
  real(8) val

  ! 作業変数
  ! Work variables
  !
  real(8) vec(3)
  integer i

  ! 実行文 ; Executable statement
  !
  vec = 0.0d0

  do i=1, pt_num
    vec(:) = vec(:) + &
       ( geodesic_arc_length(radius, r0, pts(i,:) ) - eq_d ) * calc_ei(r0(1:3), pts(i,1:3))
  end do

  val = vec_length(vec)

end function check_convergence

!
!> \~japanese
!! \brief 格子点(rcID, GP0_i, GP0_j) の近傍にある格子点 6 個に対応する, 配列のインデックスを生成する.
!!
!! @param[in] rcID
!! @param[in] GP0_i
!! @param[in] GP0_j
!! @param[in] EMin
!! @param[in] EMax
!!
!! \~english
!! \brief 格子点(rcID, GP0_i, GP0_j) の近傍にある格子点 6 個に対応する, 配列のインデックスを生成する.
!!
!! @param[in] rcID
!! @param[in] GP0_i
!! @param[in] GP0_j
!! @param[in] EMin
!! @param[in] EMax
!!
function generate_near_GP6index( &
  & rcID, GP0_i, GP0_j, EMin, EMax &  ! (in)
  & ) result(GP6Id)

  ! 宣言文 ; Declaration statements
  !

  integer, intent(in) :: rcID
  integer, intent(in) :: GP0_i
  integer, intent(in) :: GP0_j
  integer, intent(in) :: EMin
  integer, intent(in) :: EMax
  integer GP6Id(6, 2)

  ! 作業変数
  ! Work variables
  !
  integer i,j
  integer :: GP6IdShape(2) = (/ 6, 2 /)

  ! 実行文 ; Executable statements
  !

  i = GP0_i; j = GP0_j
  if ( rcID <= 5 .and. j==EMin ) then
    GP6Id(:,:) = reshape( (/ i,i-1,i-1,i-1,i,i+1, j-1,j-1,j,j+1,j+1,j /), GP6IdShape)

  else if ( rcID <= 5 .and. i==EMin ) then
    GP6Id(:,:) = reshape( (/ i,i-1,i-1,i,i+1,i+1, j-1,j-1,j,j+1,j,j-1 /), GP6IdShape)

  else if ( rcID > 5 .and. i==EMax ) then
    GP6Id(:,:) = reshape( (/ i+1,i+1,i,i-1,i-1,i, j,j+1,j+1,j+1,j,j-1 /), GP6IdShape)

  else if ( rcID > 5 .and. j==EMax ) then
    GP6Id(:,:) = reshape( (/ i+1,i+1,i+1,i,i-1,i, j-1,j,j+1,j+1,j,j-1 /), GP6IdShape)

  else
    GP6Id(:,:) = reshape( (/ i,i-1,i-1,i,i+1,i+1, j-1,j,j+1,j+1,j,j-1 /), GP6IdShape)
  end if

end function generate_near_GP6index

!
!> \~japanese
!! \brief 格子点(rcID, GP0_i, GP0_j) の近傍にある格子点 5 個に対応する, 配列のインデックスを生成する.
!!
!! @param[in,out] icgrid 構造体 IcGrid2D_FVM の変数.
!! @param[in]     rcID
!! @param[in]     GP0_i
!! @param[in]     GP0_j
!!
!! \~english
!! \brief 格子点(rcID, GP0_i, GP0_j) の近傍にある格子点 5 個に対応する, 配列のインデックスを生成する.
!!
!! @param[in,out] icgrid A variable of derived type IcGrid2D_FVM.
!! @param[in]     rcID
!! @param[in]     GP0_i
!! @param[in]     GP0_j
!!
function generate_near_GP5index( &
  & icgrid,            &  ! (inout)
  & rcID, GP0_i, GP0_j &  ! (in)
  & ) result(GP5Id)

  ! 宣言文 ; Declaration statements
  !
  type(IcGrid2D_FVM), intent(inout) :: icgrid
  integer, intent(in) :: rcID
  integer, intent(in) :: GP0_i
  integer, intent(in) :: GP0_j
  integer GP5Id(5, 2)

  ! 作業変数
  ! Work variables
  !
  integer :: tmp_GP5Id(5,3,2)

  ! 実行文 ; Executable statements
  !

  tmp_GP5Id(:,:,:) = generate_CV5_GPindex(icgrid, GP0_i, GP0_j, rcID)
  GP5Id(1:5, 1:2) = tmp_GP5Id(1:5, 2, 1:2)

end function generate_near_GP5index

!
!> \~japanese
!! \brief 固定された 6 点とバネでつながれた質点の位置の時間発展を計算する.
!!
!! 常微分方程式系の時間積分の方法は, 3 次の Runge=Kutta 法である.
!!
!! @param[in,out] r0
!! @param[in,out] w0
!! @param[in]     GP
!! @param[in]     pt_num
!! @param[in]     dt
!! @param[in]     alpha
!! @param[in]     k
!! @param[in]     eq_d
!! @param[in]     ic_radius 正二十面体を内包する球の半径.
!!
!! \~english
!! \brief 固定された 6 点とバネでつながれた質点の位置の時間発展を計算する.
!!
!! 常微分方程式系の時間積分の方法は, 3 次の Runge=Kutta 法である.
!!
!! @param[in,out] r0
!! @param[in,out] w0
!! @param[in]     GP
!! @param[in]     pt_num
!! @param[in]     dt
!! @param[in]     alpha
!! @param[in]     k
!! @param[in]     eq_d
!! @param[in]     ic_radius The radius of a sphere which an icosahedron is embedded into.
!!
subroutine integ_ode_1step_spring_dynamics( &
  & r0, w0,    &  ! (inout)
  & GP, pt_num &  ! (in)
  & )

  !
  ! 宣言文 ; Declaration statements
  !

  real(DP), intent(inout) :: r0(3)
  real(DP), intent(inout) :: w0(3)
  integer, intent(in) :: pt_num
  real(DP), intent(in) :: GP(pt_num,3)

  ! 作業変数
  ! Work variables
  !
  real(DP) r0_(3), w0_(3)
  real(DP) k1(2,3), k2(2,3), k3(2,3)

  ! 実行文 ; Executable statements
  !

  !
  ! 連立常微分方程式(Tomita etal, (2001) の式(20), (21) )を, 3 次の Runge=Kutta 法を使って解く.
  !

  ! k1 を求める.
  k1(1, 1:3) = f1(r0(:), w0(:))
  k1(2, 1:3) = f2(r0(:), w0(:), GP(:,:), pt_num, alpha, k, eq_d, radius)

  ! k2 を求める.
  r0_(:) = r0(:) + 0.5d0 * dt * k1(1,:)
  w0_(:) = w0(:) + 0.5d0 * dt * k1(2,:)
  k2(1, 1:3) = f1(r0_(:), w0_(:))
  k2(2, 1:3) = f2(r0_(:), w0_(:), GP(:,:), pt_num, alpha, k, eq_d, radius)

  ! k3 を求める.
  r0_(:) = r0(:) - dt * k1(1,:) + 2.0d0 * dt * k2(1,:)
  w0_(:) = w0(:) - dt * k1(2,:) + 2.0d0 * dt * k2(2,:)
  k3(1, 1:3) = f1(r0_(:), w0_(:))
  k3(2, 1:3) = f2(r0_(:), w0_(:), GP(:,:), pt_num, alpha, k, eq_d, radius)

  ! 質点の位置ベクトルと速度ベクトルを次のタイムレベルに更新する.
  r0(:) = r0(:) + ( k1(1,:) + 4.0d0 * k2(1,:) + k3(1,:) ) * dt / 6.0d0
  w0(:) = w0(:) + ( k1(2,:) + 4.0d0 * k2(2,:) + k3(2,:) ) * dt / 6.0d0

end subroutine integ_ode_1step_spring_dynamics

!
!> \~japanese
!! \brief 以下の常微分方程式の( Tomita etal, 2001 の式(21) ) の右辺の値を計算する.
!! \~english
!! \brief Calculates the value of right-hand side of the following ordinary differential equation described in Eq(21) in Tomita et al(2001).
!!
!! \~
!! \f[
!!  \DD{\Dvect{r_0}}{t} = f_1(\Dvect{r}_0, \Dvect{w}_0) = \Dvect{w}_0
!! \f]
!!
!! \~japanese
!! @param[in] r0 点 r0 の位置ベクトル.
!! @param[in] w0 点 r0 の速度ベクトル.
!! @return 関数 f1 の値.
!! \~english
!! @param[in] r0 The position vector of point r0.
!! @param[in] w0 The velocity vector of point r0.
!! @return The value of function f1.
!!
function f1( &
  & r0, w0 &  ! (in)
  & ) result(val)

  ! 宣言文 ; Declaration statements
  !
  real(DP), intent(in) :: r0(3)
  real(DP), intent(in) :: w0(3)
  real(DP) val(3)

  ! 実行文 ; Executable statements
  !
  val(1:3) = w0(1:3)

end function f1

!
!> \~japanese
!! \brief 以下の常微分方程式の( Tomita etal, 2001 の式(20) ) の右辺の値を計算する.
!! \~english
!! \brief Calculates the value of right-hand side of the following ordinary differential equation described in Eq(20) in Tomita et al(2001).
!!
!! \~
!! \f[
!!  M \DD{\Dvect{w_0}}{t} = f_2(\Dvect{r}_0, \Dvect{w}_0) \\
!!   = \sum_{i=1}^6 k(d_i - \bar{d}) \Dvect{e}_i - \alpha \Dvect{w}_0
!! \f]
!!
!! \~japanese
!! @param[in] r0 点 r0 の位置ベクトル.
!! @param[in] w0 点 r0 の速度ベクトル.
!! \param[in] pts
!! \param[in] pt_num
!! \param[in] alpha
!! \param[in] k
!! \param[in] eq_d
!! \param[in] radius 正二十面体を内包する球の半径.
!! \return　関数 f2 の値.
!! \~english
!! @param[in] r0 The position vector of point r0.
!! @param[in] w0 The velocity vector of point r0.
!! \param[in] pts
!! \param[in] pt_num
!! \param[in] alpha
!! \param[in] k
!! \param[in] eq_d
!! \param[in] radius The radius of a sphere which an icosahedron is embedded into.
!! \return The value of function f2.
!!
function f2( &
  & r0, w0, pts, pt_num, alpha, k, eq_d, radius &  ! (in)
  & ) result(val)

  ! 宣言文 ; Declaration statements
  !

  real(DP), intent(in) :: r0(3)
  real(DP), intent(in) :: w0(3)
  integer, intent(in) :: pt_num
  real(DP), intent(in) :: pts(pt_num,3)
  real(DP), intent(in) :: alpha
  real(DP), intent(in) :: k
  real(DP), intent(in) :: eq_d
  real(DP), intent(in) :: radius
  real(8) val(1:3)

  ! 作業変数
  ! Work variables
  !
  real(8) tmpv(3)
  integer i

  ! 実行文 ; Executable statements
  !

  tmpv(:) = 0.0d0
  do i=1, pt_num
    tmpv(:) = tmpv(:) + &
          ( geodesic_arc_length(radius, r0(1:3), pts(i,1:3)) - eq_d ) * calc_ei(r0(:), pts(i,:))
  end do

  ! M=1
  val(1:3) = k * tmpv(1:3) - alpha * w0(1:3)

end function f2

!
!> \~japanese
!! \brief 点 G_0(r0) から 点 Gi(ri) を結ぶベクトルの単位ベクトルを取得する.
!!
!! @param[in] r0 点 r0 の位置ベクトル.
!! @param[in] ri 点 ri の位置ベクトル.
!! @return \f$\Dvect{r}_i-\Dvect{r}_0\f$の単位ベクトル.
!!
!! \~english
!! \brief 点 G_0(r0) から 点 Gi(ri) を結ぶベクトルの単位ベクトルを取得する.
!!
!! @param[in] r0 The position vector of point r0.
!! @param[in] ri The position vector of point ri.
!! @return The unit vector of \f$\Dvect{r}_i-\Dvect{r}_0\f$.
!!
function calc_ei( &
  & r0, ri &  ! (in)
  & ) result(e)

  ! 宣言文 ; Declaration statements
  !
  real(8), intent(in) :: r0(3)
  real(8), intent(in) :: ri(3)
  real(8) e(3)

  ! 実行文 ; Executable statements
  !
  e(:) = cross( cross(r0(:), ri(:) ), r0)
  e(:) = vec_normalize(e(:))

end function calc_ei

end module SPRGCGrid_Builder
