LeafCoords_QuadTreeMesh2D Subroutine

  subroutine LeafCoords_QuadTreeMesh2D(this,leafIndex,geomInterp,coords)
    !! Physical geometry-node coordinates of leaf `leafIndex`, produced by repeated exact
    !! isoparametric subdivision of its root element along the quadtree path. geomInterp must be a
    !! degree-nGeo Lagrange interpolant on the mesh's geometry (quadrature) nodes - the same
    !! interpolant SELF_MeshRefinement_2D builds. Level 0 leaves return the root geometry directly.
    implicit none
    class(QuadTreeMesh2D),intent(in) :: this
    integer,intent(in) :: leafIndex
    type(Lagrange),intent(in) :: geomInterp
    real(prec),intent(out) :: coords(1:2,1:this%nGeo+1,1:this%nGeo+1)
    ! Local
    integer :: node,lvl,step
    integer,allocatable :: path(:)
    real(prec),allocatable :: cur(:,:,:),kids(:,:,:,:)

    node = this%leaf(leafIndex)
    lvl = this%level(node)

    allocate(cur(1:2,1:this%nGeo+1,1:this%nGeo+1))

    if(lvl == 0) then
      cur(1:2,:,:) = this%rootCoords(1:2,:,:,this%rootElem(node))
      coords = cur
      deallocate(cur)
      return
    endif

    ! Path of quadrant indices from root (step 1) down to the leaf (step lvl).
    allocate(path(1:lvl))
    do step = lvl,1,-1
      path(step) = this%quadrant(node)
      node = this%parent(node)
    enddo
    ! `node` is now the root.
    cur(1:2,:,:) = this%rootCoords(1:2,:,:,this%rootElem(this%leaf(leafIndex)))

    allocate(kids(1:2,1:this%nGeo+1,1:this%nGeo+1,1:4))
    do step = 1,lvl
      call SubdivideNodeCoords(geomInterp,this%nGeo,cur,kids)
      cur(1:2,:,:) = kids(1:2,:,:,path(step))
    enddo
    coords = cur

    deallocate(cur,kids,path)

  endsubroutine LeafCoords_QuadTreeMesh2D