Many elliptic partial differential equations can be solved numerically with near optimal efficiency through the uses of adaptive refinement and multigrid solution techniques. It is our goal to develop a more unified approach to the combined process of adaptive refinement and multigrid solution which can be used with high order finite elements. The basic step of the refinement process is the bisection of a pair of triangles, which corresponds 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 to determine which triangles to divide. The multigrid iteration uses a red-black Gauss-Seidel relaxation in which the black relaxations are used only locally. 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. So as to minimize the number of operations required, the duration of the refinement phase is based on increasing the dimension of the approximation space by some fixed factor which is determined to be the largest possible for the given error-reducing power of the multigrid iteration. The result is an algorithm which (i) uses only O(N) operations with a reasonable constant of proportionality, (ii) solves the discrete system to the accuracy of the discretization error, (iii) is able to achieve the optimal order of convergence of the discretization error in the presence of singularities. Numerical experiments confirm this for linear, quadratic and cubic elements. It is believed that the method can also be applied to more practical problems involving systems of PDE's, time dependence, and three spacial dimensions.