Many elliptic partial differential equations can be solved numerically with near optimal efficiency through the uses of adaptive refinement and multigrid solution techniques. This paper presents a more unified approach to the combined process of adaptive refinement and multigrid solution which can be used with high order finite elements. Refinement is achieved by the bisection of pairs of triangles, corresponding to the addition of one or more basis functions to the approximation space. An approximation of the resulting change in the solution is used as an error indicator. The multigrid iteration uses red-black Gauss-Seidel relaxation with local black relaxations. The grid transfers use the change between the nodal and hierarchical bases. This multigrid iteration requires only O(N) operations, even for highly nonuniform grids, and is defined for any finite element space. The full multigrid method is an optimal blending of the processes of adaptive refinement and multigrid iteration. To minimize the number of operations required, the duration of the refinement phase is based on increasing the dimension of the approximation space by the largest possible factor, given the error reduction of the multigrid iteration. The algorithm (i) uses only O(N) operations, (ii) solves the discrete system to the accuracy of the discretization error, and (iii) achieves optimal convergence of the discretization error in the presence of singularities. Numerical experiments confirm this for linear, quadratic, and cubic elements.