Merge Formation Names coloring work and FlowDiagnostics integration

ApplicationCode/CMakeLists.txt

@@ -31,6 +31,7 @@ include_directories(
     ${Boost_INCLUDE_DIRS}
     ${ResInsight_SOURCE_DIR}/ThirdParty/custom-opm-common/opm-common/
     ${custom-opm-parser_SOURCE_DIR}/opm-parser/
+    ${custom-opm-flowdiagnostics_SOURCE_DIR}/opm-flowdiagnostics/
 
     ${CMAKE_CURRENT_SOURCE_DIR}
     ${CMAKE_CURRENT_SOURCE_DIR}/Adm
@@ -333,7 +334,8 @@ endif()
 
 set( LINK_LIBRARIES
     custom-opm-parser
-    
+    custom-opm-flowdiagnostics
+        
     WellPathImportSsihub
 
 	ResultStatisticsCache

ApplicationCode/UnitTests/CMakeLists_files.cmake

@@ -23,6 +23,7 @@ ${CEE_CURRENT_LIST_DIR}ScalarMapper-Test.cpp
 ${CEE_CURRENT_LIST_DIR}WellPathAsciiFileReader-Test.cpp
 ${CEE_CURRENT_LIST_DIR}opm-parser-Test.cpp
 ${CEE_CURRENT_LIST_DIR}opm-parser-Performance-Test.cpp
+${CEE_CURRENT_LIST_DIR}opm-flowdiagnostics-Test.cpp
 )
 
 list(APPEND CODE_HEADER_FILES

ApplicationCode/UnitTests/RifReaderEclipseOutput-Test.cpp

@@ -40,6 +40,7 @@ TEST(RigReservoirTest, BasicTest)
 {
     RifEclipseUnifiedRestartFileAccess unrstAccess;
 
+/*
     QStringList filenames;
     //filenames << "d:/Models/Statoil/testcase_juli_2011/data/TEST10K_FLT_LGR_NNC.UNRST";
     filenames << "d:/Models/MRST/simple/SIMPLE.UNRST";
@@ -62,6 +63,7 @@ TEST(RigReservoirTest, BasicTest)
     {
         qDebug() << reportNum;
     }
+*/
 
 /*
     cvf::ref<RifReaderEclipseOutput> readerInterfaceEcl = new RifReaderEclipseOutput;

ApplicationCode/UnitTests/opm-flowdiagnostics-Test.cpp

@@ -0,0 +1,10 @@
+#include "gtest/gtest.h"
+
+#include <opm/flowdiagnostics/CellSet.hpp>
+#include <opm/utility/graph/AssembledConnections.hpp>
+
+TEST(opm_flowdiagnostics_test, basic_construction)
+{
+    auto g = Opm::AssembledConnections{};
+    auto s = Opm::FlowDiagnostics::CellSet{};
+}

CMakeLists.txt

@@ -210,6 +210,7 @@ endif(RESINSIGHT_ERT_EXTERNAL_LIB_ROOT OR RESINSIGHT_ERT_EXTERNAL_INCLUDE_ROOT)
 
 add_subdirectory(ThirdParty/custom-opm-parser)
 add_subdirectory(ThirdParty/custom-opm-parser/custom-opm-parser-tests)
+add_subdirectory(ThirdParty/custom-opm-flowdiagnostics)
 
 
 ################################################################################

Fwk/AppFwk/cafViewer/cafViewer.cpp

@@ -227,8 +227,8 @@ void caf::Viewer::setupRenderingSequence()
         quadRenderGen.addFragmentShaderCode(cvf::ShaderSourceProvider::instance()->getSourceFromRepository(cvf::ShaderSourceRepository::fs_Unlit));
         quadRenderGen.addFragmentShaderCode(cvf::ShaderSourceProvider::instance()->getSourceFromRepository(cvf::ShaderSourceRepository::src_Texture));
 
-        cvf::ref<cvf::Rendering> quadRendering = quadRenderGen.generate();
-        m_renderingSequence->addRendering(quadRendering.p());
+        m_quadRendering = quadRenderGen.generate();
+        m_renderingSequence->addRendering(m_quadRendering.p());
     }
 
     updateCamera(width(), height());
@@ -509,10 +509,7 @@ void caf::Viewer::resizeGL(int width, int height)
     {
         m_offscreenFbo->resizeAttachedBuffers(width, height);
     
-        CVF_ASSERT(m_renderingSequence->renderingCount() > 1);
-
-        cvf::ref<cvf::Rendering> quadRendering = m_renderingSequence->rendering(1);
-        quadRendering->camera()->viewport()->set(0, 0, width, height);
+        m_quadRendering->camera()->viewport()->set(0, 0, width, height);
     }
 
     updateCamera(width, height);

Fwk/AppFwk/cafViewer/cafViewer.h

@@ -231,6 +231,7 @@ private:
     // Offscreen render objects
     cvf::ref<cvf::FramebufferObject>    m_offscreenFbo;
     cvf::ref<cvf::Texture>              m_offscreenTexture;
+    cvf::ref<cvf::Rendering>            m_quadRendering;
 };
 
 } // End namespace caf

Fwk/AppFwk/cafVizExtensions/cafCategoryLegend.cpp

@@ -233,7 +233,7 @@ void CategoryLegend::setupTextDrawer(TextDrawer* textDrawer, OverlayColorLegendL
     float lastVisibleTextY = 0.0;
 
     CVF_ASSERT(m_categoryMapper.notNull());
-    size_t numLabels = m_categoryMapper->categories().size();
+    size_t numLabels = m_categoryMapper->categoryCount();
 
     float categoryHeight = static_cast<float>(layout->legendRect.height() / numLabels);
 
@@ -261,11 +261,10 @@ void CategoryLegend::setupTextDrawer(TextDrawer* textDrawer, OverlayColorLegendL
             }
         }
 
-        double tickValue = m_categoryMapper->categories()[it];
-        String valueString = String::number(tickValue);
+        String displayText = m_categoryMapper->textForCategoryIndex(it);
 
         Vec2f pos(textX, textY);
-        textDrawer->addText(valueString, pos);
+        textDrawer->addText(displayText, pos);
 
         lastVisibleTextY = textY;
         m_visibleCategoryLabels.push_back(true);

Fwk/AppFwk/cafVizExtensions/cafCategoryMapper.cpp

@@ -22,9 +22,22 @@ CategoryMapper::CategoryMapper()
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-void CategoryMapper::setCategories(const IntArray& categories)
+void CategoryMapper::setCategories(const IntArray& categoryValues)
 {
-    m_categories = categories;
+    m_categoryValues = categoryValues;
+
+    ref<Color3ubArray> colorArr = ScalarMapper::colorTableArray(ColorTable::NORMAL);
+
+    setColors(*(colorArr.p()));
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+void CategoryMapper::setCategories(const cvf::IntArray& categoryValues, const std::vector<cvf::String>& categoryNames)
+{
+    m_categoryValues = categoryValues;
+    m_categoryNames = categoryNames;
 
     ref<Color3ubArray> colorArr = ScalarMapper::colorTableArray(ColorTable::NORMAL);
 
@@ -36,9 +49,9 @@ void CategoryMapper::setCategories(const IntArray& categories)
 //--------------------------------------------------------------------------------------------------
 void CategoryMapper::setColors(const Color3ubArray& colorArray)
 {
-    m_colors.resize(m_categories.size());
+    m_colors.resize(m_categoryValues.size());
 
-    for (size_t i = 0; i < m_categories.size(); i++)
+    for (size_t i = 0; i < m_categoryValues.size(); i++)
     {
         size_t colIdx = i % colorArray.size();
         m_colors[i] = colorArray[colIdx];
@@ -48,9 +61,41 @@ void CategoryMapper::setColors(const Color3ubArray& colorArray)
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-cvf::IntArray CategoryMapper::categories() const
+void CategoryMapper::setInterpolateColors(const cvf::Color3ubArray& colorArray)
 {
-    return m_categories;
+    if (m_categoryValues.size() > 0)
+    {
+        m_colors = *interpolateColorArray(colorArray, static_cast<cvf::uint>(m_categoryValues.size()));
+    }
+    else
+    {
+        m_colors.clear();
+    }
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+size_t CategoryMapper::categoryCount() const
+{
+    return m_categoryValues.size();
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+const cvf::String CategoryMapper::textForCategoryIndex(size_t index) const
+{
+    CVF_ASSERT(index < m_categoryValues.size());
+    if (index < m_categoryNames.size())
+    {
+        return m_categoryNames[index];
+    }
+    else
+    {
+        double tickValue = m_categoryValues[index];
+        return String::number(tickValue);
+    }
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -101,9 +146,9 @@ double CategoryMapper::normalizedValue(double categoryValue) const
 
     if (catIndex != -1)
     {
-        double halfLevelHeight = 1.0 / (m_categories.size() * 2);
+        double halfLevelHeight = 1.0 / (m_categoryValues.size() * 2);
 
-        double normVal = static_cast<double>(catIndex) / static_cast<double>(m_categories.size());
+        double normVal = static_cast<double>(catIndex) / static_cast<double>(m_categoryValues.size());
 
         return normVal + halfLevelHeight;
     }
@@ -120,13 +165,13 @@ double CategoryMapper::domainValue(double normalizedValue) const
 {
     double clampedValue = cvf::Math::clamp(normalizedValue, 0.0, 1.0);
 
-    if (m_categories.size() == 0)
+    if (m_categoryValues.size() == 0)
     {
         return 0.0;
     }
 
-    size_t catIndex = static_cast<size_t>(clampedValue * m_categories.size());
-    return m_categories[catIndex];
+    size_t catIndex = static_cast<size_t>(clampedValue * m_categoryValues.size());
+    return m_categoryValues[catIndex];
 }
 
 //--------------------------------------------------------------------------------------------------
@@ -137,9 +182,9 @@ int CategoryMapper::categoryIndexForCategory(double domainValue) const
     int catIndex = -1;
 
     size_t i = 0;
-    while (i < m_categories.size() && catIndex == -1)
+    while (i < m_categoryValues.size() && catIndex == -1)
     {
-        if (m_categories[i] == domainValue)
+        if (m_categoryValues[i] == domainValue)
         {
             catIndex = static_cast<int>(i);
         }

Fwk/AppFwk/cafVizExtensions/cafCategoryMapper.h

@@ -4,6 +4,7 @@
 #include "cvfBase.h"
 #include "cvfObject.h"
 #include "cvfScalarMapper.h"
+#include "cvfString.h"
 
 namespace caf {
 
@@ -17,10 +18,14 @@ class CategoryMapper : public cvf::ScalarMapper
 public:
     CategoryMapper();
 
-    void                setCategories(const cvf::IntArray& categories);
+    void                setCategories(const cvf::IntArray& categoryValues);
+    void                setCategories(const cvf::IntArray& categoryValues, const std::vector<cvf::String>& categoryNames);
+    
     void                setColors(const cvf::Color3ubArray& colorArray);
+    void                setInterpolateColors(const cvf::Color3ubArray& colorArray);
 
-    cvf::IntArray       categories() const;
+    size_t              categoryCount() const;
+    const cvf::String   textForCategoryIndex(size_t index) const;
 
     virtual cvf::Vec2f  mapToTextureCoord(double scalarValue) const;
     virtual bool        updateTexture(cvf::TextureImage* image) const;
@@ -38,7 +43,8 @@ private:
     cvf::Color3ubArray   m_colors;
     cvf::uint            m_textureSize;      // The size of texture that updateTexture() is will produce. 
 
-    cvf::IntArray        m_categories;
+    cvf::IntArray            m_categoryValues;
+    std::vector<cvf::String> m_categoryNames;
 };
 
 }

ThirdParty/custom-opm-flowdiagnostics/CMakeLists.txt

@@ -0,0 +1,23 @@
+cmake_minimum_required (VERSION 2.8)
+
+project (custom-opm-flowdiagnostics)
+
+include_directories(
+    opm-flowdiagnostics
+)
+
+include (opm-flowdiagnostics/CMakeLists_files.cmake)
+
+set(project_source_files
+    ${MAIN_SOURCE_FILES} 
+    ${PUBLIC_HEADER_FILES}
+)
+
+foreach (file ${project_source_files} )
+    list(APPEND project_source_files_complete_path "opm-flowdiagnostics/${file}" )
+endforeach()
+
+add_library(${PROJECT_NAME} 
+    ${project_source_files_complete_path} 
+)
+

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/.gitignore

@@ -0,0 +1,38 @@
+# Editor autosave files
+*~
+*.swp
+
+# Compiled Object files
+*.slo
+*.lo
+*.o
+*.obj
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Compiled Dynamic libraries
+*.so
+*.dylib
+*.dll
+
+# Fortran module files
+*.mod
+
+# Compiled Static libraries
+*.lai
+*.la
+*.a
+*.lib
+
+# Executables
+*.exe
+*.out
+*.app
+
+# Mac OS X debug info
+*.dSYM
+
+# emacs directory setting:
+.dir-locals.el

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/.travis.yml

@@ -0,0 +1,24 @@
+language: cpp
+
+compiler:
+  - gcc
+
+addons:
+  apt:
+    sources:
+      - boost-latest
+      - ubuntu-toolchain-r-test
+    packages:
+      - libboost1.55-all-dev
+      - gcc-4.8
+      - g++-4.8
+
+install:
+    - export CXX="g++-4.8" CC="gcc-4.8"
+    - cd ..
+
+    - git clone https://github.com/OPM/opm-common.git
+
+    - opm-common/travis/build-opm-common.sh
+
+script: opm-flowdiagnostics/travis/build-and-test-opm-flowdiagnostics.sh

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/CMakeLists.txt

@@ -0,0 +1,93 @@
+# -*- mode: cmake; tab-width: 2; indent-tabs-mode: t; truncate-lines: t; compile-command: "cmake -Wdev" -*-
+# vim: set filetype=cmake autoindent tabstop=2 shiftwidth=2 noexpandtab softtabstop=2 nowrap:
+
+###########################################################################
+#                                                                         #
+# Note: The bulk of the build system is located in the cmake/ directory.  #
+#       This file only contains the specializations for this particular   #
+#       project. Most likely you are interested in editing one of these   #
+#       files instead:                                                    #
+#                                                                         #
+#       dune.module                              Name and version number  #
+#       CMakeLists_files.cmake                   Path of source files     #
+#       cmake/Modules/${project}-prereqs.cmake   Dependencies             #
+#                                                                         #
+###########################################################################
+
+cmake_minimum_required (VERSION 2.8)
+
+# additional search modules
+set( OPM_COMMON_ROOT "" CACHE PATH "Root directory containing OPM related cmake modules")
+option(SIBLING_SEARCH "Search for other modules in sibling directories?" ON)
+
+if(NOT OPM_COMMON_ROOT)
+   find_package(opm-common QUIET)
+endif()
+
+if (opm-common_FOUND)
+   include(OpmInit)
+else()
+   unset(opm-common_FOUND)
+
+   if (NOT OPM_COMMON_ROOT AND SIBLING_SEARCH)
+      set(OPM_COMMON_ROOT ${PROJECT_SOURCE_DIR}/../opm-common)
+   endif()
+   if (OPM_COMMON_ROOT)
+      list( APPEND CMAKE_MODULE_PATH "${OPM_COMMON_ROOT}/cmake/Modules")
+      include (OpmInit OPTIONAL RESULT_VARIABLE OPM_INIT)
+      set( OPM_MACROS_ROOT ${OPM_COMMON_ROOT} )
+   endif()
+
+   if (NOT OPM_INIT)
+      message( "" )
+      message( " /---------------------------------------------------------------------------------\\")
+      message( " |  Could not locate the opm build macros. The opm build macros                    |")
+      message( " |  are in a separate repository - instructions to proceed:                        |")
+      message( " |                                                                                 |")
+      message( " |    1. Clone the repository: git clone git@github.com:OPM/opm-common.git         |")
+      message( " |                                                                                 |")
+      message( " |    2. Run cmake in the current project with -DOPM_COMMON_ROOT=<path>/opm-common |")
+      message( " |                                                                                 |")
+      message( " \\---------------------------------------------------------------------------------/")
+      message( "" )
+      message( FATAL_ERROR "Could not find OPM Macros")
+   endif()
+
+endif()
+
+# not the same location as most of the other projects; this hook overrides
+macro (dir_hook)
+endmacro (dir_hook)
+
+# list of prerequisites for this particular project; this is in a
+# separate file (in cmake/Modules sub-directory) because it is shared
+# with the find module
+include ( ${project}-prereqs )
+
+# read the list of components from this file (in the project directory);
+# it should set various lists with the names of the files to include
+include (CMakeLists_files.cmake)
+
+macro (config_hook)
+endmacro (config_hook)
+
+macro (prereqs_hook)
+endmacro (prereqs_hook)
+
+macro (sources_hook)
+endmacro (sources_hook)
+
+macro (fortran_hook)
+endmacro (fortran_hook)
+
+macro (files_hook)
+endmacro (files_hook)
+
+macro (tests_hook)
+endmacro (tests_hook)
+
+macro (install_hook)
+endmacro (install_hook)
+
+# all setup common to the OPM library modules is done here
+include (OpmLibMain)

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/CMakeLists_files.cmake

@@ -0,0 +1,60 @@
+# -*- mode: cmake; tab-width: 2; indent-tabs-mode: t; truncate-lines: t; compile-command: "cmake -Wdev" -*-
+# vim: set filetype=cmake autoindent tabstop=2 shiftwidth=2 noexpandtab softtabstop=2 nowrap:
+
+# This file sets up five lists:
+#       MAIN_SOURCE_FILES     List of compilation units which will be included in
+#                             the library. If it isn't on this list, it won't be
+#                             part of the library. Please try to keep it sorted to
+#                             maintain sanity.
+#
+#       TEST_SOURCE_FILES     List of programs that will be run as unit tests.
+#
+#       TEST_DATA_FILES       Files from the source three that should be made
+#                             available in the corresponding location in the build
+#                             tree in order to run tests there.
+#
+#       EXAMPLE_SOURCE_FILES  Other programs that will be compiled as part of the
+#                             build, but which is not part of the library nor is
+#                             run as tests.
+#
+#       PUBLIC_HEADER_FILES   List of public header files that should be
+#                             distributed together with the library. The source
+#                             files can of course include other files than these;
+#                             you should only add to this list if the *user* of
+#                             the library needs it.
+
+list (APPEND MAIN_SOURCE_FILES
+        opm/flowdiagnostics/CellSet.cpp
+        opm/flowdiagnostics/CellSetValues.cpp
+        opm/flowdiagnostics/ConnectionValues.cpp
+        opm/flowdiagnostics/ConnectivityGraph.cpp
+        opm/flowdiagnostics/Solution.cpp
+        opm/flowdiagnostics/Toolbox.cpp
+        opm/flowdiagnostics/TracerTofSolver.cpp
+        opm/utility/graph/tarjan.c
+        opm/utility/graph/AssembledConnections.cpp
+        opm/utility/numeric/RandomVector.cpp
+        )
+
+list (APPEND TEST_SOURCE_FILES
+        tests/test_assembledconnections.cpp
+        tests/test_cellset.cpp
+        tests/test_cellsetvalues.cpp
+        tests/test_connectionvalues.cpp
+        tests/test_connectivitygraph.cpp
+        tests/test_flowdiagnosticstool.cpp
+        tests/test_tarjan.cpp
+        )
+
+list (APPEND PUBLIC_HEADER_FILES
+        opm/flowdiagnostics/CellSet.hpp
+        opm/flowdiagnostics/CellSetValues.hpp
+        opm/flowdiagnostics/ConnectionValues.hpp
+        opm/flowdiagnostics/ConnectivityGraph.hpp
+        opm/flowdiagnostics/Solution.hpp
+        opm/flowdiagnostics/Toolbox.hpp
+        opm/flowdiagnostics/TracerTofSolver.hpp
+        opm/utility/graph/AssembledConnections.hpp
+        opm/utility/graph/AssembledConnectionsIteration.hpp
+        opm/utility/numeric/RandomVector.hpp
+        )

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/LICENSE

@@ -0,0 +1,674 @@
+                    GNU GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+                            Preamble
+
+  The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+  The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works.  By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users.  We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors.  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+  To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights.  Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received.  You must make sure that they, too, receive
+or can get the source code.  And you must show them these terms so they
+know their rights.
+
+  Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+  For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software.  For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+  Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so.  This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software.  The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable.  Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products.  If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+  Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary.  To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+                       TERMS AND CONDITIONS
+
+  0. Definitions.
+
+  "This License" refers to version 3 of the GNU General Public License.
+
+  "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+  "The Program" refers to any copyrightable work licensed under this
+License.  Each licensee is addressed as "you".  "Licensees" and
+"recipients" may be individuals or organizations.
+
+  To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy.  The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+  A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+  To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy.  Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+  To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies.  Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+  An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License.  If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+  1. Source Code.
+
+  The "source code" for a work means the preferred form of the work
+for making modifications to it.  "Object code" means any non-source
+form of a work.
+
+  A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+  The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form.  A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+  The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities.  However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work.  For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+  The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+  The Corresponding Source for a work in source code form is that
+same work.
+
+  2. Basic Permissions.
+
+  All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met.  This License explicitly affirms your unlimited
+permission to run the unmodified Program.  The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work.  This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+  You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force.  You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright.  Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+  Conveying under any other circumstances is permitted solely under
+the conditions stated below.  Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+  3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+  No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+  When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+  4. Conveying Verbatim Copies.
+
+  You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+  You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+  5. Conveying Modified Source Versions.
+
+  You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+    a) The work must carry prominent notices stating that you modified
+    it, and giving a relevant date.
+
+    b) The work must carry prominent notices stating that it is
+    released under this License and any conditions added under section
+    7.  This requirement modifies the requirement in section 4 to
+    "keep intact all notices".
+
+    c) You must license the entire work, as a whole, under this
+    License to anyone who comes into possession of a copy.  This
+    License will therefore apply, along with any applicable section 7
+    additional terms, to the whole of the work, and all its parts,
+    regardless of how they are packaged.  This License gives no
+    permission to license the work in any other way, but it does not
+    invalidate such permission if you have separately received it.
+
+    d) If the work has interactive user interfaces, each must display
+    Appropriate Legal Notices; however, if the Program has interactive
+    interfaces that do not display Appropriate Legal Notices, your
+    work need not make them do so.
+
+  A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit.  Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+  6. Conveying Non-Source Forms.
+
+  You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+    a) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by the
+    Corresponding Source fixed on a durable physical medium
+    customarily used for software interchange.
+
+    b) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by a
+    written offer, valid for at least three years and valid for as
+    long as you offer spare parts or customer support for that product
+    model, to give anyone who possesses the object code either (1) a
+    copy of the Corresponding Source for all the software in the
+    product that is covered by this License, on a durable physical
+    medium customarily used for software interchange, for a price no
+    more than your reasonable cost of physically performing this
+    conveying of source, or (2) access to copy the
+    Corresponding Source from a network server at no charge.
+
+    c) Convey individual copies of the object code with a copy of the
+    written offer to provide the Corresponding Source.  This
+    alternative is allowed only occasionally and noncommercially, and
+    only if you received the object code with such an offer, in accord
+    with subsection 6b.
+
+    d) Convey the object code by offering access from a designated
+    place (gratis or for a charge), and offer equivalent access to the
+    Corresponding Source in the same way through the same place at no
+    further charge.  You need not require recipients to copy the
+    Corresponding Source along with the object code.  If the place to
+    copy the object code is a network server, the Corresponding Source
+    may be on a different server (operated by you or a third party)
+    that supports equivalent copying facilities, provided you maintain
+    clear directions next to the object code saying where to find the
+    Corresponding Source.  Regardless of what server hosts the
+    Corresponding Source, you remain obligated to ensure that it is
+    available for as long as needed to satisfy these requirements.
+
+    e) Convey the object code using peer-to-peer transmission, provided
+    you inform other peers where the object code and Corresponding
+    Source of the work are being offered to the general public at no
+    charge under subsection 6d.
+
+  A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+  A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling.  In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage.  For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product.  A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+  "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source.  The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+  If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information.  But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+  The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed.  Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+  Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+  7. Additional Terms.
+
+  "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law.  If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+  When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it.  (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.)  You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+  Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+    a) Disclaiming warranty or limiting liability differently from the
+    terms of sections 15 and 16 of this License; or
+
+    b) Requiring preservation of specified reasonable legal notices or
+    author attributions in that material or in the Appropriate Legal
+    Notices displayed by works containing it; or
+
+    c) Prohibiting misrepresentation of the origin of that material, or
+    requiring that modified versions of such material be marked in
+    reasonable ways as different from the original version; or
+
+    d) Limiting the use for publicity purposes of names of licensors or
+    authors of the material; or
+
+    e) Declining to grant rights under trademark law for use of some
+    trade names, trademarks, or service marks; or
+
+    f) Requiring indemnification of licensors and authors of that
+    material by anyone who conveys the material (or modified versions of
+    it) with contractual assumptions of liability to the recipient, for
+    any liability that these contractual assumptions directly impose on
+    those licensors and authors.
+
+  All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10.  If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term.  If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+  If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+  Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+  8. Termination.
+
+  You may not propagate or modify a covered work except as expressly
+provided under this License.  Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+  However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+  Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+  Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License.  If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+  9. Acceptance Not Required for Having Copies.
+
+  You are not required to accept this License in order to receive or
+run a copy of the Program.  Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance.  However,
+nothing other than this License grants you permission to propagate or
+modify any covered work.  These actions infringe copyright if you do
+not accept this License.  Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+  10. Automatic Licensing of Downstream Recipients.
+
+  Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License.  You are not responsible
+for enforcing compliance by third parties with this License.
+
+  An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations.  If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+  You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License.  For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+  11. Patents.
+
+  A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based.  The
+work thus licensed is called the contributor's "contributor version".
+
+  A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version.  For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+  Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+  In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement).  To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+  If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients.  "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+  If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+  A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License.  You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+  Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+  12. No Surrender of Others' Freedom.
+
+  If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all.  For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+  13. Use with the GNU Affero General Public License.
+
+  Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work.  The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+  14. Revised Versions of this License.
+
+  The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+  Each version is given a distinguishing version number.  If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation.  If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+  If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+  Later license versions may give you additional or different
+permissions.  However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+  15. Disclaimer of Warranty.
+
+  THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+  16. Limitation of Liability.
+
+  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+  17. Interpretation of Sections 15 and 16.
+
+  If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+                     END OF TERMS AND CONDITIONS
+
+            How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    {one line to give the program's name and a brief idea of what it does.}
+    Copyright (C) {year}  {name of author}
+
+    This program 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.
+
+    This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+
+Also add information on how to contact you by electronic and paper mail.
+
+  If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+    {project}  Copyright (C) {year}  {fullname}
+    This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+  You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+<http://www.gnu.org/licenses/>.
+
+  The GNU General Public License does not permit incorporating your program
+into proprietary programs.  If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library.  If this is what you want to do, use the GNU Lesser General
+Public License instead of this License.  But first, please read
+<http://www.gnu.org/philosophy/why-not-lgpl.html>.

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/README.md

@@ -0,0 +1,2 @@
+# opm-flowdiagnostics
+Computational Kernels for Inexpensive Flow Diagnostics Calculations

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/dune.module

@@ -0,0 +1,13 @@
+####################################################################
+# Dune module information file: This file gets parsed by dunecontrol
+# and by the CMake build scripts.
+####################################################################
+
+Module: opm-flowdiagnostics
+Description: Open Porous Media Initiative flow diagnostics
+Version: 2016.10-pre
+Label: 2016.10-pre
+Maintainer: opm@opm-project.org
+MaintainerName: OPM community
+Url: http://opm-project.org
+Depends: opm-common

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/CellSet.cpp

@@ -0,0 +1,82 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#include <opm/flowdiagnostics/CellSet.hpp>
+
+#include <utility>
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+
+
+CellSetID::CellSetID()
+{}
+
+CellSetID::CellSetID(Repr id)
+    : id_(std::move(id))
+{
+}
+
+std::string
+CellSetID::to_string() const
+{
+    return id_;
+}
+
+// =====================================================================
+
+void
+CellSet::identify(CellSetID id)
+{
+    id_ = std::move(id);
+}
+
+const CellSetID&
+CellSet::id() const
+{
+    return id_;
+}
+
+void
+CellSet::insert(const int i)
+{
+    iset_.insert(i);
+}
+
+CellSet::const_iterator
+CellSet::begin() const
+{
+    return iset_.begin();
+}
+
+CellSet::const_iterator
+CellSet::end() const
+{
+    return iset_.end();
+}
+
+} // namespace FlowDiagnostics
+} // namespace Opm

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/CellSet.hpp

@@ -0,0 +1,74 @@
+/*
+  Copyright 2016 Statoil ASA.
+  Copyright 2016 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/>.
+*/
+
+#ifndef OPM_CELLSET_HEADER_INCLUDED
+#define OPM_CELLSET_HEADER_INCLUDED
+
+#include <string>
+#include <unordered_set>
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+
+
+    class CellSetID
+    {
+    public:
+        using Repr = std::string;
+
+        CellSetID();
+
+        explicit CellSetID(Repr id);
+
+        std::string to_string() const;
+
+    private:
+        Repr id_;
+    };
+
+    class CellSet
+    {
+    private:
+        using IndexSet = std::unordered_set<int>;
+
+    public:
+        using const_iterator = IndexSet::const_iterator;
+
+        void identify(CellSetID id);
+
+        const CellSetID& id() const;
+
+        void insert(const int cell);
+
+        const_iterator begin() const;
+        const_iterator end()   const;
+
+    private:
+        CellSetID id_;
+        IndexSet  iset_;
+    };
+
+
+} // namespace FlowDiagnostics
+} // namespace Opm
+
+#endif  // OPM_CELLSET_HEADER_INCLUDED

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/CellSetValues.cpp

@@ -0,0 +1,60 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#include <opm/flowdiagnostics/CellSetValues.hpp>
+
+#include <cassert>
+#include <utility>
+
+namespace Opm {
+namespace FlowDiagnostics {
+
+CellSetValues::CellSetValues(const SizeType initialCapacity)
+{
+    assoc_.reserve(initialCapacity);
+}
+
+void
+CellSetValues::addCellValue(const int    cellIndex,
+                                 const double cellValue)
+{
+    assoc_.push_back(Association{cellIndex, cellValue});
+}
+
+CellSetValues::SizeType
+CellSetValues::cellValueCount() const
+{
+    return assoc_.size();
+}
+
+CellSetValues::CellValue
+CellSetValues::cellValue(const SizeType cellValueIndex) const
+{
+    const auto& a = assoc_[cellValueIndex];
+
+    return { a.index, a.value };
+}
+
+} // namespace FlowDiagnostics
+} // namespace Opm

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/CellSetValues.hpp

@@ -0,0 +1,76 @@
+/*
+  Copyright 2016 Statoil ASA.
+  Copyright 2016 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/>.
+*/
+
+#ifndef OPM_CELLSETVALUES_HEADER_INCLUDED
+#define OPM_CELLSETVALUES_HEADER_INCLUDED
+
+#include <utility>
+#include <vector>
+
+namespace Opm {
+namespace FlowDiagnostics {
+
+    class CellSetValues
+    {
+    public:
+        using SizeType  = std::vector<int>::size_type;
+        using CellValue = std::pair<int, double>;
+
+        /// Constructor.
+        ///
+        /// @param[in] initialCapacity Number of elements that can be stored
+        ///            in set without reallocation.
+        explicit CellSetValues(const SizeType initialCapacity = 0);
+
+        /// Associate value with specific cell, represented by its index.
+        ///
+        /// @param[in] cellIndex Index of specific cell.
+        ///
+        /// @param[in] cellValue Value associated with cell @p cellIndex.
+        void addCellValue(const int    cellIndex,
+                          const double cellValue);
+
+        /// Retrieve number of elements stored in set.
+        SizeType cellValueCount() const;
+
+        /// Retrieve value association for single set element.
+        ///
+        /// @param[in] cellValueIndex Linear ID of single cell->value
+        ///            association.  Must be in the range @code [0,
+        ///            cellValueCount()) @endcode.
+        ///
+        /// @returns Single association between cell index and numerical
+        /// value.
+        CellValue cellValue(const SizeType cellValueIndex) const;
+
+    private:
+        struct Association
+        {
+            int    index;
+            double value;
+        };
+
+        std::vector<Association> assoc_;
+    };
+
+} // namespace FlowDiagnostics
+} // namespace Opm
+
+#endif // OPM_CELLSETVALUES_HEADER_INCLUDED

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/ConnectionValues.cpp

@@ -0,0 +1,83 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#include <opm/flowdiagnostics/ConnectionValues.hpp>
+
+#include <exception>
+#include <stdexcept>
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+
+
+ConnectionValues::
+ConnectionValues(const NumConnections& nconn,
+                 const NumPhases&      nphase)
+    : nconn_ (nconn)
+    , nphase_(nphase)
+    , data_  (nconn_.total * nphase_.total, 0.0)
+{}
+
+std::vector<double>::size_type
+ConnectionValues::numConnections() const
+{
+    return nconn_.total;
+}
+
+std::size_t
+ConnectionValues::numPhases() const
+{
+    return nphase_.total;
+}
+
+double
+ConnectionValues::operator()(const ConnID&  conn,
+                                  const PhaseID& phase) const
+{
+    if ((conn .id >= nconn_ .total) ||
+        (phase.id >= nphase_.total))
+    {
+        throw std::logic_error("(Connection,Phase) pair out of range");
+    }
+
+    return data_[conn.id*nphase_.total + phase.id];
+}
+
+double&
+ConnectionValues::operator()(const ConnID&  conn,
+                                  const PhaseID& phase)
+{
+    if ((conn .id >= nconn_ .total) ||
+        (phase.id >= nphase_.total))
+    {
+        throw std::logic_error("(Connection,Phase) pair out of range");
+    }
+
+    return data_[conn.id*nphase_.total + phase.id];
+}
+
+} // namespace FlowDiagnostics
+} // namespace Opm

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/ConnectionValues.hpp

