LeafCoords_OctreeMesh3D Subroutine

  subroutine LeafCoords_OctreeMesh3D(this,leafIndex,geomInterp,coords)
    !! Physical geometry-node coordinates of leaf `leafIndex`, produced by repeated exact
    !! isoparametric subdivision of its root element along the octree path. geomInterp
    !! must be a degree-nGeo Lagrange interpolant on the mesh's geometry (quadrature)
    !! nodes - the same interpolant SELF_MeshRefinement_3D builds. Level 0 leaves return
    !! the root geometry directly. Pure function of (root coords, level, octant path), so
    !! regenerated leaf geometry is bit-identical across epochs.
    implicit none
    class(OctreeMesh3D),intent(in) :: this
    integer,intent(in) :: leafIndex
    type(Lagrange),intent(in) :: geomInterp
    real(prec),intent(out) :: coords(1:3,1:this%nGeo+1,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:3,1:this%nGeo+1,1:this%nGeo+1,1:this%nGeo+1))

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

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

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

    deallocate(cur,kids,path)

  endsubroutine LeafCoords_OctreeMesh3D