Summary: The solution accuracy based on finite element interpolation strongly depends on mesh quality. Thus, the numerical manifold method advocates that if a finite element mesh is selected as the mathematical cover, the mesh should be uniform, where all elements are congruent squares or isosceles right triangles for two-dimensional problems. However, the elements close to the points of stress concentrations or singularities must be smaller than those further away from these points. This leads to an irregular mesh in which a few nodes of certain small elements are hung on the edges of adjacent larger elements. This study implements local refinement with arbitrary \(k\)-irregular meshes by selecting piecewise polynomials as the local approximations over physical patches and enforcing appropriate constraints on the piecewise polynomials. All the elements are ordinary isoparametric elements, and elements with variable numbers of nodes on edges are not required. No linear dependency issue exists. Furthermore, refinement with any number of mesh layers is developed so that solution error can be distributed as uniformly as possible. Typical examples with different singularities are designed to illustrate the effectiveness of the proposed procedure.