From 1065dbd2f84f9198b51cd861c5a424299c021ee8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Wed, 8 May 2013 12:57:38 +0200
Subject: [PATCH] Add subset() free function.

Also add overload for UpwindSelector::select() taking constants (not AD variables).
---
 AutoDiffHelpers.hpp | 43 ++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 42 insertions(+), 1 deletion(-)

diff --git a/AutoDiffHelpers.hpp b/AutoDiffHelpers.hpp
index c0c7eaaf8..20c2f4e97 100644
--- a/AutoDiffHelpers.hpp
+++ b/AutoDiffHelpers.hpp
@@ -194,16 +194,57 @@ namespace {
             return xf;
         }
 
-        /// Apply selector to single per-cell quantity.
+        /// Apply selector to single per-cell ADB quantity.
         ADB select(const ADB& xc) const
         {
             return select_*xc;
         }
 
+        /// Apply selector to single per-cell constant quantity.
+        typename ADB::V select(const typename ADB::V& xc) const
+        {
+            return (select_*xc.matrix()).array();
+        }
+
     private:
         typename ADB::M select_;
     };
 }
 
+namespace {
+
+    template <typename Scalar, class IntVec>
+    Eigen::SparseMatrix<Scalar>
+    constructSubsetSparseMatrix(const int full_size, const IntVec& indices)
+    {
+        typedef Eigen::Triplet<Scalar> Tri;
+        const int subset_size = indices.size();
+        std::vector<Tri> triplets(subset_size);
+        for (int i = 0; i < subset_size; ++i) {
+            triplets[i] = Tri(i, indices[i], 1);
+        }
+        Eigen::SparseMatrix<Scalar> sub(subset_size, full_size);
+        sub.setFromTriplets(triplets.begin(), triplets.end());
+        return sub;
+    }
+} // anon namespace
+
+
+template <typename Scalar, class IntVec>
+AutoDiff::ForwardBlock<Scalar>
+subset(const AutoDiff::ForwardBlock<Scalar>& x,
+       const IntVec& indices)
+{
+    return ::constructSubsetSparseMatrix<Scalar>(x.value().size(), indices) * x;
+}
+
+
+template <typename Scalar, class IntVec>
+Eigen::Array<Scalar, Eigen::Dynamic, 1>
+subset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
+       const IntVec& indices)
+{
+    return (::constructSubsetSparseMatrix<Scalar>(x.size(), indices) * x.matrix()).array();
+}
 
 #endif // OPM_AUTODIFFHELPERS_HEADER_INCLUDED