Abstract: To handle the control rod cusping effect in pressurized water reactor (PWR) fuel management calculation, the variational nodal method (VNM) in the fuel management calculation code system NECP-Bamboo has been extended to tread the heterogeneous cross section distribution by expanding the volumetric cross sections into piece-wise polynomials in the early work. However, the partially inserted control rods also introduces heterogeneous discontinuity factor (DF) on nodal interface. Thus, in this paper, an ultimate solution is proposed to fully handle this problem. Firstly, the surface integral in the VNM is modified to contain the discontinuity of neutron flux, incorporating a continuous discontinuity factor in that term. Secondly, the surface DF is expanded into the sum of pieces-wise orthogonal polynomials to construct the nodal response matrixes. Comparing with current representative re-homogenization methods, the application numerical results of the BEAVRS benchmark problem demonstrate the effectiveness of the heterogeneous VNM with heterogeneous DF. It can eliminate the cusping effect by providing more accurate differential control rode worth curves and pin power distributions.