The best thing to do here is to actually reduce the size of your problem after you have assembled the full system matrices. Let's assume you've gotten to that point, so you have
$$
\mathbf{\tilde{K}d}=\mathbf{\tilde{F}}
$$
where $\mathbf{K}$ is a $N_{dof} \times N_{dof}$ matrix of stiffness coefficients, $\mathbf{d}$ is your full vector of nodal degrees of freedom (most of which are known in your case), and $\mathbf{F}$ is a $N_{dof} \times 1$ vector of "forces" which arose from integration of the source terms and Neumann boundary conditions.
Rather than solving this very large system, you want to reduce the system down to ONLY the unknowns in $\mathbf{d}$! The method I will describe here is very general, and is not limited to specially-shaped domains.
To do this, first (for the sake of simplicity) that the "top" half of your $\mathbf{d}$ vector contains your known values from the Dirichlet boundary conditions, and the bottom half contains your unknowns. Immediately, we see that we can simply ignore your first $N_{dof}/2$ equations. In practice, this means chopping off the first $N_{dof}/2$ rows of $\mathbf{\tilde{K}}$ and $\mathbf{\tilde{F}}$. Now we're left with a more manageable set of equations
$$
\mathbf{\hat{K}d}=\mathbf{\hat{F}}
$$
where $\mathbf{\hat{K}}$ is a $\frac{N_{dof}}{2} \times N_{dof}$ matrix, and $\mathbf{\hat{F}}$ is a vector of length $N_{dof}/2$. We haven't yet done anything with $\mathbf{d}$.
Now, we're going to reduce $\mathbf{d}$ down to only having the unknowns on the left side of the equation. Notice that the first $N_{dof}/2$ columns of $\mathbf{\hat{K}}$ are being multiplied with known values of $\mathbf{d}$. We will take advantage of this fact and move all of the known values to the right-hand side of the equation. That statement may have been a little cryptic, but the best way to illustrate it is with a little bit of code. I'm going to use a MATLAB-like syntax and take advantage of matrix slicing.
for ii=1:Ndof //loop through d
if(isDirichletDOF(ii)) //proceed if d(ii) is known
Fhat = Fhat - d(ii)*Khat(:,ii);
end
end
In the small snippet above, I only moved over columns associated with KNOWN values (ie, all of the dirichlet boundary values) and did not touch any of the others. Now that you've done this, you no longer need to worry about the known values of $\mathbf{d}$ any more, so just chop it down to a vector of only unknowns. Simultaneously, make a new matrix $\mathbf{K_{final}}$ from $\mathbf{\hat{K}}$ that uses only the columns that were not associated with known values. Now you finally have the system that you want to solve.
$$
\mathbf{K_{final}d_{unknowns}} = \mathbf{\hat{F}_{modified}}
$$
In your particular application, this will be MUCH MUCH faster than some other methods of applying boundary conditions. From here, you pick your favorite linear solver and go to town. I won't try to predict actual performance gains since that depends so much on the linear solver you use, but $\mathbf{K}$ is now $1/10,000$ the size that it used to be (if, for some reason you were not using a sparse matrix).
Obviously, in practice, you don't know your node numbering ahead of time (for non-trivial geometries) so you just need to do some bookkeeping to determine which rows to get rid of and which columns to move to the right hand side of the equation.
