When I introduce the values of the unknowns in the strong formulation, shouldn’t the equation be satisfied?
You don’t expect the strong form of the PDE to be satisfied exactly by the finite element solution. Considering the weighted-residual derivation, you might initially expect that the strong problem would be satisfied in an average sense over the domain, which is what you’re testing in your MWE. However, global constants are not part of the test function space, due to the DirichletBC, and the form fb1 is missing distributional contributions to the finite element solution’s second derivatives (as I discussed in more detail in my answer here).