We present a polygonal finite element method based on constrained adaptive Delaunay tessellation and conformal interpolants on arbitrary polygons. For mesh generation we use the adaptive Delaunay tessellation, an unstructured hybrid tessellation of a scattered point set that minimally covers the proximal space around each point, which is here extended to non-convex domains. Various types of polygonal interpolants are implemented. For the numerical integration of the Galerkin weak form we resort to classical Gaussian quadrature applied on triangular subdomains. The performance and efficiency of the interpolation and implementation are investigated for two dimensional elasticity in a stochastical analysis on random and regular meshes. © 2013 Elsevier Ltd. All rights reserved.