This article presents a tetrahedral mesh adaptivity algorithm for three‐dimensional elliptic partial differential equations (PDEs) using finite element methods. The main issues involved are the mesh size and mesh quality, which have great influence on the accuracy of the numerical solution and computational cost. The first issue is addressed by a posteriori error estimator based on superconvergent gradient recovery. The second issue is solved by constrained centroidal Voronoi–Delaunay tessellations (CCVDT), which guarantees good quality tetrahedrons over a large class of mesh domains even, if the grid size varies a lot at any particular refinement level. The CCVDT enjoys the energy equidistribution property so that the errors are very well equidistributed with properly chosen sizing field (density function). And with this good property, a new refinement criteria is raised which is different from the traditional bisection refinement. © 2013 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 30: 1633–1653, 2014