! 
! Copyright 2008 Shin'ya Murakami <murashin@epa.scitec.kobe-u.ac.jp>
! $Id:$
! 
! 矩形領域内で隣接格子点がグループ化された格子点リスト.
! 境界位置の格子点はどの境界か(上下左右)の情報も保持.
! 

module mod_patchedarea
  
  use mod_arearange
  use mod_pointindexlist
  use mod_boundarypointindexlist
  use mod_boundaryinfo
  use mod_integertuplelist
  
  implicit none
  
  private :: pa_id_counter
  
  type PatchedAreaNode
     type(PatchedAreaNode), pointer :: n
     type(PatchedAreaNode), pointer :: p
     type(PointIndexList), pointer :: i
     integer id
  end type PatchedAreaNode
  
  type PatchedArea
     type(PatchedAreaNode), pointer :: n
     type(BoundaryPointIndexList), pointer :: b
     type(AreaRange) r
  end type PatchedArea
    
  integer, save :: pa_id_counter = 0
  
contains
  
!!!
!!! リスト操作
!!!
  
  ! 
  ! 初期化
  ! 
  subroutine pa_init(pa)
    
    type(PatchedArea), pointer :: pa
    type(PatchedAreaNode), pointer :: node => null()
    
    nullify(pa)
    allocate(pa,node)
    pa%n => node
    nullify(pa%b, node%p, node%n, node%i)
  end subroutine pa_init
  
  ! 
  ! paに登録されているあらゆるリンクリストを解放する.
  ! 
  subroutine pa_destroy(pa)
    
    type(PatchedArea), pointer :: pa
    type(PatchedAreaNode), pointer :: node =>null(), nnode =>null()
    
    if(.not. associated(pa))then
       print *, "error: you cannot destroy null pa."
       return
    end if
    if(associated(pa%b))then
       call bpil_destroy(pa%b)
    end if
    ! pa 自身を解放.
    node => pa%n
    deallocate(pa)
    nullify(pa)
    ! 最初に node を参照するときは, 空ノードが存在するはず
    ! なので必ず associated である.
    do while(associated(node%n))
       if(associated(node%i))then
          call pil_destroy(node%i)
       end if
       nnode => node%n
       deallocate(node)
       nullify(node)
       node => nnode
    end do
    if(associated(node%i))then
       call pil_destroy(node%i)
    end if
    deallocate(node)
    nullify(node)
  end subroutine pa_destroy
  
  ! 
  ! 最初のノードを得る
  ! 
  function pa_get_first(pa) result(retval)
    
    type(PatchedArea), pointer :: pa
    type(PatchedAreaNode), pointer :: retval
    
    if(.not. associated(pa))then
       print *, "error: cannot access pa, this is null pointer."
       nullify(retval)
       return
    end if
    
    retval => pa%n%n
    
  end function pa_get_first
  
  ! 
  ! 最後のノードを得る.
  ! 
  function pa_get_last(pa) result(retval)
    
    type(PatchedArea), pointer :: pa
    type(PatchedAreaNode), pointer :: retval, node => null()
    
    node => pa_each_init(pa)
    retval => node
    do while(pa_each(node))
       retval => node
    end do
  end function pa_get_last
  
  ! 
  ! 指定したノードの次のノードを得る
  ! 
  function pan_get_next(node) result(retval)
    
    type(PatchedAreaNode), pointer :: node, retval
    
    if(associated(node))then
       retval => node%n
    else
       print *, "error: this node is null pointer"
       nullify(retval)
    end if
  end function pan_get_next
  
  !
  ! 指定したノードの次のノードを削除
  !
  function pan_remove_next(node) result(retval)
    
    type(PatchedAreaNode), pointer :: node
    type(PointIndexList), pointer :: retval
    type(PatchedAreaNode), pointer :: nnode => null()
    
    if(.not. associated(node))then
       nullify(retval)
       print *, "error: accessed null pointer"
       return
    end if
    
    if(.not. associated(node%n))then
       ! next node is null
       print *, "error: cannot remove null node."
       nullify(retval)
    else 
       nnode => node%n
       if(associated(node%n%n))then    ! 次の次のノードがある
          node%n => node%n%n
          node%n%p => node
       else
          nullify(node%n)
       end if
       retval => nnode%i
       deallocate(nnode)
    end if
  end function pan_remove_next
  
  ! 
  ! 指定したノードの次にdataを持つノードを挿入
  ! 
  subroutine pan_add(node, data)
    
    type(PatchedAreaNode), pointer :: node
    type(PointIndexList), pointer :: data
    type(PatchedAreaNode), pointer :: node_new => null()
    
    if(.not. associated(node))then
       print *, "error in pan_add: cannot access node."
       ! return
    end if
    allocate(node_new)
    node_new%i => data
    node_new%n => node%n
    node_new%p => node
    if(associated(node%n))then
       node%n%p => node_new
    end if
    node%n => node_new
  end subroutine pan_add
  
  ! 
  ! リストの末尾にdataを持つノードを追加
  ! 
  subroutine pa_add_last(pa, data)
    
    type(PatchedArea), pointer  :: pa
    type(PointIndexList), pointer :: data
    
    call pan_add(pa_get_last(pa), data)
  end subroutine pa_add_last
  
  ! 
  ! リストの先頭にdataを持つノードを追加
  ! 
  subroutine pa_add_first(pa, data)
    
    type(PatchedArea), pointer  :: pa
    type(PointIndexList), pointer :: data
    
    call pan_add(pa%n, data)
  end subroutine pa_add_first
  
  !
  ! 指定したノードのデータを取得
  ! 
  function pan_get_data(node) result(retval)
    
    type(PatchedAreaNode), pointer :: node
    type(PointIndexList), pointer :: retval
    
    retval => node%i
  end function pan_get_data
  
  ! 
  ! 指定したノードにデータをセット
  ! 
  subroutine pan_set_data(node, data)
    
    type(PatchedAreaNode), pointer :: node
    type(PointIndexList), pointer :: data
    
    if(.not. associated(node))then
       print *, "error: cannot access node."
    end if
    node%i => data
  end subroutine pan_set_data
  
