I know of at least two ways of dealing with Dirichlet boundary conditions in finite-difference schemes (and, to a lesser extent, finite-element schemes). Here I'm thinking of solving Poisson's equation $\Delta u = f$, but the question is more general.
Treat the function values at boundary nodes as further unknown variables $u_i$, but add equations of the form $u_i = g(x_i)$ to the resulting system of equations. In total we'll have $n_\text{inner} + n_\text{boundary}$ equations and unknowns, where $n_\text{inner}$ is the number of nodes in the interior of the domain and $n_\text{boundary}$ is the number of nodes on the boundary.
Don't include variables for the fixed function values at boundary nodes. In total we'll have $n_\text{inner}$ equations and unknowns.
The two approaches have the following advantages and disadvantages:
Setting up the equations in the second approach is slightly more complicated, because nodes which are neighbours of boundary nodes have to be treated slightly differently. This seems to be the main selling point of the first approach.
With the first approach, the resulting system matrix might not be symmetric. (Symmetry is very nice for iterative solvers.)
The second approach is closer to the theoretical treatment of the undiscretized case: With zero Dirichlet conditions, we look for a solution in the Sobolev space $H_0^1(\Omega)$ instead of $H^1(\Omega)$.
As soon as more than one differential operator is involved, for instance in the Stokes equation, the first approach gets more messy. One has to introduce zero rows in some of the resulting matrices.
The question: Which of the two approaches is the current recommended best practice? Are there further advantages or disadvantages which I've overlooked?