Vorticity formulations for the incompressible Navier-Stokes equations have certain advantages over primitive-variable formulations including the fact that the number of equations to be solved is reduced. However, the accurate implementation of the boundary conditions seems to continue to be an impediment to the acceptance and use of numerical methods based on vorticity formulations. Velocity boundary conditions can be implicitly satisfied by maintaining the kinematic compatibility of the velocity and vorticity fields as described by the generalized Helmholtz decomposition (GHD). This can be accomplished in one of two ways by either solving for boundary vorticity (leading to a Dirichlet boundary condition for the vorticity equation) or solving for boundary vortex sheet strengths (leading to a Neumann condition). In the past, vortex sheet strengths have often been determined by solving an over-specified set of linear equations. The over-specification arose because integral constraints were imposed on the vortex sheet strengths. These integral constraints are not necessary and typically are included to mitigate errors in determining the vortex sheet strengths themselves. Further, the constraints overspecify the linear system requiring least-squares solution techniques. To more accurately satisfy both components of the velocity boundary conditions, a Galerkin formulation is applied to the generalized Helmholtz decomposition. This formulation implicitly satisfies an integral constraint that is more general than many of the integral constraints that have been explicitly imposed. Two implementations of the Galerkin GHD are considered in the current work, one based on determining the boundary vorticity and one based on determining the boundary vortex sheet strengths. A finite element method (FEM) is implemented to solve the vorticity equation along with the boundary data generated from the GHD.