  Approximate streamsurfaces of a 3D velocity field have recently been
constructed as isosurfaces of the closest first integral of the velocity field.
Such approximate streamsurfaces enable effective and efficient visualization of
vortical regions in 3D flows. Here we propose a variational construction of
these approximate streamsurfaces to remove the limitation of Fourier series
representation of the first integral in earlier work. Specifically, we use
finite-element methods to solve a partial-differential equation that describes
the best approximate first integral for a given velocity field. We use several
examples to demonstrate the power of our approach for 3D flows in domains with
arbitrary geometries and boundary conditions. These include generalized
axisymmetric flows in the domains of a sphere (spherical vortex), a cylinder
(cylindrical vortex), and a hollow cylinder (Taylor-Couette flow) as benchmark
studies for various computational domains, non-integrable periodic flows (ABC
and Euler flows), and Rayleigh-B\'enard convection flows. We also illustrate
the use of the variational construction in extracting momentum barriers in
Rayleigh-B\'enard convection.
