using Ferrite

grid = generate_grid(Triangle, (4, 4), Vec((-2.0,-2.0)), Vec((2.0,2.0)))
addcellset!(grid, "hole", x->norm(x) ≤ 1.0)
dh = DofHandler(grid)
sdh = SubDofHandler(dh, getcellset(grid, "hole"))
add!(sdh, :v, Lagrange{RefTriangle,1}())
close!(dh)

nodes = Vec{2, Float64}.([[0.0, -1.0], [0.0, 0.0], [-1.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
ph = PointEvalHandler(grid, nodes; warn=true)
result = Ferrite.evaluate_at_points(ph, dh, zeros(ndofs(dh)))
#5-element Vector{Float64}:
# NaN
# 0.0
# NaN
# NaN
# NaN

VTKGridFile("peh-subdomain-mwe", grid) do vtk
    Ferrite.write_cellset(vtk, grid)
    Ferrite.write_solution(vtk, dh, result)
end