!!!
!!! module specificなリスト操作
!!!
  
  ! 
  ! idを持つnodeを返す.
  ! 
  function pa_get_node_by_id(pa, id) result(retval)
    
    type(PatchedArea), pointer :: pa
    integer, intent(in) :: id
    type(PatchedAreaNode), pointer :: node =>null(), retval
    
    if(.not. associated(pa))then
       print *, "error: cannot access pa."
    end if
    
    node => pa_each_init(pa)
    do while(pa_each(node))
       if(node%id == id)then
          retval => node
          return
       end if
    end do
    nullify(retval)
  end function pa_get_node_by_id
  
  ! 
  ! pa1にpa2をくっつける.
  ! 
  subroutine pa_attach_list(pa1, pa2)
    
    type(PatchedArea), pointer :: pa1, pa2
    type(PatchedAreaNode), pointer :: n1 =>null(), n2 =>null()
        
    if(.not. associated(pa1))then
       print *, "error: cannot access pa1."
    end if    
    if(.not. associated(pa2))then
       print *, "error: cannot access pa2."
    end if

    n1 => pa_get_last(pa1)
    n2 => pa_get_first(pa2)
    if(associated(n1))then
       if(associated(n2))then
          ! n1, n2 ともにある
          n1%n => n2
          n2%p => n1
          nullify(pa2%n%n)
       else
          ! n1 はある, n2は null. そのまま
       end if
    else if(associated(n2))then
       ! n1は null, n2はある
       pa1%n%n => n2
       n2%p => pa1%n
       nullify(pa2%n%n)
    else
       ! n1, n2 はともに null. そのまま
    end if
    call pa_destroy(pa2)
  end subroutine pa_attach_list

  ! 
  ! ノードのidを設定. 
  ! optional引数がない場合には, 
  ! 自動的に一意な番号を付ける. 
  ! 
  ! ノードの数がintegerの最大値を越えないことを仮定.
  ! 
  subroutine pa_set_node_id(node, id)
    
    type(PatchedAreaNode), pointer :: node
    integer, intent(in), optional :: id
    
    if(.not. associated(node))then
       print *, "error: cannot access node."       
    end if
    if(present(id))then
       node%id = id
    else
       pa_id_counter = pa_id_counter +1
       node%id = pa_id_counter
    end if
    !print *, "id is: ", node%id
  end subroutine pa_set_node_id
  
  ! 
  ! ノードに付けるidを初期化
  ! 
  subroutine pa_init_id_counter()
    
    pa_id_counter = 0
    
  end subroutine pa_init_id_counter
  