@@ -0,0 +1,78 @@
+/*
+  Copyright 2016 Statoil ASA.
+  Copyright 2016 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/>.
+*/
+
+#ifndef OPM_CONNECTIONVALUES_HEADER_INCLUDED
+#define OPM_CONNECTIONVALUES_HEADER_INCLUDED
+
+#include <cstddef>
+#include <vector>
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+
+
+    class ConnectionValues
+    {
+    public:
+        struct NumConnections
+        {
+            std::vector<double>::size_type total;
+        };
+
+        struct NumPhases
+        {
+            std::size_t total;
+        };
+
+        struct ConnID
+        {
+            std::vector<double>::size_type id;
+        };
+
+        struct PhaseID
+        {
+            std::size_t id;
+        };
+
+        ConnectionValues(const NumConnections& nconn,
+                         const NumPhases&      nphase);
+
+        std::vector<double>::size_type numConnections() const;
+        std::size_t                    numPhases()      const;
+
+        double operator()(const ConnID&  connection,
+                          const PhaseID& phase) const;
+
+        double& operator()(const ConnID&  connection,
+                           const PhaseID& phase);
+
+    private:
+        NumConnections      nconn_;
+        NumPhases           nphase_;
+
+        std::vector<double> data_;
+    };
+
+} // namespace FlowDiagnostics
+} // namespace Opm
+
+#endif // OPM_CONNECTIONVALUES_HEADER_INCLUDED

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/ConnectivityGraph.cpp

@@ -0,0 +1,80 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+    
+#include <opm/flowdiagnostics/ConnectivityGraph.hpp>
+
+#include <algorithm>
+#include <cassert>
+#include <exception>
+#include <stdexcept>
+#include <utility>
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+
+ConnectivityGraph::ConnectivityGraph(const int        num_cells,
+                                     std::vector<int> connection_to_cell)
+    : numCells_ (num_cells)
+    , connCells_(std::move(connection_to_cell))
+{
+    if ((connCells_.size() % 2) != 0) {
+        throw std::logic_error("Neighbourship must be N-by-2");
+    }
+
+    if ((! connCells_.empty()) &&
+        (*std::max_element(connCells_.begin(),
+                           connCells_.end())
+         >= num_cells))
+    {
+        throw std::logic_error("Cell indices must be in range");
+    }
+}
+
+ConnectivityGraph::SizeType
+ConnectivityGraph::numCells() const
+{
+    return numCells_;
+}
+
+ConnectivityGraph::SizeType
+ConnectivityGraph::numConnections() const
+{
+    return connCells_.size() / 2;
+}
+
+ConnectivityGraph::CellPair
+ConnectivityGraph::connection(const SizeType i) const
+{
+    assert (i < numConnections());
+
+    const auto* const start = &connCells_[2*i + 0];
+
+    return { *(start + 0), *(start + 1) };
+}
+
+} // namespace FlowDiagnostics
+} // namespace Opm
+

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/ConnectivityGraph.hpp

@@ -0,0 +1,63 @@
+/*
+  Copyright 2016 Statoil ASA.
+  Copyright 2016 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/>.
+*/
+
+#ifndef OPM_CONNECTIVITYGRAPH_HEADER_INCLUDED
+#define OPM_CONNECTIVITYGRAPH_HEADER_INCLUDED
+
+#include <vector>
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+    class ConnectivityGraph
+    {
+    public:
+        /// Construct from explicit neighbourship table.
+        ///
+        /// The @p connection_to_cell must have size equal to @code 2 * E
+        /// @endcode in which @c E is the number of connections.  Connection
+        /// @c k connects cells @code connection_to_cell[2*k + 0] @endcode
+        /// and @code connection_to_cell[2*k + 1] @endcode.
+        ConnectivityGraph(const int        num_cells,
+                          std::vector<int> connection_to_cell);
+
+        struct CellPair
+        {
+            int first;
+            int second;
+        };
+
+        using SizeType = std::vector<int>::size_type;
+
+        SizeType numCells()      const;
+        SizeType numConnections() const;
+
+        CellPair connection(const SizeType i) const;
+
+    private:
+        SizeType         numCells_;
+        std::vector<int> connCells_;
+    };
+
+} // namespace FlowDiagnostics
+} // namespace Opm
+
+#endif // OPM_CONNECTIVITYGRAPH_HEADER_INCLUDED

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/Solution.cpp

