A general high-order finite element method for solving partial differential equations, called the COnstraints Oriented Library (COOL) method, is presented. This hp approach takes into account the underlying nature of the corresponding physical problem, and thus avoids the generation of nonphysical solutions. In the COOL method, all terms in a variational form are represented by the same functional dependence and by the same regularity, thus eliminating regularity constraints imposed by standard numerical methods. External constraints, such as the incompressibility condition appearing in the Maxwell or Navier-Stokes equations, can then be satisfied identically and are eliminated algebraically. This reduces the number of variables, and leads to well-conditioned matrix problems. The consequence is that only physically relevant solutions remain. The COOL method also satisfies automatically internal constraints, such as occur for the grad(div) and curl(curl) operators, and this for any geometry. This approach can be applied to a wide range of physical problems, including fluid flows, electromagnetics, material sciences, ideal linear magnetohydrodynamic stability analysis, and Alfvén wave heating of fusion plasmas. Results obtained by applying the COOL method to the grad(div) and curl(curl) operators, the Stokes problem, and the steady and unsteady Navier-Stokes equations are presented.