BuildTransferPlan Subroutine

  subroutine BuildTransferPlan(forest,nOld,oldLeaf,plan)
    !! Build the old->new transfer plan for the adaptation epoch that took the forest from the
    !! leaf configuration (nOld, oldLeaf) - a snapshot of (forest%nLeaves, forest%leaf) taken
    !! before mutating - to its current leaf configuration. See the module documentation for the
    !! allowed mutations within one epoch. The plan is valid until the forest is mutated again.
    implicit none
    type(QuadTreeMesh2D),intent(in) :: forest
    integer,intent(in) :: nOld
    integer,intent(in) :: oldLeaf(1:nOld)
    type(TransferPlan2D),intent(out) :: plan
    ! Local
    integer :: i,li,walk,d,mx,steps
    logical :: found
    integer,allocatable :: oldElemOfNode(:)
    integer,allocatable :: famElem(:,:)
    integer,allocatable :: rev(:)

    do i = 1,nOld
      if(oldLeaf(i) < 1 .or. oldLeaf(i) > forest%nNodes) then
        print*,__FILE__,':',__LINE__, &
          ' : Error : oldLeaf snapshot contains a node id outside the forest node pool.'
        stop 1
      endif
    enddo

    ! Invert the snapshot: old element index by node id, and each old family's element indices
    ! by (child quadrant, parent node id). Coarsening zeroes the parent's child pointers but
    ! never the children's parent/quadrant entries, so famElem recovers detached families.
    allocate(oldElemOfNode(1:forest%nNodes))
    oldElemOfNode = 0
    allocate(famElem(1:4,1:forest%nNodes))
    famElem = 0
    do i = 1,nOld
      oldElemOfNode(oldLeaf(i)) = i
      if(forest%parent(oldLeaf(i)) > 0) then
        famElem(forest%quadrant(oldLeaf(i)),forest%parent(oldLeaf(i))) = i
      endif
    enddo

    plan%nOld = nOld
    plan%nNew = forest%nLeaves
    mx = max(forest%MaxLevel(),1) ! prolongation depth is bounded by the deepest leaf level
    allocate(plan%sourceKind(1:plan%nNew))
    allocate(plan%sourceElem(1:plan%nNew))
    allocate(plan%family(1:4,1:plan%nNew))
    allocate(plan%depth(1:plan%nNew))
    allocate(plan%path(1:mx,1:plan%nNew))
    plan%sourceKind = SELF_TRANSFER_COPY
    plan%sourceElem = 0
    plan%family = 0
    plan%depth = 0
    plan%path = 0
    plan%maxDepth = 0

    allocate(rev(1:mx))

    do li = 1,plan%nNew

      ! Ascend from the new leaf until reaching its data source; collect the quadrant taken at
      ! each level (bottom-up in rev, reversed into plan%path top-down). The chain has at most
      ! level+1 nodes, ending at the root.
      walk = forest%leaf(li)
      d = 0
      found = .false.
      do steps = 0,forest%level(forest%leaf(li))
        if(oldElemOfNode(walk) > 0) then
          if(d == 0) then
            plan%sourceKind(li) = SELF_TRANSFER_COPY
          else
            plan%sourceKind(li) = SELF_TRANSFER_PROLONG
          endif
          plan%sourceElem(li) = oldElemOfNode(walk)
          found = .true.
          exit
        elseif(famElem(1,walk) > 0 .and. famElem(2,walk) > 0 .and. &
               famElem(3,walk) > 0 .and. famElem(4,walk) > 0) then
          plan%sourceKind(li) = SELF_TRANSFER_RESTRICT
          plan%family(1:4,li) = famElem(1:4,walk)
          found = .true.
          exit
        endif
        if(forest%parent(walk) == 0) exit
        d = d+1
        rev(d) = forest%quadrant(walk)
        walk = forest%parent(walk)
      enddo

      if(.not. found) then
        print*,__FILE__,':',__LINE__, &
          ' : Error : new leaf has no old-leaf ancestor or coarsened old family; the snapshot'// &
          ' does not describe one adaptation epoch of this forest.'
        stop 1
      endif

      plan%depth(li) = d
      do i = 1,d
        plan%path(i,li) = rev(d-i+1)
      enddo
      plan%maxDepth = max(plan%maxDepth,d)

    enddo

    deallocate(oldElemOfNode,famElem,rev)

  endsubroutine BuildTransferPlan