@@ -0,0 +1,231 @@
+/*
+  Copyright 2016 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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#include <opm/flowdiagnostics/Solution.hpp>
+#include <map>
+#include <utility>
+
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+
+
+// ---------------------------------------------------------------------
+// Class Solution::Impl
+// ---------------------------------------------------------------------
+
+class Solution::Impl
+{
+public:
+    using GlobalToF = std::vector<double>;
+
+    void assignToF(GlobalToF&& tof);
+    void assign   (const CellSetID& i, TimeOfFlight&&  tof);
+    void assign   (const CellSetID& i, TracerConcentration&& conc);
+
+    std::vector<CellSetID> startPoints() const;
+
+    const GlobalToF& timeOfFlight() const;
+
+    CellSetValues timeOfFlight (const CellSetID& tracer) const;
+    CellSetValues concentration(const CellSetID& tracer) const;
+
+private:
+    struct CompareCellSetIDs
+    {
+        bool operator()(const CellSetID& x,
+                        const CellSetID& y) const
+        {
+            return x.to_string() < y.to_string();
+        }
+    };
+
+    using SolutionMap =
+        std::map<CellSetID, CellSetValues, CompareCellSetIDs>;
+
+    GlobalToF   tof_;
+    SolutionMap tracerToF_;
+    SolutionMap tracer_;
+
+    void assign(const CellSetID& i,
+                CellSetValues&&  x,
+                SolutionMap&     soln);
+
+    CellSetValues
+    solutionValues(const CellSetID&   i,
+                   const SolutionMap& soln) const;
+};
+
+void
+Solution::Impl::assignToF(GlobalToF&& tof)
+{
+    tof_ = std::move(tof);
+}
+
+void
+Solution::
+Impl::assign(const CellSetID& i, TimeOfFlight&& tof)
+{
+    assign(i, std::move(tof.data), tracerToF_);
+}
+
+void
+Solution::
+Impl::assign(const CellSetID& i, TracerConcentration&& conc)
+{
+    assign(i, std::move(conc.data), tracer_);
+}
+
+std::vector<CellSetID>
+Solution::Impl::startPoints() const
+{
+    auto s = std::vector<CellSetID>{};
+    s.reserve(tracer_.size());
+
+    for (const auto& t : tracer_) {
+        s.emplace_back(t.first);
+    }
+
+    return s;
+}
+
+const Solution::Impl::GlobalToF&
+Solution::Impl::timeOfFlight() const
+{
+    return tof_;
+}
+
+CellSetValues
+Solution::
+Impl::timeOfFlight(const CellSetID& tracer) const
+{
+    return solutionValues(tracer, tracerToF_);
+}
+
+CellSetValues
+Solution::
+Impl::concentration(const CellSetID& tracer) const
+{
+    return solutionValues(tracer, tracer_);
+}
+
+void
+Solution::
+Impl::assign(const CellSetID& i,
+             CellSetValues&&  x,
+             SolutionMap&     soln)
+{
+    soln[i] = std::move(x);
+}
+
+CellSetValues
+Solution::
+Impl::solutionValues(const CellSetID&   i,
+                     const SolutionMap& soln) const
+{
+    auto p = soln.find(i);
+
+    if (p == soln.end()) {
+        return CellSetValues{};
+    }
+
+    return p->second;
+}
+
+// =====================================================================
+// Implementation of public interface below separator
+// =====================================================================
+
+// ---------------------------------------------------------------------
+// Class Solution
+// ---------------------------------------------------------------------
+
+Solution::Solution()
+    : pImpl_(new Impl)
+{
+}
+
+Solution::~Solution()
+{}
+
+Solution::
+Solution(const FlowDiagnostics::Solution& rhs)
+    : pImpl_(new Impl(*rhs.pImpl_))
+{}
+
+Solution::
+Solution(FlowDiagnostics::Solution&& rhs)
+    : pImpl_(std::move(rhs.pImpl_))
+{}
+
+Solution::
+Solution(std::unique_ptr<Impl> pImpl)
+    : pImpl_(std::move(pImpl))
+{}
+
+std::vector<CellSetID>
+Solution::startPoints() const
+{
+    return pImpl_->startPoints();
+}
+
+const std::vector<double>&
+Solution::timeOfFlight() const
+{
+    return pImpl_->timeOfFlight();
+}
+
+CellSetValues
+Solution::timeOfFlight(const CellSetID& tracer) const
+{
+    return pImpl_->timeOfFlight(tracer);
+}
+
+CellSetValues
+Solution::concentration(const CellSetID& tracer) const
+{
+    return pImpl_->concentration(tracer);
+}
+
+void
+Solution::assignGlobalToF(std::vector<double>&& global_tof)
+{
+    pImpl_->assignToF(std::move(global_tof));
+}
+
+void
+Solution::assign(const CellSetID& id, TimeOfFlight&& tof)
+{
+    pImpl_->assign(id, std::move(tof));
+}
+
+void
+Solution::assign(const CellSetID& id, TracerConcentration&& conc)
+{
+    pImpl_->assign(id, std::move(conc));
+}
+
+} // namespace FlowDiagnostics
+} // namespace Opm

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/Solution.hpp

@@ -0,0 +1,96 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#ifndef OPM_SOLUTION_HEADER_INCLUDED
+#define OPM_SOLUTION_HEADER_INCLUDED
+
+
+#include <opm/flowdiagnostics/CellSet.hpp>
+#include <opm/flowdiagnostics/CellSetValues.hpp>
+#include <memory>
+#include <vector>
+
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+
+    /// Results from diagnostics computations.
+    class Solution
+    {
+    public:
+        /// Default constructor.
+        Solution();
+
+        /// Destructor.
+        ~Solution();
+
+        /// Copy constructor.
+        Solution(const Solution& rhs);
+
+        /// Move constructor.
+        Solution(Solution&& rhs);
+
+        // ------ Interface for querying of solution below ------
+
+        /// Ids of stored tracer solutions.
+        std::vector<CellSetID> startPoints() const;
+
+        /// Time-of-flight field from all start points.
+        const std::vector<double>& timeOfFlight() const;
+
+        /// Time-of-flight field restricted to single tracer region.
+        CellSetValues timeOfFlight(const CellSetID& tracer) const;
+
+        /// The computed tracer field corresponding to a single tracer.
+        ///
+        /// The \c tracer must correspond to an id passed in
+        /// computeX...Diagnostics().
+        CellSetValues concentration(const CellSetID& tracer) const;
+
+        // ------ Interface for modification of solution below ------
+
+        struct TimeOfFlight
+        {
+            CellSetValues data;
+        };
+
+        struct TracerConcentration
+        {
+            CellSetValues data;
+        };
+
+        void assignGlobalToF(std::vector<double>&& global_tof);
+        void assign(const CellSetID& id, TimeOfFlight&& tof);
+        void assign(const CellSetID& id, TracerConcentration&& conc);
+
+    private:
+        class Impl;
+
+        explicit Solution(std::unique_ptr<Impl> pImpl);
+
+        std::unique_ptr<Impl> pImpl_;
+    };
+
+} // namespace FlowDiagnostics
+} // namespace Opm
+
+#endif // OPM_SOLUTION_HEADER_INCLUDED

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/Toolbox.cpp

@@ -0,0 +1,290 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#include <opm/flowdiagnostics/Toolbox.hpp>
+
+#include <opm/flowdiagnostics/CellSet.hpp>
+#include <opm/flowdiagnostics/ConnectionValues.hpp>
+#include <opm/flowdiagnostics/ConnectivityGraph.hpp>
+#include <opm/flowdiagnostics/TracerTofSolver.hpp>
+
+#include <algorithm>
+#include <exception>
+#include <map>
+#include <stdexcept>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include <opm/utility/numeric/RandomVector.hpp>
+
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+
+// ---------------------------------------------------------------------
+// Class Toolbox::Impl
+// ---------------------------------------------------------------------
+
+class Toolbox::Impl
+{
+public:
+    explicit Impl(ConnectivityGraph g);
+
+    void assignPoreVolume(const std::vector<double>& pvol);
+    void assignConnectionFlux(const ConnectionValues& flux);
+    void assignInflowFlux(const CellSetValues& inflow_flux);
+
+    Forward injDiag (const std::vector<CellSet>& start_sets);
+    Reverse prodDiag(const std::vector<CellSet>& start_sets);
+
+private:
+    ConnectivityGraph g_;
+
+    std::vector<double> pvol_;
+    ConnectionValues    flux_;
+    CellSetValues       only_inflow_flux_;
+    CellSetValues       only_outflow_flux_;
+
+    AssembledConnections inj_conn_;
+    AssembledConnections prod_conn_;
+    bool conn_built_ = false;
+
+    void buildAssembledConnections();
+};
+
+Toolbox::Impl::Impl(ConnectivityGraph g)
+    : g_   (std::move(g))
+    , pvol_()
+    , flux_(ConnectionValues::NumConnections{ 0 },
+            ConnectionValues::NumPhases     { 0 })
+    , only_inflow_flux_()
+    , only_outflow_flux_()
+{}
+
+void
+Toolbox::Impl::assignPoreVolume(const std::vector<double>& pvol)
+{
+    if (pvol.size() != g_.numCells()) {
+        throw std::logic_error("Inconsistently sized input "
+                               "pore-volume field");
+    }
+
+    pvol_ = pvol;
+}
+
+void
+Toolbox::Impl::assignConnectionFlux(const ConnectionValues& flux)
+{
+    if (flux.numConnections() != g_.numConnections()) {
+        throw std::logic_error("Inconsistently sized input "
+                               "flux field");
+    }
+
+    flux_ = flux;
+    conn_built_ = false;
+}
+
+void
+Toolbox::Impl::assignInflowFlux(const CellSetValues& inflow_flux)
+{
+    // Count the inflow (>0) fluxes.
+    const int num_items = inflow_flux.cellValueCount();
+    int num_inflows = 0;
+    for (int item = 0; item < num_items; ++item) {
+        if (inflow_flux.cellValue(item).second > 0.0) {
+            ++num_inflows;
+        }
+    }
+
+    // Reserve memory.
+    only_inflow_flux_ = CellSetValues(num_inflows);
+    only_outflow_flux_ = CellSetValues(num_items - num_inflows);
+
+    // Build in- and out-flow structures.
+    for (int item = 0; item < num_items; ++item) {
+        auto data = inflow_flux.cellValue(item);
+        if (data.second > 0.0) {
+            only_inflow_flux_.addCellValue(data.first, data.second);
+        } else {
+            only_outflow_flux_.addCellValue(data.first, -data.second);
+        }
+    }
+}
+
+Toolbox::Forward
+Toolbox::Impl::injDiag(const std::vector<CellSet>& start_sets)
+{
+    // Check that we have specified pore volume and fluxes.
+    if (pvol_.empty() || flux_.numConnections() == 0) {
+        throw std::logic_error("Must set pore volumes and fluxes before calling diagnostics.");
+    }
+
+    if (!conn_built_) {
+        buildAssembledConnections();
+    }
+
+    Solution sol;
+    using ToF = Solution::TimeOfFlight;
+    using Conc = Solution::TracerConcentration;
+
+    TracerTofSolver solver(inj_conn_, pvol_, only_inflow_flux_);
+    sol.assignGlobalToF(solver.solveGlobal(start_sets));
+
+    for (const auto& start : start_sets) {
+        auto solution = solver.solveLocal(start);
+        sol.assign(start.id(), ToF{ solution.tof });
+        sol.assign(start.id(), Conc{ solution.concentration });
+    }
+
+    return Forward{ sol };
+}
+
+Toolbox::Reverse
+Toolbox::Impl::prodDiag(const std::vector<CellSet>& start_sets)
+{
+    // Check that we have specified pore volume and fluxes.
+    if (pvol_.empty() || flux_.numConnections() == 0) {
+        throw std::logic_error("Must set pore volumes and fluxes before calling diagnostics.");
+    }
+
+    if (!conn_built_) {
+        buildAssembledConnections();
+    }
+
+    Solution sol;
+    using ToF = Solution::TimeOfFlight;
+    using Conc = Solution::TracerConcentration;
+
+    TracerTofSolver solver(prod_conn_, pvol_, only_outflow_flux_);
+    sol.assignGlobalToF(solver.solveGlobal(start_sets));
+
+    for (const auto& start : start_sets) {
+        auto solution = solver.solveLocal(start);
+        sol.assign(start.id(), ToF{ solution.tof });
+        sol.assign(start.id(), Conc{ solution.concentration });
+    }
+
+    return Reverse{ sol };
+}
+
+void
+Toolbox::Impl::buildAssembledConnections()
+{
+    // Create the data structures needed by the tracer/tof solver.
+    const size_t num_connections = g_.numConnections();
+    const size_t num_phases = flux_.numPhases();
+    inj_conn_ = AssembledConnections();
+    prod_conn_ = AssembledConnections();
+    for (size_t conn_idx = 0; conn_idx < num_connections; ++conn_idx) {
+        auto cells = g_.connection(conn_idx);
+        using ConnID = ConnectionValues::ConnID;
+        using PhaseID = ConnectionValues::PhaseID;
+        // Adding up all phase fluxes. TODO: ensure rigor, allow phase-based calculations.
+        double connection_flux = 0.0;
+        for (size_t phase = 0; phase < num_phases; ++phase) {
+            connection_flux += flux_(ConnID{conn_idx}, PhaseID{phase});
+        }
+        if (connection_flux > 0.0) {
+            inj_conn_.addConnection(cells.first, cells.second, connection_flux);
+            prod_conn_.addConnection(cells.second, cells.first, connection_flux);
+        } else {
+            inj_conn_.addConnection(cells.second, cells.first, -connection_flux);
+            prod_conn_.addConnection(cells.first, cells.second, -connection_flux);
+        }
+    }
+    const int num_cells = g_.numCells();
+    inj_conn_.compress(num_cells);
+    prod_conn_.compress(num_cells);
+
+    // Mark as built (until flux changed).
+    conn_built_ = true;
+}
+
+// =====================================================================
+// Implementation of public interface below separator
+// =====================================================================
+
+// ---------------------------------------------------------------------
+// Class Toolbox
+// ---------------------------------------------------------------------
+
+Toolbox::
+Toolbox(const ConnectivityGraph& conn)
+    : pImpl_(new Impl(conn))
+{}
+
+Toolbox::~Toolbox()
+{}
+
+Toolbox::Toolbox(Toolbox&& rhs)
+    : pImpl_(std::move(rhs.pImpl_))
+{
+}
+
+Toolbox&
+Toolbox::operator=(Toolbox&& rhs)
+{
+    pImpl_ = std::move(rhs.pImpl_);
+
+    return *this;
+}
+
+void
+Toolbox::assignPoreVolume(const std::vector<double>& pv)
+{
+    pImpl_->assignPoreVolume(pv);
+}
+
+void
+Toolbox::assignConnectionFlux(const ConnectionValues& flux)
+{
+    pImpl_->assignConnectionFlux(flux);
+}
+
+void
+Toolbox::assignInflowFlux(const CellSetValues& inflow_flux)
+{
+    pImpl_->assignInflowFlux(inflow_flux);
+}
+
+Toolbox::Forward
+Toolbox::
+computeInjectionDiagnostics(const std::vector<CellSet>& start_sets)
+{
+    return pImpl_->injDiag(start_sets);
+}
+
+Toolbox::Reverse
+Toolbox::
+computeProductionDiagnostics(const std::vector<CellSet>& start_sets)
+{
+    return pImpl_->prodDiag(start_sets);
+}
+
+
+} // namespace FlowDiagnostics
+} // namespace Opm

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/Toolbox.hpp

@@ -0,0 +1,118 @@
+/*
+  Copyright 2016 Statoil ASA.
+  Copyright 2016 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/>.
+*/
+
+#ifndef OPM_FLOWDIAGNOSTICS_TOOLBOX_HEADER_INCLUDED
+#define OPM_FLOWDIAGNOSTICS_TOOLBOX_HEADER_INCLUDED
+
+#include <opm/flowdiagnostics/CellSet.hpp>
+#include <opm/flowdiagnostics/CellSetValues.hpp>
+#include <opm/flowdiagnostics/ConnectivityGraph.hpp>
+#include <opm/flowdiagnostics/ConnectionValues.hpp>
+#include <opm/flowdiagnostics/Solution.hpp>
+
+#include <memory>
+#include <vector>
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+
+    /// Toolbox for running flow diagnostics.
+    class Toolbox
+    {
+    public:
+        /// Construct from known neighbourship relation.
+        explicit Toolbox(const ConnectivityGraph& connectivity);
+
+        /// Destructor.
+        ~Toolbox();
+
+        /// Move constructor.
+        Toolbox(Toolbox&& rhs);
+
+        /// Move assignment.
+        Toolbox& operator=(Toolbox&& rhs);
+
+        /// Assign pore volumes associated with each active cell.
+        void assignPoreVolume(const std::vector<double>& pv);
+
+        /// Assign fluxes associated with each connection.
+        void assignConnectionFlux(const ConnectionValues& flux);
+
+        /// Assign inflow fluxes, typically from wells.
+        ///
+        /// Inflow fluxes (injection) should be positive, outflow
+        /// fluxes (production) should be negative, both should be
+        /// given in the inflow_flux argument passed to this method.
+        void assignInflowFlux(const CellSetValues& inflow_flux);
+
+        struct Forward
+        {
+            const Solution fd;
+        };
+
+        struct Reverse
+        {
+            const Solution fd;
+        };
+
+        /// Compute forward time-of-flight and tracer solutions.
+        ///
+        /// An element of \code start_sets \endcode provides a set of
+        /// starting locations for a single tracer.
+        ///
+        /// Forward time-of-flight is the time needed for a neutral fluid
+        /// particle to flow from the nearest fluid source to an arbitrary
+        /// point in the model.  The tracer solutions identify cells that
+        /// are flooded by injectors.
+        ///
+        /// You must have called assignPoreVolume() and assignConnectionFlux()
+        /// before calling this method.
+        ///
+        /// The IDs of the \code start_sets \endcode must be unique.
+        Forward computeInjectionDiagnostics(const std::vector<CellSet>& start_sets);
+
+        /// Compute reverse time-of-flight and tracer solutions.
+        ///
+        /// An element of \code start_sets \endcode provides a set of
+        /// starting locations for a single tracer.
+        ///
+        /// Reverse time-of-flight is the time needed for a neutral fluid
+        /// particle to flow from an arbitrary point to the nearest fluid
+        /// sink in the model.  The tracer solutions identify cells that are
+        /// drained by producers.
+        ///
+        /// You must have called assignPoreVolume() and assignConnectionFlux()
+        /// before calling this method.
+        ///
+        /// The IDs of the \code start_sets \endcode must be unique.
+        Reverse computeProductionDiagnostics(const std::vector<CellSet>& start_sets);
+
+    private:
+        class Impl;
+
+        std::unique_ptr<Impl> pImpl_;
+    };
+
+} // namespace FlowDiagnostics
+} // namespace Opm
+
+#endif // OPM_FLOWDIAGNOSTICS_TOOLBOX_HEADER_INCLUDED

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/TracerTofSolver.cpp

@@ -0,0 +1,330 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#include <opm/flowdiagnostics/TracerTofSolver.hpp>
+#include <opm/utility/graph/tarjan.h>
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <limits>
+#include <memory>
+#include <stdexcept>
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+
+    namespace
+    {
+        std::vector<double> expandSparse(const int n, const CellSetValues& v)
+        {
+            std::vector<double> r(n, 0.0);
+            const int num_items = v.cellValueCount();
+            for (int item = 0; item < num_items; ++item) {
+                auto data = v.cellValue(item);
+                r[data.first] = data.second;
+            }
+            return r;
+        }
+    } // anonymous namespace
+
+
+
+    // This computes both in and outflux with a single traversal of the graph.
+    struct TracerTofSolver::InOutFluxComputer
+    {
+        InOutFluxComputer(const AssembledConnections& graph)
+        {
+            const int num_cells = graph.numRows();
+            influx.resize(num_cells, 0.0);
+            outflux.resize(num_cells, 0.0);
+            for (int cell = 0; cell < num_cells; ++cell) {
+                const auto nb = graph.cellNeighbourhood(cell);
+                for (const auto& conn : nb) {
+                    influx[conn.neighbour] += conn.weight;
+                    outflux[cell] += conn.weight;
+                }
+            }
+        }
+
+        std::vector<double> influx;
+        std::vector<double> outflux;
+    };
+
+
+
+
+
+    TracerTofSolver::TracerTofSolver(const AssembledConnections& graph,
+                                     const std::vector<double>& pore_volumes,
+                                     const CellSetValues& source_inflow)
+        : TracerTofSolver(graph, pore_volumes, source_inflow, InOutFluxComputer(graph))
+    {
+    }
+
+
+
+
+
+    // The InOutFluxComputer is used so that influx_ and outflux_ can be
+    // const members of the class.
+    TracerTofSolver::TracerTofSolver(const AssembledConnections& graph,
+                                     const std::vector<double>& pore_volumes,
+                                     const CellSetValues& source_inflow,
+                                     InOutFluxComputer&& inout)
+        : g_(graph)
+        , pv_(pore_volumes)
+        , influx_(std::move(inout.influx))
+        , outflux_(std::move(inout.outflux))
+        , source_term_(expandSparse(pore_volumes.size(), source_inflow))
+    {
+    }
+
+
+
+
+
+    std::vector<double> TracerTofSolver::solveGlobal(const std::vector<CellSet>& all_startsets)
+    {
+        // Reset solver variables and set source terms.
+        prepareForSolve();
+        for (const CellSet& startset : all_startsets) {
+            setupStartArray(startset);
+        }
+
+        // Compute topological ordering and solve.
+        computeOrdering();
+        solve();
+
+        // Return computed time-of-flight.
+        return tof_;
+    }
+
+
+
+
+
+    TracerTofSolver::LocalSolution TracerTofSolver::solveLocal(const CellSet& startset)
+    {
+        // Reset solver variables and set source terms.
+        prepareForSolve();
+        setupStartArray(startset);
+
+        // Compute local topological ordering and solve.
+        computeLocalOrdering(startset);
+        solve();
+
+        // Return computed time-of-flight.
+        CellSetValues local_tof;
+        const int num_elements = component_starts_.back();
+        for (int element = 0; element < num_elements; ++element) {
+            const int cell = sequence_[element];
+            local_tof.addCellValue(cell, tof_[cell]);
+        }
+        return LocalSolution{ local_tof, CellSetValues{} }; // TODO also return tracer
+    }
+
+
+
+
+
+    void TracerTofSolver::prepareForSolve()
+    {
+        // Reset instance variables.
+        const int num_cells = pv_.size();
+        is_start_.clear();
+        is_start_.resize(num_cells, 0);
+        upwind_contrib_.clear();
+        upwind_contrib_.resize(num_cells, 0.0);
+        tof_.clear();
+        tof_.resize(num_cells, -1e100);
+        num_multicell_ = 0;
+        max_size_multicell_ = 0;
+        max_iter_multicell_ = 0;
+    }
+
+
+
+
+
+    void TracerTofSolver::setupStartArray(const CellSet& startset)
+    {
+        for (const int cell : startset) {
+            is_start_[cell] = 1;
+        }
+    }
+
+
+
+
+
+    void TracerTofSolver::computeOrdering()
+    {
+        // Compute reverse topological ordering.
+        const size_t num_cells = pv_.size();
+        assert(g_.startPointers().size() == num_cells + 1);
+        struct Deleter { void operator()(TarjanSCCResult* x) { destroy_tarjan_sccresult(x); } };
+        std::unique_ptr<TarjanSCCResult, Deleter>
+            result(tarjan(num_cells, g_.startPointers().data(), g_.neighbourhood().data()));
+
+        // Must reverse ordering, since Tarjan computes reverse ordering.
+        const int ok = tarjan_reverse_sccresult(result.get());
+        if (!ok) {
+            throw std::runtime_error("Failed to reverse topological ordering in TracerTofSolver::computeOrdering()");
+        }
+
+        // Extract data from solution.
+        sequence_.resize(num_cells);
+        const int num_comp = tarjan_get_numcomponents(result.get());
+        component_starts_.resize(num_comp + 1);
+        component_starts_[0] = 0;
+        for (int comp = 0; comp < num_comp; ++comp) {
+            const TarjanComponent tc = tarjan_get_strongcomponent(result.get(), comp);
+            std::copy(tc.vertex, tc.vertex + tc.size, sequence_.begin() + component_starts_[comp]);
+            component_starts_[comp + 1] = component_starts_[comp] + tc.size;
+        }
+        assert(component_starts_.back() == int(num_cells));
+    }
+
+
+
+
+
+    void TracerTofSolver::computeLocalOrdering(const CellSet& startset)
+    {
+        // Extract start cells.
+        std::vector<int> startcells(startset.begin(), startset.end());
+
+        // Compute reverse topological ordering.
+        const size_t num_cells = pv_.size();
+        assert(g_.startPointers().size() == num_cells + 1);
+        struct ResultDeleter { void operator()(TarjanSCCResult* x) { destroy_tarjan_sccresult(x); } };
+        std::unique_ptr<TarjanSCCResult, ResultDeleter> result;
+        {
+            struct WorkspaceDeleter { void operator()(TarjanWorkSpace* x) { destroy_tarjan_workspace(x); } };
+            std::unique_ptr<TarjanWorkSpace, WorkspaceDeleter> ws(create_tarjan_workspace(num_cells));
+            result.reset(tarjan_reachable_sccs(num_cells, g_.startPointers().data(), g_.neighbourhood().data(),
+                                               startcells.size(), startcells.data(), ws.get()));
+        }
+
+        // Must reverse ordering, since Tarjan computes reverse ordering.
+        const int ok = tarjan_reverse_sccresult(result.get());
+        if (!ok) {
+            throw std::runtime_error("Failed to reverse topological ordering in TracerTofSolver::computeOrdering()");
+        }
+
+        // Extract data from solution.
+        sequence_.resize(num_cells);
+        const int num_comp = tarjan_get_numcomponents(result.get());
+        component_starts_.resize(num_comp + 1);
+        component_starts_[0] = 0;
+        for (int comp = 0; comp < num_comp; ++comp) {
+            const TarjanComponent tc = tarjan_get_strongcomponent(result.get(), comp);
+            std::copy(tc.vertex, tc.vertex + tc.size, sequence_.begin() + component_starts_[comp]);
+            component_starts_[comp + 1] = component_starts_[comp] + tc.size;
+        }
+    }
+
+
+
+
+
+    void TracerTofSolver::solve()
+    {
+        // Solve each component.
+        const int num_components = component_starts_.size() - 1;
+        for (int comp = 0; comp < num_components; ++comp) {
+            const int comp_size = component_starts_[comp + 1] - component_starts_[comp];
+            if (comp_size == 1) {
+                solveSingleCell(sequence_[component_starts_[comp]]);
+            } else {
+                solveMultiCell(comp_size, &sequence_[component_starts_[comp]]);
+            }
+        }
+    }
+
+
+
+
+
+    void TracerTofSolver::solveSingleCell(const int cell)
+    {
+        // Compute influx (divisor of tof expression).
+        double source = 2.0 * source_term_[cell];  // Initial tof for well cell equal to half fill time.
+        if (source == 0.0 && is_start_[cell]) {
+            source = std::numeric_limits<double>::infinity(); // Gives 0 tof in start cell.
+        }
+        const double total_influx = influx_[cell] + source;
+
+        // Compute effective pv (dividend of tof expression).
+        const double eff_pv = pv_[cell] + upwind_contrib_[cell];
+
+        // Compute (capped) tof.
+        if (total_influx < eff_pv / max_tof_) {
+            tof_[cell] = max_tof_;
+        } else {
+            tof_[cell] = eff_pv / total_influx;
+        }
+
+        // Set contribution for my downwind cells (if any).
+        for (const auto& conn : g_.cellNeighbourhood(cell)) {
+            const int downwind_cell = conn.neighbour;
+            const double flux = conn.weight;
+            upwind_contrib_[downwind_cell] += tof_[cell] * flux;
+        }
+    }
+
+
+
+
+
+    void TracerTofSolver::solveMultiCell(const int num_cells, const int* cells)
+    {
+        ++num_multicell_;
+        max_size_multicell_ = std::max(max_size_multicell_, num_cells);
+        // std::cout << "Multiblock solve with " << num_cells << " cells." << std::endl;
+
+        // Using a Gauss-Seidel approach.
+        double max_delta = 1e100;
+        int num_iter = 0;
+        while (max_delta > gauss_seidel_tol_) {
+            max_delta = 0.0;
+            ++num_iter;
+            for (int ci = 0; ci < num_cells; ++ci) {
+                const int cell = cells[ci];
+                const double tof_before = tof_[cell];
+                solveSingleCell(cell);
+                max_delta = std::max(max_delta, std::fabs(tof_[cell] - tof_before));
+            }
+            // std::cout << "Max delta = " << max_delta << std::endl;
+        }
+        max_iter_multicell_ = std::max(max_iter_multicell_, num_iter);
+    }
+
+
+
+
+} // namespace FlowDiagnostics
+} // namespace Opm

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/flowdiagnostics/TracerTofSolver.hpp

