We propose in this paper an anisotropic, adaptive, finite element algorithm for steady, linear advection-diffusion-reaction problems with strong anisotropic features. The error analysis is based on the dual weighted residual methodology, allowing us to perform "goal-oriented" adaptation of a certain functional J(u) of the solution and derive an "optimal" metric tensor for local mesh adaptation with linear and quadratic finite elements. As a novelty, and to evaluate the weights of the error estimator on unstructured meshes composed of anisotropic triangles, we make use of a patchwise, higher-order interpolation recovery readily extendable to finite elements of arbitrary order. We carry out a number of numerical experiments in two dimensions so as to prove the capabilities of the goal-oriented adaptive method. We compute the convergence rate and the effectivity index for a series of output functionals of the solution. The results show the good performance of the algorithm with linear as well as quadratic finite elements.Key words. adaptive finite element algorithm, unstructured triangular anisotropic meshes, DWR method, linear and quadratic finite elements, patch-wise higher-order interpolation recovery 1. Introduction. Nowadays, the number of scientific and technological problems requiring an efficient numerical approach is arguably higher than ever before. Whether in the field of fluid mechanics, biology, rheology, or computer-generated imagery, the solution to those challenging applications involving multiple scales usually comes from a well-rounded, adaptive scheme. Our goal with this paper is to make a contribution in such a direction by providing a means to design "economical" meshes suitable for an accurate numerical discretization of linear, advection-diffusion-reaction problems, based on an a posteriori error estimate analysis.Great effort has been devoted to research in adaptive algorithms since the seminal work by Babuska and Rheinboldt in 1978 (see [7] and [8]), where they first introduced the concept of an a posteriori error estimate to design local refinement procedures. Given that meshes with equilateral shapes behave best in problems exhibiting isotropic solutions, studies in local refinement have focused mainly on isotropic mesh adaptation [3]. The former notwithstanding, a wide range of applications display "directional features" (boundary and internal layers, singularities, shock waves, jets, vortices, fracture propagation, etc.) which isotropic meshes fail to handle efficiently: the number of elements provided by standard, isotropic adaptation would be lowered by an "anisotropic" adaptive procedure making use of an optimal triangulation where not only the size of the elements, but also their shape and orientation would
