This work presents a formulation and analysis of a discontinuous least squares method (DLSM) for second order boundary value problems posed on two dimensional domains with curved boundaries. The method is developed in a broken Sobolev space framework over meshes composed of convex polygonal finite elements. In DLSM the global functional is constructed that incorporates the strong-form residuals in the domain together with boundary conditions and continuity conditions. The stationary point of that functional yields a symmetric variational problem that combines the Dirichlet and Neumann boundary conditions and controls the jumps of the solution and its normal flux on the mesh skeleton. In the proposed method, a wide range of basis functions can be employed within the finite elements. In this work, Chebyshev, Legendre, and Hermite polynomials, as well as their scaled variants, are used. The scaled bases improve numerical stability, particularly for high order approximations. The basis functions defined on each finite element are independent of the element’s shape and consequently, polygonal finite elements of arbitrary shape can be employed. In particular, elements with curved edges can be used when required. In general, curved polygonal elements are not necessary unless they are needed to accurately represent the curved outer boundary of the domain. For numerical integration, each polygonal element is subdivided into triangulars to which appropriate quadrature schemes are applied. When a triangle has only straight sides, standard numerical integration schemes are used. If a triangle has one curved side, numerical integration is performed using a scaled boundary approach [1]. In addition a collocation approach for the numerical integration is also considered [2], which naturally accommodates curved boundaries. The capabilities of the method are demonstrated through a series of numerical examples. These include benchmark Poisson problems featuring highly irregular solutions on strongly curved domains. In addition, the method is applied to the Stokes problem that combines the velocity and pressure fields. The presented examples illustrate the flexibility of the approach and its numerical robustness, even for very high order approximations.