dune/fufem/python/callable.hh

@@ -17,16 +17,12 @@
 #include <dune/common/exceptions.hh>
 #include <dune/common/shared_ptr.hh>
 #include <dune/common/fvector.hh>
-#include <dune/common/function.hh>
 #include <dune/common/parametertree.hh>
 #include <dune/common/typetraits.hh>
+#include <dune/common/indices.hh>
 
 #include <dune/grid/common/gridfactory.hh>
 
-#include <dune/fufem/functions/virtualfunctiontoboundarysegmentadapter.hh>
-#include <dune/fufem/formatstring.hh>
-
-
 #include <dune/fufem/python/reference.hh>
 
 
@@ -35,6 +31,131 @@ namespace Python
 
 // forward declarations
 void handlePythonError(const std::string&, const std::string&);
+Reference dict();
+
+namespace Imp
+{
+    // forward declarations
+    PyObject* inc(PyObject* p);
+
+
+
+    // Helper for marking a keyword argument
+    template<class Key, class Value>
+    struct KeyWordArgument
+    {
+        Key key;
+        Value value;
+    };
+
+
+
+    // Check for KeyWordArgument
+    template<class T>
+    struct IsKeyWordArgument : std::false_type {};
+
+    template<class Key, class Value>
+    struct IsKeyWordArgument<KeyWordArgument<Key, Value>> : std::true_type {};
+
+    // Count positional (non-keyword) arguments up to pos-1.
+    // The result is returned as compile time constant using index_constant.
+    template<std::size_t pos, class... Args>
+    auto positionalArgumentCount()
+    {
+        using ArgTuple = std::tuple<Args...>;
+        return Dune::unpackIntegerSequence([] (auto... i) {
+            return Dune::index_constant<(0 + ... + not(Imp::IsKeyWordArgument<std::decay_t<std::tuple_element_t<i, ArgTuple>>>::value))>();
+        }, std::make_index_sequence<pos>{});
+    }
+
+    // Helper for tagging string as keyword
+    template<class Char>
+    struct KeyWord {
+        std::string key;
+        template <typename Value>
+        KeyWordArgument<std::string, Value> operator=(Value&& value) const {
+            return {key, std::forward<Value>(value)};
+        }
+    };
+
+    // Parse argument list. The list may contain positional and keyword arguments.
+    // The parameters to this function are two callback functions, followed by
+    // the arguments to be parsed. Positional and keyword arguments may be interleaved.
+    // In any case only positional arguments are counted when determining the position
+    // of a positional argument.
+    //
+    // The first callback will be called with (pos,arg) for each positional argument arg
+    // with its position pos (as compile time value using Dune::index_constant).
+    // The second callback will be called with (keyword, value) for each
+    // keyword argument.
+    template<class PositionalCallback, class KeyWordCallBack, class... Args>
+    void parseKeywordArguments(PositionalCallback&& posCallback, KeyWordCallBack&& kwCallBack, Args&&... args)
+    {
+        auto processArgument = [&](auto i, auto&& arg)
+        {
+            if constexpr (IsKeyWordArgument<std::decay_t<decltype(arg)>>::value)
+                kwCallBack(arg.key, arg.value);
+            else
+                posCallback(positionalArgumentCount<i, Args...>(), std::forward<decltype(arg)>(arg));
+        };
+        // Unfold argument pack to process arguments. Doing this via
+        // std::initializer_list guarantees the proper evaluation order.
+        // To additionally pass the argument indices as variadic argument pack,
+        // we use unpackIntegerSequence.
+        //
+        // Here one would like to simply unpack args like this:
+        // std::initializer_list<int>{ (processArgument(i, args), 0)...};
+        // But capturing parameter packs does not work in C++17
+        // (in fact gcc allows it). Hence we have to artificially
+        // wrap in this in a tuple and access entries using std::get.
+        //
+        // Create a tuple preserving value categories.
+        // This will store appropriate l- and r-value references
+        // to the individual arguments.
+        auto argTuple = std::forward_as_tuple(std::forward<Args>(args)...);
+        Dune::unpackIntegerSequence([&](auto... i) {
+            // Extract tuple entries restoring the value category.
+            // Notice that due to reference collapsing this also returns
+            // l-value references properly. Furthermore doing multiple
+            // move's on the same tuple is OK here, because each move
+            // only changes the value category of the returned entry
+            // but does not modify the tuple itself.
+            std::initializer_list<int>{ (processArgument(i, std::get<i>(std::move(argTuple))), 0)...};
+        }, std::index_sequence_for<Args...>());
+    }
+
+} // namespace Imp
+
+
+/**
+ * \brief Create a keyword argument
+ *
+ * While the value may be stored as reference or
+ * value, if an l-value or r-value reference
+ * is passed, the key is always stored by value.
+ * As a special case, when passing a
+ * const char* it will be converted to a std::string.
+ *
+ * \param key Identifier of the argument
+ * \param value Value of the argument
+ */
+template<class Key, class Value>
+auto arg(Key key, Value&& value)
+{
+    using StoredKey = std::conditional_t<std::is_same_v<Key, const char*>, std::string, Key>;
+    return Imp::KeyWordArgument<StoredKey, Value>{key, std::forward<Value>(value)};
+}
+
+namespace Literals {
+    /**
+     * \brief User defined literal for defining string-based keyword arguments
+     *
+     * This allows to create keyword arguments using "keyword"_a=value.
+     */
+    Imp::KeyWord<char> operator "" _a(const char* keyword, std::size_t) {
+        return {keyword};
+    }
+}
 
 
 
@@ -150,46 +271,60 @@ class Callable :
         }
 
         /**
-         * \brief Call this object without arguments
+         * \brief Call this Reference with positional arguments given as tuple and keyword arguments given as dictionary
          *
-         * Shortcut for callWithArgumentTuple(makeTuple()).
-         */
-        Reference operator() () const
-        {
-            return callWithArgumentTuple(makeTuple());
-        }
-
-        /**
-         * \brief Call this object with one argument
+         * If the Reference represents a function it's called.
+         * If the Reference represents an instance its __call__ method is invoked.
+         * If the Reference represents a class a constructor is invoked.
          *
-         * Shortcut for callWithArgumentTuple(makeTuple(t0)).
-         */
-        template<class T0>
-        Reference operator() (const T0& t0) const
-        {
-            return callWithArgumentTuple(makeTuple(t0));
-        }
-
-        /**
-         * \brief Call this object with two argument
+         * Although the method is const it might change the refered object!
+         * This cannot be avoided since the python api does not know
+         * about const pointers so we use a mutable pointer internally.
+         *
+         * \param args Positional arguments represented as tuple
+         * \param keywordArgs Keyword arguments represented as dictionary
+         * \returns The result of the call
          *
-         * Shortcut for callWithArgumentTuple(makeTuple(t0,t1)).
          */
-        template<class T0, class T1>
-        Reference operator() (const T0& t0, const T1& t1) const
+        Reference callWithArgumentTupleAndKeywordArgs(const Reference& args, const Reference& keywordArgs) const
         {
-            return callWithArgumentTuple(makeTuple(t0, t1));
+            assertPyObject("Callable::callWithArgumentTupleAndKeywordArgs()");
+            PyObject* result = PyObject_Call(p_, args, keywordArgs);
+            if (not result)
+                handlePythonError("Callable::callWithArgumentTupleAndKeywordArgs()", "failed to call object");
+            return result;
         }
 
         /**
-         * \brief Call this object with three argument
+         * \brief Call this object with given arguments
          *
-         * Shortcut for callWithArgumentTuple(makeTuple(t0,t1,t2)).
+         * Convert given arguments to Python-types and pass them
+         * as arguments to Python-callable objects. Keyword arguments
+         * can be passed using either Python::arg("keyword", value)
+         * or "keyword"_a = value. For the latter you have to import
+         * the namespace Python::Literals.
+         *
+         * \returns Result of the call as Python object.
          */
-        template<class T0, class T1, class T2>
-        Reference operator() (const T0& t0, const T1& t1, const T2& t2) const
+        template<class... Args>
+        Reference operator() (const Args&... args) const
+        {
+            static constexpr std::size_t kwArgCount = (0 + ... + Imp::IsKeyWordArgument<Args>::value);
+            if constexpr (kwArgCount == 0)
+                return callWithArgumentTuple(makeTuple(args...));
+            else
             {
-            return callWithArgumentTuple(makeTuple(t0, t1, t2));
+                Reference kwArgDict = dict();
+                Reference argTuple = PyTuple_New(sizeof...(Args)-kwArgCount);
+                Imp::parseKeywordArguments(
+                    [&](auto position, auto&& arg) {
+                        PyTuple_SetItem(argTuple, position, Imp::inc(makeObject(arg)));
+                    }, [&](auto keyword, auto&& value) {
+                        PyDict_SetItem(kwArgDict, makeObject(keyword), makeObject(value));
+                        // We have to use inc() since PyTuple_SetItem STEALS a reference!
+                    }, args...);
+                return callWithArgumentTupleAndKeywordArgs(argTuple, kwArgDict);
+            }
         }
 
     protected:

dune/fufem/python/common.hh

@@ -12,9 +12,9 @@
 #include <string>
 #include <sstream>
 #include <type_traits>
+#include <tuple>
 
 #include <dune/common/exceptions.hh>
-#include <dune/common/std/apply.hh>
 
 #include <dune/functions/common/signature.hh>
 #include <dune/functions/common/differentiablefunctionfromcallables.hh>
@@ -368,6 +368,16 @@ Reference tuple(const Ts&... ts)
     return makeTuple(ts...);
 }
 