@@ -0,0 +1,124 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#ifndef OPM_TRACERTOFSOLVER_HEADER_INCLUDED
+#define OPM_TRACERTOFSOLVER_HEADER_INCLUDED
+
+#include <opm/flowdiagnostics/CellSet.hpp>
+#include <opm/flowdiagnostics/CellSetValues.hpp>
+#include <opm/utility/graph/AssembledConnections.hpp>
+#include <vector>
+
+namespace Opm
+{
+namespace FlowDiagnostics
+{
+
+
+    /// Class for solving the tracer and time-of-flight equations.
+    ///
+    /// Implements a first-order finite volume solver for
+    /// (single-phase) time-of-flight using reordering.
+    /// The equation solved is:
+    ///     \f[v \cdot \nabla\tau = \phi\f]
+    /// in which \f$ v \f$ is the fluid velocity, \f$ \tau \f$ is time-of-flight and
+    /// \f$ \phi \f$ is the porosity. This is a boundary value problem, and
+    /// \f$ \tau \f$ is considered zero on all inflow boundaries (or well inflows).
+    ///
+    /// The tracer equation is the same, except for the right hand side which is zero
+    /// instead of \f[\phi\f].
+    class TracerTofSolver
+    {
+    public:
+        /// Initialize solver with a given flow graph (a weighted,
+        /// directed asyclic graph) containing the out-fluxes from
+        /// each cell, pore volumes and inflow sources (positive).
+        TracerTofSolver(const AssembledConnections& graph,
+                        const std::vector<double>& pore_volumes,
+                        const CellSetValues& source_inflow);
+
+        /// Compute the global (combining all sources) time-of-flight of each cell.
+        ///
+        /// TODO: also compute tracer solution.
+        std::vector<double> solveGlobal(const std::vector<CellSet>& all_startsets);
+
+        /// Output data struct for solveLocal().
+        struct LocalSolution {
+            CellSetValues tof;
+            CellSetValues concentration;
+        };
+
+        /// Compute a local solution tracer and time-of-flight solution.
+        ///
+        /// Local means that only cells downwind from he startset are considered.
+        /// The solution is therefore potentially sparse.
+        /// TODO: not implemented!
+        LocalSolution solveLocal(const CellSet& startset);
+
+    private:
+
+        // --------------  Private data members --------------
+
+        const AssembledConnections& g_;
+        const std::vector<double>& pv_;
+        const std::vector<double> influx_;
+        const std::vector<double> outflux_;
+        std::vector<double> source_term_;
+        std::vector<char> is_start_; // char to avoid the nasty vector<bool> specialization
+        std::vector<int> sequence_;
+        std::vector<int> component_starts_;
+        std::vector<double> upwind_contrib_;
+        std::vector<double> tof_;
+        int num_multicell_ = 0;
+        int max_size_multicell_ = 0;
+        int max_iter_multicell_ = 0;
+        const double gauss_seidel_tol_ = 1e-3;
+        const double max_tof_ = 200.0 * 365.0 * 24.0 * 60.0 * 60.0; // 200 years.
+
+        // --------------  Private helper class --------------
+
+        struct InOutFluxComputer;
+
+        // --------------  Private methods --------------
+
+        TracerTofSolver(const AssembledConnections& graph,
+                        const std::vector<double>& pore_volumes,
+                        const CellSetValues& source_inflow,
+                        InOutFluxComputer&& inout);
+
+        void prepareForSolve();
+
+        void setupStartArray(const CellSet& startset);
+
+        void computeOrdering();
+
+        void computeLocalOrdering(const CellSet& startset);
+
+        void solve();
+
+        void solveSingleCell(const int cell);
+
+        void solveMultiCell(const int num_cells, const int* cells);
+    };
+
+} // namespace FlowDiagnostics
+} // namespace Opm
+
+#endif // OPM_TRACERTOFSOLVER_HEADER_INCLUDED

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/utility/graph/AssembledConnections.cpp

@@ -0,0 +1,487 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#include <opm/utility/graph/AssembledConnections.hpp>
+
+#include <algorithm>
+#include <cassert>
+#include <exception>
+#include <ios>
+#include <iterator>
+#include <ostream>
+#include <stdexcept>
+#include <utility>
+
+// ---------------------------------------------------------------------
+// Class Opm::AssembledConnections::Connections
+// ---------------------------------------------------------------------
+
+void
+Opm::AssembledConnections::
+Connections::add(const int i, const int j)
+{
+    i_.push_back(i);
+    j_.push_back(j);
+
+    max_i_ = std::max(max_i_, i_.back());
+    max_j_ = std::max(max_j_, j_.back());
+}
+
+void
+Opm::AssembledConnections::
+Connections::add(const int i, const int j, const double v)
+{
+    this->add(i, j);
+
+    v_.push_back(v);
+}
+
+void
+Opm::AssembledConnections::Connections::clear()
+{
+    WeightVector().swap(v_);
+    EntityVector().swap(j_);
+    EntityVector().swap(i_);
+}
+
+bool
+Opm::AssembledConnections::Connections::empty() const
+{
+    return i_.empty();
+}
+
+bool
+Opm::AssembledConnections::Connections::isValid() const
+{
+    return (i_.size() == j_.size())
+        && (v_.empty() || (v_.size() == i_.size()));
+}
+
+bool
+Opm::AssembledConnections::Connections::isWeighted() const
+{
+    return ! v_.empty();
+}
+
+int
+Opm::AssembledConnections::Connections::maxRow() const
+{
+    return max_i_;
+}
+
+int
+Opm::AssembledConnections::Connections::maxCol() const
+{
+    return max_j_;
+}
+
+Opm::AssembledConnections::Connections::EntityVector::size_type
+Opm::AssembledConnections::Connections::nnz() const
+{
+    return i_.size();
+}
+
+const Opm::AssembledConnections::Connections::EntityVector&
+Opm::AssembledConnections::Connections::i() const
+{
+    return i_;
+}
+
+const Opm::AssembledConnections::Connections::EntityVector&
+Opm::AssembledConnections::Connections::j() const
+{
+    return j_;
+}
+
+const Opm::AssembledConnections::Connections::WeightVector&
+Opm::AssembledConnections::Connections::v() const
+{
+    return v_;
+}
+
+// =====================================================================
+
+// ---------------------------------------------------------------------
+// Class Opm::AssembledConnections::CSR
+// ---------------------------------------------------------------------
+
+void
+Opm::AssembledConnections::
+CSR::create(const Connections& conns, const Offset numrows)
+{
+    if (! conns.empty() &&
+        (static_cast<Offset>(conns.maxRow()) >= numrows))
+    {
+        throw std::invalid_argument("Input graph contains more "
+                                    "source vertices than are "
+                                    "implied by explicit size of "
+                                    "adjacency matrix");
+    }
+
+    this->assemble(conns);
+
+    this->sort();
+
+    // Must be called *after* sort().
+    this->condenseDuplicates();
+
+    if (conns.isWeighted()) {
+        this->accumulateConnWeights(conns.v());
+    }
+
+    const auto nRows = this->ia().size() - 1;
+    if (nRows < numrows) {
+        this->ia_.insert(this->ia_.end(),
+                         numrows - nRows,
+                         this->ia().back());
+    }
+}
+
+const Opm::AssembledConnections::Start&
+Opm::AssembledConnections::CSR::ia() const
+{
+    return ia_;
+}
+
+const Opm::AssembledConnections::Neighbours&
+Opm::AssembledConnections::CSR::ja() const
+{
+    return ja_;
+}
+
+const Opm::AssembledConnections::ConnWeight&
+Opm::AssembledConnections::CSR::sa() const
+{
+    return sa_;
+}
+
+void
+Opm::AssembledConnections::
+CSR::assemble(const Connections& conns)
+{
+    {
+        const auto numRows = conns.maxRow() + 1;
+
+        this->accumulateRowEntries(numRows, conns.i());
+    }
+
+    this->createGraph(conns.i(), conns.j());
+
+    this->numRows_ = conns.maxRow() + 1;
+    this->numCols_ = conns.maxCol() + 1;
+}
+
+void
+Opm::AssembledConnections::CSR::sort()
+{
+    // Transposition is, effectively, a linear time bucket insert, so
+    // transposing the structure twice creates a structure with colum
+    // indices in (ascendingly) sorted order.
+
+    this->transpose();
+    this->transpose();
+}
+
+void
+Opm::AssembledConnections::CSR::condenseDuplicates()
+{
+    // Note: Must be called *after* sort().
+
+    const auto colIdx = this->ja_;
+    auto elmIdx       = this->elmIdx_;
+    auto end          = colIdx.begin();
+
+    this->ja_    .clear();
+    this->elmIdx_.clear();
+
+    for (decltype(this->ia_.size())
+             row = 0, numRows = this->ia_.size() - 1;
+         row < numRows; ++row)
+    {
+        auto begin = end;
+
+        std::advance(end, this->ia_[row + 1] - this->ia_[row + 0]);
+
+        const auto q = this->ja_.size();
+
+        this->unique(begin, end);
+
+        this->ia_[row + 0] = q;
+    }
+
+    this->remapElementIndex(std::move(elmIdx));
+
+    // Record final table sizes.
+    this->ia_.back() = this->ja_.size();
+}
+
+void
+Opm::AssembledConnections::
+CSR::accumulateConnWeights(const std::vector<double>& v)
+{
+    if (v.size() != this->elmIdx_.size()) {
+        throw std::logic_error("Connection Weights must be "
+                               "provided for each connection");
+    }
+
+    this->sa_.assign(this->ja_.size(), 0.0);
+
+    for (decltype(v.size())
+             conn = 0, numConn = v.size();
+         conn < numConn; ++conn)
+    {
+        this->sa_[ this->elmIdx_[conn] ] += v[conn];
+    }
+}
+
+void
+Opm::AssembledConnections::
+CSR::accumulateRowEntries(const int               numRows,
+                          const std::vector<int>& rowIdx)
+{
+    assert (numRows >= 0);
+
+    auto vecIdx = [](const int i)
+    {
+        return static_cast<Start::size_type>(i);
+    };
+
+    this->ia_.assign(vecIdx(numRows) + 1, vecIdx(0));
+
+    for (const auto& row : rowIdx) {
+        this->ia_[vecIdx(row) + 1] += 1;
+    }
+
+    // Note index range: 1..numRows inclusive.
+    for (Start::size_type
+             i = 1, n = vecIdx(numRows);
+         i <= n; ++i)
+    {
+        this->ia_[0] += this->ia_[i];
+        this->ia_[i]  = this->ia_[0] - this->ia_[i];
+    }
+}
+
+void
+Opm::AssembledConnections::
+CSR::createGraph(const std::vector<int>& rowIdx,
+                 const std::vector<int>& colIdx)
+{
+    assert (this->ia_[0] == rowIdx.size());
+
+    auto vecIdx = [](const int i)
+    {
+        return static_cast<Start::size_type>(i);
+    };
+
+    this->ja_.resize(rowIdx.size());
+
+    this->elmIdx_.clear();
+    this->elmIdx_.reserve(rowIdx.size());
+
+    for (decltype(rowIdx.size())
+             nz = 0, nnz = rowIdx.size();
+         nz < nnz; ++nz)
+    {
+        const auto k = ia_[vecIdx(rowIdx[nz]) + 1] ++;
+
+        this->ja_[k] = colIdx[nz];
+
+        this->elmIdx_.push_back(k);
+    }
+
+    this->ia_[0] = 0;
+}
+
+void
+Opm::AssembledConnections::CSR::transpose()
+{
+    auto elmIdx = this->elmIdx_;
+
+    {
+        const auto rowIdx = this->expandStartPointers();
+        const auto colIdx = this->ja_;
+
+        this->accumulateRowEntries(this->numCols_, colIdx);
+
+        this->createGraph(colIdx, rowIdx);
+    }
+
+    this->remapElementIndex(std::move(elmIdx));
+
+    std::swap(this->numRows_, this->numCols_);
+}
+
+std::vector<int>
+Opm::AssembledConnections::
+CSR::expandStartPointers() const
+{
+    auto rowIdx = std::vector<int>{};
+    rowIdx.reserve(this->ia_.back());
+
+    auto row = 0;
+
+    for (decltype(this->ia_.size())
+             i = 0, m = this->ia_.size() - 1; i < m; ++i, ++row)
+    {
+        const auto n = this->ia_[i + 1] - this->ia_[i + 0];
+
+        rowIdx.insert(rowIdx.end(), n, row);
+    }
+
+    return rowIdx;
+}
+
+void
+Opm::AssembledConnections::
+CSR::unique(Neighbours::const_iterator begin,
+            Neighbours::const_iterator end)
+{
+    // We assume that we're only called *after* sort() whence duplicate
+    // elements appear consecutively in [begin, end).
+    //
+    // Note: This is essentially the same as std::unique(begin, end) save
+    // for the return value.  However, we also record the 'elmIdx_' mapping
+    // to later have the ability to accumulate 'sa_'.
+
+    while (begin != end) {
+        // Note: Order of ja_ and elmIdx_ matters here.
+        this->elmIdx_.push_back(this->ja_.size());
+        this->ja_    .push_back(*begin);
+
+        while ((++begin != end) &&
+               ( *begin == this->ja_.back()))
+        {
+            this->elmIdx_.push_back(this->elmIdx_.back());
+        }
+    }
+}
+
+void
+Opm::AssembledConnections::CSR::remapElementIndex(Start&& elmIdx)
+{
+    for (auto& i : elmIdx) {
+        i = this->elmIdx_[i];
+    }
+
+    this->elmIdx_.swap(elmIdx);
+}
+
+// =====================================================================
+
+// ---------------------------------------------------------------------
+// Class Opm::AssembledConnections
+// ---------------------------------------------------------------------
+
+void
+Opm::AssembledConnections::
+addConnection(const int i,
+              const int j)
+{
+    conns_.add(i, j);
+}
+
+void
+Opm::AssembledConnections::
+addConnection(const int    i,
+              const int    j,
+              const double v)
+{
+    conns_.add(i, j, v);
+}
+
+void
+Opm::AssembledConnections::compress(const std::size_t numRows)
+{
+    if (! conns_.isValid()) {
+        throw std::logic_error("Cannot compress invalid "
+                               "connection list");
+    }
+
+    csr_.create(conns_, numRows);
+
+    conns_.clear();
+}
+
+Opm::AssembledConnections::Offset
+Opm::AssembledConnections::numRows() const
+{
+    return this->startPointers().size() - 1;
+}
+
+const Opm::AssembledConnections::Start&
+Opm::AssembledConnections::startPointers() const
+{
+    return csr_.ia();
+}
+
+const Opm::AssembledConnections::Neighbours&
+Opm::AssembledConnections::neighbourhood() const
+{
+    return csr_.ja();
+}
+
+const Opm::AssembledConnections::ConnWeight&
+Opm::AssembledConnections::connectionWeight() const
+{
+    return csr_.sa();
+}
+
+Opm::AssembledConnections::CellNeighbours
+Opm::AssembledConnections::cellNeighbourhood(const int cell) const
+{
+    const Offset beg = startPointers()[cell];
+    const Offset end = startPointers()[cell + 1];
+    const int* nb = neighbourhood().data();
+    assert(connectionWeight().size() == neighbourhood().size());
+    const double* w = connectionWeight().data();
+    return CellNeighbours{ {nb + beg, w + beg}, {nb + end, w + end} };
+}
+
+
+
+
+// Note: not a member of class AssembledConnections.
+std::ostream&
+Opm::operator<<(std::ostream& os, const Opm::AssembledConnections& ac)
+{
+    // Set output stream format, store original settings.
+    const auto oprec = os.precision(16);
+    const auto oflags = os.setf(std::ios_base::scientific);
+
+    // Write connections cell-by-cell.
+    const int num_cells = ac.numRows();
+    for (int cell = 0; cell < num_cells; ++cell) {
+        const auto nb = ac.cellNeighbourhood(cell);
+        for (const auto& conn : nb) {
+            os << cell << ' ' << conn.neighbour << ' ' << conn.weight << '\n';
+        }
+    }
+
+    // Restore original stream settings.
+    os.precision(oprec);
+    os.setf(oflags);
+
+    return os;
+}

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/utility/graph/AssembledConnections.hpp

@@ -0,0 +1,247 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#ifndef OPM_ASSEMBLEDCONNECTIONS_HEADER_INCLUDED
+#define OPM_ASSEMBLEDCONNECTIONS_HEADER_INCLUDED
+
+#include <opm/utility/graph/AssembledConnectionsIteration.hpp>
+
+#include <cstddef>
+#include <iosfwd>
+#include <vector>
+
+/// \file
+///
+/// Facility for converting collection of entity pairs or weighted entity
+/// triplets into a sparse (CSR) adjacency matrix representation of a graph.
+/// Supports O(nnz) compression and, if applicable, accumulation of weight
+/// values for repeated entity pairs.
+
+namespace Opm {
+
+    /// Form CSR adjacency matrix representation of graph provided as a list
+    /// of connections between entities, possibly including edge weights.
+    class AssembledConnections
+    {
+    public:
+        /// Add connection between entities.
+        ///
+        /// \param[in] i Primary (source) entity.  Used as row index.
+        ///
+        /// \param[in] j Secondary (sink) entity.  Used as column index.
+        void addConnection(const int i, const int j);
+
+        /// Add weighted connection between entities.
+        ///
+        /// \param[in] i Primary (source) entity.  Used as row index.
+        ///
+        /// \param[in] j Secondary (sink) entity.  Used as column index.
+        ///
+        /// \param[in] v Edge weight associated to connection.
+        void addConnection(const int i, const int j, const double v);
+
+        /// Form CSR adjacency matrix representation of input graph from
+        /// connections established in previous calls to addConnection().
+        ///
+        /// A call to function compress() will fail unless all previously
+        /// established connections have an explicit edge weight (\code
+        /// addConnection(i,j,v) \endcode) or none of those connections have
+        /// an explicit edge weight (\code addConnection(i,j) \endcode).
+        ///
+        /// This method destroys the connection list so if there are
+        /// subsequent calls to method addConnection() then those will
+        /// effectively create a new graph.
+        ///
+        /// \param[in] numRows Number of rows in resulting CSR matrix.  If
+        ///     prior calls to addConnection() supply source entity IDs (row
+        ///     indices) greater than or equal to \p numRows, then method
+        ///     compress() will throw \code std::invalid_argument \endcode.
+        void compress(const std::size_t numRows);
+
+        /// Representation of neighbouring entities.
+        using Neighbours = std::vector<int>;
+
+        /// Offset into neighbour array.
+        using Offset = Neighbours::size_type;
+
+        /// CSR start pointers.
+        using Start = std::vector<Offset>;
+
+        /// Aggregate connection weights
+        using ConnWeight = std::vector<double>;
+
+        /// Range of neighbours connected to a particular entity.  Includes
+        /// edge weights.
+        using CellNeighbours = SimpleIteratorRange<NeighbourhoodIterator>;
+
+        /// Retrieve number of rows (source entities) in input graph.
+        /// Corresponds to value of argument passed to compress().  Valid
+        /// only after calling compress().
+        Offset numRows() const;
+
+        /// Retrieve raw CSR start pointers pertaining to current input
+        /// graph.  Only valid after compress() is called.
+        const Start& startPointers() const;
+
+        /// Retrieve raw CSR column indices pertaining to current input
+        /// graph.  Only valid after compress() is called.  The neighbours
+        /// of entity \c i are
+        ///
+        ///   \code
+        ///     neighbourhood()[startPointers()[i + 0]] ...
+        ///     neighbourhood()[startPointers()[i + 1] - 1]
+        ///   \endcode
+        const Neighbours& neighbourhood() const;
+
+        /// Retrieve accumulated connection weights for each aggregate
+        /// neighbour relation in the adjacency matrix.  Only valid if the
+        /// input graph is specified in terms of explicit edge weights
+        /// (\code addConnection(i,j,v) \endcode).
+        ///
+        /// Weights have the same ordering as the neighbours represented by
+        /// neighbourhood().
+        const ConnWeight& connectionWeight() const;
+
+        /// Provide a range over a cell neighbourhood.
+        ///
+        /// Example usage:
+        ///
+        /// \code
+        ///    for (const auto& connection : cellNeighbourhood(cell) {
+        ///        // connection.neigbour is the neigbouring cell
+        ///        // connection.weight   is the corresponding connection weight
+        ///     }
+        /// \endcode
+        ///
+        /// This method can only be called if the weight-providing
+        /// overload of addConnection() was used to build the instance.
+        CellNeighbours cellNeighbourhood(const int cell) const;
+
+    private:
+        class Connections
+        {
+        public:
+            using EntityVector = std::vector<int>;
+            using WeightVector = std::vector<double>;
+
+            void add(const int i, const int j);
+            void add(const int i, const int j, const double v);
+
+            void clear();
+
+            bool empty() const;
+
+            bool isValid() const;
+
+            bool isWeighted() const;
+
+            int maxRow() const;
+            int maxCol() const;
+
+            EntityVector::size_type nnz() const;
+
+            const EntityVector& i() const;
+            const EntityVector& j() const;
+            const WeightVector& v() const;
+
+        private:
+            EntityVector i_;
+            EntityVector j_;
+            WeightVector v_;
+
+            int max_i_{ -1 };
+            int max_j_{ -1 };
+        };
+
+        class CSR
+        {
+        public:
+            void create(const Connections& conns,
+                        const Offset       numRows);
+
+            const Start&      ia() const;
+            const Neighbours& ja() const;
+            const ConnWeight& sa() const;
+
+        private:
+            Start      ia_;
+            Neighbours ja_;
+            ConnWeight sa_;
+
+            Start      elmIdx_;
+
+            int        numRows_{ 0 };
+            int        numCols_{ 0 };
+
+            // ---------------------------------------------------------
+            // Implementation of create()
+            // ---------------------------------------------------------
+
+            void assemble(const Connections& conns);
+
+            void sort();
+
+            void condenseDuplicates();
+
+            void accumulateConnWeights(const std::vector<double>& v);
+
+            // ---------------------------------------------------------
+            // Implementation of assemble()
+            // ---------------------------------------------------------
+
+            void accumulateRowEntries(const int               numRows,
+                                      const std::vector<int>& rowIdx);
+
+            void createGraph(const std::vector<int>& rowIdx,
+                             const std::vector<int>& colIdx);
+
+            // ---------------------------------------------------------
+            // General utilities
+            // ---------------------------------------------------------
+
+            void transpose();
+
+            std::vector<int> expandStartPointers() const;
+
+            void unique(Neighbours::const_iterator begin,
+                        Neighbours::const_iterator end);
+
+            void remapElementIndex(Start&& elmIdx);
+        };
+
+        Connections conns_;
+        CSR         csr_;
+    };
+
+
+
+
+
+    /// Output an AssembledConnections object to a stream in text format.
+    ///
+    /// Assumes compress() has been called, and that the
+    /// weight-providing overload of addConnection() was used to build
+    /// the ac instance, as for AssembledConnections::cellNeighbourhood().
+    std::ostream& operator<<(std::ostream& os, const AssembledConnections& ac);
+
+
+} // namespace Opm
+
+#endif // OPM_ASSEMBLEDCONNECTIONS_HEADER_INCLUDED

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/utility/graph/AssembledConnectionsIteration.hpp

