From c5205498cb7f7fdbf1231be17c516a582ae269b9 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Wed, 14 Nov 2018 11:50:30 +0100
Subject: [PATCH] Also implement the boundarydofs method for dune-functions
 bases

With this patch, the method works both for old-style dune-fufem bases
and for dune-functions bases directly.  The distinction is done
using some minor SFINAE trickery.
---
 dune/fufem/functiontools/boundarydofs.hh | 49 +++++++++++++++++++++++-
 1 file changed, 48 insertions(+), 1 deletion(-)

diff --git a/dune/fufem/functiontools/boundarydofs.hh b/dune/fufem/functiontools/boundarydofs.hh
index 57dd6b94..b156363c 100644
--- a/dune/fufem/functiontools/boundarydofs.hh
+++ b/dune/fufem/functiontools/boundarydofs.hh
@@ -15,11 +15,58 @@
     of freedom on the patch.
 
     \param boundaryDofs Bit is set if corresponding dofs belongs to the boundary patch
+    \param sfinae Dummy parameter to instantiate this method only for dune-functions bases
 */
 template <class GridView, class Basis, int blocksize>
 void constructBoundaryDofs(const BoundaryPatch<GridView>& boundaryPatch,
                            const Basis& basis,
-                           Dune::BitSetVector<blocksize>& boundaryDofs)
+                           Dune::BitSetVector<blocksize>& 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!");
+
+    //////////////////////////////////////////////////////////
+    //   Init output bitfield
+    //////////////////////////////////////////////////////////
+
+    boundaryDofs.resize(basis.size({}));
+    boundaryDofs.unsetAll();
+
+    auto localView = basis.localView();
+
+    for (auto it = boundaryPatch.begin(); it != boundaryPatch.end(); ++it)
+    {
+        localView.bind(it->inside());
+        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))
+                boundaryDofs[localView.index(i)] = true;
+        }
+    }
+}
+
+
+/** \brief For a given basis and boundary patch, determine all degrees
+    of freedom on the patch.
+
+    Same method as above, but this time for dune-fufem bases.
+*/
+template <class GridView, class Basis, int blocksize>
+void constructBoundaryDofs(const BoundaryPatch<GridView>& boundaryPatch,
+                           const Basis& basis,
+                           Dune::BitSetVector<blocksize>& boundaryDofs,
+                           typename Basis::LocalFiniteElement* sfinae=nullptr)
 {
     // Check consistency of the input
     static_assert((std::is_same<GridView, typename Basis::GridView>::value),
-- 
GitLab