+/**
+ * \brief Create an empty Python tuple.
+ *
+ * Entries can be inserted using setItem()
+ */
+Reference dict()
+{
+    return PyDict_New();
+}
+
 /**
  * \brief Get size of a python sequence
  *
@@ -560,7 +570,7 @@ auto makeDifferentiableFunction(const R&... f)
     auto signatureTag = Dune::Functions::SignatureTag<Signature, DerivativeTraits>();
     auto derivativeSignatureTags = Dune::Functions::derivativeSignatureTags<sizeof...(f)-1>(signatureTag);
 
-    return Dune::Std::apply([&f..., &signatureTag](auto... tag) {
+    return std::apply([&f..., &signatureTag](auto... tag) {
             return Dune::Functions::makeDifferentiableFunctionFromCallables(signatureTag,
                 makeFunction<typename decltype(tag)::Signature>(f)...);
         }, derivativeSignatureTags);

dune/fufem/python/conversion.hh

@@ -20,6 +20,10 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/parametertree.hh>
 #include <dune/common/typetraits.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/tuplevector.hh>
+#include <dune/common/indices.hh>
+#include <dune/common/rangeutilities.hh>
 
 #include <dune/geometry/utility/typefromvertexcount.hh>
 
@@ -29,12 +33,35 @@
 
 #include <dune/geometry/utility/typefromvertexcount.hh>
 
-#include <dune/fufem/functions/virtualfunctiontoboundarysegmentadapter.hh>
-#include <dune/fufem/formatstring.hh>
-
 #include <dune/fufem/python/reference.hh>
 #include <dune/fufem/python/common.hh>
+#include <dune/fufem/python/callable.hh>
+
+namespace Impl
+{
 
+template<int dim, int dimworld = dim, class ctype = double>
+class BoundarySegmentFromCallable
+    : public Dune::BoundarySegment<dim, dimworld, ctype>
+{
+    using Domain = Dune::FieldVector<ctype, dim-1>;
+    using Range = Dune::FieldVector<ctype, dimworld>;
+
+public:
+    BoundarySegmentFromCallable(std::function<Range(Domain)> f)
+        : f_(f)
+    {}
+
+    Range operator()(const Domain& x) const
+    {
+        return f_(x);
+    }
+
+private:
+    std::function<Range(Domain)> f_;
+};
+
+}
 
 namespace Python
 {
@@ -373,6 +400,45 @@ struct Conversion<Dune::BlockVector<T> >
 };
 
 
+// conversion of FieldVector
+template<class... T>
+struct Conversion<Dune::TupleVector<T...> >
+{
+    enum {useDefaultConstructorConversion=true};
+    static void toC(PyObject* list, Dune::TupleVector<T...>& v)
+    {
+        static constexpr auto n = Dune::index_constant<sizeof...(T)>{};
+        if (not PySequence_Check(list))
+            DUNE_THROW(Dune::Exception, "cannot convert a non-sequence to a Dune::TupleVector<T...>");
+        else
+        {
+            if ((std::size_t)PySequence_Size(list)!=v.size())
+                DUNE_THROW(Dune::Exception, "cannot convert a sequence of size " << PySequence_Size(list) << " to a Dune::TupleVector of size " << v.size());
+
+            Dune::Hybrid::forEach(Dune::range(n), [&](auto i) {
+                Reference(PySequence_GetItem(list, i)).toC(v[i]);
+            });
+            // The above is OK since PySequence_GetItem returns a new reference.
+            // With PyTuple_GetItem we would need one of the following since it
+            // returns a borrowed reference:
+            // Conversion<T>::toC(PySequence_GetItem(list, i), v[i]);
+            // Reference(Imp::inc(PySequence_GetItem(list, i))).toC(v[i]);
+        }
+    }
+
+    static PyObject* toPy(const Dune::TupleVector<T...>& v)
+    {
+        static constexpr auto n = Dune::index_constant<sizeof...(T)>{};
+        PyObject* tuple = PyTuple_New(n);
+        Dune::Hybrid::forEach(Dune::range(n), [&](auto i) {
+            PyTuple_SetItem(tuple, i, Imp::inc(makeObject(v[i])));
+        });
+        // We have to use inc() since PyTuple_SetItem STEALS a reference!
+        return tuple;
+    }
+};
+
+
 // conversion to dict to ParameterTree
 template<>
 struct Conversion<Dune::ParameterTree>
@@ -403,6 +469,31 @@ struct Conversion<Dune::ParameterTree>
     }
 };
 
+template<class Domain, class Range>
+struct Conversion<std::function<Range(Domain)>>
+{
+    enum {useDefaultConstructorConversion=true};
+    static void toC(PyObject* pyF, std::function<Range(Domain)>& f)
+    {
+        f = Python::makeFunction<Range(Domain)>(Imp::inc(pyF));
+    }
+};
+
+template<int dim, int dimworld, class ctype>
+struct Conversion<Dune::BoundarySegment<dim, dimworld, ctype>*>
+{
+    enum {useDefaultConstructorConversion=true};
+    static void toC(PyObject* pyF, Dune::BoundarySegment<dim, dimworld, ctype>*& segmentPointer)
+    {
+        using Range = Dune::FieldVector<double,dimworld>;
+
+        using Segment = Impl::BoundarySegmentFromCallable<dim, dimworld, ctype>;
+
+        segmentPointer = new Segment(make_function<Range>(Imp::inc(pyF)));
+    }
+};
+
+
 /**
  * \brief Conversion of python object to GridFactory
  *

dune/fufem/python/function.hh

-// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set ts=8 sw=4 et sts=4:
-#ifndef DUNE_FUFEM_PYTHON_FUNCTION_HH
-#define DUNE_FUFEM_PYTHON_FUNCTION_HH
-
-// Only introduce the dune-python interface if python
-// was found and enabled.
-#if HAVE_PYTHON || DOXYGEN
-
-#include <Python.h>
-
-#include <string>
-#include <sstream>
-
-#include <dune/common/exceptions.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/function.hh>
-#include <dune/common/typetraits.hh>
-
-#include <dune/fufem/functions/virtualfunctiontoboundarysegmentadapter.hh>
-#include <dune/fufem/functions/virtualdifferentiablefunction.hh>
-
-#include <dune/fufem/python/reference.hh>
-#include <dune/fufem/python/callable.hh>
-#include <dune/fufem/python/common.hh>
-#include <dune/fufem/python/conversion.hh>
-
-
-
-/**
- * \brief Function using embedded python for evaluation.
- *
- * \tparam DT Domain type
- * \tparam RT Range type
- */
-template<class DT, class RT>
-class PythonFunction :
-    public Dune::VirtualFunction<DT, RT>
-{
-    public:
-
-        /**
-         * \brief Create PythonFunction wrapping a callable python object
-         */
-        PythonFunction(const Python::Callable callable)
-        {
-            if (not callable)
-                DUNE_THROW(Dune::Exception, "Trying to construct PythonFunction from NULL");
-            callable_ = callable;
-        }
-
-        /**
-         * \brief Create PythonFunction by name of a python function
-         *
-         * This creates a PythonFunction that
-         * wraps an already existing function
-         * with the given name.
-         */
-        PythonFunction(const std::string& name) DUNE_DEPRECATED_MSG("Use PythonFunction(Python::main().get(name) instead")
-        {
-            // get reference for python function
-            Python::Callable callable = Python::main().get(name);
-            if (not callable)
-                DUNE_THROW(Dune::Exception, "Trying to construct PythonFunction from NULL");
-            callable_ = callable;
-        }
-
-        /**
-         * \brief Create PythonFunction from string expression.
-         *
-         * This will first create a python function with the given name
-         * that evaluates the given expression and than wrap
-         * this as PythonFunction.
-         *
-         * If expression contains no newline it is directly
-         * evaluated and returned. If it contains a newline
-         * it is used as function body and expected to call
-         * return itself.
-         */
-        PythonFunction(const std::string& name, const std::string& expression)
-        {
-            // create a python function
-            std::stringstream s;
-            std::string::size_type newline = expression.find("\n");
-            if (newline != std::string::npos)
-                s << "def " << name << "(x):\n" << expression;
-            else
-                s << "def " << name << "(x):\n    return " << expression;
-
-            // Execute code in __main__ and get reference for python function
-            Python::Module main = Python::main();
-            main.run(s.str());
-            callable_ = main.get(name);
-        }
-
-        /**
-         * \copydoc VirtualFunction::evaluate
-         */
-        void evaluate(const DT& x, RT& y) const
-        {
-            callable_(x).toC(y);
-        }
-
-        /**
-         * \brief Evaluate function
-         */
-        RT operator()(const DT& x) const
-        {
-            return callable_(x).template toC<RT>();
-        }
-
-    protected:
-
-        Python::Callable callable_;
-};
-
-
-/**
- * \brief Function using embedded python for evaluation of values and derivatives.
- *
- * \tparam DT Domain type
- * \tparam RT Range type
- */
-template<class DT, class RT>
-class DifferentiablePythonFunction :
-    public VirtualDifferentiableFunction<DT, RT>
-{
-        typedef VirtualDifferentiableFunction<DT, RT> Base;
-    public:
-        typedef typename Base::DomainType DomainType;
-        typedef typename Base::RangeType RangeType;
-        typedef typename Base::DerivativeType DerivativeType;
-
-        /**
-         * \brief Create DifferentiablePythonFunction wrapping two callable python object
-         */
-        DifferentiablePythonFunction(const Python::Callable value, const Python::Callable derivative)
-        {
-            if (not value)
-                DUNE_THROW(Dune::Exception, "Trying to construct DifferentiablePythonFunction with NULL as value argument");
-            if (not derivative)
-                DUNE_THROW(Dune::Exception, "Trying to construct DifferentiablePythonFunction with NULL as derivative argument");
-            value_ = value;
-            derivative_ = derivative;
-        }
-
-        /**
-         * \brief Create DifferentiablePythonFunction from sequence of callables
-         */
-        DifferentiablePythonFunction(const Python::Reference fdf)
-        {
-            if (not Python::isSequence(fdf))
-                DUNE_THROW(Dune::Exception, "Trying to construct DifferentiablePythonFunction from nonsequence python object");
-            int size = Python::size(fdf);
-            if (size<2)
-                DUNE_THROW(Dune::Exception, "Trying to construct DifferentiablePythonFunction from sequence with size<2");
-            Python::Callable value(Python::getItem(fdf, 0));
-            Python::Callable derivative(Python::getItem(fdf, 1));
-            if (not value)
-                DUNE_THROW(Dune::Exception, "Trying to construct DifferentiablePythonFunction with NULL as value argument");
-            if (not derivative)
-                DUNE_THROW(Dune::Exception, "Trying to construct DifferentiablePythonFunction with NULL as derivative argument");
-            value_ = value;
-            derivative_ = derivative;
-        }
-
-        /**
-         * \copydoc VirtualFunction::evaluate
-         */
-        void evaluate(const DomainType& x, RangeType& y) const
-        {
-            value_(x).toC(y);
-        }
-
-        /**
-         * \copydoc VirtualDifferentiableFunction::evaluateDerivative
-         */
-        void evaluateDerivative(const DomainType& x, DerivativeType& y) const
-        {
-            derivative_(x).toC(y);
-        }
-
-        /**
-         * \brief Evaluate function
-         */
-        RT operator()(const DT& x) const
-        {
-            return callable_(x).template toC<RT>();
-        }
-
-    protected:
-
-        Python::Callable value_;
-        Python::Callable derivative_;
-};
-
-
-namespace Python
-{
-
-// conversion of callable to PythonFunction*
-template<class DT, class RT>
-struct Conversion<Dune::VirtualFunction<DT, RT>*>
-{
-    enum {useDefaultConstructorConversion=true};
-    static void toC(PyObject* pyF, Dune::VirtualFunction<DT, RT>*& f)
-    {
-        // We don't want to steal a reference - hence we must call inc().
-        f = new PythonFunction<DT, RT>(Reference(Imp::inc(pyF)));
-    }
-};
-
-// conversion of callable to PythonFunction*
-template<class DT, class RT>
-struct Conversion<PythonFunction<DT, RT>*>
-{
-    enum {useDefaultConstructorConversion=true};
-    static void toC(PyObject* pyF, PythonFunction<DT, RT>*& f)
-    {
-        // We don't want to steal a reference - hence we must call inc().
-        f = new PythonFunction<DT, RT>(Reference(Imp::inc(pyF)));
-    }
-};
-
-// conversion of callable to BoundarySegment*
-template<int dim, int dimworld>
-struct Conversion<Dune::BoundarySegment<dim,dimworld>*>
-{
-    enum {useDefaultConstructorConversion=true};
-    static void toC(PyObject* pyF, Dune::BoundarySegment<dim,dimworld>*& segmentPointer)
-    {
-        typedef Dune::FieldVector<double,dim-1> DomainType;
-        typedef Dune::FieldVector<double,dimworld> RangeType;
-
-        // Wrap raw PyObject without stealing a reference
-        Reference f = Reference(Imp::inc(pyF));
-
-        typedef VirtualFunctionToBoundarySegmentAdapter<dim, dimworld> Adapter;
-        segmentPointer = new Adapter(typename Adapter::FunctionSharedPtr(new PythonFunction<DomainType, RangeType>(f)));
-    }
-};
-
-// conversion of sequence of callables to VirtualDifferentiableFunction*
-template<class DT, class RT>
-struct Conversion<VirtualDifferentiableFunction<DT, RT>*>
-{
-    enum {useDefaultConstructorConversion=true};
-    static void toC(PyObject* pyFDF, VirtualDifferentiableFunction<DT, RT>*& f)
-    {
-        // We don't want to steal a reference - hence we must call inc().
-        f = new DifferentiablePythonFunction<DT, RT>(Reference(Imp::inc(pyFDF)));
-    }
-};
-
-// conversion of sequence of callables to DifferentiablePythonFunction*
-template<class DT, class RT>
-struct Conversion<DifferentiablePythonFunction<DT, RT>*>
-{
-    enum {useDefaultConstructorConversion=true};
-    static void toC(PyObject* pyFDF, DifferentiablePythonFunction<DT, RT>*& f)
-    {
-        // We don't want to steal a reference - hence we must call inc().
-        f = new DifferentiablePythonFunction<DT, RT>(Reference(Imp::inc(pyFDF)));
-    }
-};
-
-} // end of namespace Python
-
-
-#else
-    #warning dunepython.hh was included but python was not found or enabled!
-#endif // DUNE_FUFEM_PYTHON_FUNCTION_HH
-
-
-#endif
-