!!!
!!! ユーティリティ
!!!
  
  ! 
  ! 要素数を数える
  ! 
  integer function pa_size(pa) result(retval)
    
    type(PatchedArea), pointer :: pa    
    type(PatchedAreaNode), pointer :: node => null()
    integer n
    
    if(.not. associated(pa))then
       print *, "error: cannot access pa."
    end if
    
    n = 0
    node => pa_each_init(pa)
    do while(pa_each(node))
       n = n + 1
    end do
    retval = n
  end function pa_size
  
  ! 
  ! ポインタのスワップ
  ! 
  subroutine pa_swap_pointer(pa1, pa2)
    
    type(PatchedArea), pointer :: pa1, pa2
    type(PatchedArea), pointer :: tmp => null()
    
    tmp => pa1
    pa1 => pa2
    pa2 => tmp
    
  end subroutine pa_swap_pointer
  
  ! 
  ! while loop用初期化
  ! 
  function pa_each_init(pa) result(retval)
    
    type(PatchedArea), pointer :: pa
    type(PatchedAreaNode), pointer :: retval
    
    retval => pa%n
  end function pa_each_init
  
  ! 
  ! while ループ用
  ! 
  logical function pa_each(node) result(retval)
    
    type(PatchedAreaNode), pointer :: node
    
    node => node%n
    retval = associated(node)
  end function pa_each
  
!!! 
!!! リストの作成, リストから配列への書き戻し
!!!
  
  ! 
  ! 配列から, グループ化された配列の格子点列を持つPatchedAreaリストを構成
  ! 
  function pa_array2list(lary, nx, ny, ix0, ix1, iy0, iy1) result(retval)
    
    logical, intent(in) :: lary(0:nx-1, 0:ny-1)
    integer, intent(in) :: nx, ny, ix0, ix1, iy0, iy1
    type(PatchedArea), pointer :: pa => null(), retval
    type(PointIndexList), pointer :: pil => null()
    
    !print *, "nx, ny, ix0, ix1, iy0, iy1", nx, ny, ix0, ix1, iy0, iy1
    call pa_init(pa)
    pa%r = AreaRange(ix0, ix1, iy0, iy1)
    pil => pil_parray2list(lary, nx, ny, ix0, ix1, iy0, iy1)
    call pa_grouping(pa, pil) ! 隣接チェックしてグループ化
    call pa_generate_bpil(pa)
    retval => pa
  end function pa_array2list
  
  !
  ! AreaRange版
  ! 
  function pa_array2list2(lary, nx, ny, r) result(retval)
    
    logical, intent(in) :: lary(0:nx-1, 0:ny-1)
    integer, intent(in) :: nx, ny
    type(AreaRange), intent(in) :: r
    type(PatchedArea), pointer :: retval
    
    retval => pa_array2list(lary, nx, ny, r%ix0, r%ix1, r%iy0, r%iy1)
  end function pa_array2list2
  
  ! 
  ! リストから, 保持している格子点を論理配列に書き出す.
  ! 
  subroutine pa_list2array(pa, lary, nx, ny)
    
    type(PatchedArea), pointer :: pa    
    integer, intent(in) :: nx, ny
    logical, intent(out) :: lary(0:nx-1, 0:ny-1)
    type(PatchedAreaNode), pointer :: node => null()
    
    lary = .false.
    node => pa_each_init(pa)
    do while(pa_each(node))
       call pil_list2array(lary, nx, ny, node%i)
    end do
  end subroutine pa_list2array
  
  ! 
  ! 領域の範囲情報と格子点列から境界位置格子点を取り出し, 
  ! 境界位置にある, 境界情報付き格子点リストを作成する
  !
  ! pa%rが適切に設定されていることを仮定.
  ! pa%bが設定されていない(nullである)ことを仮定
  ! 
  subroutine pa_generate_bpil(pa)
    
    type(PatchedArea), pointer :: pa
    type(PatchedAreaNode), pointer :: pa_node => null()
    type(BoundaryPointIndexList), pointer :: bpil => null()
    type(PointIndexListNode), pointer :: pil_node => null()
    type(BoundaryInfo) b
    type(BoundaryPointIndex) bpi
    
    call bpil_init(bpil)
    pa_node => pa_each_init(pa)
    do while(pa_each(pa_node))
       pil_node => pil_each_init(pa_node%i)
       do while(pil_each(pil_node))
          if(ar_is_at_boundary(pa%r, pil_node%data%ix, &
               &                     pil_node%data%iy, b))then
             call bpi_set_data(bpi, b, pil_node%data, pa_node%id)
             call bpil_add_first(bpil, bpi)
          end if
       end do
    end do
    pa%b => bpil
  end subroutine pa_generate_bpil
  
  !
  ! 最小分割領域内で並んだ格子点リスト pool をグループ分けする.
  ! 分けたらそれぞれに領域のidをふる
  !
  subroutine pa_grouping(pa, pool)
    
    type(PatchedArea), pointer :: pa
    type(PointIndexList), pointer :: pool
    type(PointIndexList), pointer :: pil => null()
    type(PointIndexList), pointer :: pil_waitings => null()
    type(PointIndexListNode), pointer :: pool_node => null()
    type(PointIndexListNode), pointer :: piln_tmp => null()
    type(PointIndex) pi, tmp
    
    call pil_init(pil_waitings)
    outer: do while(pil_size(pool) > 0) ! プールに要素が1こ以上ある