@@ -0,0 +1,100 @@
+/*
+  Copyright 2016 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/>.
+*/
+
+#ifndef OPM_ASSEMBLEDCONNECTIONSITERATION_HEADER_INCLUDED
+#define OPM_ASSEMBLEDCONNECTIONSITERATION_HEADER_INCLUDED
+
+#include <cassert>
+
+namespace Opm {
+
+    /// Holds data for a single connection.
+    /// The neighbour will typically be a cell
+    /// index, and the weight will typically be
+    /// a flux.
+    struct ConnectionData
+    {
+        int neighbour;
+        double weight;
+    };
+
+
+
+
+    /// Iterator over the graph neighbourhood
+    /// of a cell, typically used from AssembledConnections.
+    class NeighbourhoodIterator
+    {
+    public:
+        NeighbourhoodIterator(const int* neighbour_iter,
+                              const double* weight_iter)
+            : neighbour_iter_(neighbour_iter),
+              weight_iter_(weight_iter)
+        {
+        }
+
+        ConnectionData operator*()
+        {
+            return { *neighbour_iter_, *weight_iter_ };
+        }
+
+        bool operator!=(const NeighbourhoodIterator& other) const
+        {
+            assert((neighbour_iter_ != other.neighbour_iter_) == (weight_iter_ != other.weight_iter_));
+            return neighbour_iter_ != other.neighbour_iter_;
+        }
+
+        NeighbourhoodIterator& operator++()
+        {
+            ++neighbour_iter_;
+            ++weight_iter_;
+            return *this;
+        }
+
+    private:
+        const int* neighbour_iter_;
+        const double* weight_iter_;
+    };
+
+
+
+
+    /// A straightforward range class whose only purpose is
+    /// to easily allow range-for loops over [begin_, end_).
+    template <class Iterator>
+    struct SimpleIteratorRange
+    {
+        Iterator begin() const
+        {
+            return begin_;
+        }
+
+        Iterator end() const
+        {
+            return end_;
+        }
+
+        Iterator begin_;
+        Iterator end_;
+    };
+
+
+} // namespace Opm
+
+#endif // OPM_ASSEMBLEDCONNECTIONSITERATION_HEADER_INCLUDED

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/utility/graph/tarjan.c

@@ -0,0 +1,546 @@
+/*
+Copyright (C) 2012 (c) Jostein R. Natvig <jostein natvig at gmail.com>
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+of the Software, and to permit persons to whom the Software is furnished to do
+so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif /* HAVE_CONFIG_H */
+
+#include <opm/utility/graph/tarjan.h>
+
+#include <assert.h>
+#include <stddef.h>
+#include <stdlib.h>
+#include <string.h>
+
+enum VertexMark { DONE = -2, REMAINING = -1 };
+
+struct TarjanWorkSpace
+{
+    size_t nvert;               /**< Number of vertices */
+
+    int *status;                /**< Vertex processing status */
+    int *link;                  /**< Vertex low-link */
+    int *time;                  /**< Vertex DFS discovery time */
+};
+
+struct TarjanSCCResult
+{
+    size_t  ncomp;              /**< Running number of SCCs */
+    size_t *start;              /**< SCC start pointers (CSR) */
+    int    *vert;               /**< SCC contents (CSR format) */
+
+    size_t *vstack, *vstart;    /**< Vertex stack (reuse 'start') */
+    int    *cstack, *cstart;    /**< Component stack (reuse 'vert') */
+};
+
+static void
+initialise_stacks(const size_t nvert, struct TarjanSCCResult *scc)
+{
+    /* Vertex and component stacks grow from high end downwards */
+
+    scc->vstack = scc->vstart = scc->start + (nvert + 0);
+    scc->cstack = scc->cstart = scc->vert  + (nvert - 1);
+}
+
+static void
+assign_int_vector(size_t n, const int val, int *v)
+{
+    size_t i;
+
+    for (i = 0; i < n; i++) { v[i] = val; }
+}
+
+static void
+copy_int_vector(const size_t n, const int* src, int *dst)
+{
+    size_t i, thres;
+
+    /* Cut-off point from thin air (i.e., no measurement) */
+    thres = (size_t) (1LU << 10);
+
+    if (n < thres) {
+        for (i = 0; i < n; i++) {
+            dst[i] = src[i];
+        }
+    }
+    else {
+        memcpy(dst, src, n * sizeof *dst);
+    }
+}
+
+static struct TarjanSCCResult *
+allocate_scc_result(const size_t nvert)
+{
+    struct TarjanSCCResult *scc, scc0 = { 0 };
+
+    scc = malloc(1 * sizeof *scc);
+
+    if (scc != NULL) {
+        *scc = scc0;
+
+        scc->start = malloc((nvert + 1) * sizeof *scc->start);
+        scc->vert  = malloc(nvert       * sizeof *scc->vert);
+
+        if ((scc->start == NULL) || (scc->vert == NULL)) {
+            destroy_tarjan_sccresult(scc);
+
+            scc = NULL;
+        }
+    }
+
+    return scc;
+}
+
+#define MIN(a,b) (((a) < (b)) ? (a) : (b))
+
+#define SWAP(x,y,tmp)                           \
+    do {                                        \
+        (tmp) = (x);                            \
+        (x)   = (y);                            \
+        (y)   = (tmp);                          \
+    } while (0)
+
+/* Stack grows to lower addresses */
+#define peek(stack) (*((stack) + 1))
+#define push(stack) *(stack)--
+
+static void
+capture_scc(const size_t            c,
+            struct TarjanWorkSpace *work,
+            struct TarjanSCCResult *scc)
+{
+    int    vertex;
+    size_t v;
+
+    /* Initialise strong component end pointer.
+     *
+     * Supports the equivalent of push_back() on the
+     * component. */
+    scc->start[ scc->ncomp + 1 ] =
+        scc->start[ scc->ncomp + 0 ];
+
+    do {
+        assert (scc->cstack != scc->cstart);
+
+        /* Pop strong component stack */
+        v = vertex = *++scc->cstack;
+        work->status[v] = DONE;
+
+        /* Capture component vertex in VERT while
+         * advancing component end pointer */
+        scc->vert[ scc->start[ scc->ncomp + 1 ] ++ ] = vertex;
+    } while (v != c);
+
+    /* Record completion of strong component.
+     * Prepare for next. */
+    scc->ncomp += 1;
+}
+
+static void
+complete_dfs_from_vertex(const size_t            c,
+                         struct TarjanWorkSpace *work,
+                         struct TarjanSCCResult *scc)
+{
+    size_t v;
+
+    /* Record strong component if 'c' is root */
+    if (work->link[c] == work->time[c]) {
+        capture_scc(c, work, scc);
+    }
+
+    /* Pop 'c' from DFS (vertex) stack */
+    ++scc->vstack;
+
+    if (scc->vstack != scc->vstart) {
+        v = peek(scc->vstack);
+
+        work->link[v] = MIN(work->link[v], work->link[c]);
+    }
+}
+
+static void
+dfs_from_descendant(const size_t            c,
+                    const size_t           *ia,
+                    const int              *ja,
+                    struct TarjanWorkSpace *work,
+                    struct TarjanSCCResult *scc)
+{
+    size_t child;
+
+    assert (work->status[c] > 0);
+
+    child = ja[ia[c] + (work->status[c] - 1)];
+
+    /* decrement descendant count of c*/
+    --work->status[c];
+
+    if (work->status[child] == REMAINING) {
+        /* push child */
+        push(scc->vstack) = child;
+    }
+    else if (work->status[child] >= 0) {
+        work->link[c] = MIN(work->link[c], work->time[child]);
+    }
+    else {
+        assert (work->status[child] == DONE);
+    }
+}
+
+static void
+discover_vertex(const size_t            c,
+                const size_t           *ia,
+                int                    *time,
+                struct TarjanWorkSpace *work,
+                struct TarjanSCCResult *scc)
+{
+    /* number of descendants of c */
+    work->status[c] = (int) (ia[c + 1] - ia[c]);
+    work->time[c]   = work->link[c] = (*time)++;
+
+    push(scc->cstack) = (int) c;
+}
+
+/*
+ * Note: 'work->status' serves dual purpose during processing.
+ *
+ *   status[c] = DONE      -> Vertex 'c' fully processed.
+ *               REMAINING -> Vertex 'c' undiscovered.
+ *               0         -> Vertex 'c' not fully classified--there are no
+ *                            remaining descendants for this vertex but we
+ *                            do not yet know if the vertex is an SCC all of
+ *                            itself or if it is part of a larger component.
+ *
+ *   status[c] > 0 -> Vertex 'c' has 'status[c]' remaining
+ *                    descendants (i.e., successors/children).
+ */
+
+static void
+tarjan_initialise_csr(struct TarjanSCCResult *scc)
+{
+    scc->ncomp                   = 0;
+    scc->start[ scc->ncomp + 0 ] = 0;
+}
+
+static void
+tarjan_initialise(struct TarjanWorkSpace *work,
+                  struct TarjanSCCResult *scc)
+{
+    initialise_stacks(work->nvert, scc);
+
+    /* Initialise status for all vertices */
+    assign_int_vector(work->nvert, REMAINING, work->status);
+
+    tarjan_initialise_csr(scc);
+}
+
+static void
+tarjan_local(const size_t            start,
+             const size_t           *ia,
+             const int              *ja,
+             struct TarjanWorkSpace *work,
+             struct TarjanSCCResult *scc)
+{
+    int    time;
+    size_t c;
+
+    push(scc->vstack) = start;
+
+    time = 0;
+
+    while (scc->vstack != scc->vstart) {
+        c = peek(scc->vstack);
+
+        assert (work->status[c] != DONE);
+        assert (work->status[c] >= -2);
+
+        if (work->status[c] == REMAINING) {
+            /* Discover() writes to work->status[c] */
+
+            discover_vertex(c, ia, &time, work, scc);
+        }
+
+        if (work->status[c] == 0) {
+            /* All descendants of 'c' processed */
+
+            complete_dfs_from_vertex(c, work, scc);
+        }
+        else {
+            /* Process next descendant of 'c' */
+
+            dfs_from_descendant(c, ia, ja, work, scc);
+        }
+    }
+
+    assert (scc->cstack == scc->cstart);
+}
+
+static size_t
+tarjan_find_start(const size_t                  seed,
+                  const struct TarjanWorkSpace *work)
+{
+    size_t start = seed;
+
+    while ((start < work->nvert) &&
+           (work->status[start] == DONE))
+    {
+        ++start;
+    }
+
+    return start;
+}
+
+static void
+tarjan_global(const size_t           *ia,
+              const int              *ja,
+              struct TarjanWorkSpace *work,
+              struct TarjanSCCResult *scc)
+{
+    size_t start;
+
+    assert ((work != NULL) && "Work array must be non-NULL");
+    assert ((scc  != NULL) && "Result array must be non-NULL");
+
+    tarjan_initialise(work, scc);
+
+    for (start = tarjan_find_start(  0  , work);
+         start < work->nvert;
+         start = tarjan_find_start(start, work))
+    {
+        tarjan_local(start, ia, ja, work, scc);
+    }
+}
+
+static void
+reverse_scc_core(const struct TarjanSCCResult *src,
+                 struct TarjanSCCResult       *dst)
+{
+    size_t comp, n;
+
+    tarjan_initialise_csr(dst);
+
+    for (comp = src->ncomp; comp > 0; --comp, ++dst->ncomp) {
+        n = src->start[ comp - 0 ] - src->start[ comp - 1 ];
+
+        dst->start[ dst->ncomp + 1 ] =
+            dst->start[ dst->ncomp + 0 ];
+
+        copy_int_vector(n, src->vert + src->start[ comp - 1 ],
+                        dst->vert + dst->start[ dst->ncomp + 1 ]);
+
+        dst->start[ dst->ncomp + 1 ] += n;
+    }
+}
+
+static void
+swap_scc(struct TarjanSCCResult *scc1,
+         struct TarjanSCCResult *scc2)
+{
+    SWAP(scc1->start, scc2->start, scc1->vstart);
+    SWAP(scc1->vert , scc2->vert , scc1->cstart);
+}
+
+/* ======================================================================
+ * Public interface below separator
+ * ====================================================================== */
+
+struct TarjanWorkSpace *
+create_tarjan_workspace(const size_t nvert)
+{
+    struct TarjanWorkSpace *ws, ws0 = { 0 };
+
+    ws = malloc(1 * sizeof *ws);
+
+    if (ws != NULL) {
+        *ws = ws0;
+
+        ws->status = malloc(3 * nvert * sizeof *ws->status);
+
+        if (ws->status == NULL) {
+            destroy_tarjan_workspace(ws);
+
+            ws = NULL;
+        }
+        else {
+            ws->nvert = nvert;
+
+            ws->link = ws->status + ws->nvert;
+            ws->time = ws->link   + ws->nvert;
+        }
+    }
+
+    return ws;
+}
+
+void
+destroy_tarjan_workspace(struct TarjanWorkSpace *ws)
+{
+    if (ws != NULL) {
+        free(ws->status);
+    }
+
+    free(ws);
+}
+
+void
+destroy_tarjan_sccresult(struct TarjanSCCResult *scc)
+{
+    if (scc != NULL) {
+        /* Reverse order of acquisition. */
+        free(scc->vert);
+        free(scc->start);
+    }
+
+    free(scc);
+}
+
+size_t
+tarjan_get_numcomponents(const struct TarjanSCCResult *scc)
+{
+    return scc->ncomp;
+}
+
+struct TarjanComponent
+tarjan_get_strongcomponent(const struct TarjanSCCResult *scc,
+                           const size_t                  compID)
+{
+    size_t start;
+
+    assert ((compID < tarjan_get_numcomponents(scc)) &&
+            "Component ID out of bounds");
+
+    struct TarjanComponent c = { 0 };
+
+    start = scc->start[compID];
+
+    c.size   = scc->start[compID + 1] - start;
+    c.vertex = &scc->vert[start];
+
+    return c;
+}
+
+/*
+  Compute the strong components of directed graph G(edges, vertices),
+  return components in reverse topological sorted sequence.
+  Complexity O(|vertices|+|edges|). See "http://en.wikipedia.org/wiki/
+  Tarjan's_strongly_connected_components_algorithm".
+
+  nv    - number of vertices
+
+  ia,ja - adjacency matrix for directed graph in compressed sparse row
+          format: vertex i has directed edges to vertices ja[ia[i]],
+          ..., ja[ia[i+1]-1].
+ */
+
+/*--------------------------------------------------------------------*/
+struct TarjanSCCResult *
+tarjan(const size_t  nv,
+       const size_t *ia,
+       const int    *ja)
+/*--------------------------------------------------------------------*/
+{
+    struct TarjanWorkSpace *work = NULL;
+    struct TarjanSCCResult *scc  = NULL;
+
+    work = create_tarjan_workspace(nv);
+    scc  = allocate_scc_result(nv);
+
+    if ((work == NULL) || (scc == NULL)) {
+        destroy_tarjan_sccresult(scc);
+        destroy_tarjan_workspace(work);
+
+        work = NULL;
+        scc  = NULL;
+    }
+
+    if (scc != NULL) {
+        tarjan_global(ia, ja, work, scc);
+    }
+
+    destroy_tarjan_workspace(work);
+
+    return scc;
+}
+
+struct TarjanSCCResult *
+tarjan_reachable_sccs(const size_t            nv,
+                      const size_t           *ia,
+                      const int              *ja,
+                      const size_t            nstart,
+                      const int              *start_pts,
+                      struct TarjanWorkSpace *work)
+{
+    size_t                  i, start;
+    struct TarjanSCCResult *scc;
+
+    assert ((work != NULL) && "Work-space must be non-NULL");
+    assert ((work->nvert >= nv) &&
+            "Work-space must be large enough to accommodate graph");
+
+#if defined(NDEBUG)
+    /* Suppress 'unused argument' diagnostic ('nv' only in assert()) */
+    (void) nv;
+#endif  /* defined(NDEBUG) */
+
+    scc = allocate_scc_result(work->nvert);
+
+    if (scc != NULL) {
+        tarjan_initialise(work, scc);
+
+        for (i = 0; i < nstart; i++) {
+            start = start_pts[i];
+
+            if (work->status[start] == DONE) {
+                continue;
+            }
+
+            tarjan_local(start, ia, ja, work, scc);
+        }
+    }
+
+    return scc;
+}
+
+int
+tarjan_reverse_sccresult(struct TarjanSCCResult *scc)
+{
+    int ok;
+
+    struct TarjanSCCResult *rev;
+
+    rev = allocate_scc_result(scc->start[ scc->ncomp ]);
+
+    ok  = rev != NULL;
+
+    if (ok) {
+        reverse_scc_core(scc, rev);
+
+        swap_scc(rev, scc);
+
+        destroy_tarjan_sccresult(rev);
+    }
+
+    return ok;
+}
+
+/* Local Variables:    */
+/* c-basic-offset:4    */
+/* End:                */

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/utility/graph/tarjan.h

