I am implementing with finite element method the equation of elastostatic. The weak form gives a matrix to be computed in the interior of the elements, and a right hand side vector to be computed along the boundary of the domain. This vector should be split into two parts: one for the Neumann boundary conditions, and one for the Dirichlet boundary conditions. I do not have any Neumann boundary conditions. In the process of assembly of the matrix and rhs, I fix the degrees of freedom subject to Dirichlet boundary conditions setting 1 in the corresponding element of the diagonal of the matrix, and setting the Dirichlet value in the corresponding row of the rhs. Since I do not have any traction applied on the boundary, our rhs is made only of the fixed Dirichlet values and zeroes otherwise. The result is that the solution of the assembled system is always zero for all the interior degrees of freedom. What am I doing wrong?