Multilevel representations and mesh reduction techniques have been used for accelerating the processing and the rendering of large datasets representing scalar or vector valued functions defined on complex 2 or 3 dimensional meshes. We present a method based on finite element approximations which combines these two approaches in a new and unique way that is conceptually simple and theoretically sound. The main idea is to consider mesh reduction as an approximation problem in appropriate finite element spaces. Starting with a very coarse triangulation of the functional domain a hierarchy of highly non-uniform tetrahedral (or triangular in 2D) meshes is generated adaptively by local refinement. This process is driven by controlling the local error of the piecewise linear finite element approximation of the function on each mesh element. A reliable and efficient computation of the global approximation error combined with a multilevel preconditioned conjugate gradient solver are the key components of the implementation. In order to analyze the properties and advantages of the adaptively generated tetrahedral meshes we implemented two volume visualization algorithms: an iso-surface extractor and a ray-easter. Both algorithms, while conceptually simple, show significant speedups over conventional methods delivering comparable rendering quality from adaptively compressed datasets.