! 
! Copyright 2008 Shin'ya Murakami <murashin@epa.scitec.kobe-u.ac.jp>
! $Id:$
! 
! 2次元平面に等間隔で分布する格子点を, 
! 指定したlogicalな配列がtrueである位置ごとにグループ化する. 
! データに格子点インデックスを持つ連結リストを要素とする連結リストを返す. 
! 

module Divide2DField
  
  use mod_arearange
  use mod_pointindexlist
  use mod_boundarypointindexlist
  use mod_boundaryinfo
  use mod_integertuplelist
  use mod_patchedarea
  
  implicit none
  
  type AreaQuadTreeBranch
     type(AreaRange) r                 ! ツリーが指す範囲
     type(AreaQuadTree), pointer :: ul ! upper left
     type(AreaQuadTree), pointer :: ur ! upper right
     type(AreaQuadTree), pointer :: ll ! lower left
     type(AreaQuadTree), pointer :: lr ! lower right
  end type AreaQuadTreeBranch
  
  type AreaQuadTree
     logical is_leaf   ! leaf or branch
     type(AreaQuadTreeBranch) b ! branch
     type(PatchedArea), pointer :: l
  end type AreaQuadTree
  
  private
  public :: d2df_array2patchedarea

contains
  
  !
  ! 配列からPatchedAreaデータを作成
  !
  function d2df_array2patchedarea(lary, nx, ny, x_periodic, y_periodic) result(retval)
    
    logical, intent(in) :: lary(0:nx-1, 0:ny-1), x_periodic, y_periodic
    integer, intent(in) :: nx, ny
    type(PatchedArea), pointer :: retval
    type(AreaQuadTree), pointer :: aqt => null()
    integer n
    
    n = d2df_determine_dividing_times(nx) ! 分割回数を決定.
    call aqt_init(aqt, nx, ny)            ! ツリー全体を初期化
    call pa_init_id_counter()             ! paのノードにつけるエリアid初期化
    call aqt_make_tree(aqt, lary, nx, ny, n) ! 場を再帰的に分割し, ツリー作成. 
    call d2df_merge_recursive(aqt, n-1) ! ツリーと領域を再帰的に統合
    ! 最後に周期境界条件に伴った分の隣接領域を統合.
    if(x_periodic)then
       call pa_merge_periodic(aqt%l, .true.)
    end if
    if(y_periodic)then
       call pa_merge_periodic(aqt%l, .false.)
    end if
    ! おわり
    retval => aqt%l
    nullify(aqt%l)
    deallocate(aqt)
  end function d2df_array2patchedarea
  
  ! 
  ! 分割回数計算
  ! 
  integer function d2df_determine_dividing_times(n) result(retval)
    
    integer, intent(in) :: n
    
    integer i
    
    i = floor(log(dble(n))/log(2.0d0))
    if(i > 3) then
       i = i - 3
    end if
    !print *, "number of division: ", i
    retval = i
  end function d2df_determine_dividing_times
  
  ! 
  ! ツリーを再帰的に統合
  ! 
  recursive subroutine d2df_merge_recursive(aqt, depth)
    
    type(AreaQuadTree), pointer :: aqt
    integer, intent(in) :: depth
    integer n
    
    n = depth-1
    if(n > 0)then
       call d2df_merge_recursive(aqt%b%ll, n)
       call d2df_merge_recursive(aqt%b%ul, n)
       call d2df_merge_recursive(aqt%b%lr, n)
       call d2df_merge_recursive(aqt%b%ur, n)
    end if
    ! print *, "--- n-th stage merging: ", n, " ---"
    call pa_merge(aqt%b%ll%l, aqt%b%lr%l)
    call pa_merge(aqt%b%ul%l, aqt%b%ur%l)
    call pa_merge(aqt%b%ll%l, aqt%b%ul%l)
    aqt%l => aqt%b%ll%l
    nullify(aqt%b%ll%l)
    call aqt_branch_destroy(aqt)
    
  end subroutine d2df_merge_recursive
  
  ! 
  ! ツリー初期化. 
  ! 
  subroutine aqt_init(aqt, nx, ny)
    
    type(AreaQuadTree), pointer :: aqt
    integer, intent(in) :: nx, ny
    
    call aqt_branch_init(aqt, 0, nx-1, 0, ny-1)
  end subroutine aqt_init
  
  ! 
  ! ツリーの解放
  ! 
  subroutine aqt_branch_destroy(aqt)
    
    type(AreaQuadTree), pointer :: aqt
    
    deallocate(aqt%b%ul, aqt%b%ur, aqt%b%ll, aqt%b%lr)
    aqt%is_leaf = .true.
  end subroutine aqt_branch_destroy
  
  ! 
  ! ブランチ初期化
  ! 
  subroutine aqt_branch_init(aqt, ixmin, ixmax, iymin, iymax)
    
    type(AreaQuadTree), pointer :: aqt
    integer, intent(in) :: ixmin, ixmax, iymin, iymax
    
    allocate(aqt)
    aqt%is_leaf = .false.
    call ar_set(aqt%b%r, ixmin, ixmax, iymin, iymax)
  end subroutine aqt_branch_init
  
  ! 
  ! リーフ初期化
  ! 
  subroutine aqt_leaf_init(aqt, pa)
    
    type(AreaQuadTree), pointer :: aqt
    type(PatchedArea), pointer :: pa
    
    aqt%is_leaf = .true.
    aqt%l => pa
  end subroutine aqt_leaf_init
  
  ! 
  ! qtを4つに分解
  ! 
  recursive subroutine aqt_make_tree(aqt, lary, nx, ny, depth)
    
    type(AreaQuadTree), pointer :: aqt
    logical, intent(in) :: lary(0:nx-1, 0:ny-1)
    integer, intent(in) :: nx, ny, depth
    
    type(PatchedArea), pointer :: pa => null()
    type(AreaRange) r
    integer ixc,iyc
    
    if(depth-1 > 0)then
       r = aqt%b%r
       ixc = (r%ix1 - r%ix0 + 1)/2 + r%ix0
       iyc = (r%iy1 - r%iy0 + 1)/2 + r%iy0
       call aqt_branch_init(aqt%b%ll, r%ix0,   ixc-1, r%iy0,   iyc-1) !左下
       call aqt_branch_init(aqt%b%ul, r%ix0,   ixc-1,   iyc, r%iy1)   !左上
       call aqt_branch_init(aqt%b%lr,   ixc, r%ix1, r%iy0,   iyc-1)   !右下
       call aqt_branch_init(aqt%b%ur,   ixc, r%ix1,   iyc, r%iy1)     !右上
       call aqt_make_tree(aqt%b%ll, lary, nx, ny, depth-1)
       call aqt_make_tree(aqt%b%ul, lary, nx, ny, depth-1)
       call aqt_make_tree(aqt%b%lr, lary, nx, ny, depth-1)
       call aqt_make_tree(aqt%b%ur, lary, nx, ny, depth-1)
    else
       ! print *, "dividing: ", r%ix0,ixc-1,r%iy0,iyc-1
       ! print *, "array2list: ", aqt%b%r
       pa => pa_array2list2(lary, nx, ny, aqt%b%r) !配列からリスト生成
       call aqt_leaf_init(aqt, pa) ! 登録する
    end if
  end subroutine aqt_make_tree
  
end module Divide2DField