!    outer: do while(.not. pil_is_empty(pool)) ! プールに要素が1こ以上ある
       ! プールの先頭から1こpiを取り出して, pilを作成. 
       pi = piln_remove_next(pool%n)
       call pil_init(pil)
       call pil_add_first(pil, pi)
       ! paの先頭に新しいノードを登録. 
       call pa_add_first(pa, pil)
       call pa_set_node_id(pa_get_first(pa))
       ! プールの残りをチェックして, pilの隣にあったら pil_waitings に登録.
       inner: do while(.true.)
          pool_node => pil_each_init(pool)
          do while(pil_each(pool_node))
             if(pi_is_adjacent(pi, pool_node%data))then
                ! pil_waitingsに隣接しているデータを登録
                piln_tmp => pool_node%p
                tmp = piln_remove_next(pool_node%p)
                call pil_add_first(pil_waitings, tmp)
                pool_node => piln_tmp
             end if
          end do
          if(pil_size(pil_waitings) == 0)then
!          if(pil_is_empty(pil_waitings))then
             exit inner
          else ! pil_waitingsの先頭から1こ取り出して, pilに登録.
             pi = piln_remove_next(pil_waitings%n)
             call pil_add_last(pil, pi)
          end if
       end do inner
    end do outer
    call pil_destroy(pil_waitings)
    call pil_destroy(pool)
  end subroutine pa_grouping
  
