SubdivideNodeCoords Subroutine

  subroutine SubdivideNodeCoords(geomInterp,nGeo,parentCoords,childCoords)
    !! Isoparametric subdivision of one element's geometry node coordinates into its eight
    !! children. The parent geometry is the degree-nGeo Lagrange interpolant through
    !! parentCoords (sampled at geomInterp's control points, i.e. the mesh geometry nodes);
    !! each child node coordinate is that interpolant evaluated at the corresponding point
    !! of the parent reference cube. Exact for any polynomial geometry of degree <= nGeo
    !! and for any control-node type, so straight-sided and curved (isoparametric) elements
    !! are both handled without approximation.
    !!
    !!   childCoords(d,i,j,k,c) =
    !!     sum_{ii,jj,kk} H(ii,i,ax) H(jj,j,ay) H(kk,k,az) parentCoords(d,ii,jj,kk)
    !!
    !! where H(ii,i,a) = l_ii( halfmap(s_i,a) ), halfmap(s,a) = (s + (2a-1))/2 maps the
    !! child node reference coordinate s_i into the parent half [a-1,a], and l_ii is the
    !! parent's Lagrange basis. d = 1,2,3 are the physical coordinate components.
    implicit none
    type(Lagrange),intent(in) :: geomInterp
    integer,intent(in) :: nGeo
    real(prec),intent(in) :: parentCoords(1:3,1:nGeo+1,1:nGeo+1,1:nGeo+1)
    real(prec),intent(out) :: childCoords(1:3,1:nGeo+1,1:nGeo+1,1:nGeo+1,1:8)
    ! Local
    integer :: i,ii,j,jj,k,kk,c,a,ax,ay,az
    real(prec) :: s
    real(prec) :: H(1:nGeo+1,1:nGeo+1,0:1)
    real(prec) :: rowsum(1:3),colsum(1:3),acc(1:3)

    ! Precompute the half-interval interpolation weights H(ii,i,a) for the lower (a=0) and
    ! upper (a=1) halves. H(:,i,a) are the parent basis values at the mapped position of
    ! child node i.
    do a = 0,1
      do i = 1,nGeo+1
        s = geomInterp%controlPoints(i)
        H(1:nGeo+1,i,a) = geomInterp%CalculateLagrangePolynomials(0.5_prec*(s+real(2*a-1,prec)))
      enddo
    enddo

    do c = 1,8
      ax = octAxc(c)
      ay = octAyc(c)
      az = octAzc(c)
      do k = 1,nGeo+1
        do j = 1,nGeo+1
          do i = 1,nGeo+1
            acc(1:3) = 0.0_prec
            do kk = 1,nGeo+1
              colsum(1:3) = 0.0_prec
              do jj = 1,nGeo+1
                rowsum(1:3) = 0.0_prec
                do ii = 1,nGeo+1
                  rowsum(1:3) = rowsum(1:3)+H(ii,i,ax)*parentCoords(1:3,ii,jj,kk)
                enddo
                colsum(1:3) = colsum(1:3)+H(jj,j,ay)*rowsum(1:3)
              enddo
              acc(1:3) = acc(1:3)+H(kk,k,az)*colsum(1:3)
            enddo
            childCoords(1:3,i,j,k,c) = acc(1:3)
          enddo
        enddo
      enddo
    enddo

  endsubroutine SubdivideNodeCoords