diff --git a/NumLib/NamedFunction.cpp b/NumLib/NamedFunction.cpp
deleted file mode 100644
index 2ce999199763a87509d242f6a63a12ea2e4c0a46..0000000000000000000000000000000000000000
--- a/NumLib/NamedFunction.cpp
+++ /dev/null
@@ -1,199 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "NamedFunction.h"
-#include "BaseLib/Functional.h"
-
-//! Helper struct used in conjunction with std::integer_sequence<int, ...> to
-//! get a sequence of types, where each type is double.
-template <int>
-struct Double
-{
-    using type = double;
-};
-
-/*! Calls the given \c function with the given  \c arguments.
- *
- * \tparam Indices sequence of integers used to expand the \c arguments vector
- * to individual arguments.
- */
-template <int... Indices>
-double call_(void* function, std::vector<double> const& arguments)
-{
-    assert(arguments.size() == sizeof...(Indices));
-    auto* fct = reinterpret_cast<
-        std::function<double(typename Double<Indices>::type...)>*>(function);
-    auto* args = arguments.data();
-    return (*fct)(args[Indices]...);
-}
-
-using CallerFunction = double (*)(void*, const std::vector<double>&);
-
-//! Helps instantiating the call_() function.
-template <int... Indices>
-CallerFunction generateCallerFromIntegerSequence(
-    std::integer_sequence<int, Indices...> /*unused*/)
-{
-    return call_<Indices...>;
-}
-
-//! Instantiates the call_() function for the provided number of arguments.
-template <int NumArguments>
-CallerFunction generateCaller()
-{
-    return generateCallerFromIntegerSequence(
-        std::make_integer_sequence<int, NumArguments>{});
-}
-
-//! Holds instantiations of the call_() function for various numbers of
-//! arguments.
-static const CallerFunction callers[] = {
-    generateCaller<0>(),  generateCaller<1>(),  generateCaller<2>(),
-    generateCaller<3>(),  generateCaller<4>(),  generateCaller<5>(),
-    generateCaller<6>(),  generateCaller<7>(),  generateCaller<8>(),
-    generateCaller<9>(),  generateCaller<10>(), generateCaller<11>(),
-    generateCaller<12>(), generateCaller<13>(), generateCaller<14>(),
-    generateCaller<15>(), generateCaller<16>(), generateCaller<17>(),
-    generateCaller<18>(), generateCaller<19>(), generateCaller<20>(),
-    generateCaller<21>(), generateCaller<22>(), generateCaller<23>(),
-    generateCaller<24>(), generateCaller<25>(), generateCaller<26>(),
-    generateCaller<27>(), generateCaller<28>(), generateCaller<29>(),
-    generateCaller<30>(), generateCaller<31>(), generateCaller<32>()};
-static_assert(sizeof(callers) / sizeof(CallerFunction) ==
-                  NumLib::NamedFunction::MAX_FUNCTION_ARGS + 1,
-              "You did not instantiate the right number of callers.");
-
-/*! Deletes the given \c function.
- *
- * \tparam Indices sequence of integers used to cast the \c function to the
- * correct type.
- */
-template <int... Indices>
-void delete_(void* function)
-{
-    auto* fct = reinterpret_cast<
-        std::function<double(typename Double<Indices>::type...)>*>(function);
-    delete fct;
-}
-
-using DeleterFunction = void (*)(void*);
-
-//! Helps instantiating the delete_() function.
-template <int... Indices>
-DeleterFunction generateDeleterFromIntegerSequence(
-    std::integer_sequence<int, Indices...> /*unused*/)
-{
-    return delete_<Indices...>;
-}
-
-//! Instantiates the delete_() function for the provided number of arguments.
-template <int NumArguments>
-DeleterFunction generateDeleter()
-{
-    return generateDeleterFromIntegerSequence(
-        std::make_integer_sequence<int, NumArguments>{});
-}
-
-//! Holds instantiations of the delete_() function for various numbers of
-//! arguments.
-static const DeleterFunction deleters[] = {
-    generateDeleter<0>(),  generateDeleter<1>(),  generateDeleter<2>(),
-    generateDeleter<3>(),  generateDeleter<4>(),  generateDeleter<5>(),
-    generateDeleter<6>(),  generateDeleter<7>(),  generateDeleter<8>(),
-    generateDeleter<9>(),  generateDeleter<10>(), generateDeleter<11>(),
-    generateDeleter<12>(), generateDeleter<13>(), generateDeleter<14>(),
-    generateDeleter<15>(), generateDeleter<16>(), generateDeleter<17>(),
-    generateDeleter<18>(), generateDeleter<19>(), generateDeleter<20>(),
-    generateDeleter<21>(), generateDeleter<22>(), generateDeleter<23>(),
-    generateDeleter<24>(), generateDeleter<25>(), generateDeleter<26>(),
-    generateDeleter<27>(), generateDeleter<28>(), generateDeleter<29>(),
-    generateDeleter<30>(), generateDeleter<31>(), generateDeleter<32>()};
-static_assert(sizeof(deleters) / sizeof(DeleterFunction) ==
-                  NumLib::NamedFunction::MAX_FUNCTION_ARGS + 1,
-              "You did not instantiate the right number of deleters.");
-
-/*! Copies the given \c function.
- *
- * \tparam Indices sequence of integers used to cast the \c function to the
- * correct type.
- */
-template <int... Indices>
-void* copy_(void* function)
-{
-    auto* fct = reinterpret_cast<
-        std::function<double(typename Double<Indices>::type...)>*>(function);
-    return new std::function<double(typename Double<Indices>::type...)>(*fct);
-}
-
-using CopierFunction = void* (*)(void*);
-
-//! Helps instantiating the copy_() function.
-template <int... Indices>
-CopierFunction generateCopierFromIntegerSequence(
-    std::integer_sequence<int, Indices...> /*unused*/)
-{
-    return copy_<Indices...>;
-}
-
-//! Instantiates the copy_() function for the provided number of arguments.
-template <int NumArguments>
-CopierFunction generateCopier()
-{
-    return generateCopierFromIntegerSequence(
-        std::make_integer_sequence<int, NumArguments>{});
-}
-
-//! Holds instantiations of the copy_() function for various numbers of
-//! arguments.
-static const CopierFunction copiers[] = {
-    generateCopier<0>(),  generateCopier<1>(),  generateCopier<2>(),
-    generateCopier<3>(),  generateCopier<4>(),  generateCopier<5>(),
-    generateCopier<6>(),  generateCopier<7>(),  generateCopier<8>(),
-    generateCopier<9>(),  generateCopier<10>(), generateCopier<11>(),
-    generateCopier<12>(), generateCopier<13>(), generateCopier<14>(),
-    generateCopier<15>(), generateCopier<16>(), generateCopier<17>(),
-    generateCopier<18>(), generateCopier<19>(), generateCopier<20>(),
-    generateCopier<21>(), generateCopier<22>(), generateCopier<23>(),
-    generateCopier<24>(), generateCopier<25>(), generateCopier<26>(),
-    generateCopier<27>(), generateCopier<28>(), generateCopier<29>(),
-    generateCopier<30>(), generateCopier<31>(), generateCopier<32>()};
-static_assert(sizeof(copiers) / sizeof(CopierFunction) ==
-                  NumLib::NamedFunction::MAX_FUNCTION_ARGS + 1,
-              "You did not instantiate the right number of deleters.");
-
-namespace NumLib
-{
-NamedFunction::NamedFunction(NamedFunction&& other)
-    : _name(std::move(other._name)),
-      _argument_names(std::move(other._argument_names)),
-      _function(other._function)
-{
-    other._function = nullptr;
-}
-
-NamedFunction::NamedFunction(NamedFunction const& other)
-    : _name(other._name),
-      _argument_names(other._argument_names),
-      _function(copiers[_argument_names.size()](other._function))
-{
-}
-
-NamedFunction::~NamedFunction()
-{
-    deleters[_argument_names.size()](_function);
-}
-
-double NamedFunction::call(const std::vector<double>& arguments) const
-{
-    assert(arguments.size() == _argument_names.size());
-    return callers[_argument_names.size()](_function, arguments);
-}
-
-}  // namespace NumLib
diff --git a/NumLib/NamedFunction.h b/NumLib/NamedFunction.h
deleted file mode 100644
index 11e2851346957a6760e92278d51287f6fbfe6c86..0000000000000000000000000000000000000000
--- a/NumLib/NamedFunction.h
+++ /dev/null
@@ -1,115 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include <cassert>
-#include <functional>
-#include <memory>
-#include <string>
-#include <type_traits>
-#include <utility>
-#include <vector>
-
-namespace detail
-{
-//! Checks if all types \c Ts are the same as \c TRef.
-template<typename TRef, typename... Ts>
-struct AllTypesSameAs;
-
-template<typename TRef, typename T1, typename... Ts>
-struct AllTypesSameAs<TRef, T1, Ts...>
-{
-    static const bool value = std::is_same<TRef, T1>::value
-        && AllTypesSameAs<TRef, Ts...>::value;
-};
-
-template<typename TRef>
-struct AllTypesSameAs<TRef>
-{
-    static const bool value = true;
-};
-
-template<typename TRef, typename T1, typename... Ts>
-const bool AllTypesSameAs<TRef, T1, Ts...>::value;
-
-template<typename TRef>
-const bool AllTypesSameAs<TRef>::value;
-
-} // namespace detail
-
-namespace NumLib
-{
-//! Stores a function object along with a name for it and information about its
-//! arguments.
-class NamedFunction final
-{
-public:
-    /*! Constructs a new named function.
-     *
-     * \param name the function's name
-     * \param argument_names names  of arguments of the function
-     * \param function the actual function object
-     */
-    template <typename ReturnType, typename... Arguments>
-    NamedFunction(std::string name,
-                  std::vector<std::string>&& argument_names,
-                  std::function<ReturnType(Arguments...)>&& function);
-
-    NamedFunction(NamedFunction&& other);
-    NamedFunction(NamedFunction const& other);
-
-    ~NamedFunction();
-
-    //! Returns the function's name.
-    std::string const& getName() const { return _name; }
-    //! Returns the names of the function's arguments.
-    std::vector<std::string> const& getArgumentNames() const
-    {
-        return _argument_names;
-    }
-
-    //! Call the function with the supplied arguments.
-    double call(std::vector<double> const& arguments) const;
-
-    //! Maximum number of function arguments supported by NamedFunction.
-    static const int MAX_FUNCTION_ARGS = 32;
-
-private:
-    //! The function's name.
-    std::string _name;
-
-    //! Information about the function's arguments.
-    std::vector<std::string> _argument_names;
-
-    //! The function handle.
-    void* _function;
-};
-
-template <typename ReturnType, typename... Arguments>
-NamedFunction::NamedFunction(std::string name,
-                             std::vector<std::string>&& argument_names,
-                             std::function<ReturnType(Arguments...)>&& function)
-    : _name(std::move(name)),
-      _argument_names(std::move(argument_names)),
-      _function(
-          new std::function<ReturnType(Arguments...)>(std::move(function)))
-{
-    static_assert(
-        ::detail::AllTypesSameAs<double, ReturnType, Arguments...>::value,
-        "Some of the arguments or the return type of the passed function are "
-        "not of the type double.");
-    static_assert(sizeof...(Arguments) <= MAX_FUNCTION_ARGS,
-                  "The function you passed has too many arguments.");
-
-    assert(sizeof...(Arguments) == _argument_names.size());
-}
-
-}  // namespace NumLib
diff --git a/NumLib/NamedFunctionCaller.cpp b/NumLib/NamedFunctionCaller.cpp
deleted file mode 100644
index 1340ac1cee88ec06ede059caa0bb77f1d74894fa..0000000000000000000000000000000000000000
--- a/NumLib/NamedFunctionCaller.cpp
+++ /dev/null
@@ -1,314 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "NamedFunctionCaller.h"
-
-#include <algorithm>
-#include <limits>
-
-#include "BaseLib/Algorithm.h"
-
-/// To understand the working of the algorithm the following two lemmata are
-/// useful:
-/// If a directed graph has a topological order, then it is a Directed Acyclic
-/// Graph (DAG).
-/// If a graph is a DAG, then it has a node with no incoming edges.
-///
-/// \note Negative indices in the graph adjacency list are leaves of the graph.
-///
-/// \param graph represents directed graph given as an adjacency list.
-/// \return true if the given directed graph is acyclic.
-bool hasTopologicalOrdering(std::vector<std::vector<int>> const& graph)
-{
-    // Number of incoming edges for each node of the graph
-    std::vector<unsigned> number_incoming_edges(graph.size());
-
-    // Count all incoming edges (i,j) for each node (i).
-    for (auto const& node_i_adjacencies : graph)
-    {
-        for (int node_j : node_i_adjacencies)
-        {
-            if (node_j >= 0)
-            {  // ignore negative indices.
-                ++number_incoming_edges[node_j];
-            }
-        }
-    }
-
-    // Working queue: a set of nodes with no incoming edges.
-    std::vector<std::size_t> q;
-    for (std::size_t node_i = 0; node_i < number_incoming_edges.size();
-         ++node_i)
-    {
-        if (number_incoming_edges[node_i] == 0)
-        {
-            q.push_back(node_i);
-        }
-    }
-
-
-    auto num_dependent = number_incoming_edges.size() - q.size();
-
-    while (!q.empty())
-    {
-        auto const node_i = q.back();
-        q.pop_back();
-
-        // Decrement counts for all edges (i,j).
-        for (int node_j : graph[node_i])
-        {
-            if (node_j < 0)
-            {
-                continue;  // ignore negative indices
-            }
-            if (--number_incoming_edges[node_j] == 0)
-            {  // Add a node without incoming edges to the queue.
-                q.push_back(node_j);
-                --num_dependent;
-            }
-        }
-    }
-
-    return num_dependent == 0;
-}
-
-enum class TraversePosition { StartNode, BetweenChildren, EndNode };
-
-/*! Traverses the graph given by the adjacency list \c map_sink_source in a
- * depth-first manner starting at the node \c sink_fct.
- *
- * At the beginning and end of each node as well as between every two child
- * nodes the given \c callback is called.
- */
-template <typename Callback>
-void traverse(std::vector<std::vector<int>> const& map_sink_source,
-              int sink_fct, Callback&& callback)
-{
-    assert(sink_fct < static_cast<int>(map_sink_source.size()));
-    callback(sink_fct, TraversePosition::StartNode);
-
-    if (sink_fct < 0)
-    {
-        return;
-    }
-
-    auto const& si_so = map_sink_source[sink_fct];
-    std::size_t const num_args = si_so.size();
-    for (std::size_t sink_arg = 0; sink_arg != num_args; ++sink_arg) {
-        if (sink_arg != 0)
-        {
-            callback(sink_fct, TraversePosition::BetweenChildren);
-        }
-        traverse(map_sink_source, si_so[sink_arg], callback);
-    }
-
-    callback(sink_fct, TraversePosition::EndNode);
-}
-
-namespace NumLib
-{
-NamedFunctionCaller::NamedFunctionCaller(
-    std::initializer_list<std::string> unbound_argument_names)
-    : _uninitialized(-1 - static_cast<int>(unbound_argument_names.size()))
-{
-    assert(unbound_argument_names.size() <
-           static_cast<std::size_t>(std::numeric_limits<int>::max()));
-    int idx = -1;
-    for (auto arg : unbound_argument_names) {
-        BaseLib::insertIfKeyUniqueElseError(
-            _map_name_idx, arg, idx,
-            "The name of the unbound argument is not unique.");
-        --idx;
-    }
-}
-
-void NamedFunctionCaller::addNamedFunction(NamedFunction&& fct)
-{
-    DBUG("Adding named function `%s'", fct.getName().c_str());
-
-    BaseLib::insertIfKeyUniqueElseError(
-        _map_name_idx, fct.getName(), _named_functions.size(),
-        "The name of the function is not unique.");
-
-    _map_sink_source.emplace_back(fct.getArgumentNames().size(),
-                                  _uninitialized);
-    _named_functions.push_back(std::move(fct));
-}
-
-void NamedFunctionCaller::plug(const std::string& sink_fct,
-                               const std::string& sink_arg,
-                               const std::string& source_fct)
-{
-    _deferred_plugs.push_back({sink_fct, sink_arg, source_fct});
-}
-
-void NamedFunctionCaller::applyPlugs()
-{
-    while (!_deferred_plugs.empty())
-    {
-        auto const& plug = _deferred_plugs.back();
-        auto const& sink_fct = plug.sink_fct;
-        auto const& sink_arg = plug.sink_arg;
-        auto const& source = plug.source;
-
-        auto const source_it = _map_name_idx.find(source);
-        if (source_it == _map_name_idx.end())
-        {
-            OGS_FATAL("A function with the name `%s' has not been found.",
-                      source.c_str());
-        }
-        auto const source_idx = source_it->second;
-
-        auto const sink_it = _map_name_idx.find(sink_fct);
-        if (sink_it == _map_name_idx.end())
-        {
-            OGS_FATAL("A function with the name `%s' has not been found.",
-                      sink_fct.c_str());
-        }
-        auto const sink_fct_idx = sink_it->second;
-
-        auto const& sink_args =
-            _named_functions[sink_it->second].getArgumentNames();
-        auto const sink_arg_it =
-            std::find(sink_args.begin(), sink_args.end(), sink_arg);
-        if (sink_arg_it == sink_args.end())
-        {
-            OGS_FATAL(
-                "An argument with the name `%s' has not been found for the "
-                "function `%s'.",
-                sink_arg.c_str(), sink_fct.c_str());
-        }
-        std::size_t const sink_arg_idx =
-            std::distance(sink_args.begin(), sink_arg_it);
-
-        auto& sis_sos = _map_sink_source[sink_fct_idx];
-        if (sis_sos[sink_arg_idx] != _uninitialized)
-        {
-            OGS_FATAL("A dependency for `%s'.`%s' has already been introduced.",
-                      sink_fct.c_str(), sink_arg.c_str());
-        }
-        sis_sos[sink_arg_idx] = source_idx;
-        if (!hasTopologicalOrdering(_map_sink_source))
-        {
-            OGS_FATAL(
-                "The call graph being plugged together must be an acyclic "
-                "graph. The added dependency for `%s'.`%s' introduces a cycle "
-                "into the graph.",
-                sink_fct.c_str(), sink_arg.c_str());
-        }
-
-        _deferred_plugs.pop_back();
-    }
-}
-
-double NamedFunctionCaller::call(
-    std::size_t function_idx,
-    const std::vector<double>& unbound_arguments) const
-{
-    assert(_deferred_plugs.empty() &&
-           "You must call applyPlugs() before this method!");
-
-    auto const& sis_sos = _map_sink_source[function_idx];
-    assert(sis_sos.size() ==
-           _named_functions[function_idx].getArgumentNames().size());
-    std::vector<double> fct_args(sis_sos.size());
-
-    for (std::size_t sink=0; sink<sis_sos.size(); ++sink)
-    {
-        auto const source = sis_sos[sink];
-
-        if (source >= 0) {
-            fct_args[sink] = call(source, unbound_arguments);
-        } else {
-            assert(source != _uninitialized);
-            fct_args[sink] = unbound_arguments[-source-1];
-        }
-    }
-
-    return _named_functions[function_idx].call(fct_args);
-}
-
-std::string NamedFunctionCaller::getCallExpression(
-    std::string const& function_name) const
-{
-    auto const fct_it = _map_name_idx.find(function_name);
-    if (fct_it == _map_name_idx.end()) {
-        OGS_FATAL("A function with the name `%s' has not been found.",
-                  function_name.c_str());
-    }
-
-    std::string expr;
-    auto callback = [&](int fct_idx, TraversePosition pos)
-    {
-        switch (pos) {
-        case TraversePosition::StartNode:
-        {
-            if (fct_idx < 0) {
-                auto it = std::find_if(
-                    _map_name_idx.begin(), _map_name_idx.end(),
-                    [fct_idx](std::pair<std::string, int> const& e) {
-                        return e.second == fct_idx;
-                    });
-                if (it == _map_name_idx.end()) {
-                    OGS_FATAL("The function index %i has not been found.", fct_idx);
-                }
-                expr += it->first;
-            } else {
-                expr += _named_functions[fct_idx].getName() + "(";
-            }
-            break;
-        }
-        case TraversePosition::BetweenChildren:
-            expr += ", ";
-            break;
-        case TraversePosition::EndNode:
-            expr += ")";
-        }
-    };
-
-    traverse(_map_sink_source, fct_it->second, callback);
-    DBUG("expression: %s", expr.c_str());
-    return expr;
-}
-
-SpecificFunctionCaller
-NamedFunctionCaller::getSpecificFunctionCaller(const std::string &function_name)
-{
-    auto const fct_it = _map_name_idx.find(function_name);
-    if (fct_it == _map_name_idx.end()) {
-        OGS_FATAL("A function with the name `%s' has not been found.",
-                  function_name.c_str());
-    }
-    return SpecificFunctionCaller(fct_it->second, *this);
-}
-
-std::size_t NamedFunctionCaller::getNumberOfUnboundArguments() const
-{
-    return -_uninitialized - 1;
-}
-
-SpecificFunctionCaller::SpecificFunctionCaller(const std::size_t function_idx,
-                                             const NamedFunctionCaller& caller)
-    : _function_idx(function_idx), _caller(caller)
-{
-}
-
-double SpecificFunctionCaller::call(
-    const std::vector<double>& unbound_arguments) const
-{
-    return _caller.call(_function_idx, unbound_arguments);
-}
-
-std::size_t SpecificFunctionCaller::getNumberOfUnboundArguments() const
-{
-    return _caller.getNumberOfUnboundArguments();
-}
-
-} // namespace NumLib
diff --git a/NumLib/NamedFunctionCaller.h b/NumLib/NamedFunctionCaller.h
deleted file mode 100644
index 276867c77eaadaeae48943a1176cada5c9322c3a..0000000000000000000000000000000000000000
--- a/NumLib/NamedFunctionCaller.h
+++ /dev/null
@@ -1,126 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include <map>
-#include <vector>
-#include "NamedFunction.h"
-
-namespace NumLib
-{
-class SpecificFunctionCaller;
-
-//! Builds expression trees of named functions dynamically at runtime.
-class NamedFunctionCaller final
-{
-public:
-    //! Constructs an instance whose unbound arguments have the given names.
-    NamedFunctionCaller(
-        std::initializer_list<std::string> unbound_argument_names);
-
-    //! Adds the given named function
-    void addNamedFunction(NamedFunction&& fct);
-
-    //! Returns all named functions associated with the caller instance.
-    std::vector<NamedFunction> const& getNamedFunctions() const
-    {
-        return _named_functions;
-    }
-
-    //! Declares that the argument with name \c sink_arg of the function \c
-    //! sink_fct is being computed by the function \c source_fct.
-    //!
-    //! The functions involved need not already be known to the
-    //! NamedFunctionCaller.
-    void plug(std::string const& sink_fct, std::string const& sink_arg,
-              std::string const& source_fct);
-
-    //! Actually plug all the plugs previously declared.
-    //!
-    //! \pre All functions involved must have been added.
-    void applyPlugs();
-
-    //! Creates a function caller that is able to call the function with the
-    //! given name.
-    SpecificFunctionCaller getSpecificFunctionCaller(
-        std::string const& function_name);
-
-    //! Returns a string representing the expression graph of the given
-    //! function.
-    //!
-    //! \pre applyPlugs() must have been called before.
-    std::string getCallExpression(std::string const& function_name) const;
-
-    //! Returns the number of unbound arguments.
-    std::size_t getNumberOfUnboundArguments() const;
-
-private:
-    //! Calls the function with the given index with the given unbound
-    //! arguments.
-    double call(std::size_t function_idx,
-                std::vector<double> const& unbound_arguments) const;
-
-    //! Maps function names to indices.
-    //! Negative indices refer to unbound arguments.
-    std::map<std::string, int> _map_name_idx;
-
-    //! Contains all named functions.
-    std::vector<NamedFunction> _named_functions;
-
-    //! The expression graph.
-    //! Contains for each named function (outer vector) a vector which maps each
-    //! function argument to the source function index that computes this
-    //! argument.
-    std::vector<std::vector<int>> _map_sink_source;
-
-    //! Magic number used to mark function arguments in \c _map_sink_source
-    //! whose source functions have not yet been set up.
-    const int _uninitialized;
-
-    struct SinkSource
-    {
-        std::string const sink_fct;
-        std::string const sink_arg;
-        std::string const source;
-    };
-
-    //! Saves plugs declared by plug().
-    std::vector<SinkSource> _deferred_plugs;
-
-    friend class SpecificFunctionCaller;
-};
-
-//! A function caller that can call one specific function.
-//!
-//! \todo Use this class to provide some optimizations of the expression
-//! evaluation.
-class SpecificFunctionCaller final
-{
-public:
-    //! Constructs a new instance.
-    SpecificFunctionCaller(std::size_t const function_idx,
-                          NamedFunctionCaller const& caller);
-
-    //! Call the function set up with the given unbound arguments.
-    double call(std::vector<double> const& unbound_arguments) const;
-
-    //! Returns the number of unbound arguments.
-    std::size_t getNumberOfUnboundArguments() const;
-
-private:
-    //! Index of the referenced function.
-    std::size_t const _function_idx;
-
-    //! The named function caller used for the evaluation.
-    NamedFunctionCaller const& _caller;
-};
-
-} // namespace NumLib
diff --git a/ProcessLib/Output/GlobalVectorFromNamedFunction.cpp b/ProcessLib/Output/GlobalVectorFromNamedFunction.cpp
deleted file mode 100644
index d6c73154bb9510f209a37013c94279f9e1329ad1..0000000000000000000000000000000000000000
--- a/ProcessLib/Output/GlobalVectorFromNamedFunction.cpp
+++ /dev/null
@@ -1,68 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "GlobalVectorFromNamedFunction.h"
-#include "MathLib/LinAlg/FinalizeVectorAssembly.h"
-#include "MathLib/LinAlg/MatrixVectorTraits.h"
-#include "NumLib/DOF/DOFTableUtil.h"
-
-namespace ProcessLib
-{
-GlobalVectorFromNamedFunction::GlobalVectorFromNamedFunction(
-    NumLib::SpecificFunctionCaller&& function_caller,
-    MeshLib::Mesh const& mesh,
-    NumLib::LocalToGlobalIndexMap const& dof_table_single,
-    SecondaryVariableContext& context)
-    : _function_caller(std::move(function_caller)),
-      _mesh(mesh),
-      _dof_table_single(dof_table_single),
-      _context(context)
-{
-    assert(dof_table_single.getNumberOfComponents() == 1);
-}
-
-GlobalVector const& GlobalVectorFromNamedFunction::call(
-    const double /*t*/,
-    GlobalVector const& x,
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::unique_ptr<GlobalVector>& result)
-{
-    result = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
-        {_dof_table_single.dofSizeWithoutGhosts(),
-         _dof_table_single.dofSizeWithoutGhosts(),
-         &_dof_table_single.getGhostIndices(), nullptr});
-
-    GlobalIndexType nnodes = _mesh.getNumberOfNodes();
-
-    auto const n_args = _function_caller.getNumberOfUnboundArguments();
-    assert(dof_table.getNumberOfComponents() == static_cast<int>(n_args));
-    std::vector<double> args(n_args);
-
-    for (GlobalIndexType node_id = 0; node_id < nnodes; ++node_id)
-    {
-        // TODO maybe fill args via callback mechanism or remove this class
-        // entirely. Caution: The order of args will be the same as the order of
-        // the components in the global vector!
-        for (std::size_t i = 0; i < n_args; ++i)
-        {
-            args[i] = NumLib::getNodalValue(x, _mesh, dof_table, node_id, i);
-        }
-
-        _context.index = node_id;
-        auto const value = _function_caller.call(args);
-
-        result->set(node_id, value);
-    }
-
-    MathLib::finalizeVectorAssembly(*result);
-    return *result;
-}
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/Output/GlobalVectorFromNamedFunction.h b/ProcessLib/Output/GlobalVectorFromNamedFunction.h
deleted file mode 100644
index 656d73e2225ea5f33fbc27b423cd80538aa38922..0000000000000000000000000000000000000000
--- a/ProcessLib/Output/GlobalVectorFromNamedFunction.h
+++ /dev/null
@@ -1,56 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include "MeshLib/Mesh.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "NumLib/NamedFunctionCaller.h"
-#include "NumLib/NumericsConfig.h"
-#include "SecondaryVariableContext.h"
-
-namespace ProcessLib
-{
-//! Computes a global vector from a NamedFunction (which can only compute double
-//! values).
-class GlobalVectorFromNamedFunction final
-{
-public:
-    /*! Constructs a new instance.
-     *
-     * \param function_caller will provide the individual entries of the
-     * GlobalVector to be computed.
-     * \param mesh to which the \c dof_table_single is associated
-     * \param dof_table_single used for constructing the GlobalVector
-     * \param context used by the \c function_caller to access "auxiliary
-     * unbound arguments".
-     */
-    GlobalVectorFromNamedFunction(
-        NumLib::SpecificFunctionCaller&& function_caller,
-        MeshLib::Mesh const& mesh,
-        NumLib::LocalToGlobalIndexMap const& dof_table_single,
-        SecondaryVariableContext& context);
-
-    //! Computes the GlobalVector.
-    //!
-    //! The signature of this method matches
-    //! SecondaryVariableFunctions::Function, i.e., this method can be used to
-    //! compute a secondary variable.
-    GlobalVector const& call(const double t, GlobalVector const& x,
-                             NumLib::LocalToGlobalIndexMap const& dof_table,
-                             std::unique_ptr<GlobalVector>& result);
-
-private:
-    NumLib::SpecificFunctionCaller _function_caller;
-    MeshLib::Mesh const& _mesh;
-    NumLib::LocalToGlobalIndexMap const& _dof_table_single;
-    SecondaryVariableContext& _context;
-};
-}  // namespace ProcessLib
diff --git a/ProcessLib/Output/SecondaryVariableContext.h b/ProcessLib/Output/SecondaryVariableContext.h
deleted file mode 100644
index 6f4c072fa94af70d346eeee6fe58718598147528..0000000000000000000000000000000000000000
--- a/ProcessLib/Output/SecondaryVariableContext.h
+++ /dev/null
@@ -1,31 +0,0 @@
-/**
- * \file
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
-
-namespace ProcessLib
-{
-/*! Provides a "global" context, e.g., for NamedFunction's.
- *
- * Currently the context comprises an index. It can be used, e.g. to compute a
- * NamedFunction which needs additional external data apart from the unbound
- * variables that are configured.
- */
-struct SecondaryVariableContext
-{
-public:
-    //! Points to the position in a GlobalVector being read or written right
-    //! now.
-    GlobalIndexType index = 0;
-};
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 048a872ee70212118b6025903d73c4b83e5e8bfe..dac131b0eb2e7be072def79edc08674bd5aea615 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -17,7 +17,6 @@
 #include "NumLib/ODESolver/TimeDiscretization.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/BoundaryCondition/BoundaryConditionCollection.h"
