A grid generation problem in two-dimensional domains is considered by using a quasi-conformal mapping of the parametric domain with a given square mesh onto the physical domain where the grid is required. To this end, a harmonic mapping is first applied, which, by the Radó theorem, is a diffeomorphism subject to some known conditions. However, the discrete harmonic mapping may produce folded meshes on a nonconvex domain with a strongly bent boundary. We demonstrate that it is caused by the truncation error. With the aim of controlling grid node location, an additional mapping is used. The Dirichlet problem for the universal elliptic partial differential equations is solved to construct the mapping. This allows to produce unfolded grids with a prescribed cell shape.