  Numerical studies of the dynamics of gravitational systems, e.g., black
hole-neutron star systems, require physical and constraint-satisfying initial
data. In this article, we present the newly developed pseudo-spectral code
Elliptica, an infrastructure for construction of initial data for various
binary and single gravitational systems of all kinds. The elliptic equations
under consideration are solved on a single spatial hypersurface of the
spacetime manifold. Using coordinate maps, the hypersurface is covered by
patches whose boundaries can adapt to the surface of the compact objects. To
solve elliptic equations with arbitrary boundary condition, Elliptica deploys a
Schur complement domain decomposition method with a direct solver. In this
version, we use cubed sphere coordinate maps and the fields are expanded using
Chebyshev polynomials of the first kind. Here, we explain the building blocks
of Elliptica and the initial data construction algorithm for a black
hole-neutron star binary system. We perform convergence tests and evolve the
data to validate our results. Within our framework, the neutron star can reach
spin values close to breakup with arbitrary direction, while the black hole can
have arbitrary spin with dimensionless spin magnitude $\sim 0.8$.
