diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp
index 046dae38..ce04aa59 100644
--- a/examples/compute_tof.cpp
+++ b/examples/compute_tof.cpp
@@ -225,7 +225,8 @@ main(int argc, char** argv)
     // Solve time-of-flight.
     std::vector<double> tof;
     if (use_dg) {
-        Opm::TransportModelTracerTofDiscGal tofsolver(*grid->c_grid());
+        bool use_cvi = param.getDefault("use_cvi", false);
+        Opm::TransportModelTracerTofDiscGal tofsolver(*grid->c_grid(), use_cvi);
         transport_timer.start();
         tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof);
         transport_timer.stop();
diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
index 2d58ec91..718f2cf6 100644
--- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
+++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
@@ -20,10 +20,11 @@
 #include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/VelocityInterpolation.hpp>
 #include <opm/core/linalg/blas_lapack.h>
 #include <algorithm>
-#include <numeric>
 #include <cmath>
+#include <numeric>
 
 namespace Opm
 {
@@ -403,32 +404,6 @@ namespace Opm
 
 
 
-    // Initial version: only a constant interpolation.
-    static void interpolateVelocity(const UnstructuredGrid& grid,
-                                    const int cell,
-                                    const double* darcyflux,
-                                    const double* /*x*/,
-                                    double* v)
-    {
-        const int dim = grid.dimensions;
-        std::fill(v, v + dim, 0.0);
-        const double* cc = grid.cell_centroids + cell*dim;
-        for (int hface = grid.cell_facepos[cell]; hface < grid.cell_facepos[cell+1]; ++hface) {
-            const int face = grid.cell_faces[hface];
-            const double* fc = grid.face_centroids + face*dim;
-            double flux = 0.0;
-            if (cell == grid.face_cells[2*face]) {
-                flux = darcyflux[face];
-            } else {
-                flux = -darcyflux[face];
-            }
-            for (int dd = 0; dd < dim; ++dd) {
-                v[dd] += flux * (fc[dd] - cc[dd]) / grid.cell_volumes[cell];
-            }
-        }
-    }
-
-
 
     // ---------------   Methods of TransportModelTracerTofDiscGal ---------------
 
@@ -436,11 +411,19 @@ namespace Opm
 
     /// Construct solver.
     /// \param[in] grid      A 2d or 3d grid.
-    TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid)
+    /// \param[in] use_cvi   If true, use corner point velocity interpolation.
+    ///                      Otherwise, use the basic constant interpolation.
+    TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
+                                                                   const bool use_cvi)
         : grid_(grid),
           coord_(grid.dimensions),
           velocity_(grid.dimensions)
     {
+        if (use_cvi) {
+            velocity_interpolation_.reset(new VelocityInterpolationECVI(grid));
+        } else {
+            velocity_interpolation_.reset(new VelocityInterpolationConstant(grid));
+        }
     }
 
 
@@ -485,6 +468,7 @@ namespace Opm
         basis_.resize(num_basis);
         basis_nb_.resize(num_basis);
         grad_basis_.resize(num_basis*grid_.dimensions);
+        velocity_interpolation_->setupFluxes(darcyflux);
         reorderAndTransport(grid_, darcyflux);
     }
 
@@ -560,7 +544,7 @@ namespace Opm
                 quad.quadPtCoord(quad_pt, &coord_[0]);
                 DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
                 DGBasis::evalGrad(grid_, cell, degree_, &coord_[0], &grad_basis_[0]);
-                interpolateVelocity(grid_, cell, darcyflux_, &coord_[0], &velocity_[0]);
+                velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]);
                 const double w = quad.quadPtWeight(quad_pt);
                 for (int j = 0; j < num_basis; ++j) {
                     for (int i = 0; i < num_basis; ++i) {
diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp
index 0285e1f0..68307e1c 100644
--- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp
+++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp
@@ -21,15 +21,18 @@
 #define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED
 
 #include <opm/core/transport/reorder/TransportModelInterface.hpp>
+#include <boost/shared_ptr.hpp>
 #include <vector>
 #include <map>
 #include <ostream>
+
 struct UnstructuredGrid;
 
 namespace Opm
 {
 
     class IncompPropertiesInterface;
+    class VelocityInterpolationInterface;
 
     /// Implements a discontinuous Galerkin solver for
     /// (single-phase) time-of-flight using reordering.
@@ -45,7 +48,10 @@ namespace Opm
     public:
         /// Construct solver.
         /// \param[in] grid      A 2d or 3d grid.
-        TransportModelTracerTofDiscGal(const UnstructuredGrid& grid);
+        /// \param[in] use_cvi   If true, use corner point velocity interpolation.
+        ///                      Otherwise, use the basic constant interpolation.
+        TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
+                                       const bool use_cvi);
 
 
         /// Solve for time-of-flight.
@@ -72,7 +78,12 @@ namespace Opm
         virtual void solveMultiCell(const int num_cells, const int* cells);
 
     private:
+        // Disable copying and assignment.
+        TransportModelTracerTofDiscGal(const TransportModelTracerTofDiscGal&);
+        TransportModelTracerTofDiscGal& operator=(const TransportModelTracerTofDiscGal&);
+
         const UnstructuredGrid& grid_;
+        boost::shared_ptr<VelocityInterpolationInterface> velocity_interpolation_;
         const double* darcyflux_;   // one flux per grid face
         const double* porevolume_;  // one volume per cell
         const double* source_;      // one volumetric source term per cell