Based on the hierarchical structure, two types of the p-version of the finite element code are developed for solving the two-dimensional neutron diffusion equations. One is the conventional p-type FEM, and the reactor domain is discretized into finite elements, and each element matrices are assembled and solved. In the other, with the domain decomposition approach, the reactor domain is decomposed into subdomains and each subdomain is solved independently. They are coupled by use of incoming and outgoing partial currents along the interfaces, which are obtained analytically from fluxes. Based on the multi-p V-cycle method, multilevel acceleration is implemented into both of the methods. The methods it ere tested on the IAEA-2D benchmark problem, the initial core of UIckin Unit 1, and a MOX-fuel loaded core. It turns out that the methods can generate accurate core multiplication factor and assembly average quantities.