We propose a methodology aimed at the generation and adaptation of two dimensional curvilinear anisotropic meshes. The metric tensor field used to generate the $P^2$ triangular elements is defined based on an interpolation error estimate, using recovered derivatives of order up to three of a numerical solution computed e.g. with the finite element method. A quasi-unit straight-sided mesh is first generated by spawning points in the computational domain in such a way that two neighbouring points lie within the range of $[1/\sqrt{2},\sqrt{2}]$ when computing distances with respect to the metric field, then by connecting them using a standard anisotropic mesh generator. Straight-sided edges are then curved by moving the mid-edge nodes to approximate the geodesic between the edges's extremities while ensuring the validity of the elements at all time. Topological operations such as curvilinear edge swaps and curvilinear small polygon reconnection (CSPR) are finally applied to the curved triangulation to further increase mesh quality, yielding a quasi-unit curved mesh with only valid elements. For a comparable complexity, adapted meshes exhibit a reduced interpolation error compared to straight-sided anisotropic meshes. Numerical applications, such as tracer advection and Von Karman streets, are presented to illustrate the method.