There is mathematical justification for setting Dirichlet boundary degrees of freedom to a value. However, you should adjust your variational form accordingly. If you are looking at a general problem, say:
Find $u\in\mathcal{U}$ such that
$a(u,w)=l(w) \ \ \forall w\in\mathcal{V}$
where
$\mathcal{U}=\{u:\int \nabla u^2 < \infty, u=g\text{ on }\Gamma_D\}$
$\mathcal{V}=\{u:\int \nabla u^2 < \infty, u=0\text{ on }\Gamma_D\}$
Instead we can write $u = v + g$ where $v\in\mathcal{V}$ and $g$ is the Dirichlet condition. Then the variational form becomes
$a(v+g,w)=l(w)$
or by using the linearity of $a(.,.)$
$a(v,w)=l(w)-a(g,w)$
In a finite element code, you can form your element stiffness matrix as if there were no boundary conditions. Then you take the column of the local matrix which corresponds to the Dirichlet boundary condition, scale it by the coefficient you want to enforce, and subtract it from the right-hand-side. This is the discrete form of what I wrote above, $-a(g,w)$. Then you zero out that column and the corresponding Dirichlet row, placing a 1 in the diagonal and the coefficient you wish to enforce. This decouples the equation from the system and yet sets the value you wish to enforce.
I recommend The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, by Tom Hughes. He has an expanded discussion of this issue starting on page 8.
