# Turn a region and points inside that region into a mesh

I have a region defined by a polygon and a few points that are inside that region. I'd like to create a mesh where its boundary is the boundary of the region and the vertices interior to the mesh are the points that I created (or more if needed for compatibility). How can I create a mesh like this?

region = Polygon[{{0, 0}, {1, 0}, {1, -1}, {2, -1}, {2, 0}, {3,
0}, {3, -2}, {0, -2}}];
points = RandomPoint[region, 1000];
mesh = createMeshFromRegionAndPoints[region, points];  (* This is the function I need to create *)


I know I can create a mesh from the region with DiscretizeRegion, but this won't capture the interior vertices. I can also create a DelaunayMesh of the interior points, but given my region is non-convex, the boundary will be exceeded. If I could somehow "trim" the edges of the DelaunayMesh using the region, then that would get me most of what I want (I'm trying to make a ListDensityPlot that draws on my region).

You could use IncludePoints from the Finite Element Package so that the required interior vertices are in the mesh (along with other vertices to create a valid mesh). Here is an example:

region = Polygon[{{0, 0}, {1, 0}, {1, -1}, {2, -1}, {2, 0}, {3,
0}, {3, -2}, {0, -2}}];
points = RandomPoint[region, 1000];
Needs["NDSolveFEM"];