In this contribution, we present a novel polygonal finite element method applied to hyperelastic analysis. For generating polygonal meshes in a bounded period of time, we use the adaptive Delaunay tessellation (ADT) proposed by Constantinu et al. [1]. ADT is an unstructured hybrid tessellation of a scattered point set that minimally covers the proximal space around each point. In this work, we have extended the ADT to nonconvex domains using concepts from constrained Delaunay triangulation (CDT). The proposed method is thus based on a constrained adaptive Delaunay tessellation (CADT) for the discretization of domains into polygonal regions. We involve the metric coordinate (Malsch) method for obtaining the interpolation over convex and nonconvex domains. For the numerical integration of the Galerkin weak form, we resort to classical Gaussian quadrature based on triangles. Numerical examples of two-dimensional hyperelasticity are considered to demonstrate the advantages of the polygonal finite element method. © 2018, © 2018 Taylor & Francis Group, LLC.