@@ -0,0 +1,235 @@
+/*
+Copyright (C) 2012 (c) Jostein R. Natvig <jostein natvig at gmail.com>
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+of the Software, and to permit persons to whom the Software is furnished to do
+so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
+
+/**
+ * \file
+ *
+ * Simple implementation of of Tarjan's algorithm for computing the
+ * strongly connected components of a directed graph, \f$G(V,E)\f$.
+ * Run-time complexity is \f$O(|V| + |E|)\f$.
+ *
+ * The implementation is based on
+ * http://en.wikipedia.org/wiki/Tarjan's_strongly_connected_components_algorithm
+ */
+
+#ifndef TARJAN_H_INCLUDED
+#define TARJAN_H_INCLUDED
+
+#include <stddef.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif  /* __cplusplus */
+
+/**
+ * Abstract datatype providing backing store for working dataset in SCC
+ * processing.
+ *
+ * Intentionally incomplete.
+ */
+struct TarjanWorkSpace;
+
+/**
+ * Abstract datatype providing backing store for SCC result data.
+ *
+ * Intentionally incomplete.
+ */
+struct TarjanSCCResult;
+
+/**
+ * Contents of single strong component.
+ *
+ * Essentially a reference type.
+ */
+struct TarjanComponent
+{
+    /**
+     * Number of vertices in component.
+     */
+    size_t  size;
+
+    /**
+     * IDs of vertices in strong component.  Non-owning pointer.  The
+     * vertices of this component are
+     * \code
+     *   vertex[0 .. size-1]
+     * \endcode
+     */
+    const int *vertex;
+};
+
+/**
+ * Create work-space for Tarjan algorithm on graph with specified number of
+ * nodes (vertices).
+ *
+ * \param[in] nvert Number of graph vertices.
+ *
+ * \return Backing store for intermediate data created during SCC
+ * processing.  \c NULL in case of allocation failure.  Dispose of
+ * work-space using function destroy_tarjan_workspace().
+ */
+struct TarjanWorkSpace *
+create_tarjan_workspace(const size_t nvert);
+
+/**
+ * Dispose of backing store for intermediate SCC/Tarjan data.
+ *
+ * \param[in,out] ws Work-space structure procured in a previous call to
+ * function create_tarjan_workspace().  Invalid upon return.
+ */
+void
+destroy_tarjan_workspace(struct TarjanWorkSpace *ws);
+
+/**
+ * Dispose of backing store for SCC result data.
+ *
+ * \param[in,out] scc SCC result data from a previous call to function
+ * tarjan() or tarjan_reachable_sccs().  Invalid upon return.
+ */
+void
+destroy_tarjan_sccresult(struct TarjanSCCResult *scc);
+
+/**
+ * Retrieve number of strong components in result dataset.
+ *
+ * \param[in] scc Collection of strongly connected components obtained from
+ * a previous call to function tarjan() or tarjan_reachable_sccs().
+ *
+ * \return Number of strong components in result dataset.
+ */
+size_t
+tarjan_get_numcomponents(const struct TarjanSCCResult *scc);
+
+/**
+ * Get access to single strong component from SCC result dataset.
+ *
+ * \param[in] scc Collection of strongly connected components obtained from
+ * a previous call to function tarjan() or tarjan_reachable_sccs().
+ *
+ * \param[in] compID Linear ID of single strong component.  Must be in the
+ * range \code [0 .. tarjan_get_numcomponents(scc) - 1] \endcode.
+ *
+ * \return Single strong component corresponding to explicit linear ID.
+ */
+struct TarjanComponent
+tarjan_get_strongcomponent(const struct TarjanSCCResult *scc,
+                           const size_t                  compID);
+
+/**
+ * Compute the strongly connected components of a directed graph,
+ * \f$G(V,E)\f$.
+ *
+ * The components are returned in reverse topological sorted sequence.
+ *
+ * \param[in] nv Number of graph vertices.
+ *
+ * \param[in] ia CSR sparse matrix start pointers corresponding to
+ *               downstream vertices.
+ *
+ * \param[in] ja CSR sparse matrix representation of out-neighbours in a
+ *               directed graph: vertex \c i has directed edges to vertices
+ *               \code ja[ia[i]], ..., ja[ia[i + 1] - 1] \endcode.
+ *
+ * \return Strong component result dataset.  Owning pointer.  Dispose of
+ * associate memory by calling the destructor destroy_tarjan_sccresult().
+ * Returns \c NULL in case of failure to allocate the result set or internal
+ * work-space.
+ */
+struct TarjanSCCResult *
+tarjan(const size_t  nv,
+       const size_t *ia,
+       const int    *ja);
+
+/**
+ * Compute the strongly connected components of reachable set in a directed
+ * graph \f$G(V,E)\$ when starting from a sparse collection of start points.
+ *
+ * The components are returned in reverse topological order.  Use function
+ * tarjan_reverse_sccresult() to access components in topological order.
+ *
+ * \param[in] nv Number of graph vertices.
+ *
+ * \param[in] ia CSR sparse matrix start pointers corresponding to
+ *               downstream vertices.
+ *
+ * \param[in] ja CSR sparse matrix representation of out-neighbours in a
+ *               directed graph: vertex \c i has directed edges to vertices
+ *               \code ja[ia[i]], ..., ja[ia[i + 1] - 1] \endcode.
+ *
+ * \param[in] nstart Number of start points from which to initiate reachable
+ *               set calculations.
+ *
+ * \param[in] start_pts Vertices from which to initiate reachable set
+ *               calculations.  Array of size \p nstart.
+ *
+ * \param[in,out] ws Backing-store for quantities needed during SCC
+ *               processing.  Obtained from a previous call to function
+ *               create_tarjan_workspace().
+ *
+ * \return Strong component result dataset.  Owning pointer.  Dispose of
+ * associate memory by calling the destructor destroy_tarjan_sccresult().
+ * Returns \c NULL in case of failure to allocate the result set or internal
+ * work-space.
+ */
+struct TarjanSCCResult *
+tarjan_reachable_sccs(const size_t            nv,
+                      const size_t           *ia,
+                      const int              *ja,
+                      const size_t            nstart,
+                      const int              *start_pts,
+                      struct TarjanWorkSpace *ws);
+
+/**
+ * Reverse order of SCCs represented by result set.
+ *
+ * Maintain order of vertices within each strong component.
+ *
+ * If successful, invalidates any TarjanComponent data held by caller.
+ *
+ * If unsuccessful, the original component order is maintained in the input
+ * argument.  In this case, the requested order can be obtained by
+ * traversing the linear component IDs in reverse using accessor function
+ * tarjan_get_strongcomponent().
+ *
+ * \param[in,out] scc On input, collection of strongly connected components
+ *                    obtained from a previous call to function tarjan() or
+ *                    function tarjan_reachable_sccs().  On output,
+ *                    reordered collection of SCCs such that linear access
+ *                    in ascending order of component IDs accesses the
+ *                    original strong components in reverse order.
+ *
+ * \return Whether or not the components could be reversed.  One (true) if
+ * order reversal successful and zero (false) otherwise--typically due to
+ * failure to allocate internal resources for operation.
+ */
+int
+tarjan_reverse_sccresult(struct TarjanSCCResult *scc);
+
+#ifdef __cplusplus
+}
+#endif  /* __cplusplus */
+
+#endif /* TARJAN_H_INCLUDED */
+
+/* Local Variables:    */
+/* c-basic-offset:4    */
+/* End:                */

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/utility/numeric/RandomVector.cpp

@@ -0,0 +1,100 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#include <opm/utility/numeric/RandomVector.hpp>
+
+#include <random>
+
+class Opm::RandomVector::Impl
+{
+public:
+    Sample normal(const Size   n,
+                  const double mean,
+                  const double stdev);
+
+    std::vector<int> integer(const Size n,
+                             const int  min,
+                             const int  max);
+
+private:
+    using BitGenerator = std::mt19937;
+
+    BitGenerator gen_;
+};
+
+Opm::RandomVector::Sample
+Opm::RandomVector::
+Impl::normal(const Size n, const double mean, const double stdev)
+{
+    auto distr = std::normal_distribution<>{ mean, stdev };
+
+    auto sample = Sample{};
+    sample.reserve(n);
+
+    for (auto i = 0*n; i < n; ++i) {
+        sample.push_back(distr(gen_));
+    }
+
+    return sample;
+}
+
+std::vector<int>
+Opm::RandomVector::
+Impl::integer(const Size n, const int min, const int max
+)
+{
+    auto distr = std::uniform_int_distribution<>{ min, max };
+
+    auto idx = std::vector<int>{};
+    idx.reserve(n);
+
+    for (auto i = 0*n; i < n; ++i) {
+        idx.push_back(distr(gen_));
+    }
+
+    return idx;
+}
+
+// =====================================================================
+
+Opm::RandomVector::RandomVector()
+    : pImpl_(new Impl())
+{}
+
+Opm::RandomVector::~RandomVector()
+{}
+
+Opm::RandomVector::Sample
+Opm::RandomVector::normal(const Size   n,
+                          const double mean,
+                          const double stdev)
+{
+    return pImpl_->normal(n, mean, stdev);
+}
+
+std::vector<int>
+Opm::RandomVector::index(const Size n, const int maxIdx)
+{
+    return pImpl_->integer(n, 0, maxIdx);
+}

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/opm/utility/numeric/RandomVector.hpp