-#include "ProcessLib/Output/SecondaryVariableContext.h"
 #include "ProcessLib/Output/ExtrapolatorData.h"
 #include "ProcessLib/Output/IntegrationPointWriter.h"
 #include "ProcessLib/Output/SecondaryVariable.h"
diff --git a/Tests/NumLib/TestNamedFunction.cpp b/Tests/NumLib/TestNamedFunction.cpp
deleted file mode 100644
index 2fd5a402b9bfb52b41cfaa174c395b64e1a2bd5d..0000000000000000000000000000000000000000
--- a/Tests/NumLib/TestNamedFunction.cpp
+++ /dev/null
@@ -1,196 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include <gtest/gtest.h>
-#include <logog/include/logog.hpp>
-
-#include "BaseLib/Functional.h"
-#include "NumLib/NamedFunctionProvider.h"
-#include "NumLib/NamedFunctionCaller.h"
-#include "Tests/InstanceCounter.h"
-
-class F : public NumLib::NamedFunctionProvider
-{
-public:
-    double f(double arg_g, double arg_y) const
-    {
-        return arg_g + arg_y;
-    }
-
-    std::vector<NumLib::NamedFunction>
-    getNamedFunctions() const override
-    {
-        return {{"f", {"g_arg", "y"}, BaseLib::easyBind(&F::f, this)}};
-    }
-};
-
-class G : public NumLib::NamedFunctionProvider
-{
-public:
-    double g(double arg_x) const
-    {
-        return -arg_x;
-    }
-
-    std::vector<NumLib::NamedFunction>
-    getNamedFunctions() const override
-    {
-        return {{"g", {"x"}, BaseLib::easyBind(&G::g, this)}};
-    }
-};
-
-class H : public InstanceCounter<H>
-{
-public:
-    explicit H(double const z) : _z(z) {}
-    double h(double arg_x, double arg_y) { return arg_x * arg_y - _z; }
-    double setZ(double const z)
-    {
-        _z = z;
-        return z;
-    }
-
-private:
-    double _z;
-};
-
-class I : public NumLib::NamedFunctionProvider
-{
-public:
-    double i(double arg_x) const
-    {
-        return -arg_x;
-    }
-
-    std::vector<NumLib::NamedFunction>
-    getNamedFunctions() const override
-    {
-        return {{"i", {"x"}, BaseLib::easyBind(&I::i, this)}};
-    }
-};
-
-TEST(NumLib, NamedFunctionCaller)
-{
-    F f_inst;
-    G g_inst;
-
-    NumLib::NamedFunctionCaller caller{ "x", "y" };
-
-    for (auto&& f_named : f_inst.getNamedFunctions()) {
-        caller.addNamedFunction(std::move(f_named));
-    }
-
-    caller.plug("f", "g_arg", "g");
-    caller.plug("f", "y", "y");
-    caller.plug("g", "x", "x");
-
-    // test if adding function after plug works
-    for (auto&& g_named : g_inst.getNamedFunctions()) {
-        caller.addNamedFunction(std::move(g_named));
-    }
-
-    caller.applyPlugs();
-
-    double x = 1.0;
-    double y = 2.0;
-
-    auto const g_caller = caller.getSpecificFunctionCaller("g");
-    DBUG("calling %s", caller.getCallExpression("g").c_str());
-    EXPECT_EQ(g_inst.g(x), g_caller.call({x, y}));
-
-    auto const f_caller = caller.getSpecificFunctionCaller("f");
-    DBUG("calling %s", caller.getCallExpression("f").c_str());
-    EXPECT_EQ(f_inst.f(g_inst.g(x), y), f_caller.call({x, y}));
-}
-
-TEST(NumLib, NamedFunctionCallerCyclicGraph)
-{
-    // Construct a cyclic case with f(g(i(f(...), y)))
-    F f_inst;
-    G g_inst;
-    I i_inst;
-
-    NumLib::NamedFunctionCaller caller{ "x", "y" };
-
-    for (auto&& f_named : f_inst.getNamedFunctions()) {
-        caller.addNamedFunction(std::move(f_named));
-    }
-    for (auto&& g_named : g_inst.getNamedFunctions()) {
-        caller.addNamedFunction(std::move(g_named));
-    }
-    for (auto&& i_named : i_inst.getNamedFunctions()) {
-        caller.addNamedFunction(std::move(i_named));
-    }
-
-    caller.plug("f", "g_arg", "g");
-    caller.plug("f", "y", "y");
-    caller.plug("g", "x", "i");
-    caller.plug("i", "x", "f");
-
-    ASSERT_ANY_THROW(caller.applyPlugs());
-}
-
-TEST(NumLib, NamedFunctionNoLeaks)
-{
-    auto num_const = InstanceCounter<H>::getNumberOfConstructions();
-    auto num_move = InstanceCounter<H>::getNumberOfMoves();
-    auto num_copy = InstanceCounter<H>::getNumberOfCopies();
-    auto num_inst = InstanceCounter<H>::getNumberOfInstances();
-
-    {
-        H h_inst(1.0);
-        EXPECT_EQ(num_const+1, InstanceCounter<H>::getNumberOfConstructions());
-        EXPECT_EQ(num_inst+1, InstanceCounter<H>::getNumberOfInstances());
-        InstanceCounter<H>::update(num_const, num_move, num_copy, num_inst);
-
-        auto h_fct = NumLib::NamedFunction("h", {"x", "y"},
-                                           BaseLib::easyBind(&H::h, h_inst));
-        EXPECT_EQ(num_const, InstanceCounter<H>::getNumberOfConstructions());
-        EXPECT_EQ(num_inst, InstanceCounter<H>::getNumberOfInstances());
-
-        EXPECT_EQ(5.0, h_fct.call({ 2.0, 3.0 }));
-        h_inst.setZ(2.0);
-        EXPECT_EQ(4.0, h_fct.call({ 2.0, 3.0 }));
-
-        auto h_bind = BaseLib::easyBind(&H::h, H{3.0});
-        EXPECT_EQ(num_const+1, InstanceCounter<H>::getNumberOfConstructions());
-        EXPECT_EQ(num_inst+1, InstanceCounter<H>::getNumberOfInstances());
-        InstanceCounter<H>::update(num_const, num_move, num_copy, num_inst);
-
-        // Move an object. NamedFunction will implicitly do memory management.
-        auto h_fct2 = NumLib::NamedFunction("h", {"x", "y"}, std::move(h_bind));
-        EXPECT_EQ(num_copy, InstanceCounter<H>::getNumberOfCopies());
-        EXPECT_EQ(num_const, InstanceCounter<H>::getNumberOfConstructions());
-        EXPECT_EQ(num_inst, InstanceCounter<H>::getNumberOfInstances());
-
-        EXPECT_EQ(5.0, h_fct2.call({ 2.0, 4.0 }));
-
-        // copy
-        NumLib::NamedFunction h_fct3 = h_fct2;
-        EXPECT_EQ(num_const, InstanceCounter<H>::getNumberOfConstructions());
-        EXPECT_EQ(num_copy+1, InstanceCounter<H>::getNumberOfCopies());
-        EXPECT_EQ(num_inst+1, InstanceCounter<H>::getNumberOfInstances());
-        InstanceCounter<H>::update(num_const, num_move, num_copy, num_inst);
-
-        EXPECT_EQ(-1.0, h_fct3.call({ 2.0, 1.0 }));
-
-        // move
-        NumLib::NamedFunction h_fct4 = std::move(h_fct3);
-        EXPECT_INSTANCES(H, num_const, num_move, num_copy, num_inst);
-        InstanceCounter<H>::update(num_const, num_move, num_copy, num_inst);
-
-        EXPECT_EQ(3.0, h_fct4.call({ 3.0, 2.0 }));
-    }
-    EXPECT_EQ(num_const, InstanceCounter<H>::getNumberOfConstructions());
-    EXPECT_EQ(num_copy, InstanceCounter<H>::getNumberOfCopies());
-    // If zero instances are left, the destructor has been called the right
-    // number of times, i.e., all internal casts in NamedFunction have been
-    // successful.
-    EXPECT_EQ(0, InstanceCounter<H>::getNumberOfInstances());
-}