!!!
!!! 併合に関する手続き
!!!
  
  ! 
  ! 領域の併合手順.
  ! 
  ! 1. 「領域の境界情報のピックアップ」
  !    左(下)領域から右(上)側境界をもつ領域の境界情報を取り出す.
  !    右(上)領域から左(下)側境界をもつ領域の境界情報を取り出す.
  ! 
  ! 2. 「隣接する塊領域の列挙」
  !    領域の境界情報から, 
  !    領域 A に属する格子点と領域 B に属する格子点で隣接したものを探す.
  !    隣接したものが見付かったら, 属するエリアidの組を覚えておく.
  ! 
  ! 3. 「隣接する塊領域 に関する情報を整理」
  !    エリアidの組をソートする. 重複分を消す.
  !    エリア境界情報から隣接情報を消す. 境界でなくなったら消去.
  ! 
  ! 4. 「単純併合作業」
  !    領域を単純にくっつける. つまり, エリアの範囲をくっつけ, 
  !    境界情報をくっつけ, 塊領域のリストもくっつける.
  !    また, 3で使ったエリア境界情報を戻しておく.
  ! 
  ! 5. 「塊領域の併合」
  !    エリアidの組を使って, 塊領域を併合する.
  !    このとき, 併合する二つの領域のうち, 片方はエリアidが消滅するので, 
  !    エリアidの組に現れる消滅したエリアidは, 
  !    併合先の新しいエリアidにしなければいけない.
  !    これをエリアidの組がなくなるまで続ける.
  !
  subroutine pa_merge(pa1, pa2)
    
    type(PatchedArea), pointer :: pa1, pa2
    type(BoundaryPointIndexList), pointer :: b1 =>null(), b2 =>null()
    
    type(IntegerTupleList), pointer :: area_id_tuple =>null()
    logical dir, need_swap
    
    !隣接境界位置を調べる
    call ar_check_adjacency_direction(dir, need_swap, pa1%r, pa2%r)
    if(need_swap)then
       print *, "need_swap called"
       call pa_swap_pointer(pa1, pa2)
    end if
    ! 以降ではpa1,pa2の順に並んでいると仮定.
    
    ! 1.「領域の境界情報のピックアップ」    
    call bpil_pickup_at_boundary_in_adjacent_area(pa1%b, pa2%b, b1, b2, dir)
    
    ! 2. 「隣接する塊領域の列挙」
    area_id_tuple => bpil_listup_in_adjacent(b1, b2, dir)
    
    ! 3. 「隣接する塊領域 に関する情報を整理」
    call itl_sort(area_id_tuple)
    call itl_uniq(area_id_tuple)
    call itl_delete_if_equal(area_id_tuple)
    call bpil_delete_at_boundary_in_adjacent_area(b1, b2, dir) !境界情報消去
    
    ! 4. 「単純併合作業」
    call bpil_attach_list(pa1%b, pa2%b) ! エリア境界情報をくっつける
    call bpil_attach_list(pa1%b, b1)    !
    call bpil_attach_list(pa1%b, b2)    !
    pa1%r = ar_merge(pa1%r, pa2%r, dir) ! エリア範囲設定
    call pa_attach_list(pa1, pa2) ! pa1にpa2をくっつける.
    
    ! 5. 「塊領域の併合」
    call pa_merge_patched_area(pa1, area_id_tuple)
    
  end subroutine pa_merge
  
  ! 
  ! 周期境界の場合に, 周期境界をまたいだ隣接領域を併合
  ! 
  ! periodic_in_x_or_y = .true. のとき x 方向に周期境界として併合
  !                      .false.のとき y 方向に周期境界として併合
  ! このとき, 境界情報は消さない. 
  ! 
  subroutine pa_merge_periodic(pa, periodic_in_x_or_y)
    
    type(PatchedArea), pointer :: pa
    type(BoundaryPointIndexList), pointer :: b1 => null(), b2 => null()
    logical, intent(in) :: periodic_in_x_or_y
    
    type(IntegerTupleList), pointer :: area_id_tuple => null()
    logical dir
    
    ! 与えられた周期境界位置から, 適当な境界情報を取り出す.
    if(periodic_in_x_or_y)then
       call bpil_pickup_at_right(pa%b, b1)
       call bpil_pickup_at_left(pa%b, b2)
       dir = .false.
    else
       call bpil_pickup_at_top(pa%b, b1)
       call bpil_pickup_at_bottom(pa%b, b2)
       dir = .true.
    end if
    area_id_tuple => bpil_listup_in_adjacent(b1, b2, dir)
    ! 3. 「隣接する塊領域 に関する情報を整理」
    call itl_sort(area_id_tuple)
    call itl_uniq(area_id_tuple)
    call itl_delete_if_equal(area_id_tuple)
    ! 4. 「単純併合作業」
    !call bpil_delete_at_boundary_in_adjacent_area(b1, b2, dir) !境界情報消去
    call bpil_attach_list(pa%b, b1) ! エリア境界情報をくっつける
    call bpil_attach_list(pa%b, b2) ! エリア境界情報をくっつける
    
    ! 5. 「塊領域の併合」
    call pa_merge_patched_area(pa, area_id_tuple)
    
  end subroutine pa_merge_periodic
  
  ! 名前がよくないが, 保留.
  ! 
  ! 領域内の塊領域を併合する.
  ! area_id_tupleに合併を行うべき領域idの組を置く
  ! 
  subroutine pa_merge_patched_area(pa, area_id_tuple)
    
    type(PatchedArea), pointer :: pa
    type(IntegerTupleList), pointer :: area_id_tuple
    type(IntegerTuple) tpl
    
    do while(itl_exist(area_id_tuple))
       tpl = itl_remove_first(area_id_tuple)
       !tplに挙げられたエリアidを持つnodeのpilどうしをくっつける.
       call pa_merge_nodes_by_id(pa, tpl)
       call itl_replace(area_id_tuple, tpl%j, tpl%i) ! tpl%j -> tpl%i
       call itl_sort(area_id_tuple)
       call itl_uniq(area_id_tuple)
       call itl_delete_if_equal(area_id_tuple)
       call bpil_replace_id(pa%b, tpl%j, tpl%i)
    end do
    call itl_destroy(area_id_tuple)
  end subroutine pa_merge_patched_area
  
  ! 
  ! idを指定して領域を併合する.
  ! 
  subroutine pa_merge_nodes_by_id(pa, tpl)
    
    type(PatchedArea), pointer :: pa
    type(IntegerTuple), intent(in) :: tpl
    type(PatchedAreaNode), pointer :: n1 => null(), n2 => null()
    type(PointIndexList), pointer :: tmp => null()
    
    n1 => pa_get_node_by_id(pa, tpl%i)
    n2 => pa_get_node_by_id(pa, tpl%j)
    if((.not. associated(n1)) .or. (.not. associated(n2)))then !
       print *, "error: no such id in pa" !
       print *, "status:"
       print *, "   tuple:", tpl%i, tpl%j
       print *, "   associated:", associated(n1), associated(n2)
       call pa_dump(pa)
       call exit(1) !
    end if !
    tmp => pan_remove_next(n2%p)
    call pil_attach_list(n1%i, tmp)
  end subroutine pa_merge_nodes_by_id
  
  ! 
  ! PatchedAreaのダンプ
  ! 
  subroutine pa_dump(pa)
    
    type(PatchedArea), pointer :: pa    
    type(PatchedAreaNode), pointer :: node => null()
    integer i 
    
    i = 0
    print *, "dumping patched area...."
    node => pa_each_init(pa)
    do while(pa_each(node))
       i = i + 1
       if(associated(node%i))then
          print *, i, node%id, associated(node%i), pil_size(node%i)
       else
          print *, i, node%id, associated(node%i)
       end if
    end do
    print *, "end of dumping"
  end subroutine pa_dump
  
  ! 
  ! idがゼロのノードがあるかどうかを返す.
  ! 
  logical function pa_id_zero_exist(pa) result(retval)
    
    type(PatchedArea), pointer :: pa
    type(PatchedAreaNode), pointer :: node => null()
    
    node => pa_each_init(pa)
    do while(pa_each(node))
       if(node%id == 0)then
          retval = .true.
          return
       end if
    end do
    retval = .false.
  end function pa_id_zero_exist
  
  ! 
  ! PatchedAreaの整合性チェック
  ! いったん最後のノードまで行って, 最初の空ノードまで戻る.
  ! 
  subroutine pa_check_consistency(pa)
    
    type(PatchedArea), pointer :: pa    
    type(PatchedAreaNode), pointer :: node => null(), tmp => null()
    integer i
    
    i = 0
    node => pa_each_init(pa)
    do while(pa_each(node))
       tmp => node
       i = i + 1
    end do
    node => tmp
    do while(associated(node%p))
       node => node%p
       i = i - 1
    end do
    if(i == 0)then
       !
    else
       print *, "warning: inconsistent of pa have found."
       print *, "status: i=", i
    end if
  end subroutine pa_check_consistency
  
end module mod_patchedarea