@@ -0,0 +1,59 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#ifndef OPM_RANDOMVECTOR_HEADER_INCLUDED
+#define OPM_RANDOMVECTOR_HEADER_INCLUDED
+
+#include <memory>
+#include <vector>
+
+namespace Opm
+{
+    class RandomVector
+    {
+    public:
+        using Sample = std::vector<double>;
+        using Size   = Sample::size_type;
+
+        RandomVector();
+        ~RandomVector();
+
+        RandomVector(const RandomVector& rhs) = delete;
+        RandomVector(RandomVector&&      rhs) = delete;
+
+        RandomVector& operator=(const RandomVector& rhs) = delete;
+        RandomVector& operator=(RandomVector&&      rhs) = delete;
+
+        Sample normal(const Size   n,
+                      const double mean  = 0.0,
+                      const double stdev = 1.0);
+
+        std::vector<int> index(const Size n,
+                               const int  maxIdx);
+
+    private:
+        class Impl;
+
+        std::unique_ptr<Impl> pImpl_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_RANDOMVECTOR_HEADER_INCLUDED

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_assembledconnections.cpp

@@ -0,0 +1,479 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE
+
+#define BOOST_TEST_MODULE TEST_ASSEMBLED_CONNECTIONS
+
+#include <opm/common/utility/platform_dependent/disable_warnings.h>
+#include <boost/test/unit_test.hpp>
+#include <opm/common/utility/platform_dependent/reenable_warnings.h>
+
+#include <opm/utility/graph/AssembledConnections.hpp>
+
+#include <exception>
+#include <stdexcept>
+
+namespace {
+
+    template <class Collection1, class Collection2>
+    void check_is_close(const Collection1& c1, const Collection2& c2)
+    {
+        BOOST_REQUIRE_EQUAL(c1.size(), c2.size());
+
+        if (! c1.empty()) {
+            auto i1 = c1.begin(), e1 = c1.end();
+            auto i2 = c2.begin();
+
+            for (; i1 != e1; ++i1, ++i2) {
+                BOOST_CHECK_CLOSE(*i1, *i2, 1.0e-10);
+            }
+        }
+    }
+
+} // Namespace Anonymous
+
+BOOST_AUTO_TEST_SUITE(Two_By_Two)
+
+BOOST_AUTO_TEST_CASE (Constructor)
+{
+    auto g = Opm::AssembledConnections{};
+}
+
+BOOST_AUTO_TEST_CASE (Zero_To_One)
+{
+    auto g = Opm::AssembledConnections{};
+
+    g.addConnection(0, 1);
+
+    g.compress(4);
+
+    BOOST_CHECK_EQUAL(g.numRows(), 4);
+
+    {
+        const auto start = g.startPointers();
+
+        const auto expect_ia = std::vector<int>{ 0, 1, 1, 1, 1 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(start    .begin(), start    .end(),
+                                      expect_ia.begin(), expect_ia.end());
+    }
+
+    {
+        const auto neigh = g.neighbourhood();
+
+        const auto expect_ja = std::vector<int>{ 1 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(neigh    .begin(), neigh    .end(),
+                                      expect_ja.begin(), expect_ja.end());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (Zero_To_One_Two)
+{
+    auto g = Opm::AssembledConnections{};
+
+    g.addConnection(0, 2);
+    g.addConnection(0, 1);
+
+    g.compress(4);
+
+    BOOST_CHECK_EQUAL(g.numRows(), 4);
+
+    {
+        const auto start = g.startPointers();
+
+        const auto expect_ia = std::vector<int>{ 0, 2, 2, 2, 2 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(start    .begin(), start    .end(),
+                                      expect_ia.begin(), expect_ia.end());
+    }
+
+    {
+        const auto neigh = g.neighbourhood();
+
+        const auto expect_ja = std::vector<int>{ 1, 2 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(neigh    .begin(), neigh    .end(),
+                                      expect_ja.begin(), expect_ja.end());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (Zero_To_One_Two_Duplicate)
+{
+    auto g = Opm::AssembledConnections{};
+
+    for (auto i = 0, n = 3; i < n; ++i) {
+        g.addConnection(0, 2);
+    }
+
+    g.addConnection(0, 1);
+
+    for (auto i = 0, n = 3; i < n; ++i) {
+        g.addConnection(0, 2);
+    }
+
+    for (auto i = 0, n = 3; i < n; ++i) {
+        g.addConnection(0, 1);
+    }
+
+    g.compress(4);
+
+    BOOST_CHECK_EQUAL(g.numRows(), 4);
+
+    {
+        const auto start = g.startPointers();
+
+        const auto expect_ia = std::vector<int>{ 0, 2, 2, 2, 2 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(start    .begin(), start    .end(),
+                                      expect_ia.begin(), expect_ia.end());
+    }
+
+    {
+        const auto neigh = g.neighbourhood();
+
+        const auto expect_ja = std::vector<int>{ 1, 2 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(neigh    .begin(), neigh    .end(),
+                                      expect_ja.begin(), expect_ja.end());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (No_Out_Edge_From_High_Vertex)
+{
+    // Vertex of highest numerical ID not referenced as source vertex in
+    // connection list.  We must still produce a complete adjacency
+    // representation that is aware of all vertices when the expected size
+    // is provided.
+    auto g = Opm::AssembledConnections{};
+
+    g.addConnection(0, 1, 1.0);
+    g.addConnection(0, 2, 1.0);
+    g.addConnection(1, 3, 1.0);
+    g.addConnection(2, 3, 1.0);
+
+    g.compress(4);
+
+    BOOST_CHECK_EQUAL(g.numRows(), 4);
+
+    {
+        const auto start = g.startPointers();
+
+        const auto expect_ia = std::vector<int>{ 0, 2, 3, 4, 4 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(start    .begin(), start    .end(),
+                                      expect_ia.begin(), expect_ia.end());
+    }
+
+    {
+        const auto neigh = g.neighbourhood();
+
+        const auto expect_ja = std::vector<int>{ 1, 2, 3, 3 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(neigh    .begin(), neigh    .end(),
+                                      expect_ja.begin(), expect_ja.end());
+    }
+
+    {
+        const auto w = g.connectionWeight();
+
+        const auto expect_w = std::vector<double>{
+            1.0, 1.0, 1.0, 1.0
+        };
+
+        check_is_close(w, expect_w);
+    }
+}
+
+BOOST_AUTO_TEST_CASE (Unexpected_Vertex_ID)
+{
+    // We expect a lower number of vertices than are actually present.  This
+    // is an error.
+
+    auto g = Opm::AssembledConnections{};
+
+    g.addConnection(0, 1);
+    g.addConnection(0, 2);
+    g.addConnection(1, 3);
+    g.addConnection(3, 2);
+
+    BOOST_CHECK_THROW(g.compress(3), std::invalid_argument);
+}
+
+BOOST_AUTO_TEST_CASE (Isolated_Node)
+{
+    auto g = Opm::AssembledConnections{};
+
+    g.addConnection(2, 0);
+    g.addConnection(0, 2);
+    g.addConnection(0, 0);
+
+    g.compress(3);
+
+    BOOST_CHECK_EQUAL(g.numRows(), 3);
+
+    {
+        const auto start = g.startPointers();
+
+        const auto expect_ia = std::vector<int>{ 0, 2, 2, 3 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(start    .begin(), start    .end(),
+                                      expect_ia.begin(), expect_ia.end());
+    }
+
+    {
+        const auto neigh = g.neighbourhood();
+
+        const auto expect_ja = std::vector<int>{ 0, 2, 0 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(neigh    .begin(), neigh    .end(),
+                                      expect_ja.begin(), expect_ja.end());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (All_To_All)
+{
+    auto g = Opm::AssembledConnections{};
+
+    const auto n = 4;
+    for (auto i = 0*n; i < n; ++i) {
+        for (auto j = 0*n; j < n; ++j) {
+            g.addConnection(j, i);
+        }
+    }
+
+    g.compress(n);
+
+    BOOST_CHECK_EQUAL(g.numRows(), n);
+
+    {
+        const auto start = g.startPointers();
+
+        const auto expect_ia = std::vector<int>{ 0, 4, 8, 12, 16 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(start    .begin(), start    .end(),
+                                      expect_ia.begin(), expect_ia.end());
+    }
+
+    {
+        const auto neigh = g.neighbourhood();
+
+        const auto expect_ja = std::vector<int>{
+            0, 1, 2, 3,
+            0, 1, 2, 3,
+            0, 1, 2, 3,
+            0, 1, 2, 3,
+        };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(neigh    .begin(), neigh    .end(),
+                                      expect_ja.begin(), expect_ja.end());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (All_To_All_Duplicate)
+{
+    auto g = Opm::AssembledConnections{};
+
+    const auto n = 4;
+
+    for (auto k = 0*n; k < n; ++k) {
+        for (auto i = 0*n; i < n; ++i) {
+            for (auto j = 0*n; j < n; ++j) {
+                g.addConnection(i, i);
+                g.addConnection(i, j);
+                g.addConnection(j, i);
+                g.addConnection(j, j);
+            }
+        }
+    }
+
+    g.compress(n);
+
+    BOOST_CHECK_EQUAL(g.numRows(), n);
+
+    {
+        const auto start = g.startPointers();
+
+        const auto expect_ia = std::vector<int>{ 0, 4, 8, 12, 16 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(start    .begin(), start    .end(),
+                                      expect_ia.begin(), expect_ia.end());
+    }
+
+    {
+        const auto neigh = g.neighbourhood();
+
+        const auto expect_ja = std::vector<int>{
+            0, 1, 2, 3,
+            0, 1, 2, 3,
+            0, 1, 2, 3,
+            0, 1, 2, 3,
+        };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(neigh    .begin(), neigh    .end(),
+                                      expect_ja.begin(), expect_ja.end());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (Weighted_Graph_Single)
+{
+    auto g = Opm::AssembledConnections{};
+
+    g.addConnection(0, 2,   0.2);
+    g.addConnection(0, 1, - 0.1);
+
+    g.addConnection(1, 3,  13.0);
+
+    g.addConnection(2, 3, -23.0);
+
+    g.compress(4);
+
+    BOOST_CHECK_EQUAL(g.numRows(), 4);
+
+    {
+        const auto start = g.startPointers();
+
+        const auto expect_ia = std::vector<int>{ 0, 2, 3, 4, 4 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(start    .begin(), start    .end(),
+                                      expect_ia.begin(), expect_ia.end());
+    }
+
+    {
+        const auto neigh = g.neighbourhood();
+
+        const auto expect_ja = std::vector<int>{ 1, 2, 3, 3 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(neigh    .begin(), neigh    .end(),
+                                      expect_ja.begin(), expect_ja.end());
+    }
+
+    {
+        const auto w = g.connectionWeight();
+
+        const auto expect_w = std::vector<double>{
+            - 0.1, 0.2,
+             13.0,
+            -23.0,
+        };
+
+        check_is_close(w, expect_w);
+    }
+
+    {
+        const int num_cells = g.startPointers().size() - 1;
+        int num_conn = 0;
+        double weight_sum = 0.0;
+        for (int cell = 0; cell < num_cells; ++cell) {
+            for (const auto& conn : g.cellNeighbourhood(cell)) {
+                ++num_conn;
+                weight_sum += conn.weight;
+            }
+        }
+        BOOST_CHECK_EQUAL(num_conn, 4);
+        BOOST_CHECK_CLOSE(weight_sum, -9.9, 1e-10);
+    }
+}
+
+BOOST_AUTO_TEST_CASE (Weighted_Graph_Multiple)
+{
+    auto g = Opm::AssembledConnections{};
+
+    const auto count = 4;
+
+    for (auto i = 0*count; i < count; ++i) {
+        g.addConnection(0, 2,   0.2);
+        g.addConnection(0, 1, - 0.1);
+
+        g.addConnection(1, 3,  13.0);
+
+        g.addConnection(2, 3, -23.0);
+    }
+
+    g.compress(4);
+
+    BOOST_CHECK_EQUAL(g.numRows(), 4);
+
+    {
+        const auto start = g.startPointers();
+
+        const auto expect_ia = std::vector<int>{ 0, 2, 3, 4, 4 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(start    .begin(), start    .end(),
+                                      expect_ia.begin(), expect_ia.end());
+    }
+
+    {
+        const auto neigh = g.neighbourhood();
+
+        const auto expect_ja = std::vector<int>{ 1, 2, 3, 3 };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(neigh    .begin(), neigh    .end(),
+                                      expect_ja.begin(), expect_ja.end());
+    }
+
+    {
+        const auto w = g.connectionWeight();
+
+        const auto expect_w = std::vector<double>{
+            count * (- 0.1), count * 0.2,
+            count *   13.0 ,
+            count * (-23.0),
+        };
+
+        check_is_close(w, expect_w);
+    }
+
+    {
+        const int num_cells = g.startPointers().size() - 1;
+        int num_conn = 0;
+        double weight_sum = 0.0;
+        for (int cell = 0; cell < num_cells; ++cell) {
+            for (const auto& conn : g.cellNeighbourhood(cell)) {
+                ++num_conn;
+                weight_sum += conn.weight;
+            }
+        }
+        BOOST_CHECK_EQUAL(num_conn, 4);
+        BOOST_CHECK_CLOSE(weight_sum, count*(-9.9), 1e-10);
+    }
+}
+
+BOOST_AUTO_TEST_CASE (Compress_Invalid_Throw)
+{
+    auto g = Opm::AssembledConnections{};
+
+    g.addConnection(0, 1);
+    g.addConnection(0, 2, 1.0);
+
+    // Can't mix weighted and unweighted edges.
+    BOOST_CHECK_THROW(g.compress(3), std::logic_error);
+}
+
+BOOST_AUTO_TEST_SUITE_END()

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_cellset.cpp

@@ -0,0 +1,128 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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 <config.h>
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE
+
+#define BOOST_TEST_MODULE TEST_CELLSET
+
+#include <boost/test/unit_test.hpp>
+
+#include <opm/flowdiagnostics/CellSet.hpp>
+
+#include <algorithm>
+
+
+using Opm::FlowDiagnostics::CellSet;
+using Opm::FlowDiagnostics::CellSetID;
+
+
+BOOST_AUTO_TEST_SUITE(CellSetIDTest)
+
+BOOST_AUTO_TEST_CASE (Construct)
+{
+    {
+        const auto i = CellSetID{};
+
+        BOOST_CHECK_EQUAL(i.to_string(), "");
+    }
+
+    {
+        const auto name = std::string("Injection");
+
+        const auto i = CellSetID(name);
+
+        BOOST_CHECK_EQUAL(i.to_string(), name);
+    }
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+BOOST_AUTO_TEST_SUITE(CellSetTest)
+
+BOOST_AUTO_TEST_CASE (Constructor)
+{
+    {
+        auto s = CellSet{};
+
+        BOOST_CHECK_EQUAL(s.id().to_string(), "");
+    }
+
+    {
+        const auto name = std::string("Test-Ctor");
+
+        auto s = CellSet{};
+        {
+            s.identify(CellSetID(name));
+        }
+
+        BOOST_CHECK_EQUAL(s.id().to_string(), name);
+    }
+}
+
+BOOST_AUTO_TEST_CASE (AssignCells)
+{
+    auto s = CellSet{};
+
+    const auto cells = std::vector<int>
+        { 0, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000 };
+
+    for (const auto& cell : cells) {
+        s.insert(cell);
+    }
+
+    auto out = std::vector<int>(s.begin(), s.end());
+    {
+        std::sort(out.begin(), out.end());
+    }
+
+    BOOST_CHECK_EQUAL_COLLECTIONS(out  .begin(), out  .end(),
+                                  cells.begin(), cells.end());
+}
+
+BOOST_AUTO_TEST_CASE (Duplicates)
+{
+    auto s = CellSet{};
+
+    const auto cells = std::vector<int>
+        { 0, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000 };
+
+    for (auto i = 0; i < 2; ++i) {
+        for (const auto& cell : cells) {
+            s.insert(cell);
+        }
+    }
+
+    auto out = std::vector<int>(s.begin(), s.end());
+    {
+        std::sort(out.begin(), out.end());
+    }
+
+    BOOST_CHECK_EQUAL_COLLECTIONS(out  .begin(), out  .end(),
+                                  cells.begin(), cells.end());
+}
+
+BOOST_AUTO_TEST_SUITE_END()

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_cellsetvalues.cpp

@@ -0,0 +1,89 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE
+
+#define BOOST_TEST_MODULE TEST_CELLSETVALUES
+
+#include <opm/common/utility/platform_dependent/disable_warnings.h>
+#include <boost/test/unit_test.hpp>
+#include <opm/common/utility/platform_dependent/reenable_warnings.h>
+
+#include <opm/flowdiagnostics/CellSetValues.hpp>
+
+using Opm::FlowDiagnostics::CellSetValues;
+
+BOOST_AUTO_TEST_SUITE(CellSet_Values)
+
+BOOST_AUTO_TEST_CASE (Constructor)
+{
+    {
+        CellSetValues s{};
+    }
+
+    {
+        auto s = CellSetValues{ 100 };
+    }
+}
+
+BOOST_AUTO_TEST_CASE (AssignValues)
+{
+    auto s = CellSetValues{ 100 };
+
+    for (decltype(s.cellValueCount())
+             i = 0, n = 100;
+         i < n; ++i)
+    {
+        s.addCellValue(100 - i, i * 10.0);
+    }
+
+    BOOST_CHECK_EQUAL(s.cellValueCount(), 100);
+
+    {
+        const auto a = s.cellValue(0);
+
+        BOOST_CHECK_EQUAL(a.first , 100);
+        BOOST_CHECK_CLOSE(a.second, 0.0, 1.0e-10);
+    }
+
+    {
+        const auto a = s.cellValue(s.cellValueCount() - 1);
+
+        BOOST_CHECK_EQUAL(a.first , 1);
+        BOOST_CHECK_CLOSE(a.second, 990.0, 1.0e-10);
+    }
+
+    {
+        const auto a = s.cellValue(50);
+
+        BOOST_CHECK_EQUAL(a.first , 50);
+        BOOST_CHECK_CLOSE(a.second, 500.0, 1.0e-10);
+    }
+}
+
+BOOST_AUTO_TEST_SUITE_END()

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_connectionvalues.cpp

@@ -0,0 +1,97 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE
+
+#define BOOST_TEST_MODULE TEST_CONNECTIONVALUES
+
+#include <opm/common/utility/platform_dependent/disable_warnings.h>
+#include <boost/test/unit_test.hpp>
+#include <opm/common/utility/platform_dependent/reenable_warnings.h>
+
+#include <opm/flowdiagnostics/ConnectionValues.hpp>
+
+BOOST_AUTO_TEST_SUITE(Connection_Values)
+
+
+using Opm::FlowDiagnostics::ConnectionValues;
+
+
+BOOST_AUTO_TEST_CASE (Constructor)
+{
+    using NConn = ConnectionValues::NumConnections;
+    using NPhas = ConnectionValues::NumPhases;
+
+    const auto nconn = NConn{4};
+    const auto nphas = NPhas{2};
+
+    auto v = ConnectionValues(nconn, nphas);
+
+    BOOST_CHECK_EQUAL(v.numPhases()     , nphas.total);
+    BOOST_CHECK_EQUAL(v.numConnections(), nconn.total);
+}
+
+BOOST_AUTO_TEST_CASE (AssignValues)
+{
+    using NConn = ConnectionValues::NumConnections;
+    using NPhas = ConnectionValues::NumPhases;
+
+    const auto nconn = NConn{4};
+    const auto nphas = NPhas{2};
+
+    using ConnID = ConnectionValues::ConnID;
+    using PhasID = ConnectionValues::PhaseID;
+
+    auto v = ConnectionValues(nconn, nphas);
+
+    {
+        for (decltype(v.numConnections())
+                 conn = 0, numconn = v.numConnections();
+             conn < numconn; ++conn)
+        {
+            for (decltype(v.numPhases())
+                     phas = 0, numphas = v.numPhases();
+                 phas < numphas; ++phas)
+            {
+                v(ConnID{conn}, PhasID{phas}) =
+                    conn*numphas + phas;
+            }
+        }
+    }
+
+    {
+        const auto w = v;
+
+        BOOST_CHECK_CLOSE(w(ConnID{0}, PhasID{0}), 0.0, 1.0e-10);
+
+        BOOST_CHECK_CLOSE(w(ConnID{ nconn.total - 1 },
+                            PhasID{ 0 }), 6.0, 1.0e-10);
+    }
+}
+
+BOOST_AUTO_TEST_SUITE_END()

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_connectivitygraph.cpp

@@ -0,0 +1,94 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE
+
+#define BOOST_TEST_MODULE TEST_CONNECTIVITY_GRAPH
+
+#include <opm/common/utility/platform_dependent/disable_warnings.h>
+#include <boost/test/unit_test.hpp>
+#include <opm/common/utility/platform_dependent/reenable_warnings.h>
+
+#include <opm/flowdiagnostics/ConnectivityGraph.hpp>
+
+using Opm::FlowDiagnostics::ConnectivityGraph;
+
+BOOST_AUTO_TEST_SUITE(Two_By_Two)
+
+BOOST_AUTO_TEST_CASE (Constructor)
+{
+    const auto g =
+        ConnectivityGraph(4,
+                          { 0 , 1 ,
+                            0 , 2 ,
+                            1 , 3 ,
+                            2 , 3 });
+
+    BOOST_CHECK_EQUAL(g.numCells(),       4);
+    BOOST_CHECK_EQUAL(g.numConnections(), 4);
+}
+
+BOOST_AUTO_TEST_CASE (ConnectionList)
+{
+    const auto g =
+        ConnectivityGraph(4,
+                          { 0, 1 ,
+                            0, 2 ,
+                            1, 3 ,
+                            2, 3 });
+
+    {
+        const auto c = g.connection(0);
+
+        BOOST_CHECK_EQUAL(c.first , 0);
+        BOOST_CHECK_EQUAL(c.second, 1);
+    }
+
+    {
+        const auto c = g.connection(1);
+
+        BOOST_CHECK_EQUAL(c.first , 0);
+        BOOST_CHECK_EQUAL(c.second, 2);
+    }
+
+    {
+        const auto c = g.connection(2);
+
+        BOOST_CHECK_EQUAL(c.first , 1);
+        BOOST_CHECK_EQUAL(c.second, 3);
+    }
+
+    {
+        const auto c = g.connection(3);
+
+        BOOST_CHECK_EQUAL(c.first , 2);
+        BOOST_CHECK_EQUAL(c.second, 3);
+    }
+}
+
+BOOST_AUTO_TEST_SUITE_END()

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_flowdiagnosticstool.cpp

@@ -0,0 +1,417 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif // HAVE_CONFIG_H
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif // HAVE_DYNAMIC_BOOST_TEST
+
+#define NVERBOSE
+
+#define BOOST_TEST_MODULE TEST_FLOWDIAGNOSTICSTOOL
+
+#include <boost/test/unit_test.hpp>
+
+#include <opm/flowdiagnostics/Toolbox.hpp>
+
+#include <opm/flowdiagnostics/CellSet.hpp>
+#include <opm/flowdiagnostics/ConnectionValues.hpp>
+#include <opm/flowdiagnostics/ConnectivityGraph.hpp>
+#include <opm/utility/numeric/RandomVector.hpp>
+
+#include <algorithm>
+
+using namespace Opm::FlowDiagnostics;
+
+namespace
+{
+    std::size_t
+    numIntConn(const std::size_t nx,
+               const std::size_t ny)
+    {
+        return (nx - 1)*ny + nx*(ny - 1);
+    }
+
+    std::vector<int>
+    internalConnections(const std::size_t nx,
+                        const std::size_t ny)
+    {
+        auto cellID = [](const std::size_t start,
+                         const std::size_t off)
+        {
+            return static_cast<int>(start + off);
+        };
+
+        auto neighbours = std::vector<int>{};
+        neighbours.reserve(2 * numIntConn(nx, ny));
+
+        // I connections
+        {
+            for (auto j = 0*ny; j < ny; ++j) {
+                const auto start = j * nx;
+
+                for (auto i = 0*nx + 1; i < nx; ++i) {
+                    neighbours.push_back(cellID(start, i - 1));
+                    neighbours.push_back(cellID(start, i - 0));
+                }
+            }
+        }
+
+        // J connections
+        {
+            for (auto j = 0*ny + 1; j < ny; ++j) {
+                const auto start = (j - 1)*nx;
+
+                for (auto i = 0*nx; i < nx; ++i) {
+                    neighbours.push_back(cellID(start, i + 0 ));
+                    neighbours.push_back(cellID(start, i + nx));
+                }
+            }
+        }
+
+        return neighbours;
+    }
+
+    std::vector<double>
+    flowField(const std::vector<double>::size_type n)
+    {
+        static Opm::RandomVector genRandom{};
+
+        return genRandom.normal(n);
+    }
+
+} // Namespace anonymous
+
+class Setup
+{
+public:
+    Setup(const std::size_t nx,
+          const std::size_t ny);
+
+    const ConnectivityGraph& connectivity() const;
+    const std::vector<double>&    poreVolume()   const;
+    const ConnectionValues&  flux()         const;
+
+private:
+    ConnectivityGraph g_;
+    std::vector<double>    pvol_;
+    ConnectionValues  flux_;
+};
+
+Setup::Setup(const std::size_t nx,
+             const std::size_t ny)
+    : g_   (nx * ny, internalConnections(nx, ny))
+    , pvol_(g_.numCells(), 0.3)
+    , flux_(ConnectionValues::NumConnections{ g_.numConnections() },
+            ConnectionValues::NumPhases     { 1 })
+{
+    const auto flux = flowField(g_.numConnections());
+
+    using ConnID = ConnectionValues::ConnID;
+
+    const auto phaseID =
+        ConnectionValues::PhaseID{ 0 };
+
+    for (decltype(flux_.numConnections())
+             conn = 0, nconn = flux_.numConnections();
+         conn < nconn; ++conn)
+    {
+        flux_(ConnID{conn}, phaseID) = flux[conn];
+    }
+}
+
+const ConnectivityGraph&
+Setup::connectivity() const
+{
+    return g_;
+}
+
+const std::vector<double>&
+Setup::poreVolume() const
+{
+    return pvol_;
+}
+
+const ConnectionValues&
+Setup::flux() const
+{
+    return flux_;
+}
+
+BOOST_AUTO_TEST_SUITE(FlowDiagnostics_Toolbox)
+
+BOOST_AUTO_TEST_CASE (Constructor)
+{
+    const auto cas = Setup(2, 2);
+
+    Toolbox diagTool(cas.connectivity());
+
+    diagTool.assignPoreVolume(cas.poreVolume());
+    diagTool.assignConnectionFlux(cas.flux());
+}
+
+BOOST_AUTO_TEST_CASE (InjectionDiagnostics)
+{
+    const auto cas = Setup(2, 2);
+
+    Toolbox diagTool(cas.connectivity());
+
+    diagTool.assignPoreVolume(cas.poreVolume());
+    diagTool.assignConnectionFlux(cas.flux());
+
+    auto start = std::vector<CellSet>{};
+    {
+        start.emplace_back();
+
+        auto& s = start.back();
+
+        s.identify(CellSetID("I-1"));
+        s.insert(0);
+    }
+
+    {
+        start.emplace_back();
+
+        auto& s = start.back();
+
+        s.identify(CellSetID("I-2"));
+        s.insert(cas.connectivity().numCells() - 1);
+    }
+
+    const auto fwd = diagTool
+        .computeInjectionDiagnostics(start);
+
+    // Global ToF field (accumulated from all injectors)
+    {
+        const auto tof = fwd.fd.timeOfFlight();
+
+        BOOST_CHECK_EQUAL(tof.size(), cas.connectivity().numCells());
+    }
+
+    // Verify set of start points.
+    {
+        const auto startpts = fwd.fd.startPoints();
+
+        BOOST_CHECK_EQUAL(startpts.size(), start.size());
+
+        for (const auto& pt : startpts) {
+            auto pos =
+                std::find_if(start.begin(), start.end(),
+                    [&pt](const CellSet& s)
+                    {
+                        return s.id().to_string() == pt.to_string();
+                    });
+
+            // ID of 'pt' *MUST* be in set of identified start points.
+            BOOST_CHECK(pos != start.end());
+        }
+    }
+
+    // Tracer-ToF
+    {
+        const auto tof = fwd.fd
+            .timeOfFlight(CellSetID("I-1"));
+
+        for (decltype(tof.cellValueCount())
+                 i = 0, n = tof.cellValueCount();
+             i < n; ++i)
+        {
+            const auto v = tof.cellValue(i);
+
+            BOOST_TEST_MESSAGE("[" << i << "] -> ToF["
+                               << v.first << "] = "
+                               << v.second);
+        }
+    }
+
+    // Tracer Concentration
+    {
+        const auto conc = fwd.fd
+            .concentration(CellSetID("I-2"));
+
+        BOOST_TEST_MESSAGE("conc.cellValueCount() = " <<
+                           conc.cellValueCount());
+
+        for (decltype(conc.cellValueCount())
+                 i = 0, n = conc.cellValueCount();
+             i < n; ++i)
+        {
+            const auto v = conc.cellValue(i);
+
+            BOOST_TEST_MESSAGE("[" << i << "] -> Conc["
+                               << v.first << "] = "
+                               << v.second);
+        }
+    }
+}
+
+
+
+
+namespace {
+
+    template <class Collection1, class Collection2>
+    void check_is_close(const Collection1& c1, const Collection2& c2)
+    {
+        BOOST_REQUIRE_EQUAL(c1.size(), c2.size());
+
+        if (! c1.empty()) {
+            auto i1 = c1.begin(), e1 = c1.end();
+            auto i2 = c2.begin();
+
+            for (; i1 != e1; ++i1, ++i2) {
+                BOOST_CHECK_CLOSE(*i1, *i2, 1.0e-10);
+            }
+        }
+    }
+
+} // Namespace Anonymous
+
+
+
+
+
+
+BOOST_AUTO_TEST_CASE (OneDimCase)
+{
+    using namespace Opm::FlowDiagnostics;
+
+    const auto cas = Setup(5, 1);
+    const auto& graph = cas.connectivity();
+
+    // Create fluxes.
+    ConnectionValues flux(ConnectionValues::NumConnections{ graph.numConnections() },
+                          ConnectionValues::NumPhases     { 1 });
+    const size_t nconn = cas.connectivity().numConnections();
+    for (size_t conn = 0; conn < nconn; ++conn) {
+        flux(ConnectionValues::ConnID{conn}, ConnectionValues::PhaseID{0}) = 0.3;
+    }
+
+    // Create well in/out flows.
+    CellSetValues wellflow;
+    wellflow.addCellValue(0, 0.3);
+    wellflow.addCellValue(4, -0.3);
+
+    Toolbox diagTool(graph);
+    diagTool.assignPoreVolume(cas.poreVolume());
+    diagTool.assignConnectionFlux(flux);
+    diagTool.assignInflowFlux(wellflow);
+
+    auto start = std::vector<CellSet>{};
+    {
+        start.emplace_back();
+
+        auto& s = start.back();
+
+        s.identify(CellSetID("I-1"));
+        s.insert(0);
+    }
+
+    {
+        start.emplace_back();
+
+        auto& s = start.back();
+
+        s.identify(CellSetID("I-2"));
+        s.insert(cas.connectivity().numCells() - 1);
+    }
+
+    const auto fwd = diagTool.computeInjectionDiagnostics(start);
+    const auto rev = diagTool.computeProductionDiagnostics(start);
+
+    // Global ToF field (accumulated from all injectors)
+    {
+        const auto tof = fwd.fd.timeOfFlight();
+
+        BOOST_REQUIRE_EQUAL(tof.size(), cas.connectivity().numCells());
+        std::vector<double> expected = { 0.5, 1.5, 2.5, 3.5, 0.0 };
+        check_is_close(tof, expected);
+    }
+
+    // Global ToF field (accumulated from all producers)
+    {
+        const auto tof = rev.fd.timeOfFlight();
+
+        BOOST_REQUIRE_EQUAL(tof.size(), cas.connectivity().numCells());
+        std::vector<double> expected = { 0.0, 3.5, 2.5, 1.5, 0.5 };
+        check_is_close(tof, expected);
+    }
+
+    // Verify set of start points.
+    {
+        const auto startpts = fwd.fd.startPoints();
+
+        BOOST_CHECK_EQUAL(startpts.size(), start.size());
+
+        for (const auto& pt : startpts) {
+            auto pos =
+                std::find_if(start.begin(), start.end(),
+                    [&pt](const CellSet& s)
+                    {
+                        return s.id().to_string() == pt.to_string();
+                    });
+
+            // ID of 'pt' *MUST* be in set of identified start points.
+            BOOST_CHECK(pos != start.end());
+        }
+    }
+
+    // Tracer-ToF
+    {
+        const auto tof = fwd.fd
+            .timeOfFlight(CellSetID("I-2"));
+
+        for (decltype(tof.cellValueCount())
+                 i = 0, n = tof.cellValueCount();
+             i < n; ++i)
+        {
+            const auto v = tof.cellValue(i);
+
+            BOOST_TEST_MESSAGE("[" << i << "] -> ToF["
+                               << v.first << "] = "
+                               << v.second);
+        }
+    }
+
+    // Tracer Concentration
+    {
+        const auto conc = fwd.fd
+            .concentration(CellSetID("I-2"));
+
+        BOOST_TEST_MESSAGE("conc.cellValueCount() = " <<
+                           conc.cellValueCount());
+
+        for (decltype(conc.cellValueCount())
+                 i = 0, n = conc.cellValueCount();
+             i < n; ++i)
+        {
+            const auto v = conc.cellValue(i);
+
+            BOOST_TEST_MESSAGE("[" << i << "] -> Conc["
+                               << v.first << "] = "
+                               << v.second);
+        }
+    }
+}
+
+BOOST_AUTO_TEST_SUITE_END()

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/tests/test_tarjan.cpp

@@ -0,0 +1,722 @@
+/*
+  Copyright 2016 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 Statoil ASA.
+
+  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 <config.h>
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE
+
+#define BOOST_TEST_MODULE TarjanImplementationTest
+
+#include <opm/common/utility/platform_dependent/disable_warnings.h>
+#include <boost/test/unit_test.hpp>
+#include <opm/common/utility/platform_dependent/reenable_warnings.h>
+
+#include <opm/utility/graph/tarjan.h>
+
+#include <memory>
+
+namespace {
+    struct DestroyWorkSpace
+    {
+        void operator()(TarjanWorkSpace* ws);
+    };
+
+    struct DestroySCCResult
+    {
+        void operator()(TarjanSCCResult *scc);
+    };
+
+    void DestroyWorkSpace::operator()(TarjanWorkSpace* ws)
+    {
+        destroy_tarjan_workspace(ws);
+    }
+
+    void DestroySCCResult::operator()(TarjanSCCResult* scc)
+    {
+        destroy_tarjan_sccresult(scc);
+    }
+
+    using WorkSpace =
+        std::unique_ptr<TarjanWorkSpace, DestroyWorkSpace>;
+
+    using SCCResult =
+        std::unique_ptr<TarjanSCCResult, DestroySCCResult>;
+
+    void check_scc(const std::size_t*     expect_size,
+                   const int*             expect_vert,
+                   const TarjanSCCResult* scc)
+    {
+        const auto ncomp = tarjan_get_numcomponents(scc);
+
+        for (auto comp = 0*ncomp, k = 0*ncomp; comp < ncomp; ++comp) {
+            const auto c = tarjan_get_strongcomponent(scc, comp);
+
+            BOOST_CHECK_EQUAL(c.size, expect_size[comp]);
+
+            for (auto i = 0*c.size; i < c.size; ++i, ++k) {
+                BOOST_CHECK_EQUAL(c.vertex[i], expect_vert[k]);
+            }
+        }
+    }
+} // Anonymous
+
+BOOST_AUTO_TEST_SUITE(Two_By_Two)
+
+// +-----+-----+
+// |  2  |  3  |
+// +-----+-----+
+// |  0  |  1  |
+// +-----+-----+
+
+BOOST_AUTO_TEST_CASE (FullySeparable)
+{
+    // Quarter five-spot pattern:
+    //   0 -> 1
+    //   0 -> 2
+    //   1 -> 3
+    //   2 -> 3
+    //
+    // Note: (ia,ja) is INFLOW graph whence tarjan() returns SCCs in
+    // topological order from sources to sinks.
+    const std::size_t ia[] = { 0, 0, 1, 2, 4 };
+    const int         ja[] = { 0, 0, 1, 2 };
+
+    const std::size_t nv           = (sizeof ia) / (sizeof ia[0]) - 1;
+    const std::size_t expect_ncomp = 4;
+
+    auto scc = SCCResult{ tarjan(nv, ia, ja) };
+
+    const auto ncomp = tarjan_get_numcomponents(scc.get());
+
+    BOOST_CHECK_EQUAL(ncomp, expect_ncomp);
+
+    {
+        const std::size_t expect_size[] = { 1, 1, 1, 1 };
+        const int         expect_vert[] = { 0, 1, 2, 3 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+
+    if (tarjan_reverse_sccresult(scc.get())) {
+        const std::size_t expect_size[] = { 1, 1, 1, 1 };
+        const int         expect_vert[] = { 3, 2, 1, 0 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (FullySeparableSparse)
+{
+    // Quarter five-spot pattern:
+    //   0 -> 1
+    //   0 -> 2
+    //   1 -> 3
+    //   2 -> 3
+    //
+    // Note: (ia,ja) is OUTFLOW graph.  We use tarjan_reverse_sccresult() to
+    // access the SCCs in topological order.
+    const std::size_t ia[] = { 0, 2, 3, 4, 4 };
+    const int         ja[] = { 1, 2, 3, 3 };
+
+    const int         start_pts[] = { 2 };
+    const std::size_t nv          = (sizeof ia) / (sizeof ia[0]) - 1;
+    const std::size_t nstart      = (sizeof start_pts) / (sizeof start_pts[0]);
+
+    auto work = WorkSpace{ create_tarjan_workspace(nv) };
+
+    auto scc = SCCResult{
+        tarjan_reachable_sccs(nv, ia, ja, nstart,
+                              start_pts, work.get())
+    };
+
+    BOOST_CHECK_EQUAL(tarjan_get_numcomponents(scc.get()), 2);
+
+    {
+        const std::size_t expect_size[] = { 1, 1 };
+        const int         expect_vert[] = { 3, 2 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+
+    if (tarjan_reverse_sccresult(scc.get())) {
+        const std::size_t expect_size[] = { 1, 1 };
+        const int         expect_vert[] = { 2, 3 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (Loop)
+{
+    // Circulation:
+    //   0 -> 1
+    //   1 -> 3
+    //   3 -> 2
+    //   2 -> 0
+    //
+    // Note: (ia,ja) is INFLOW graph whence tarjan() returns SCCs in
+    // topological order from sources to sinks.
+    const std::size_t ia[] = { 0, 1, 2, 3, 4 };
+    const int         ja[] = { 2, 0, 3, 1 };
+
+    const std::size_t nv           = (sizeof ia) / (sizeof ia[0]) - 1;
+    const std::size_t expect_ncomp = 1;
+
+    auto scc = SCCResult{ tarjan(nv, ia, ja) };
+
+    const auto ncomp = tarjan_get_numcomponents(scc.get());
+
+    BOOST_CHECK_EQUAL(ncomp, expect_ncomp);
+
+    // Cell indices within component returned in (essentially) arbitrary
+    // order.  This particular order happened to be correct at the time the
+    // test was implemented so the assertion on 'vertex' is only usable as a
+    // regression test.
+    {
+        const std::size_t expect_size[] = { nv };
+        const int         expect_vert[] = { 1, 3, 2, 0 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+
+    if (tarjan_reverse_sccresult(scc.get())) {
+        // Reversing order of components maintains internal order of
+        // vertices within each component.
+
+        const std::size_t expect_size[] = { nv };
+        const int         expect_vert[] = { 1, 3, 2, 0 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (LoopSparse)
+{
+    // Circulation:
+    //   0 -> 1
+    //   1 -> 3
+    //   3 -> 2
+    //   2 -> 0
+    //
+    // Note: (ia,ja) is OUTFLOW graph.  We use tarjan_reverse_sccresult() to
+    // access the SCCs in topological order.
+    const std::size_t ia[] = { 0, 1, 2, 3, 4 };
+    const int         ja[] = { 1, 3, 0, 2 };
+
+    const int         start_pts[] = { 2 };
+    const std::size_t nv          = (sizeof ia) / (sizeof ia[0]) - 1;
+    const std::size_t nstart      = (sizeof start_pts) / (sizeof start_pts[0]);
+
+    const std::size_t expect_ncomp = 1;
+
+    auto work = WorkSpace{ create_tarjan_workspace(nv) };
+
+    auto scc = SCCResult{
+        tarjan_reachable_sccs(nv, ia, ja, nstart,
+                              start_pts, work.get())
+    };
+
+    BOOST_CHECK_EQUAL(tarjan_get_numcomponents(scc.get()),
+                      expect_ncomp);
+
+    {
+        const std::size_t expect_size[] = { nv };
+        const int         expect_vert[] = { 3, 1, 0, 2 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+
+    if (tarjan_reverse_sccresult(scc.get())) {
+        // Reversing order of components maintains internal order of
+        // vertices within each component.
+
+        const std::size_t expect_size[] = { nv };
+        const int         expect_vert[] = { 3, 1, 0, 2 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (DualPath)
+{
+    // Two flow paths from cell 0 to cell 2:
+    //   0 -> 1
+    //   1 -> 3
+    //   3 -> 2
+    //   0 -> 2
+    //
+    // Note: (ia,ja) is INFLOW graph whence tarjan() returns SCCs in
+    // topological order from sources to sinks.
+    const std::size_t ia[] = { 0, 0, 2, 3, 4 };
+    const int         ja[] = { 0, 0, 3, 1 };
+
+    const std::size_t nv           = (sizeof ia) / (sizeof ia[0]) - 1;
+    const std::size_t expect_ncomp = 4;
+
+    auto scc = SCCResult{ tarjan(nv, ia, ja) };
+
+    const auto ncomp = tarjan_get_numcomponents(scc.get());
+
+    BOOST_CHECK_EQUAL(ncomp, expect_ncomp);
+
+    // Cell 0 is a source and cell 2 is a sink so first and last cells must
+    // be 0 and 2 respectively.  The order of cells 1 and 3 is determined by
+    // the flow path in which 1 precedes 3.
+
+    {
+        const std::size_t expect_size[] = { 1, 1, 1, 1 };
+        const int         expect_vert[] = { 0, 1, 3, 2 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+
+    if (tarjan_reverse_sccresult(scc.get())) {
+        const std::size_t expect_size[] = { 1, 1, 1, 1 };
+        const int         expect_vert[] = { 2, 3, 1, 0 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (DualPathSparse)
+{
+    // Two flow paths from cell 0 to cell 2:
+    //   0 -> 1
+    //   1 -> 3
+    //   3 -> 2
+    //   0 -> 2
+    //
+    // Note: (ia,ja) is OUTFLOW graph.  We use tarjan_reverse_sccresult() to
+    // access the SCCs in topological order.
+    const std::size_t ia[] = { 0, 2, 3, 3, 4 };
+    const int         ja[] = { 1, 2, 3, 2 };
+
+    const int         start_pts[] = { 1 };
+    const std::size_t nv          = (sizeof ia) / (sizeof ia[0]) - 1;
+    const std::size_t nstart      = (sizeof start_pts) / (sizeof start_pts[0]);
+
+    const std::size_t expect_ncomp = 3;
+
+    auto work = WorkSpace{ create_tarjan_workspace(nv) };
+
+    auto scc = SCCResult{
+        tarjan_reachable_sccs(nv, ia, ja, nstart,
+                              start_pts, work.get())
+    };
+
+    BOOST_CHECK_EQUAL(tarjan_get_numcomponents(scc.get()),
+                      expect_ncomp);
+
+    {
+        const std::size_t expect_size[] = { 1, 1, 1 };
+        const int         expect_vert[] = { 2, 3, 1 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+
+    if (tarjan_reverse_sccresult(scc.get())) {
+        const std::size_t expect_size[] = { 1, 1, 1 };
+        const int         expect_vert[] = { 1, 3, 2 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (IsolatedFlows)
+{
+    // Compartmentalised reservoir with source and sink in each compartment.
+    //   0 -> 2
+    //   3 -> 1
+    //
+    // Note: (ia,ja) is INFLOW graph whence tarjan() returns SCCs in
+    // topological order from sources to sinks.
+    const std::size_t ia[] = { 0, 0, 1, 2, 2 };
+    const int         ja[] = { 3, 0 };
+
+    const std::size_t nv           = (sizeof ia) / (sizeof ia[0]) - 1;
+    const std::size_t expect_ncomp = 4;
+
+    auto scc = SCCResult{ tarjan(nv, ia, ja) };
+
+    const auto ncomp = tarjan_get_numcomponents(scc.get());
+
+    BOOST_CHECK_EQUAL(ncomp, expect_ncomp);
+
+    // Sources before sinks, but no a priori ordering between sources or
+    // between sinks.  This particular order happened to be correct at the
+    // time the test was implemented so the assertion on 'vert' is only
+    // usable as a regression test.
+
+    {
+        const std::size_t expect_size[] = { 1, 1, 1, 1 };
+        const int         expect_vert[] = { 0, 3, 1, 2 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+
+    if (tarjan_reverse_sccresult(scc.get())) {
+        const std::size_t expect_size[] = { 1, 1, 1, 1 };
+        const int         expect_vert[] = { 2, 1, 3, 0 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (IsolatedFlowsSparse)
+{
+    // Compartmentalised reservoir with source and sink in each compartment.
+    //   0 -> 2
+    //   3 -> 1
+    //
+    // Note: (ia,ja) is OUTFLOW graph.  We use tarjan_reverse_sccresult() to
+    // access the SCCs in topological order.
+    const std::size_t ia[] = { 0, 1, 1, 1, 2 };
+    const int         ja[] = { 2, 1 };
+
+    const int         start_pts[] = { 3 };
+    const std::size_t nv          = (sizeof ia) / (sizeof ia[0]) - 1;
+    const std::size_t nstart      = (sizeof start_pts) / (sizeof start_pts[0]);
+
+    const std::size_t expect_ncomp = 2;
+
+    auto work = WorkSpace{ create_tarjan_workspace(nv) };
+
+    auto scc = SCCResult{
+        tarjan_reachable_sccs(nv, ia, ja, nstart,
+                              start_pts, work.get())
+    };
+
+    BOOST_CHECK_EQUAL(tarjan_get_numcomponents(scc.get()),
+                      expect_ncomp);
+
+    {
+        const std::size_t expect_size[] = { 1, 1 };
+        const int         expect_vert[] = { 1, 3 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+
+    if (tarjan_reverse_sccresult(scc.get())) {
+        const std::size_t expect_size[] = { 1, 1 };
+        const int         expect_vert[] = { 3, 1 };
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+BOOST_AUTO_TEST_SUITE(Four_By_Four)
+
+// +-----+-----+-----+-----+
+// | 12  | 13  | 14  | 15  |
+// +-----+-----+-----+-----+
+// |  8  |  9  | 10  | 11  |
+// +-----+-----+-----+-----+
+// |  4  |  5  |  6  |  7  |
+// +-----+-----+-----+-----+
+// |  0  |  1  |  2  |  3  |
+// +-----+-----+-----+-----+
+
+BOOST_AUTO_TEST_CASE (CentreLoop)
+{
+    //  From -> To
+    //     0 -> [  1, 4 ],
+    //     1 -> [  2, 5 ],
+    //     2 -> [  3, 6 ],
+    //     3 -> [  7 ],
+    //     4 -> [  5, 8 ],
+    //     5 -> [  6 ],
+    //     6 -> [  7, 10 ],
+    //     7 -> [ 11 ],
+    //     8 -> [  9, 12 ],
+    //     9 -> [  5 ],
+    //    10 -> [  9, 11 ],
+    //    11 -> [ 15 ],
+    //    12 -> [ 13 ],
+    //    13 -> [  9, 14 ],
+    //    14 -> [ 10, 15 ],
+    //    15 -> Void (sink cell)
+    //
+    // Note: (ia,ja) is OUTFLOW graph.  We use tarjan_reverse_sccresult() to
+    // access the SCCs in topological order.
+
+    const std::size_t ia[] = {0, 2, 4, 6, 7, 9, 10, 12,
+                              13, 15, 16, 18, 19, 20, 22, 24, 24};
+    const int         ja[] = {1, 4, 2, 5, 3, 6, 7, 5, 8, 6, 7, 10, 11,
+                              9, 12, 5, 9, 11, 15, 13, 9, 14, 10, 15};
+
+    const std::size_t nv           = (sizeof ia) / (sizeof ia[0]) - 1;
+    const std::size_t expect_ncomp = 13;
+
+    auto scc = SCCResult{ tarjan(nv, ia, ja) };
+
+    const auto ncomp = tarjan_get_numcomponents(scc.get());
+
+    BOOST_CHECK_EQUAL(ncomp, expect_ncomp);
+
+    {
+        const std::size_t expect_size[] =
+            { 1, 1, 1, 4,      // 0 ..  3
+              1, 1, 1, 1,      // 4 ..  7
+              1, 1, 1, 1, 1 }; // 8 .. 12
+
+        const int expect_vert[] =
+            { 15, 11, 7,        // 0 ..  2
+              6, 5, 9, 10,      // 3
+              14, 13, 12, 8,    // 4 ..  7
+              4, 3, 2, 1, 0 };  // 8 .. 12
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+
+    if (tarjan_reverse_sccresult(scc.get())) {
+        const std::size_t expect_size[] =
+            { 1, 1, 1, 1, 1,    // 0 ..  4
+              1, 1, 1, 1,       // 5 ..  8
+              4, 1, 1, 1 };     // 9 .. 12
+
+        const int expect_vert[] =
+            { 0, 1, 2, 3, 4, //  0 ..  4
+              8, 12, 13, 14, //  5 ..  8
+              6, 5, 9, 10,   //  9
+              7, 11, 15 };   // 10 .. 12
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (CentreLoopSparse)
+{
+    //  From -> To
+    //     0 -> [  1, 4 ],
+    //     1 -> [  2, 5 ],
+    //     2 -> [  3, 6 ],
+    //     3 -> [  7 ],
+    //     4 -> [  5, 8 ],
+    //     5 -> [  6 ],
+    //     6 -> [  7, 10 ],
+    //     7 -> [ 11 ],
+    //     8 -> [  9, 12 ],
+    //     9 -> [  5 ],
+    //    10 -> [  9, 11 ],
+    //    11 -> [ 15 ],
+    //    12 -> [ 13 ],
+    //    13 -> [  9, 14 ],
+    //    14 -> [ 10, 15 ],
+    //    15 -> Void (sink cell)
+    //
+    // Note: (ia,ja) is OUTFLOW graph.  We use tarjan_reverse_sccresult() to
+    // access the SCCs in topological order.
+
+    const std::size_t ia[] = {0, 2, 4, 6, 7, 9, 10, 12,
+                              13, 15, 16, 18, 19, 20, 22, 24, 24};
+    const int         ja[] = {1, 4, 2, 5, 3, 6, 7, 5, 8, 6, 7, 10, 11,
+                              9, 12, 5, 9, 11, 15, 13, 9, 14, 10, 15};
+
+    const int         start_pts[] = { 5, 12 };
+    const std::size_t nv          = (sizeof ia) / (sizeof ia[0]) - 1;
+    const std::size_t nstart      = (sizeof start_pts) / (sizeof start_pts[0]);
+
+    const std::size_t expect_ncomp = 7;
+
+    auto work = WorkSpace{ create_tarjan_workspace(nv) };
+
+    auto scc = SCCResult{
+        tarjan_reachable_sccs(nv, ia, ja, nstart,
+                              start_pts, work.get())
+    };
+
+    BOOST_CHECK_EQUAL(tarjan_get_numcomponents(scc.get()),
+                      expect_ncomp);
+
+    {
+        const std::size_t expect_size[] =
+            { 1, 1, 1, 4,       // 0 .. 3
+              1, 1, 1 };        // 4 .. 6
+
+        const int expect_vert[] =
+            { 15, 11, 7,        // 0 .. 2
+              9, 10, 6, 5,      // 3
+              14, 13, 12 };     // 4 .. 6
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+
+    if (tarjan_reverse_sccresult(scc.get())) {
+        // Order of vertices preserved within strong component.
+
+        const std::size_t expect_size[] =
+            { 1, 1, 1, 4,       // 0 .. 3
+              1, 1, 1 };        // 4 .. 6
+
+        const int expect_vert[] =
+            { 12, 13, 14,       // 0 .. 2
+              9, 10, 6, 5,      // 3
+              7, 11, 15 };      // 4 .. 6
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (CentreLoopSource12)
+{
+    //  From -> To
+    //     0 -> [  1 ],
+    //     1 -> [  2 ],
+    //     2 -> [  3, 6 ],
+    //     3 -> [  7 ],
+    //     4 -> [  0, 5 ],
+    //     5 -> [  6 ],
+    //     6 -> [  7, 10 ],
+    //     7 -> [ 11 ],
+    //     8 -> [  4, 9 ],
+    //     9 -> [  5 ],
+    //    10 -> [  9, 11 ],
+    //    11 -> [ 15 ],
+    //    12 -> [  8, 13 ],
+    //    13 -> [ 14 ],
+    //    14 -> [ 10, 15 ],
+    //    15 -> Void (sink cell)
+
+    const std::size_t ia[] = {0, 1, 2, 4, 5, 7, 8, 10, 11,
+                              13, 14, 16, 17, 19, 20, 22, 22};
+
+    const int         ja[] = {1, 2, 3, 6, 7, 0, 5, 6, 7, 10, 11,
+                              4, 9, 5, 9, 11, 15, 8, 13, 14, 10, 15};
+
+    const std::size_t nv           = (sizeof ia) / (sizeof ia[0]) - 1;
+    const std::size_t expect_ncomp = 13;
+
+    auto scc = SCCResult{ tarjan(nv, ia, ja) };
+
+    const auto ncomp = tarjan_get_numcomponents(scc.get());
+
+    BOOST_CHECK_EQUAL(ncomp, expect_ncomp);
+
+    {
+        const std::size_t expect_size[] =
+            { 1, 1, 1, 4,      // 0 ..  3
+              1, 1, 1, 1,      // 4 ..  7
+              1, 1, 1, 1, 1 }; // 8 .. 12
+
+        const int expect_vert[] =
+            { 15, 11, 7,           // 0 ..  2
+              5, 9, 10, 6,         // 3
+              3, 2, 1, 0,          // 4 ..  7
+              4, 8, 14, 13, 12 };  // 8 .. 12
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+
+    if (tarjan_reverse_sccresult(scc.get())) {
+        const std::size_t expect_size[] =
+            { 1, 1, 1, 1, 1,    // 0 ..  4
+              1, 1, 1, 1,       // 5 ..  8
+              4, 1, 1, 1 };     // 9 .. 12
+
+        const int expect_vert[] =
+            { 12, 13, 14, 8, 4, //  0 ..  4
+              0, 1, 2, 3,       //  5 ..  8
+              5, 9, 10, 6,      //  9
+              7, 11, 15 };      // 10 .. 12
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+}
+
+BOOST_AUTO_TEST_CASE (CentreLoopSource12Sparse)
+{
+    //  From -> To
+    //     0 -> [  1 ],
+    //     1 -> [  2 ],
+    //     2 -> [  3, 6 ],
+    //     3 -> [  7 ],
+    //     4 -> [  0, 5 ],
+    //     5 -> [  6 ],
+    //     6 -> [  7, 10 ],
+    //     7 -> [ 11 ],
+    //     8 -> [  4, 9 ],
+    //     9 -> [  5 ],
+    //    10 -> [  9, 11 ],
+    //    11 -> [ 15 ],
+    //    12 -> [  8, 13 ],
+    //    13 -> [ 14 ],
+    //    14 -> [ 10, 15 ],
+    //    15 -> Void (sink cell)
+    //
+    // Note: (ia,ja) is OUTFLOW graph.  We use tarjan_reverse_sccresult() to
+    // access the SCCs in topological order.
+
+    const std::size_t ia[] = {0, 1, 2, 4, 5, 7, 8, 10, 11,
+                              13, 14, 16, 17, 19, 20, 22, 22};
+
+    const int         ja[] = {1, 2, 3, 6, 7, 0, 5, 6, 7, 10, 11,
+                              4, 9, 5, 9, 11, 15, 8, 13, 14, 10, 15};
+
+    const int         start_pts[] = { 7, 14 };
+    const std::size_t nv          = (sizeof ia) / (sizeof ia[0]) - 1;
+    const std::size_t nstart      = (sizeof start_pts) / (sizeof start_pts[0]);
+
+    const std::size_t expect_ncomp = 5;
+
+    auto work = WorkSpace{ create_tarjan_workspace(nv) };
+
+    auto scc = SCCResult{
+        tarjan_reachable_sccs(nv, ia, ja, nstart,
+                              start_pts, work.get())
+    };
+
+    BOOST_CHECK_EQUAL(tarjan_get_numcomponents(scc.get()),
+                      expect_ncomp);
+
+    {
+        const std::size_t expect_size[] =
+            { 1, 1, 1, 4, 1 };  // 0 .. 4
+
+        const int expect_vert[] =
+            { 15, 11, 7,           // 0 .. 2
+              6, 5, 9, 10,         // 3
+              14 };                // 4
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+
+    if (tarjan_reverse_sccresult(scc.get())) {
+        // Order of vertices preserved within each strong component.
+
+        const std::size_t expect_size[] =
+            { 1, 4, 1, 1, 1 };  // 0 .. 4
+
+        const int expect_vert[] =
+            { 14,               // 0
+              6, 5, 9, 10,      // 1
+              7, 11, 15 };      // 2 .. 4
+
+        check_scc(expect_size, expect_vert, scc.get());
+    }
+}
+
+BOOST_AUTO_TEST_SUITE_END()

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/travis/build-and-test-opm-flowdiagnostics.sh

@@ -0,0 +1,7 @@
+#!/bin/sh
+set -e
+
+sh ./opm-flowdiagnostics/travis/build-opm-flowdiagnostics.sh
+
+(cd opm-flowdiagnostics/build
+ ctest --output-on-failure)

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/travis/build-opm-common-shared.sh

@@ -0,0 +1,10 @@
+#!/bin/sh
+
+set -e
+
+mkdir -p opm-flowdiagnostics/build
+
+(cd opm-flowdiagnostics/build
+ cmake ../ -DBUILD_SHARED_LIBS=ON
+
+ make)

ThirdParty/custom-opm-flowdiagnostics/opm-flowdiagnostics/travis/build-opm-flowdiagnostics.sh

@@ -0,0 +1,10 @@
+#!/bin/sh
+
+set -e
+
+mkdir -p opm-flowdiagnostics/build
+
+(cd opm-flowdiagnostics/build
+ cmake ../
+
+ make)