In the second form to get the values at specific points you can write sol/.{x->3,y->3}, or else redefine sol with arguments:
sol[x_,y_] = NDSolveValue[{op == 0, pbc, Subscript[\[CapitalGamma], Dc]},
u[x, y], {x, y} \[Element] \[CapitalOmega]]
It is similar to the difference between f = Function[x,x^2] and f=x^2. With the former you can call f[1], but not with the latter.
