From 29ba859a805a67a08d7b6a57e51e3d58a644c32d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= <Bard.Skaflestad@sintef.no>
Date: Mon, 12 Mar 2012 16:18:02 +0100
Subject: [PATCH] Re-implement computeTotalMobilityOmega() in terms of
 computePhaseMobilities().

---
 opm/core/utility/miscUtilities.cpp | 40 ++++++++++++++----------------
 1 file changed, 19 insertions(+), 21 deletions(-)

diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp
index cd840e267..f99f8f331 100644
--- a/opm/core/utility/miscUtilities.cpp
+++ b/opm/core/utility/miscUtilities.cpp
@@ -119,27 +119,25 @@ namespace Opm
 				   std::vector<double>& totmob,
 				   std::vector<double>& omega)
     {
-	int num_cells = cells.size();
-	int num_phases = props.numPhases();
-	totmob.resize(num_cells);
-	omega.resize(num_cells);
-	ASSERT(int(s.size()) == num_cells*num_phases);
-	std::vector<double> kr(num_cells*num_phases);
-	props.relperm(num_cells, &s[0], &cells[0], &kr[0], 0);
-	const double* mu = props.viscosity();
-	for (int cell = 0; cell < num_cells; ++cell) {
-	    totmob[cell] = 0.0;
-	    for (int phase = 0; phase < num_phases; ++phase) {	
-		totmob[cell] += kr[num_phases*cell + phase]/mu[phase];
-	    }
-	}
-	const double* rho = props.density();
-	for (int cell = 0; cell < num_cells; ++cell) {
-	    omega[cell] = 0.0;
-	    for (int phase = 0; phase < num_phases; ++phase) {	
-		omega[cell] += rho[phase]*(kr[num_phases*cell + phase]/mu[phase])/totmob[cell];
-	    }
-	}
+        std::vector<double> pmobc;
+
+        computePhaseMobilities(props, cells, s, pmobc);
+
+        const std::size_t                 np = props.numPhases();
+        const std::vector<int>::size_type nc = cells.size();
+
+        std::vector<double>(cells.size(), 0.0).swap(totmob);
+        std::vector<double>(cells.size(), 0.0).swap(omega );
+
+        const double* rho = props.density();
+        for (std::vector<int>::size_type c = 0; c < nc; ++c) {
+            for (std::size_t p = 0; p < np; ++p) {
+                totmob[ c ] += pmobc[c*np + p];
+                omega [ c ] += pmobc[c*np + p] * rho[ p ];
+            }
+
+            omega[ c ] /= totmob[ c ];
+        }
     }
 
     void computePhaseMobilities(const Opm::IncompPropertiesInterface& props,