TransportModel* classes are now expected to have a custom solve() method. More:

- Using new solve() method in spu_2p.
 - solve() implemented in terms of protected superclass method reorderAndTransport().
 - Removed unused code being replaced by solve().

opm/core/transport/reorder/TransportModelInterface.cpp

@@ -0,0 +1,49 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/core/transport/reorder/TransportModelInterface.hpp>
+#include <opm/core/transport/reorder/twophasetransport.hpp>
+#include <opm/core/transport/reorder/reordersequence.h>
+#include <opm/core/grid.h>
+
+#include <vector>
+#include <cassert>
+
+
+void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux)
+{
+    // Compute sequence of single-cell problems
+    std::vector<int> sequence(grid.number_of_cells);
+    std::vector<int> components(grid.number_of_cells + 1);
+    int ncomponents;
+    compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents);
+
+    // Assume all strong components are single-cell domains.
+    assert(ncomponents == grid.number_of_cells);
+    for (int i = 0; i < grid.number_of_cells; ++i) {
+#ifdef MATLAB_MEX_FILE
+        if (interrupt_signal) {
+            mexPrintf("Reorder loop interrupted by user: %d of %d "
+                      "cells finished.\n", i, grid.number_of_cells);
+            break;
+        }
+#endif
+        solveSingleCell(sequence[i]);
+    }
+}

opm/core/transport/reorder/TransportModelInterface.hpp

@@ -20,16 +20,28 @@
 #ifndef OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED
 #define OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED
 
+class UnstructuredGrid;
+
 namespace Opm
 {
 
+    /// Interface for reordering transport models.
+    /// A transport model must provide the solveSingleCell()
+    /// method, and is expected to implement a solve() method
+    /// that will have an interface geared to the model's needs.
+    /// (The solve() method is therefore not virtual in this class).
+    /// The reorderAndTransport() method is provided as an
+    /// aid to implementing solve() in subclasses.
     class TransportModelInterface
     {
     public:
 	virtual ~TransportModelInterface() {}
 	virtual void solveSingleCell(int cell) = 0;
+    protected:
+	void reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux);
     };
 
+
 } // namespace Opm
 
 #endif // OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED

opm/core/transport/reorder/TransportModelTwophase.cpp

@@ -27,28 +27,35 @@ namespace Opm
 {
 
 
-    TransportModelTwophase::TransportModelTwophase(const UnstructuredGrid* grid,
-						   const Opm::IncompPropertiesInterface* props,
-						   const double* darcyflux,
+    TransportModelTwophase::TransportModelTwophase(const UnstructuredGrid& grid,
 						   const double* porevolume,
-						   const double* source,
-						   const double dt,
-						   double* saturation)
+						   const Opm::IncompPropertiesInterface& props)
 	: grid_(grid),
-	  props_(props),
-	  darcyflux_(darcyflux),
 	  porevolume_(porevolume),
-	  source_(source),
-	  dt_(dt),
-	  saturation_(saturation),
-	  fractionalflow_(grid->number_of_cells, 0.0)
+	  props_(props),
+	  darcyflux_(0),
+	  source_(0),
+	  dt_(0.0),
+	  saturation_(0),
+	  fractionalflow_(grid.number_of_cells, -1.0)
     {
-	if (props->numPhases() != 2) {
+	if (props.numPhases() != 2) {
 	    THROW("Property object must have 2 phases");
 	}
-	visc_ = props->viscosity();
+	visc_ = props.viscosity();
     }
 
+    void TransportModelTwophase::solve(const double* darcyflux,
+				       const double* source,
+				       const double dt,
+				       double* saturation)
+    {
+	darcyflux_ = darcyflux;
+	source_ = source;
+	dt_ = dt;
+	saturation_ = saturation;
+	reorderAndTransport(grid_, darcyflux);
+    }
 
     // Residual function r(s) for a single-cell implicit Euler transport
     //
@@ -73,17 +80,17 @@ namespace Opm
 	    outflux = tm.source_[cell] <= 0 ? -tm.source_[cell] : 0.0;
 	    dtpv    = tm.dt_/tm.porevolume_[cell];
 
-	    for (int i = tm.grid_->cell_facepos[cell]; i < tm.grid_->cell_facepos[cell+1]; ++i) {
-		int f = tm.grid_->cell_faces[i];
+	    for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
+		int f = tm.grid_.cell_faces[i];
 		double flux;
 		int other;
 		// Compute cell flux
-		if (cell == tm.grid_->face_cells[2*f]) {
+		if (cell == tm.grid_.face_cells[2*f]) {
 		    flux  = tm.darcyflux_[f];
-		    other = tm.grid_->face_cells[2*f+1];
+		    other = tm.grid_.face_cells[2*f+1];
 		} else {
 		    flux  =-tm.darcyflux_[f];
-		    other = tm.grid_->face_cells[2*f];
+		    other = tm.grid_.face_cells[2*f];
 		}
 		// Add flux to influx or outflux, if interiour.
 		if (other != -1) {
@@ -118,7 +125,7 @@ namespace Opm
     {
 	double sat[2] = { s, 1.0 - s };
 	double mob[2];
-	props_->relperm(1, sat, &cell, mob, 0);
+	props_.relperm(1, sat, &cell, mob, 0);
 	mob[0] /= visc_[0];
 	mob[1] /= visc_[1];
 	return mob[0]/(mob[0] + mob[1]);

opm/core/transport/reorder/TransportModelTwophase.hpp

@@ -33,21 +33,23 @@ namespace Opm
     class TransportModelTwophase : public TransportModelInterface
     {
     public:
-	TransportModelTwophase(const UnstructuredGrid* grid,
-			       const Opm::IncompPropertiesInterface* props,
-			       const double* darcyflux,
+	TransportModelTwophase(const UnstructuredGrid& grid,
 			       const double* porevolume,
-			       const double* source,
-			       const double dt,
-			       double* saturation);
+			       const Opm::IncompPropertiesInterface& props);
+
+	void solve(const double* darcyflux,
+		   const double* source,
+		   const double dt,
+		   double* saturation);
 
 	virtual void solveSingleCell(int cell);
 
     private:
-	const UnstructuredGrid* grid_;
-	const IncompPropertiesInterface* props_;
-	const double* darcyflux_;   // one flux per grid face
+	const UnstructuredGrid& grid_;
 	const double* porevolume_;  // one volume per cell
+	const IncompPropertiesInterface& props_;
+
+	const double* darcyflux_;   // one flux per grid face
 	const double* source_;      // one source per cell
 	double dt_;
 	double* saturation_;        // one per cell