A new methodology is presented for model reduction of linear parabolic partial differential equations (PDEs) on general geometries using multivariate splines on triangulations. State-space descriptions are derived that can be used for control design. This method uses Galerkin projection with B-splines to derive a finite set of ordinary differential equations (ODEs). Any desired smoothness conditions between elements as well as the boundary conditions are flexibly imposed as a system of side constraints on the set of ODEs. Projection of the set of ODEs on the null space of the system of side constraints naturally produces a reduced-order model that satisfies these constraints. This method can be applied for both in-domain control and boundary control of parabolic PDEs with spatially varying coefficients on general geometries. The reduction method is applied to design and implement feedback controllers for stabilisation of a 1-D unstable heat equation and a more challenging 2-D reaction-convection-diffusion equation on an irregular domain. It is shown that effective feedback stabilisation can be achieved using low-order control models.