Solution on the boundary points

Hello everyone,
I am working on an elasticity problem where I am calculating the anti-plane displacement. I want to calculate the energy for the points on the vertical line just before the PML on the left and right vertical boundaries. I wanted to know how to access those points on the mesh and find the values of the solution on those points. Please help. The code for the domain is attached.

geo = geom2d.SplineGeometry()
p1,p2,p3,p4,p5,p6,p7,p8 = [geo.AppendPoint(x,y) for x,y in [(-5,-1.5),(5,-1.5),(5,0),(4.5,0),(4.5,-1),(-4.5,-1),(-4.5,0),(-5,0)] ]
origin = geo.AppendPoint(0,0)

order = 6
maxh = 0.15

geo.Append (["line", p4, origin],bc="bc1",leftdomain=1,rightdomain=0)
geo.Append (["line", origin, p7],bc="bc2",leftdomain=1,rightdomain=0)

geo.Append (["line", p1, p2],leftdomain=2,rightdomain=0)
geo.Append (["line", p2, p3],leftdomain=2,rightdomain=0)
geo.Append (["line", p3, p4],leftdomain=2,rightdomain=0)
geo.Append (["line", p4, p5],leftdomain=2,rightdomain=1)
geo.Append (["line", p5, p6],leftdomain=2,rightdomain=1)
geo.Append (["line", p6, p7],leftdomain=2,rightdomain=1)
geo.Append (["line", p7, p8],leftdomain=2,rightdomain=0)
geo.Append (["line", p8, p1],leftdomain=2,rightdomain=0)

geo.SetMaterial(1, "inner")
geo.SetMaterial(2, "pmlregion")

mesh=Mesh(geo.GenerateMesh (maxh=maxh))

There is another doubt. I want to add the solution coming out of the NGSolve to a function defined on the domain. I have used the following piece of the code:


gfu = GridFunction(fes, name="u")
ui = gfu.vec = ui+a.mat.Inverse() * b.vec

Please let me know if this is correct. From the figures I am plotting using this approach, I think I am on the right track.

Thank you.