Boundary dofs: extent to generic dune-functions basis

This change provides the usage of `constructBoundaryDofs`
with any dune-funcions compatible basis.

The user has to choose the correct data type for the BitVector.

For backward-compatibility an ISTLVectorBackend wrapper is included.

CHANGELOG.md

 # Master (will become release 2.8)
 
-- ...
+- constructBoundaryDofs:
+    - Small interface change in the template parameters: drop `blocksize` and replace it by `BitSetVector`
+    - The method can now handle generic `dune-functions` basis types, as long as we have consistency in the data types
 
 ## Deprecations and removals
 

dune/fufem/functiontools/boundarydofs.hh

@@ -10,59 +10,56 @@
 
 #include <dune/fufem/boundarypatch.hh>
 
+#include <dune/functions/backends/concepts.hh>
+#include <dune/functions/backends/istlvectorbackend.hh>
+#include <dune/functions/functionspacebases/subentitydofs.hh>
+
 
 /** \brief For a given basis and boundary patch, determine all degrees
     of freedom on the patch.
 
-    \param boundaryDofs Bit is set if corresponding dofs belongs to the boundary patch
+    A bit of boundaryDofs is set if corresponding dof belongs to the boundary patch.
+
+    \param boundaryDofs BitSetVector corresponding to the basis DOFs
     \param sfinae Dummy parameter to instantiate this method only for dune-functions bases
 */
-template <class GridView, class Basis, int blocksize>
+template <class GridView, class Basis, class BitSetVector>
 void constructBoundaryDofs(const BoundaryPatch<GridView>& boundaryPatch,
                            const Basis& basis,
-                           Dune::BitSetVector<blocksize>& boundaryDofs,
+                           BitSetVector& boundaryDofs,
                            typename Basis::LocalView* sfinae = nullptr)
 {
     // Check consistency of the input
     static_assert((std::is_same<GridView, typename Basis::GridView>::value),
                        "BoundaryPatch and global basis must be for the same grid view!");
 
+    // Small helper function to wrap vectors using istlVectorBackend
+    // if they do not already satisfy the VectorBackend interface.
+    auto toVectorBackend = [&](auto& v) -> decltype(auto) {
+        if constexpr (Dune::models<Dune::Functions::Concept::VectorBackend<Basis>, decltype(v)>()) {
+            return v;
+        } else {
+            return Dune::Functions::istlVectorBackend(v);
+        }
+    };
+
     //////////////////////////////////////////////////////////
     //   Init output bitfield
     //////////////////////////////////////////////////////////
 
-    boundaryDofs.resize(basis.size({}));
-    boundaryDofs.unsetAll();
+    // ensure IstlVectorBackend interface
+    auto boundaryDofsBackend = toVectorBackend(boundaryDofs);
 
-    auto localView = basis.localView();
-#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,6)
-    auto localIndexSet = basis.localIndexSet();
-#endif
+    boundaryDofsBackend.resize(basis);
+    boundaryDofsBackend = false;
 
-    for (auto it = boundaryPatch.begin(); it != boundaryPatch.end(); ++it)
+    auto localView = basis.localView();
+    auto seDOFs = Dune::Functions::subEntityDOFs(basis);
+    for(const auto& intersection : boundaryPatch)
     {
-        localView.bind(it->inside());
-#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,6)
-        localIndexSet.bind(localView);
-#endif
-        const auto& localCoefficients = localView.tree().finiteElement().localCoefficients();
-
-        for (size_t i=0; i<localCoefficients.size(); i++)
-        {
-            ////////////////////////////////////////////////////
-            //   Test whether dof is on the boundary face
-            ////////////////////////////////////////////////////
-
-            unsigned int entity = localCoefficients.localKey(i).subEntity();
-            unsigned int codim  = localCoefficients.localKey(i).codim();
-
-            if (it.containsInsideSubentity(entity, codim))
-#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,6)
-                boundaryDofs[localIndexSet.index(i)] = true;
-#else
-                boundaryDofs[localView.index(i)] = true;
-#endif
-        }
+        localView.bind(intersection.inside());
+        for(auto localIndex: seDOFs.bind(localView, intersection))
+            boundaryDofsBackend[localView.index(localIndex)] = true;
     }
 }
 

dune/fufem/test/CMakeLists.txt

@@ -8,6 +8,7 @@ set(GRID_BASED_TESTS
   boundarypatchprolongatortest
   coarsegridfunctionwrappertest
   composedfunctiontest
+  constructboundarydofstest
   dunefunctionsassemblertest
   dunefunctionsipdgassemblertest
   functionintegratortest

dune/fufem/test/constructboundarydofstest.cc

+#include <config.h>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/test/testsuite.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/functiontools/boundarydofs.hh>
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/test/common.hh>
+
+#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <dune/istl/bvector.hh>
+
+using namespace Dune;
+
+template<int order, class GW>
+bool checkConstructBoundaryDofs(const GW& gw)
+{
+    // create active vertices
+    const auto dim = GW::dimension;
+    const auto size = gw.size(dim);
+    BitSetVector<dim> vertices(size);
+    // set some values
+    vertices[0].set();
+    vertices[size / 2].set();
+    vertices[size-1].set();
+
+    // create BoundaryPatch
+    BoundaryPatch bp(gw,vertices);
+
+    // create scalar and powerBasis
+    using namespace Functions::BasisFactory;
+    auto lagrangeBasis = makeBasis(gw, lagrange<order>());
+    auto powerBasis = makeBasis(gw, power<3>(lagrange<order>()));
+
+    // create common bitsetvector
+    BitSetVector<1> lagrangeDOFs;
+    BitSetVector<3> powerDOFs;
+
+    // ... and functions-compatible bitVectors
+    Dune::BlockVector<Dune::FieldVector<char,1>> lagrangeBits;
+    Dune::BlockVector<Dune::FieldVector<char,3>> powerBits;
+    auto lagrangeBackend = Dune::Functions::istlVectorBackend(lagrangeBits);
+    auto powerBackend = Dune::Functions::istlVectorBackend(powerBits);
+
+
+    // check whether this compiles
+    constructBoundaryDofs(bp,lagrangeBasis,lagrangeDOFs);
+    constructBoundaryDofs(bp,powerBasis,powerDOFs);
+    constructBoundaryDofs(bp,lagrangeBasis,lagrangeBackend);
+    constructBoundaryDofs(bp,powerBasis,powerBackend);
+
+    return true;
+}
+
+
+struct Suite
+{
+    template<class GridType>
+    bool check(const GridType& grid)
+    {
+        auto maxLevel = grid.maxLevel();
+        auto gridView = grid.levelGridView(maxLevel);
+
+        bool passed = checkConstructBoundaryDofs<0>(gridView);
+        passed = passed and checkConstructBoundaryDofs<1>(gridView);
+        passed = passed and checkConstructBoundaryDofs<2>(gridView);
+        return passed;
+    }
+};
+
+int main(int argc, char** argv)
+{
+    Dune::MPIHelper::instance(argc, argv);
+
+    Suite tests;
+
+    bool passed = checkWithStructuredGrid<2>(tests, 3);
+    passed = passed and checkWithStructuredGrid<3>(tests, 3);
+
+    return passed ? 0 : 1;
+}