dune/fufem/python/module.hh

@@ -17,16 +17,12 @@
 #include <dune/common/exceptions.hh>
 #include <dune/common/shared_ptr.hh>
 #include <dune/common/fvector.hh>
-#include <dune/common/function.hh>
 #include <dune/common/parametertree.hh>
 #include <dune/common/typetraits.hh>
+#include <dune/common/stringutility.hh>
 
 #include <dune/grid/common/gridfactory.hh>
 
-#include <dune/fufem/functions/virtualfunctiontoboundarysegmentadapter.hh>
-#include <dune/fufem/formatstring.hh>
-
-
 #include <dune/fufem/python/reference.hh>
 
 
@@ -198,7 +194,7 @@ class Module :
         {
             assertPyObject("Module::run()");
             if (not PyRun_String(code.c_str(), Py_file_input, dict_, dict_))
-                handlePythonError("Module::run()", Dune::Fufem::formatString("failed to execute the following code '%s'", code.c_str()));
+                handlePythonError("Module::run()", Dune::formatString("failed to execute the following code '%s'", code.c_str()));
         }
 
         /**

dune/fufem/python/reference.hh

@@ -14,8 +14,8 @@
 
 #include <dune/common/exceptions.hh>
 #include <dune/common/typetraits.hh>
+#include <dune/common/stringutility.hh>
 
-#include <dune/fufem/formatstring.hh>
 
 
 namespace Python
@@ -237,7 +237,7 @@ class Reference
             // This returns a new reference, i.e., it increments the reference count.
             PyObject* value = PyObject_GetAttrString(p_, name.c_str());
             if (not value)
-                handlePythonError("get()", Dune::Fufem::formatString("failed to get attribute '%s'", name.c_str()));
+                handlePythonError("get()", Dune::formatString("failed to get attribute '%s'", name.c_str()));
             return value;
         }
 
@@ -257,7 +257,7 @@ class Reference
         {
             assertPyObject("set()");
             if (PyObject_SetAttrString(p_, name.c_str(), makeObject(value))==-1)
-                handlePythonError("set()", Dune::Fufem::formatString("failed to set attribute '%s'", name.c_str()));
+                handlePythonError("set()", Dune::formatString("failed to set attribute '%s'", name.c_str()));
         }
 
         /**
@@ -277,7 +277,8 @@ class Reference
          * \returns The result of the call
          *
          */
-        Reference callWithArgumentTuple(const Reference& args) const DUNE_DEPRECATED_MSG("Replace ref.callWithArgumentTuple(args) by Callable(ref).callWithArgumentTuple(args).")
+        [[deprecated("Replace ref.callWithArgumentTuple(args) by Callable(ref).callWithArgumentTuple(args).")]]
+        Reference callWithArgumentTuple(const Reference& args) const
         {
             return callWithArgumentTupleImp(args);
         }
@@ -287,7 +288,8 @@ class Reference
          *
          * Shortcut for callWithArgumentTuple(makeTuple()).
          */
-        Reference operator() () const DUNE_DEPRECATED_MSG("Replace ref() by Callable(ref)().")
+        [[deprecated("Replace ref() by Callable(ref)().")]]
+        Reference operator() () const
         {
             return callWithArgumentTupleImp(makeTuple());
         }
@@ -298,7 +300,7 @@ class Reference
          * Shortcut for callWithArgumentTuple(makeTuple(t0)).
          */
         template<class T0>
