Finite element methods play a vital role in computational materials science. Moreover, developing efficient high order finite element solvers, and visualizing the finite element solution effectively are crucial in modern scientific applications. In this talk, we discuss recent developments in p-hierarchical basis finite element methods. We first present a newly developed multigrid solver that uses space decomposition as a smoother. The proposed solver combines the features of p-hierarchical bases, multigrid and space decomposition effectively. Unlike the existing p-version multigrid solvers that perform smoothing only on a subspace of the multigrid space, we adopt a global smoothing strategy at each multigrid level via space decomposition, in a multiplicative Schwarz manner. As the function spaces, space decomposition and the multigrid levels are p-hierarchical, this solver is particularly advantageous on p-nonconforming meshes. We will analyse the error operator and will establish the p-version strengthened Cauchy Buniakowskii Schwarz inequality constants yielding the first convergence estimates. Furthermore, we present numerical results to demonstrate the effectiveness of this solver. The numerical experiments show that the proposed solver outperforms the benchmark p-version space decomposition solver in terms of the convergence rates as well as the computational complexity. We finally discuss the visualization of the p-hierarchical basis finite element solution, which is challenging with the existing software such as Paraview. Visualization on ParaView with a VTK file format is done by a linear or quadratic interpolation of the solution values at specific nodal locations. However, the p-hierarchical basis finite element solution does not support it directly. In this talk, we present a novel strategy to overcome this barrier. This new strategy comprises of two key steps: “p-hierarchical to nodal projection” and “higher order to lower order projection. We will discuss the method and will present results obtained solving Poisson problems in two and three dimensions.