Interpolation on an unstructured mesh

Here is how to do it, just let ToElementMesh create the mesh:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[pts];
values = {#[[1]], #[[2]], #[[1]] #[[2]]^2} & /@ pts;
int = ElementMeshInterpolation[{mesh}, values[[All, 3]]];
{Length[pts], Length[mesh["Coordinates"]]}

{1001, 1001}

Plot3D[int[x, y], {x, y} \[Element] mesh]

Generallly speaking, given a set of points ToBoundayMesh will return a convex hull and ToElementMesh will return a Delaunay triangulation.

Try

mesh = ToElementMesh[DelaunayMesh@pts, MeshQualityGoal -> 0 ,"MeshOrder" -> 1 ];

This is a simple triangle mesh (no additional points!)

Length[pts]==Length[mesh["Coordinates"]]
(*True*) 

values = #[[1]] #[[2]]^2 & /@ mesh["Coordinates"];
fFE = ElementMeshInterpolation[{mesh}, values];
Plot3D[fFE[x, y], Element[{x, y}, mesh] ]

I think Interpolation does what you want under the hood (basically what @MarcoB said):

ifn = Interpolation[Transpose@{pts, values[[All, 3]]}, 
   InterpolationOrder -> 1];

emesh = ifn@"ElementMesh";
emesh["Wireframe"]

Note that it controlled the construction of the mesh in the way you wanted.

Update. You can construct the mesh this way:

mymesh = ToElementMesh[ConvexHullMesh@pts, MeshQualityGoal -> 0, 
   MaxCellMeasure -> Infinity, "IncludePoints" -> pts, 
   "MeshOrder" -> 1];
Normal@mymesh["Wireframe"] === Normal@emesh["Wireframe"]
(* True  *)

The coordinates are ordered differently from pts, so we need to permute the values accordingly to construct the ElementMeshInterpolation:

myvalues = 
  values[[All, 3]][[
    Ordering@ pts]][[
    InversePermutation@ Ordering@ mymesh@"Coordinates"]];
myifn = ElementMeshInterpolation[{mymesh}, myvalues];

Plot3D[myifn[x, y], {x, y} ∈ mymesh]

Check equivalence:

myifn[##] === ifn[##] & @@
 Transpose@RandomPoint[MeshRegion@emesh, 10000]
(*  True  *)