Divergence_TwoPointVector3D_t Subroutine

  subroutine Divergence_TwoPointVector3D_t(this,df)
    !! Computes the split-form (two-point) divergence of a 3-D vector field
    !! in the reference element (computational coordinates).
    !!
    !! The split-form divergence at node (i,j,k) is
    !!
    !!   (nabla.F)_{i,j,k} = 2 sum_n [ D_{n,i} F^1(n,i,j,k)
    !!                                + D_{n,j} F^2(n,i,j,k)
    !!                                + D_{n,k} F^3(n,i,j,k) ]
    !!
    !! where D is the standard derivative matrix (dMatrix) and
    !! F^idir(n,i,j,k,...) = interior(n,i,j,k,iEl,iVar,idir) stores the
    !! two-point flux between nodes i (j or k) and n in the idir-th direction.
    !!
    !! The interior array is assumed to hold contravariant (J-scaled) two-point
    !! fluxes.  To obtain the physical divergence, divide the result by the
    !! element Jacobian J(i,j,k,iEl).
    implicit none
    class(TwoPointVector3D_t),intent(in) :: this
    real(prec),intent(out) :: df(1:this%N+1,1:this%N+1,1:this%N+1,1:this%nElem,1:this%nVar)
    ! Local
    integer :: i,j,k,nn,iEl,iVar
    real(prec) :: dfLoc

    do concurrent(i=1:this%N+1,j=1:this%N+1,k=1:this%N+1, &
                  iEl=1:this%nElem,iVar=1:this%nVar)

      dfLoc = 0.0_prec
      do nn = 1,this%N+1
        dfLoc = dfLoc+ &
                this%interp%dSplitMatrix(nn,i)*this%interior(nn,i,j,k,iEl,iVar,1)+ &
                this%interp%dSplitMatrix(nn,j)*this%interior(nn,i,j,k,iEl,iVar,2)+ &
                this%interp%dSplitMatrix(nn,k)*this%interior(nn,i,j,k,iEl,iVar,3)
      enddo
      df(i,j,k,iEl,iVar) = 2.0_prec*dfLoc

    enddo

  endsubroutine Divergence_TwoPointVector3D_t