EvalScalar_3D_Points_t Subroutine

  subroutine EvalScalar_3D_Points_t(this,scalar,values)
    !! Evaluate a 3D MappedScalar at all located points by tensor-product
    !! Lagrange interpolation. Points with elements(p) == 0 receive zero.
    !! When LocatePoints has cached the basis at the matching degree, this
    !! routine reuses it.
    implicit none
    class(Points_t),intent(in) :: this
    class(MappedScalar3D_t),intent(in) :: scalar
    real(prec),intent(out) :: values(1:this%nPoints,1:scalar%nVar)
    ! Local
    integer :: p,iEl,iVar,i,j,k,N
    real(prec) :: fijk,fij,fi
    real(prec),allocatable :: lS(:),lT(:),lU(:)
    logical :: useCache

    if(this%nDim /= 3) then
      print*,"SELF_Points_t::EvaluateScalar (3D): nDim must be 3"
      stop 1
    endif

    N = scalar%interp%N
    useCache = (this%nCached == N) .and. associated(this%lS_cache) .and. &
               associated(this%lT_cache) .and. associated(this%lU_cache)

    allocate(lS(0:N),lT(0:N),lU(0:N))
    values = 0.0_prec

    do p = 1,this%nPoints
      iEl = this%elements(p)
      if(iEl <= 0) cycle

      if(useCache) then
        lS = this%lS_cache(:,p)
        lT = this%lT_cache(:,p)
        lU = this%lU_cache(:,p)
      else
        lS = scalar%interp%CalculateLagrangePolynomials(this%coordinates(p,1))
        lT = scalar%interp%CalculateLagrangePolynomials(this%coordinates(p,2))
        lU = scalar%interp%CalculateLagrangePolynomials(this%coordinates(p,3))
      endif

      do iVar = 1,scalar%nVar
        fijk = 0.0_prec
        do k = 1,N+1
          fij = 0.0_prec
          do j = 1,N+1
            fi = 0.0_prec
            do i = 1,N+1
              fi = fi+lS(i-1)*scalar%interior(i,j,k,iEl,iVar)
            enddo
            fij = fij+lT(j-1)*fi
          enddo
          fijk = fijk+lU(k-1)*fij
        enddo
        values(p,iVar) = fijk
      enddo
    enddo

    deallocate(lS,lT,lU)

  endsubroutine EvalScalar_3D_Points_t