We present a novel approach for the construction of basis functions to be employed in selective or adaptive h-refined finite element applications with arbitrary-level hanging node configurations. Our analysis is not restricted to 1-irregular meshes, as it is usually done in the literature, allowing our results to be applicable to a broader class of local refinement strategies. The proposed method does not require the solution of any linear system to obtain the constraints necessary to enforce continuity of the basis functions and it can be easily implemented. A mathematical analysis is carried out to prove that the proposed basis functions are continuous and linearly independent. Finite element spaces are then defined as the spanning sets of such functions, and the implementation of a multigrid algorithm built on these spaces is discussed. A spectral analysis of the multigrid algorithm highlights superior convergence properties of the proposed method over existing strategies based on a local smoothing procedure. Finally, linear and nonlinear numerical examples are tested to show the robustness and versatility of the multigrid algorithm.restrict the residual to V J r i J = R J r i J ; 6. Let D J be an "easy to invert" approximation of the matrix A J , then solve D J w i J = r i J , and set u i+1 J = u i J + w i J ; 7. if w i J ≤ ε exit, otherwise set i = i + 1 and go back to step 3. A few remarks on the above algorithm:We construct the matrix A J and the residual r J in V J , since the assembling is done element-wise following the standard finite element approach for unstructured grids. Moreover, let x i , x j in X l J be two nodal points belonging to a generic element E of the triangulation T l J . The two pairs of test functions ϕ J,i , ϕ J,j and ϕ l J,i , ϕ l J,j EUGENIO AULISA, GIACOMO CAPODAGLIO, GUOYI KE PETSc [4] as a black box, providing as inputs the matrices A k for k = J, . . . , 0, the matrices Q k for k = J, . . . , 1 and the residual vector r i J , and letting PETSc compute w i J . This algebraic approach is convenient and elegant since all the burden of dealing with the hanging nodes at different level triangulations T k has been incorporated in the prolongation operator Q k and no ad-hoc multigrid solver or preconditioner needs to be implemented.The above algorithm can be easily extended to nonlinear problems. In this case the residual r i J at step 4 is a more complex function of u i J and D J is an "easy to invert" approximation of the tangent matrixThis corresponds to a Newton-Raphson iterative scheme [24]. Such a scheme is used in the last example where a Navier-Stokes flow is considered.