  Polynomial spectral methods produce fast, accurate, and flexible solvers for
broad ranges of PDEs with one bounded dimension, where the incorporation of
general boundary conditions is well understood. However, automating extensions
to domains with multiple bounded dimensions is challenging because of
difficulties in imposing boundary conditions at shared edges and corners. Past
work has included various workarounds, such as the anisotropic inclusion of
partial boundary data at shared edges or approaches that only work for specific
boundary conditions. Here we present a general system for imposing boundary
conditions for elliptic equations on hypercubes. We take an approach based on
the generalized tau method, which allows for a wide range of boundary
conditions for many different spectral schemes. The generalized tau method has
the distinct advantage that the specified polynomial residual determines the
exact algebraic solution; afterwards, any stable numerical scheme will find the
same result. We can, therefore, provide one-to-one comparisons to traditional
collocation and Galerkin methods within the tau framework. As an essential
requirement, we add specific tau corrections to the boundary conditions, in
addition to the bulk PDE, which produce a unique set of compatible boundary
data at shared subsurfaces. Our approach works with general boundary conditions
that commute on intersecting subsurfaces, including Dirichlet, Neumann, Robin,
and any combination of these on all boundaries. The boundary tau corrections
can be made hyperoctahedrally symmetric and easily incorporated into existing
solvers. We present the method explicitly for the Poisson equation in two and
three dimensions and describe its extension to arbitrary elliptic equations
(e.g. biharmonic) in any dimension.
