Many problems in geometric modeling can be described using variational formulations that define the smoothness of the shape and its behavior w.r.t. the posed modeling constraints. For example, high-quality C2 surfaces that obey boundary conditions on positions, tangents and curvatures can be conveniently defined as solutions of high-order geometric PDEs; the advantage of such a formulation is its conceptual representation-independence. In practice, solving high-order problems efficiently and accurately for surfaces approximated by meshes is notoriously difficult. Classical FEM approaches require high-order elements which are complex to construct and expensive to compute. Recent discrete geometric schemes are more efficient, but their convergence properties are hard to analyze, and they often lack a systematic way to impose boundary conditions. In this paper, we present an approach to discretizing common PDEs on meshes using mixed finite elements, where additional variables for the der...