EvalScalar_2D_Points_t Subroutine

  subroutine EvalScalar_2D_Points_t(this,scalar,values)
    !! Evaluate a 2D MappedScalar at all located points by tensor-product
    !! Lagrange interpolation at the stored reference coordinates. Points with
    !! elements(p) == 0 receive a value of zero. When LocatePoints has cached
    !! the per-point basis at the matching polynomial degree, this routine
    !! reuses it and skips Lagrange-polynomial evaluation altogether.
    implicit none
    class(Points_t),intent(in) :: this
    class(MappedScalar2D_t),intent(in) :: scalar
    real(prec),intent(out) :: values(1:this%nPoints,1:scalar%nVar)
    ! Local
    integer :: p,iEl,iVar,i,j,N
    real(prec) :: fij,fi
    real(prec),allocatable :: lS(:),lT(:)
    logical :: useCache

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

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

    allocate(lS(0:N),lT(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)
      else
        lS = scalar%interp%CalculateLagrangePolynomials(this%coordinates(p,1))
        lT = scalar%interp%CalculateLagrangePolynomials(this%coordinates(p,2))
      endif

      do iVar = 1,scalar%nVar
        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,iEl,iVar)
          enddo
          fij = fij+lT(j-1)*fi
        enddo
        values(p,iVar) = fij
      enddo
    enddo

    deallocate(lS,lT)

  endsubroutine EvalScalar_2D_Points_t