-        DUNE_DEPRECATED_MSG("Replace ref(t0) by Callable(ref)(t0).")
+        [[deprecated("Replace ref(t0) by Callable(ref)(t0).")]]
         Reference operator() (const T0& t0) const
         {
             return callWithArgumentTupleImp(makeTuple(t0));
@@ -310,7 +312,7 @@ class Reference
          * Shortcut for callWithArgumentTuple(makeTuple(t0,t1)).
          */
         template<class T0, class T1>
-        DUNE_DEPRECATED_MSG("Replace ref(t0,t1) by Callable(ref)(t0,t1).")
+        [[deprecated("Replace ref(t0,t1) by Callable(ref)(t0,t1).")]]
         Reference operator() (const T0& t0, const T1& t1) const
         {
             return callWithArgumentTupleImp(makeTuple(t0, t1));
@@ -322,7 +324,7 @@ class Reference
          * Shortcut for callWithArgumentTuple(makeTuple(t0,t1,t2)).
          */
         template<class T0, class T1, class T2>
-        DUNE_DEPRECATED_MSG("Replace ref(t0,t1,t2) by Callable(ref)(t0,t1,t2).")
+        [[deprecated("Replace ref(t0,t1,t2) by Callable(ref)(t0,t1,t2).")]]
         Reference operator() (const T0& t0, const T1& t1, const T2& t2) const
         {
             return callWithArgumentTupleImp(makeTuple(t0, t1, t2));

dune/fufem/quadraturerules/compositequadraturerule.hh

@@ -23,9 +23,9 @@ class CompositeQuadratureRule:
         CompositeQuadratureRule(const Dune::QuadratureRule<ct,dim>& quad, int refinement)
         {
             
-            typedef Dune::StaticRefinement<Dune::Impl::SimplexTopology<dim>::type::id, 
+            typedef Dune::StaticRefinement<Dune::GeometryTypes::simplex(dim).id(),
                                            ct,
-                                           Dune::Impl::SimplexTopology<dim>::type::id,
+                                           Dune::GeometryTypes::simplex(dim).id(),
                                            dim> Refinement;
 
             int numberOfSubelements = (1<<(dim*refinement));

dune/fufem/quadraturerules/quadraturerulecache.hh

@@ -21,6 +21,15 @@ class QuadratureRuleKey
 {
     public:
 
+        /** \brief Default constructor, same as QuadratureRuleKey(0,0)
+         */
+        QuadratureRuleKey() :
+            gt_(Dune::GeometryTypes::none(0)),
+            order_(0),
+            refinement_(0),
+            lumping_(false)
+        {}
+
         /** \brief Create a key for a quadrature rule that can integrate a given local basis exactly
          */
         template<class LocalFiniteElement>
@@ -109,7 +118,7 @@ class QuadratureRuleKey
             return refinement_;
         }
 
-        DUNE_DEPRECATED_MSG("Please use setRefinement() as a replacement.")
+        [[deprecated("Please use setRefinement() as a replacement.")]]
         void refinement(int r)
         {
             setRefinement(r);

dune/fufem/readbitfield.hh

@@ -67,7 +67,7 @@ inline void readBitField(Dune::BitSetVector<ncomp>& field, int size, const std::
 
 //! Read level boundarypatch from an AmiraMesh file. */
 template <class GridType>
-DUNE_DEPRECATED_MSG("This method is deprecated. Please use the version with GridView as template.")
+[[deprecated("This method is deprecated. Please use the version with GridView as template.")]]
 void readBoundaryPatch(BoundaryPatch<typename GridType::LevelGridView>& patch,
                        const std::string& filename) {
     /** \todo Not very memory efficient! */

dune/fufem/referenceelementhelper.hh

@@ -28,7 +28,7 @@ struct ReferenceElementHelper
      * \param j      Number of possibly contained subentity
      * \param jCodim Codim of possibly contained subentity
      */
-    DUNE_DEPRECATED_MSG("Please referenceElement.subEntities(i,iCodim,jCodim).contains(j) instead")
+    [[deprecated("Please referenceElement.subEntities(i,iCodim,jCodim).contains(j) instead")]]
     static bool subEntityContainsSubEntity(const typename Dune::GeometryType& gt, int i, int iCodim, int j, int jCodim)
     {
         if (jCodim < iCodim)

dune/fufem/symmetricmatrix.hh

@@ -42,7 +42,7 @@ public:
   static constexpr SymmetricMatrix<T,N> identityMatrix()
   {
     SymmetricMatrix<T,N> id;
-    for( size_t i=0; i<N*(N+1)/2; ++i)
+    for( size_t i=0; i<N; ++i)
       id.setEntry(i,i,1);
 
     return id;

dune/fufem/test/CMakeLists.txt

-# Put your test in here if it needs access to external grids
-set(GRID_BASED_TESTS
-  basisgridfunctiontest
-  basisinterpolatortest
-  basisinterpolationmatrixtest
-  assembletransferoperatortest
-  boundarypatchtest
-  boundarypatchprolongatortest
-  coarsegridfunctionwrappertest
-  composedfunctiontest
-  constructboundarydofstest
-  dunefunctionsassemblertest
-  dunefunctionsipdgassemblertest
-  functionintegratortest
-  functionspacebasistest
-  generalizedlaplaceassemblertest
-  gradientassemblertest
-  gridconstructiontest
-  gridfunctiontest
-  gridfunctionadaptortest
-  h1functionalassemblertest
-  integraloperatorassemblertest
-  istlbackendtest
-  laplaceassemblertest
-  polynomialtest
-  secondorderassemblertest
-  subgridxyfunctionalassemblertest
-  sumfunctiontest
-  tensortest
-  transferoperatorassemblertest
-  )
-
-if (ADOLC_FOUND)
-    set(GRID_BASED_TESTS ${GRID_BASED_TESTS} adolctest)
-endif(ADOLC_FOUND)
-
-dune_add_test(SOURCES constantfunctiontest.cc)
+# Tests that should be run unconditionally
+dune_add_test(SOURCES assembletransferoperatortest.cc)
+dune_add_test(SOURCES basisgridfunctiontest.cc)
+dune_add_test(SOURCES basisinterpolatortest.cc)
+dune_add_test(SOURCES boundarypatchprolongatortest.cc)
+dune_add_test(SOURCES boundarypatchtest.cc)
+dune_add_test(SOURCES coarsegridfunctionwrappertest.cc)
+dune_add_test(SOURCES constructboundarydofstest.cc)
+dune_add_test(SOURCES dunefunctionsipdgassemblertest.cc)
+dune_add_test(SOURCES functionintegratortest.cc)
+dune_add_test(SOURCES functionspacebasistest.cc)
+dune_add_test(SOURCES generalizedlaplaceassemblertest.cc)
+dune_add_test(SOURCES gradientassemblertest.cc)
+dune_add_test(SOURCES gridconstructiontest.cc)
+dune_add_test(SOURCES gridfunctionadaptortest.cc)
+dune_add_test(SOURCES gridfunctiontest.cc)
+dune_add_test(SOURCES h1functionalassemblertest.cc)
+dune_add_test(SOURCES integraloperatorassemblertest.cc)
+dune_add_test(SOURCES istlbackendtest.cc)
+dune_add_test(SOURCES laplaceassemblertest.cc)
+dune_add_test(SOURCES localassemblertest.cc)
 dune_add_test(SOURCES makerefinedsimplexgeometrytest.cc)
+dune_add_test(SOURCES mappedmatrixtest.cc)
 dune_add_test(SOURCES newpfeassemblertest.cc)
 dune_add_test(SOURCES pgmtest.cc)
 dune_add_test(SOURCES ppmtest.cc)
 dune_add_test(SOURCES refinedsimplexgeometrytest.cc)
-dune_add_test(SOURCES test-polyhedral-minimisation.cc)
+dune_add_test(SOURCES secondorderassemblertest.cc)
+dune_add_test(SOURCES subgridxyfunctionalassemblertest.cc)
 dune_add_test(SOURCES symmetricmatrixtest.cc)
+dune_add_test(SOURCES tensortest.cc)
+dune_add_test(SOURCES test-polyhedral-minimisation.cc)
+dune_add_test(SOURCES transferoperatorassemblertest.cc)
+dune_add_test(SOURCES vintagebasisgridfunctiontest.cc)
 
-set(TESTS ${GRID_BASED_TESTS})
-
-# Setup targets, register tests, and add dune flags
-foreach(_test ${TESTS})
-  dune_add_test(SOURCES ${_test}.cc)
-endforeach()
-
-# Add external grid manager flags
-foreach(_test ${GRID_BASED_TESTS})
-  if(HAVE_UG)
-    add_dune_ug_flags(${_test})
-  endif(HAVE_UG)
-endforeach()
-
-if(PYTHONLIBS_FOUND)
+# PYTHONLIBS_FOUND is just placed for backward compatibility with 2.7 Core tests
+# and can be removed once tests against 2.7 Core are disabled
+if (Python3_FOUND OR PYTHONLIBS_FOUND)
   dune_add_test(SOURCES dunepythontest.cc)
 
   # Add python flags to corresponding test
@@ -68,18 +46,14 @@ if(PYTHONLIBS_FOUND)
 endif()
 
 if (ADOLC_FOUND)
+    dune_add_test(SOURCES adolctest.cc)
     add_dune_adolc_flags(adolctest)
 endif()
 
-# disable test involving boos for clang with c++17
-# In the standard, auto_ptr got removed and while there is a flag
-# to tell boost not to use auto_ptr anymore, this leads to undefined reference linker error
-if (NOT "${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" OR CMAKE_CXX_COMPILER_VERSION VERSION_LESS 6)
 if(Boost_SERIALIZATION_FOUND)
   dune_add_test(SOURCES serializationtest.cc)
   set_property(TARGET serializationtest APPEND PROPERTY INCLUDE_DIRECTORIES ${Boost_INCLUDE_DIRS})
-        target_link_libraries(serializationtest ${Boost_SERIALIZATION_LIBRARY})
-      endif()
+  target_link_libraries(serializationtest PUBLIC ${Boost_SERIALIZATION_LIBRARY})
 endif()
 
 if (HDF5_FOUND)

dune/fufem/test/assembletransferoperatortest.cc

@@ -7,15 +7,36 @@
 #include <dune/functions/functionspacebases/defaultglobalbasis.hh>
 #include <dune/functions/functionspacebases/lagrangebasis.hh>
 #include <dune/functions/functionspacebases/raviartthomasbasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
 #include <dune/functions/functionspacebases/interpolate.hh>
 #include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
 
 #include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
-
-#include "common.hh"
+#include "dune/fufem/test/common.hh"
 
 using namespace Dune;
 
+template<class Matrix, class Function, class FineBasis, class FineCoeff, class CoarseBasis, class CoarseCoeff>
+double checkInterpolation(Matrix& matrix, const Function& f,
+                        const FineBasis& fineBasis, FineCoeff& fineCoeffs,
+                        const CoarseBasis& coarseBasis, CoarseCoeff& coarseCoeffs)
+{
+  // assemble the matrix
+  assembleGlobalBasisTransferMatrix(matrix,coarseBasis,fineBasis);
+
+  // interpolate f on the coarse grid
+  Functions::interpolate(coarseBasis,coarseCoeffs,f);
+
+  // interpolate f on the fine grid
+  Functions::interpolate(fineBasis,fineCoeffs,f);
+
+  // now we should have
+  // fineCoeffs - matrix * coarseCoeffs == 0
+  matrix.mmv(coarseCoeffs,fineCoeffs);
+
+  return fineCoeffs.two_norm2();
+}
+
 template<int order, class GW0, class GW1>
 bool testLagrange(const GW0& gw0, const GW1& gw1)
 {
@@ -39,22 +60,114 @@ bool testLagrange(const GW0& gw0, const GW1& gw1)
   BCRSMatrix<FieldMatrix<double,1,1>> matrix;
   BlockVector<FieldVector<double,1>> coarseCoeffs, fineCoeffs;
 
-  // assemble the matrix
-  assembleGlobalBasisTransferMatrix(matrix,coarseLagrange,fineLagrange);
+  const auto diff = checkInterpolation(matrix, f, fineLagrange, fineCoeffs, coarseLagrange, coarseCoeffs);
 
-  // interpolate f on the coarse grid
-  Functions::interpolate(coarseLagrange,coarseCoeffs,f);
+  std::cout << " test Lagrange " << GW0::dimension << "D with order " << order  << " : error = " << diff << std::endl;
 
-  // interpolate f on the fine grid
-  Functions::interpolate(fineLagrange,fineCoeffs,f);
+  return diff < 1e-12;
+}
 
-  // now we should have
-  // fineCoeffs - matrix * coarseCoeffs == 0
-  matrix.mmv(coarseCoeffs,fineCoeffs);
 
-  std::cout << " test Lagrange " << GW0::dimension << "D with order " << order  << " : error = " << fineCoeffs.two_norm2()<< std::endl;
+template<int order, class GW0, class GW1>
+bool testPowerBasis(const GW0& gw0, const GW1& gw1)
+{
+  // our test funcion for interpolation
+  auto f = [](const auto& x)
+  {
+    // const function for order 0
+    if ( order == 0 )
+      return FieldVector<double,3>{1.,1.,1.};
 
-  return fineCoeffs.two_norm2() < 1e-12;
+    // for the other cases affine function
+    FieldVector<double,3> y;
+    y = 1.;
+    for( auto xs : x )
+    {
+      y[0] += xs;
+      y[1] += 2*xs;
+      y[2] += 3*xs;
+    }
+    return y;
+  };
+
+  using namespace Functions::BasisFactory;
+  bool passed = true;
+  {
+    auto coarseLagrange = makeBasis(gw0, power<3>(lagrange<order>(),blockedInterleaved()));
+    auto fineLagrange   = makeBasis(gw1, power<3>(lagrange<order>(),blockedInterleaved()));
+
+    BCRSMatrix<FieldMatrix<double,3,3>> matrix;
+    BlockVector<FieldVector<double,3>> coarseCoeffs, fineCoeffs;
+
+    const auto diff = checkInterpolation(matrix, f, fineLagrange, fineCoeffs, coarseLagrange, coarseCoeffs);
+
+    std::cout << " test PowerBasis blockedInterleaved " << GW0::dimension << "D with order " << order  << " : error = " << diff << std::endl;
+
+    passed = passed and diff < 1e-12;
+  }
+  {
+    auto coarseLagrange = makeBasis(gw0, power<3>(lagrange<order>(),flatInterleaved()));
+    auto fineLagrange   = makeBasis(gw1, power<3>(lagrange<order>(),flatInterleaved()));
+
+    BCRSMatrix<FieldMatrix<double,1,1>> matrix;
+    BlockVector<FieldVector<double,1>> coarseCoeffs, fineCoeffs;
+
+    const auto diff = checkInterpolation(matrix, f, fineLagrange, fineCoeffs, coarseLagrange, coarseCoeffs);
+
+    std::cout << " test PowerBasis flatInterleaved " << GW0::dimension << "D with order " << order  << " : error = " << diff << std::endl;
+
+    passed = passed and diff < 1e-12;
+  }
+  {
+    auto coarseLagrange = makeBasis(gw0, power<3>(lagrange<order>(),flatLexicographic()));
+    auto fineLagrange   = makeBasis(gw1, power<3>(lagrange<order>(),flatLexicographic()));
+
+    BCRSMatrix<FieldMatrix<double,1,1>> matrix;
+    BlockVector<FieldVector<double,1>> coarseCoeffs, fineCoeffs;
+
+    const auto diff = checkInterpolation(matrix, f, fineLagrange, fineCoeffs, coarseLagrange, coarseCoeffs);
+
+    std::cout << " test PowerBasis flatLexicographic " << GW0::dimension << "D with order " << order  << " : error = " << diff << std::endl;
+
+    passed = passed and diff < 1e-12;
+  }
+
+  return passed;
+}
+
+
+
+
+template<int order0, int order1, class GW0, class GW1>
+bool testOrder(const GW0& gw0, const GW1& gw1)
+{
+  // our test funcion for interpolation
+  auto f = [](const auto& x)
+  {
+    FieldVector<double,3> y;
+    y = 1.;
+    for( auto xs : x )
+    {
+      y[0] += std::sin(xs);
+      y[1] += 2*xs +5.;
+      y[2] += 3*xs*xs;
+    }
+    return y;
+  };
+
+  using namespace Functions::BasisFactory;
+
+  auto lagrange0 = makeBasis(gw0, power<3>(lagrange<order0>(),blockedInterleaved()));
+  auto lagrange1 = makeBasis(gw1, power<3>(lagrange<order1>(),blockedInterleaved()));
+
+  BCRSMatrix<FieldMatrix<double,3,3>> matrix;
+  BlockVector<FieldVector<double,3>> coeffs0, coeffs1;
+
+  const auto diff = checkInterpolation(matrix, f, lagrange0, coeffs0, lagrange1, coeffs1);
+
+  std::cout << " test interpolation P" << order0 << " -> P" << order1 << " : error = " << diff << std::endl;
+
+  return diff < 1e-12;
 }
 
 
@@ -99,9 +212,7 @@ bool testRaviartThomas(const GW0& gw0, const GW1& gw1)
     FieldVector<double,2> testPoint(0);
     auto testPointInFather = element.geometryInFather().global(testPoint);
 
-    // TODO: The run-time test is disabled for the time being, because it requires
-    // https://gitlab.dune-project.org/staging/dune-functions/-/merge_requests/241
-    if ( false && (coarseLocalFunction(testPointInFather) - fineLocalFunction(testPoint) ).two_norm() > 1e-10 )
+    if ( (coarseLocalFunction(testPointInFather) - fineLocalFunction(testPoint) ).two_norm() > 1e-10 )
     {
       std::cerr << "Raviart-Thomas: transfer matrix is not the natural injection!" << std::endl;
       return false;
@@ -156,17 +267,33 @@ bool checkTransferMatrix(const GW0& gw0, const GW1& gw1)
   bool passed = true;
 
   // test Lagrange elements for orders 0-3
+  // for scalar bases ...
   passed = passed and testLagrange<0>(gw0,gw1);
   passed = passed and testLagrange<1>(gw0,gw1);
   passed = passed and testLagrange<2>(gw0,gw1);
+  // ... and powerBases
+  passed = passed and testPowerBasis<0>(gw0,gw1);
+  passed = passed and testPowerBasis<1>(gw0,gw1);
+  passed = passed and testPowerBasis<2>(gw0,gw1);
   // order 3 only for 2D grids
-  if ( GW0::dimension < 3 )
+  if constexpr ( GW0::dimension < 3 )
+  {
     passed = passed and testLagrange<3>(gw0,gw1);
+    passed = passed and testPowerBasis<3>(gw0,gw1);
+  }
 
   // test transfer operators for Raviart-Thomas bases
   if constexpr ( GW0::dimension==2 )
     passed = passed and testRaviartThomas<0>(gw0,gw1);
 
+    // 2D: check lagrange interpolation for different orders on the same grid
+  if constexpr ( GW0::dimension < 3 )
+  {
+    passed = passed and testOrder<1,2>(gw1,gw1);
+    passed = passed and testOrder<2,3>(gw1,gw1);
+    passed = passed and testOrder<1,3>(gw1,gw1);
+  }
+
   return passed;
 }
 

dune/fufem/test/basisinterpolationmatrixtest.cc

-#include <config.h>
-
-#include <dune/common/parallel/mpihelper.hh>
-
-#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
-
-#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/functiontools/basisinterpolator.hh>
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/functionspacebases/p2nodalbasis.hh>
-#include <dune/fufem/functionspacebases/p3nodalbasis.hh>
-#include <dune/fufem/functionspacebases/refinedp1nodalbasis.hh>
-#include <dune/fufem/functionspacebases/conformingbasis.hh>
-#include "common.hh"
-
-template<class DT, class RT>
-class F:
-    public Dune::VirtualFunction<DT, RT>
-{
-    public:
-        void evaluate(const DT& x, RT& y) const
-        {
-            y = 1;
-        }
-};
-
-
-struct Suite
-{
-    template<class GridType>
-    bool check(const GridType& grid)
-    {
-        bool passed = true;
-
-        typedef P1NodalBasis<typename GridType::LeafGridView> P1Basis;
-        P1Basis p1Basis(grid.leafGridView());
-        typedef P2NodalBasis<typename GridType::LeafGridView> P2Basis;
-        P2Basis p2Basis(grid.leafGridView());
-        typedef P3NodalBasis<typename GridType::LeafGridView> P3Basis;
-        P3Basis p3Basis(grid.leafGridView());
-        typedef RefinedP1NodalBasis<typename GridType::LeafGridView> RefinedP1Basis;
-        RefinedP1Basis refP1Basis(grid.leafGridView());
-
-        typedef ConformingBasis<P1Basis> ConfP1Basis;
-        ConfP1Basis confP1Basis(p1Basis);
-        typedef ConformingBasis<RefinedP1Basis> ConfRefinedP1Basis;
-        ConfRefinedP1Basis confRefP1Basis(refP1Basis);
-
-        std::cout<<"Checking Refined P1 -> P1 \n";
-        passed = passed and checkInterpolation<P1Basis, RefinedP1Basis,3>(p1Basis,refP1Basis);
-
-        if (!passed)
-            return false;
-
-        std::cout<<"Checking P2 -> P1 \n";
-        passed = passed and checkInterpolation<P1Basis, P2Basis,2>(p1Basis,p2Basis);
-
-        if (!passed)
-            return false;
-
-        std::cout<<"Checking P3 -> P1 \n";
-        passed = passed and checkInterpolation<P1Basis, P3Basis,1>(p1Basis,p3Basis);
-
-        if (!passed)
-            return false;
-
-        std::cout<<"Checking Conforming<P1> -> P1 \n";
-        passed = passed and checkInterpolation<P1Basis,ConfP1Basis,3> (p1Basis, confP1Basis);
-
-        if (!passed)
-            return false;
-
-        std::cout<<"Checking P1 -> Conforming<P1> \n";
-        passed = passed and checkInterpolation<ConfP1Basis,P1Basis,3> (confP1Basis, p1Basis);
-
-        if (!passed)
-            return false;
-
-        std::cout<<"Checking Conforming<RefinedP1> -> Conforming<P1> \n";
-        passed = passed and checkInterpolation<ConfP1Basis,ConfRefinedP1Basis,1> (confP1Basis, confRefP1Basis);
-
-        if (!passed)
-            return false;
-
-        std::cout<<"Checking P3 -> P2 \n";
-        passed = passed and checkInterpolation<P2Basis, P3Basis,2>(p2Basis,p3Basis);
-
-        return passed;
-    }
-
-    template<class BasisType0, class BasisType1, int ncomp>
-    bool checkInterpolation(const BasisType0& coarseBase, const BasisType1& fineBase)
-    {
-        typedef Dune::FieldVector<double, BasisType0::GridView::dimensionworld> DT;
-        typedef Dune::FieldVector<double, ncomp> RT;
-
-        typedef Dune::BlockVector<RT> VectorType;
-        typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,ncomp,ncomp> > MatrixType;
-
-        // Constant 1 function
-        F<DT, RT> f;
-
-        // interpolate in coarse basis
-        VectorType v0;
-        Functions::interpolate(coarseBase, v0, f);
-
-        // interpolate from coarse to fine
-        VectorType v1;
-        Functions::interpolate(fineBase, v1, Functions::makeFunction(coarseBase, v0));
-
-        // setup interpolation matrix
-        MatrixType interpolMat;
-        assembleBasisInterpolationMatrix(interpolMat, coarseBase, fineBase);
-
-        // interpolate again
-        VectorType v2(interpolMat.N());
-        interpolMat.mv(v0,v2);
-
-
-        for (size_t i=0; i<v1.size(); ++i)
-            if ((v1[i]-v2[i]).two_norm()>1e-14)
-                return false;
-
-        return (v1.size() == v2.size());
-    }
-
-};
-
-int main(int argc, char** argv)
-{
-    Dune::MPIHelper::instance(argc, argv);
-
-    Suite tests;
-
-    bool passed = checkWithStandardAdaptiveGrids(tests);
-
-    return passed ? 0 : 1;
-}

dune/fufem/test/boundarypatchtest.cc

@@ -47,7 +47,7 @@ struct BoundaryPatchTestSuite
 
         // Test copy construction
         BP fooBoundary = boundary;
-        int c DUNE_UNUSED;
+        [[maybe_unused]] int c;
         c = fooBoundary.gridView().indexSet().size(0);   // make the copy do something useful
 
         // Test assignment
@@ -70,7 +70,7 @@ struct BoundaryPatchTestSuite
         // Test leaf copy construction
         BoundaryPatch<typename GridType::LeafGridView> leafBoundary(grid.leafGridView(), true);
         BoundaryPatch<typename GridType::LeafGridView> foo = leafBoundary;
-        int c DUNE_UNUSED;
+        [[maybe_unused]] int c;
         c = foo.gridView().indexSet().size(0);   // make the copy do something useful
 
 

dune/fufem/test/common.hh

@@ -12,7 +12,7 @@
 #include <dune/grid/yaspgrid.hh>
 #include <dune/grid/common/gridfactory.hh>
 
-#if HAVE_UG
+#if HAVE_DUNE_UGGRID
 #include <dune/grid/uggrid.hh>
 #include <dune/grid/uggrid/uggridfactory.hh>
 #endif
@@ -70,6 +70,37 @@ bool isCloseAbs(FT t1, Other t2)
   return std::abs<FT>(t1 - t2) < ToleranceTraits<FT>::tol;
 }
 
+template<class GridView, class FunctionA, class SignatureB>
+bool compareEvaluateByGridViewQuadrature(const FunctionA& fA, const std::function<SignatureB>& fB, const GridView& gridView, int order)
+{
+    typedef typename FunctionA::RangeType RangeTypeA;
+
+    constexpr int dim = GridView::dimension;
+
+    for(const auto& it : elements(gridView))
+    {
+        const auto& quad = Dune::QuadratureRules<typename GridView::ctype, dim>::rule(it.type(), order);
+        const auto& geometry = it.geometry();
+        for (const auto& pt : quad)
+        {
+            auto x = geometry.global(pt.position());
+            RangeTypeA yA;
+
+            fA.evaluate(x,yA);
+            auto yB = fB(x);
+            yA -= yB;
+
+            if (yA.infinity_norm() > 1e-12)
+            {
+                std::cout << "Result of evaluate() differs at global coordinate " << x << std::endl;
+                return false;
+            }
+        }
+    }
+    return true;
+}
+
+
 template<class GridView, class FunctionA, class FunctionB>
 bool compareEvaluateByGridViewQuadrature(const FunctionA& fA, const FunctionB& fB, const GridView& gridView, int order)
 {
@@ -103,7 +134,6 @@ bool compareEvaluateByGridViewQuadrature(const FunctionA& fA, const FunctionB& f
 }
 
 
-
 template<class GridView, class FunctionA, class FunctionB>
 bool compareEvaluateDerivativeByGridViewQuadrature(const FunctionA& fA, const FunctionB& fB, const GridView& gridView, int order)
 {
@@ -251,7 +281,7 @@ std::unique_ptr<GridType> constructCoarseGrid()
 }
 
 template<int dim>
-Dune::YaspGrid<dim>* constructCoarseYaspGrid()
+auto constructCoarseYaspGrid()
 {
     typedef Dune::YaspGrid<dim> GridType;
 
@@ -259,7 +289,7 @@ Dune::YaspGrid<dim>* constructCoarseYaspGrid()
     std::array<int,dim> s;
     std::fill(s.begin(), s.end(), 1);
 
-    GridType* grid = new GridType(L,s);
+    auto grid = std::make_unique<GridType>(L,s);
 
     return grid;
 }
@@ -306,8 +336,7 @@ bool checkWithAdaptiveGrid(TestSuite& suite, int uniformRefinements, int localRe
 template<int dim, class TestSuite>
 bool checkWithStructuredGrid(TestSuite& suite, int uniformRefinements)
 {
-    typedef Dune::YaspGrid<dim> GridType;
-    GridType* grid(constructCoarseYaspGrid<dim>());
+    auto  grid = constructCoarseYaspGrid<dim>();
 
     grid->globalRefine(uniformRefinements);
 
@@ -317,7 +346,6 @@ bool checkWithStructuredGrid(TestSuite& suite, int uniformRefinements)
     if (passed)
         std::cout << "All tests passed with " << Dune::className(grid) << "." << std::endl;
 
-    delete grid;
     return passed;
 }
 
@@ -379,13 +407,9 @@ bool checkWithStandardAdaptiveGrids(TestSuite& suite)
 {
     bool passed = true;
 
-#if HAVE_UG
+#if HAVE_DUNE_UGGRID
     passed = passed and checkWithAdaptiveGrid<Dune::UGGrid<2> >(suite, 3, 3);
-// This is disabled due to a bug (https://gitlab.dune-project.org/core/dune-grid/issues/27)
-// in parallel UGGrid. Once this is fixed we can reenable this in favour of the check
-// with less refinement below.
-//    passed = passed and checkWithAdaptiveGrid<Dune::UGGrid<3> >(suite, 1, 2);
-    passed = passed and checkWithAdaptiveGrid<Dune::UGGrid<3> >(suite, 1, 1);
+    passed = passed and checkWithAdaptiveGrid<Dune::UGGrid<3> >(suite, 1, 2);
 #if HAVE_DUNE_SUBGRID
     passed = passed and checkWithSubGrid<Dune::UGGrid<2> >(suite, 3, 3, "SubGrid< 2,UGGrid<2> >");
     passed = passed and checkWithSubGrid<Dune::UGGrid<3> >(suite, 1, 2, "SubGrid< 3,UGGrid<3> >");

dune/fufem/test/composedfunctiontest.cc

-#include <config.h>
-
-#include <cmath>
-#include <cstdio>
-
-#include <dune/common/parallel/mpihelper.hh>
-#include <dune/common/exceptions.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/function.hh>
-
-
-#include <dune/fufem/functions/composedfunction.hh>
-#include <dune/fufem/functions/composedgridfunction.hh>
-#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/functions/virtualgridfunction.hh>
-#include <dune/fufem/functions/virtualdifferentiablefunction.hh>
-
-#include <dune/fufem/functiontools/basisinterpolator.hh>
-
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-
-#include "common.hh"
-
-/** \brief inflates by zeros or cuts off trailing components according to Domain and Range dimensions. for domaindim==rangedim it's simply the identity
- *
- */
-template <class DT, class RT>
-class InflateOrCut_Base:
-    public Dune::VirtualFunction<DT,RT>
-{
-        typedef Dune::VirtualFunction<DT,RT> BaseType;
-
-    public:
-        typedef typename BaseType::DomainType DomainType;
-        typedef typename BaseType::RangeType RangeType;
-
-        enum {DomainDim = DomainType::dimension,
-              RangeDim  = RangeType::dimension};
-
-        InflateOrCut_Base()
-        {}
-
-        virtual void evaluate(const DomainType& x, RangeType& y) const
-        {
-            y = 0.0;
-
-            for (int i=0; i<DomainDim and i<RangeDim; ++i)
-            {
-                y[i]=x[i];
-            }
-
-            return;
-        }
-};
-
-/** \brief inflates by zeros or cuts off trailing components according to Domain and Range dimensions. for domaindim==rangedim it's simply the identity
- *
- */
-template <class GridType, class RT>
-class InflateOrCut_GridF:
-    public VirtualGridFunction<GridType,RT>
-{
-        typedef VirtualGridFunction<GridType,RT> BaseType;
-
-    public:
-        typedef typename BaseType::LocalDomainType LocalDomainType;
-        typedef typename BaseType::DomainType DomainType;
-        typedef typename BaseType::RangeType RangeType;
-
-        typedef typename BaseType::Element Element;
-
-        enum {DomainDim = DomainType::dimension,
-              RangeDim  = RangeType::dimension};
-
-        InflateOrCut_GridF(GridType& grid):
-            BaseType(grid)
-        {}
-
-        virtual void evaluateLocal(const Element& e, const LocalDomainType& x, RangeType& y) const
-        {
-            y = 0.0;
-
-            DomainType xglobal = e.geometry().global(x);
-
-            for (int i=0; i<DomainDim and i<RangeDim; ++i)
-            {
-                y[i]=xglobal[i];
-            }
-
-            return;
-        }
-
-        virtual bool isDefinedOn(const Element&) const
-        {
-            return true;
-        }
-
-};
-/** \brief inflates by zeros or cuts off trailing components according to Domain and Range dimensions. for domaindim==rangedim it's simply the identity
- *
- */
-template <class GridViewType, class RT>
-class InflateOrCut_GridVF:
-    public VirtualGridViewFunction<GridViewType,RT>
-{
-        typedef VirtualGridViewFunction<GridViewType,RT> BaseType;
-
-    public:
-        typedef typename BaseType::LocalDomainType LocalDomainType;
-        typedef typename BaseType::DomainType DomainType;
-        typedef typename BaseType::RangeType RangeType;
-
-        typedef typename BaseType::Element Element;
-
-        enum {DomainDim = DomainType::dimension,
-              RangeDim  = RangeType::dimension};
-
-        InflateOrCut_GridVF(const GridViewType& gridview):
-            BaseType(gridview)
-        {}
-
-        virtual void evaluateLocal(const Element& e, const LocalDomainType& x, RangeType& y) const
-        {
-            y = 0.0;
-
-            DomainType xglobal = e.geometry().global(x);
-
-            for (int i=0; i<DomainDim and i<RangeDim; ++i)
-            {
-                y[i]=xglobal[i];
-            }
-
-            return;
-        }
-};
-
-/** \brief TestSuite for ComposedFunction
- *
- *  Takes InflateOrCut_Base functions f,g such that \f$f\circ g:\mathbb R^n\longrightarrow\mathbb R^m\longrightarrow \mathbb R^d\f$ with \f$m > n \geq d\f$. So on the first d components
- *  the composition should be the identity.
- */
-struct ComposedFunctionTestSuite
-{
-    bool check()
-    {
-        const std::size_t n=3,
-                  m=5,
-                  d=2;
-
-        std::size_t n_testpoints=100;
-
-        typedef Dune::FieldVector<double,n> DomainType;
-        typedef Dune::FieldVector<double,m> IntermediateRangeType;
-        typedef Dune::FieldVector<double,d> RangeType;
-
-        typedef InflateOrCut_Base<DomainType,IntermediateRangeType> InnerFunctionType;
-        typedef InflateOrCut_Base<IntermediateRangeType,RangeType> OuterFunctionType;
-
-        InnerFunctionType g;
-        OuterFunctionType f;
-
-        ComposedFunction<DomainType,IntermediateRangeType,RangeType> f_after_g(f,g);
-
-        DomainType x(0.0);
-        RangeType  y(0.0);
-
-        for (std::size_t run=0; run<n_testpoints; ++run)
-        {
-            for (std::size_t dcomp=0; dcomp<n; ++dcomp)
-                x[dcomp] = (5.0*rand())/RAND_MAX - 2.5;
-
-            f_after_g.evaluate(x,y);
-
-            for (std::size_t rcomp=0; rcomp<d; ++rcomp)
-                if (std::abs(y[rcomp]-x[rcomp])>1e-15)
-                {
-                    std::cout << "y[" << rcomp << "] = " << y[rcomp] << ", x[" << rcomp << "] = " << x[rcomp] << std::endl;
-
-                    return false;
-                }
-        }
-
-        return true;
-    }
-};
-/** \brief TestSuite for ComposedFunction
- *
- *  Takes InflateOrCut functions f,g such that \f$f\circ g:\mathbb R^n\longrightarrow\mathbb R^m\longrightarrow \mathbb R^d\f$ with \f$m > n,\ m > d\f$. So on the first min(n,d) components
- *  the composition should be the identity.
- */
-struct ComposedGridFunctionTestSuite
-{
-    template <class GridType>
-    bool check(GridType& grid)
-    {
-        const unsigned int n=GridType::dimension,
-                  m=5,
-                  d=4;
-
-        const std::size_t n_testpoints=100;
-
-        typedef Dune::FieldVector<typename GridType::ctype,n> DomainType;
-        typedef Dune::FieldVector<double,m> IntermediateRangeType;
-        typedef Dune::FieldVector<double,d> RangeType;
-
-        typedef InflateOrCut_GridF<GridType,IntermediateRangeType> InnerFunctionType;
-        typedef InflateOrCut_GridVF<typename GridType::LeafGridView,IntermediateRangeType> InnerFunctionType2;
-        typedef InflateOrCut_Base<IntermediateRangeType,RangeType> OuterFunctionType;
-
-        InnerFunctionType g(grid);
-        InnerFunctionType2 gg(grid.leafGridView());
-        OuterFunctionType f;
-
-        ComposedGridFunction<GridType,IntermediateRangeType,RangeType> f_after_g(f,g);
-        ComposedGridFunction<GridType,IntermediateRangeType,RangeType> f_after_gg(f,gg);
-
-        DomainType x(0.0);
-        RangeType  y(0.0),yy(0.0);
-
-        for (std::size_t run=0; run<n_testpoints; ++run)
-        {
-            double sum=0;
-            for (std::size_t dcomp=0; dcomp<n; ++dcomp)
-            {
-                x[dcomp] = std::min((1.0*rand())/(n-1)/RAND_MAX, 1-sum);
-                sum += x[dcomp];
-            }
-
-            f_after_g.evaluate(x,y); // let's just use evaluate. It's the default implementation using evaluateLocal internally so the GridFunction aspect should be checked well enough
-            f_after_gg.evaluate(x,yy);
-
-            std::size_t kmin = std::min(d,n);
-
-            for (std::size_t rcomp=0; rcomp<kmin; ++rcomp)
-                if (std::abs(y[rcomp]-x[rcomp])>1e-15 or std::abs(yy[rcomp]-x[rcomp])>1e-15)
-                {
-                    std::cout << "y[" << rcomp << "] = " << y[rcomp] << ", x[" << rcomp << "] = " << x[rcomp] << std::endl;
-
-                    return false;
-                }
-
-            for (std::size_t rcomp=kmin; rcomp<d; ++rcomp)
-                if (y[rcomp]!=0.0 or yy[rcomp]!=0.0)
-                    return false;
-        }
-
-        return true;
-    }
-};
-
-int main(int argc, char** argv)
-{
-    Dune::MPIHelper::instance(argc, argv);
-
-    std::cout << "This is the ComposedFunctionTest" << std::endl;
-
-    bool passed = true;
-
-    {
-        ComposedFunctionTestSuite cftests;
-
-        passed = cftests.check();
-        if (passed)
-            std::cout << "ComposedFunction test passed" << std::endl;
-        else
-            std::cout << "ComposedFunction test failed" << std::endl;
-    }
-
-    {
-        ComposedGridFunctionTestSuite cftests;
-        passed = checkWithStandardAdaptiveGrids(cftests);
-    }
-
-    return passed ? 0 : 1;
-
-}

dune/fufem/test/constantfunctiontest.cc

-#include <config.h>
-
-#include <type_traits>
-
-#include <dune/common/parallel/mpihelper.hh>
-#include <dune/common/exceptions.hh>
-
-#include <dune/fufem/functions/constantfunction.hh>
-
-using namespace Dune;
-
-template<class DT, class RT>
-struct ConstantFunctionTestSuite
-{
-    private:
-        typedef typename ConstantFunction<DT, RT>::DerivativeType D;
-
-        static bool checkEvaluate(const ConstantFunction<DT,RT>& constantFunction, const RT& constantFunctionValue)
-        {
-            DT x(5);
-            RT y(0);
-
-            constantFunction.evaluate(x,y);
-
-            return (y == constantFunctionValue);
-        }
-
-        static bool checkEvaluateDerivative(const ConstantFunction<DT,RT>& constantFunction)
-        {
-            DT x(5);
-            D d;
-            constantFunction.evaluateDerivative(x,d);
-
-            return checkDerivativeEqualsZero(d);
-        }
-
-        static bool checkDerivativeEqualsZero(DerivativeTypeNotImplemented& d)
-        {
-            return true;
-        }
-
-        template<class T, typename std::enable_if<not(std::is_same<T, DerivativeTypeNotImplemented>::value),int >::type = 0>
-        static bool checkDerivativeEqualsZero(T& d)
-        {
-            return d == D(0);
-        }
-
-    public:
-        static bool check()
-        {
-            bool passed = true;
-
-            const RT constantFunctionValue(RT(1));
-            const ConstantFunction<DT, RT> constantFunction(constantFunctionValue);
-
-            passed = passed and checkEvaluate(constantFunction,constantFunctionValue);
-            passed = passed and checkEvaluateDerivative(constantFunction);
-
-            return passed;
-        }
-};
-
-int main (int argc, char *argv[])
-{
-    Dune::MPIHelper::instance(argc, argv);
-
-    bool passed = true;
-
-    passed = passed and ConstantFunctionTestSuite< double, double >::check();
-    passed = passed and ConstantFunctionTestSuite< float, double >::check();
-    passed = passed and ConstantFunctionTestSuite< double, int >::check();
-
-    passed = passed and ConstantFunctionTestSuite< Dune::FieldVector<double,2>, Dune::FieldVector<double,2> >::check();
-    passed = passed and ConstantFunctionTestSuite< Dune::FieldVector<double,3>, Dune::FieldVector<double,2> >::check();
-    passed = passed and ConstantFunctionTestSuite< Dune::FieldVector<double,2>, Dune::FieldVector<double,3> >::check();
-    passed = passed and ConstantFunctionTestSuite< Dune::FieldMatrix<double,2,2>, Dune::FieldVector<double,2> >::check();
-    passed = passed and ConstantFunctionTestSuite< Dune::FieldVector<double,2>, Dune::FieldMatrix<double,1,1> >::check();
-    passed = passed and ConstantFunctionTestSuite< Dune::FieldMatrix<double,1,1>, Dune::FieldMatrix<double,1,1> >::check();
-
-    passed = passed and ConstantFunctionTestSuite< Dune::FieldVector<float,3>, Dune::FieldVector<double,2> >::check();
-    passed = passed and ConstantFunctionTestSuite< Dune::FieldVector<float,3>, Dune::FieldVector<float,2> >::check();
-    passed = passed and ConstantFunctionTestSuite< Dune::FieldVector<double,3>, Dune::FieldVector<float,2> >::check();
-    passed = passed and ConstantFunctionTestSuite< Dune::FieldVector<double,3>, Dune::FieldVector<int,2> >::check();
-
-    if (not passed)
-        std::cout << "ConstantFunctionTest failed." << std::endl;
-
-    return passed ? 0 : 1;
-}

dune/fufem/test/dunepythontest.cc

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
 #include <config.h>
 
 #ifdef HAVE_PYTHON
@@ -10,10 +12,11 @@
 
 #include <dune/common/parallel/mpihelper.hh>
 #include <dune/common/fvector.hh>
-#include <dune/common/function.hh>
 #include <dune/common/parametertree.hh>
+#include <dune/common/tuplevector.hh>
+#include <dune/common/indices.hh>
 
-#if HAVE_UG
+#if HAVE_DUNE_UGGRID
 #include <dune/grid/uggrid.hh>
 #include <dune/grid/uggrid/uggridfactory.hh>
 #include <dune/grid/common/gridfactory.hh>
@@ -25,6 +28,7 @@
 #include <dune/fufem/dunepython.hh>
 
 
+using namespace Dune::Indices;
 
 int main(int argc, char** argv)
 {
@@ -373,6 +377,15 @@ int main(int argc, char** argv)
         // You can also call a lambda function
         auto lambda = Python::Callable(Python::evaluate("lambda x,y: (x+y)"));
         std::cout << "Result of call to lambda function: " << lambda(foo, bar).str() << std::endl;
+
+        // Named arguments can also be passed and mixed with positional ones
+        using namespace Python::Literals;
+        using Python::arg;
+        auto checkIncreasing = Python::Callable(module.get("checkIncreasing"));
+        std::cout << "Result of checkIncreasing(0,1): " << checkIncreasing(0,1).str() << std::endl;
+        std::cout << "Result of checkIncreasing(0,1,reverse=False): " << checkIncreasing(0, 1, arg("reverse", false)).str() << std::endl;
+        std::cout << "Result of checkIncreasing(0,1,reverse=True): " << checkIncreasing(0, 1, "reverse"_a=true).str() << std::endl;
+        std::cout << "Result of checkIncreasing(1,reverse=True,0): " << checkIncreasing(1, "reverse"_a=true,0).str() << std::endl;
     }
 
 
@@ -401,6 +414,7 @@ int main(int argc, char** argv)
         auto int_list = module.get("int_list");
         auto str_list = module.get("str_list");
         auto int_tuple = module.get("int_tuple");
+        auto mixed_tuple = module.get("mixed_tuple");
 
         // Modify first argument of int_list.
         // Similar to calling a Callable you can either hand
@@ -431,6 +445,15 @@ int main(int argc, char** argv)
         std::cout << "python integer tuple converted to vector<int>:" << std::endl;
         for(size_t i=0; i<int_vec2.size(); ++i)
             std::cout << int_vec2[i] << std::endl;
+
+        // Convert python mixed tuple to Dune::TupleVector and print results:
+        Dune::TupleVector<std::vector<int>,std::string> mixed_tuple_c;
+        mixed_tuple.toC(mixed_tuple_c);
+        std::cout << "python mixed tuple (string representation):" << mixed_tuple.str() << std::endl;
+        std::cout << "python mixed tuple converted to Dune::TupleVector:" << std::endl;
+        for(size_t i=0; i<mixed_tuple_c[_0].size(); ++i)
+            std::cout << mixed_tuple_c[_0][i] << std::endl;
+        std::cout << mixed_tuple_c[_1] << std::endl;
     }
 
 
@@ -472,19 +495,16 @@ int main(int argc, char** argv)
 
 
     std::cout << std::endl;
-    std::cout << "Example 12: Using python to define Dune::VirtualFunction's *********************" << std::endl;
+    std::cout << "Example 12: Using python to define std::functions's *********************" << std::endl;
     {
         // There are various ways to wrap callable objects from python into a
-        // PythonFunction that implements the Dune::VirtualFunction interface.
-        // The basic requirement is that there is a C->python conversion
+        // std::function.  The basic requirement is that there is a C->python conversion
         // for the domain type and a python->C conversion for the range type.
-        // PythonFunction also implements classic function call syntax via
-        // operator().
 
         int x=2;
         double y;
 
-        // We can construct a PythonFunction from a callable object in python.
+        // We can construct a std::function from a callable object in python.
 
         // This can be a function, ...
         pyMain.runStream()
@@ -492,39 +512,21 @@ int main(int argc, char** argv)
             << std::endl << "    return x + 3.141593 + 1"
             << std::endl;
 
-        PythonFunction<int, double> f(pyMain.get("f"));
-        f.evaluate(x, y);
-        std::cout << "Result of evaluating the wrapped function f through VirtualFunction interface: " << y << std::endl;
+        auto f = pyMain.get("f").toC<std::function<double(int)>>();
         std::cout << "Result of evaluating the wrapped function f: " << f(x) << std::endl;
 
         // ... a lambda function, ...
-        PythonFunction<int, double> lambda(Python::evaluate("lambda x: (x+2)"));
-        lambda.evaluate(x, y);
+        auto lambda = Python::evaluate("lambda x: (x+2)").toC<std::function<double(int)>>();
         std::cout << "Result of evaluating a wrapped lambda function: " << lambda(x) << std::endl;
 
         // ... a functor (instance with __call__ method), ...
         Python::run("c = C(1)");
-        PythonFunction<int, double> c(pyMain.get("c"));
+        auto c = pyMain.get("c").toC<std::function<double(int)>>();
         std::cout << "Result of evaluating the wrapped functor c: " << c(x) << std::endl;
 
         // ... or a method of an instance.
-        PythonFunction<int, double> c_f(pyMain.get("c").get("f"));
+        auto c_f = pyMain.get("c").get("f").toC<std::function<double(int)>>();
         std::cout << "Result of evaluating the wrapped method c.f: " << c_f(x) << std::endl;
-
-        // We can also construct a PythonFunction from a python expression
-        // using the argument named x. This will create a python
-        // function with the given name, taking the argument x,
-        // and returning the result of the expression evaluation.
-        PythonFunction<int, double> g("g", "x + 3.14159");
-        std::cout << "Result of evaluating the function g created from an expression: " << g(x) << std::endl;
-
-        // You can now also call the function in python directly...
-        std::cout << "Result of evaluating the created function g in python directly: ";
-        pyMain.run("print(g(2))");
-
-        // ... or by calling it manually.
-        std::cout << "Result of evaluating the created function g in python manually through the C++ interface: ";
-        std::cout << Python::Callable(pyMain.get("g"))(2).str() << std::endl;
     }
 
 
@@ -560,26 +562,26 @@ int main(int argc, char** argv)
     std::cout << std::endl;
     std::cout << "Example 14: Functions acting on dune types *************************************" << std::endl;
     {
-        // Conversion is also implementing for FieldVector so we can use this in PythonFunction
+        // Conversion is also implementing for FieldVector so we can use this in std::function
 
         // We can use vector valued functions
         {
             typedef Dune::FieldVector<double, 3> DT;
             typedef Dune::FieldVector<double, 2> RT;
+            using F = std::function<RT(DT)>;
             pyMain.runStream()
                 << std::endl << "def h(x):"
                 << std::endl << "    return (x[0] + x[1], x[1] + x[2])"
                 << std::endl;
-            PythonFunction<DT, RT> h(pyMain.get("h"));
-            std::cout << "Result of calling a PythonFunction<FV<double,3> , FV<double,2> > : " << h({1, 2, 3}) << std::endl;
+            std::function<RT(DT)> h = pyMain.get("h").toC<F>();
+            std::cout << "Result of calling a std::function< FV<double,2>(FV<double,3>) > : " << h({1, 2, 3}) << std::endl;
         }
 
         // Python scalars are treated like sequences with size one
         {
             typedef Dune::FieldVector<double, 2> DT;
             typedef Dune::FieldVector<double, 1> RT;
-            typedef Dune::VirtualFunction<DT, RT> FBase;
-            typedef PythonFunction<DT, RT> F;
+            typedef std::function<RT(DT)> F;
 
             DT x = {1, 2};
             RT y;
@@ -589,77 +591,8 @@ int main(int argc, char** argv)
                 << std::endl << "def h3(x):"
                 << std::endl << "    return (x[0] + x[1],)"  // Result is a list with size one
                 << std::endl;
-
-            // Call constructor directly
-            F h2(pyMain.get("h2"));
-
-            // Convert callable to shared_ptr of interface class
-            auto h3 = pyMain.get("h3").toC<std::shared_ptr<FBase>>();
-
-            std::cout << "Result of calling a PythonFunction<FV<double,3> , FV<double,1> > (python result is scalar): " << h2(x) << std::endl;
-            h3->evaluate(x,y);
-            std::cout << "Result of calling a PythonFunction<FV<double,3> , FV<double,1> > (python result is sigleton list): " << y << std::endl;
         }
 
-        // We can also construct scalar differentiable functions...
-        {
-            typedef Dune::FieldVector<double, 2> DT;
-            typedef Dune::FieldVector<double, 1> RT;
-            typedef DifferentiablePythonFunction<DT, RT> F;
-            typedef F::DerivativeType DRT;
-
-            DT x = {1, 2};
-            RT y;
-            DRT z;
-
-            pyMain.import("math");
-            pyMain.runStream()
-                << std::endl << "def f(x):"
-                << std::endl << "    return math.sin(x[0]) + math.cos(x[1])"     // Result is a scalar
-                << std::endl << "def df(x):"
-                << std::endl << "    return ((math.cos(x[0]),-math.sin(x[1])),)"  // Single row FieldMatrix can only be constructed from sequence of sequences
-                << std::endl;
-
-            // Call constructor directly
-            F f1(pyMain.get("f"), pyMain.get("df"));
-
-            f1.evaluate(x,y);
-            std::cout << "Result of calling a DifferentiablePythonFunction<FV<double,3> , FV<double,1> >::evaluate: " << y << std::endl;
-
-            f1.evaluateDerivative(x,z);
-            std::cout << "Result of calling a DifferentiablePythonFunction<FV<double,3> , FV<double,1> >::evaluateDerivative: " << z << std::endl;
-        }
-
-        // .. and vector valued differentiable functions...
-        {
-            typedef Dune::FieldVector<double, 3> DT;
-            typedef Dune::FieldVector<double, 2> RT;
-            typedef DifferentiablePythonFunction<DT, RT> F;
-            typedef VirtualDifferentiableFunction<DT, RT> FBase;
-            typedef F::DerivativeType DRT;
-
-            DT x = {1, 2, 3};
-            RT y;
-            DRT z;
-
-            pyMain.import("math");
-            pyMain.runStream()
-                << std::endl << "def f(x):"
-                << std::endl << "    return (x[0]**2+x[1], x[1]**2+x[2])"     // Result is a list
-                << std::endl << "def df(x):"
-                << std::endl << "    return ((2*x[0], 1, 0),(0, 2*x[1], 1))"  // Nested tuple as matrix
-                << std::endl << "fdf = (f, df)"
-                << std::endl;
-
-//            F f(pyMain.get("f"), pyMain.get("df"));
-            auto f = pyMain.get("fdf").toC<std::shared_ptr<FBase>>();
-
-            f->evaluate(x,y);
-            std::cout << "Result of calling a DifferentiablePythonFunction<FV<double,3> , FV<double,2> >::evaluate: " << y << std::endl;
-
-            f->evaluateDerivative(x,z);
-            std::cout << "Result of calling a DifferentiablePythonFunction<FV<double,3> , FV<double,2> >::evaluateDerivative: " << z << std::endl;
-        }
     }
 
 
@@ -671,8 +604,7 @@ int main(int argc, char** argv)
         //
         //   std::map,
         //   shared_ptr<T>,
-        //   PythonFunction*,
-        //   VirtualFunction*
+        //   std::function
         //
         // it is also possible to import dictionaries of callable
         // python objects into a map of functions.
@@ -705,9 +637,8 @@ int main(int argc, char** argv)
 
         typedef Dune::FieldVector<double, 3> DT;
         typedef Dune::FieldVector<double, 2> RT;
-        typedef Dune::VirtualFunction<DT, RT> Function;
-        typedef std::shared_ptr<Function> FunctionPointer;
-        typedef std::map<std::string, FunctionPointer> FunctionMap;
+        typedef std::function<RT(DT)> Function;
+        typedef std::map<std::string, Function> FunctionMap;
 
         // Import dictionaries of callables to map of functions.
         FunctionMap functions;
@@ -721,11 +652,8 @@ int main(int argc, char** argv)
 
         FunctionMap::iterator it = functions.begin();
         FunctionMap::iterator end = functions.end();
-        for(; it != end; ++it)
-        {
-            it->second->evaluate(x,y);
-            std::cout << it->first << "(x) = " << y << std::endl;
-        }
+        for(auto&& [key,f] : functions)
+            std::cout << key << "(x) = " << f(x) << std::endl;
     }
 
     std::cout << std::endl;
@@ -772,7 +700,7 @@ int main(int argc, char** argv)
     }
 
 
-#if HAVE_UG
+#if HAVE_DUNE_UGGRID
     std::cout << std::endl;
     std::cout << "Example 17: Filling a GridFactory from python **********************************" << std::endl;
     {

dune/fufem/test/dunepythontest.py

@@ -10,7 +10,13 @@ def add(x,y):
 int_list=[11,12,13]
 str_list=["string1","string2","string3"]
 int_tuple=(1,2,3)
+mixed_tuple=((1,2),"three")
 
+def checkIncreasing(a, b, reverse=False):
+    if (not reverse):
+        return a<=b
+    elif reverse:
+        return b<=a
 
 
 # some utilities for creating parametrized boundaries