dune/fufem/assemblers/localassemblers/adolchessianassembler.hh

@@ -33,8 +33,8 @@ class AdolcHessianAssembler :
 
     typedef VirtualGridFunction<GridType, Dune::FieldVector<field_type,blocksize> > GridFunction;
 
-    typedef typename Dune::LocalFiniteElementFunctionBase<TrialLocalFE>::type FunctionBaseClass;
-    typedef LocalFunctionComponentWrapper<GridFunction, GridType, FunctionBaseClass> LocalWrapper;
+    typedef typename TrialLocalFE::Traits::LocalBasisType::Traits FunctionTraits;
+    typedef LocalFunctionComponentWrapper<GridFunction, GridType, FunctionTraits> LocalWrapper;
 
   public:
     typedef typename Base::Element Element;

dune/fufem/assemblers/localassemblers/adolclinearizationassembler.hh

@@ -29,8 +29,8 @@ class AdolcLinearizationAssembler :
 
     typedef VirtualGridFunction<GridType, T> GridFunction;
 
-    typedef typename Dune::LocalFiniteElementFunctionBase<TrialLocalFE>::type FunctionBaseClass;
-    typedef LocalFunctionComponentWrapper<GridFunction, GridType, FunctionBaseClass> LocalWrapper;
+    typedef typename TrialLocalFE::Traits::LocalBasisType::Traits FunctionTraits;
+    typedef LocalFunctionComponentWrapper<GridFunction, GridType, FunctionTraits> LocalWrapper;
 
   public:
     typedef typename Base::Element Element;

dune/fufem/assemblers/localassemblers/convolutionassembler.hh

@@ -3,13 +3,15 @@
 
 #include <utility>
 
-#include <dune/common/function.hh>
-
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
 
 #include <dune/istl/matrix.hh>
 
+// Borrow internal helper from dune-localfunctions
+// for transition of function interface
+#include <dune/localfunctions/common/localinterpolation.hh>
+
 #include <dune/matrix-vector/addtodiagonal.hh>
 
 #include "dune/fufem/staticmatrixtools.hh"
@@ -37,10 +39,11 @@ class ConvolutionAssembler
         typedef typename Dune::template FieldVector<double,AnsatzGridView::dimensionworld> AnsatzGlobalCoordinate;
         typedef typename std::template pair<TrialGlobalCoordinate, AnsatzGlobalCoordinate> KernelArgument;
 
-        typedef typename Dune::VirtualFunction<KernelArgument, double> KernelFunction;
+        typedef typename std::function<double(KernelArgument)> KernelFunction;
 
-        ConvolutionAssembler(const KernelFunction& kernel, int tQuadOrder=2, int aQuadOrder=2):
-            kernel_(kernel),
+        template<class KFunc>
+        ConvolutionAssembler(const KFunc& kernel, int tQuadOrder=2, int aQuadOrder=2):
+            kernel_(Dune::Impl::makeFunctionWithCallOperator<KernelArgument>(kernel)),
             tQuadOrder_(tQuadOrder),
             aQuadOrder_(aQuadOrder)
         {}
@@ -112,7 +115,7 @@ class ConvolutionAssembler
 
 
                     // evauluate integration kernel at quadrature points
-                    kernel_.evaluate(std::make_pair(tGlobalQuadPos, aGlobalQuadPos), kernelValue);
+                    kernelValue = kernel_(std::make_pair(tGlobalQuadPos, aGlobalQuadPos));
 
                     for (size_t i=0; i<tFE.localBasis().size(); ++i)
                     {
@@ -131,7 +134,7 @@ class ConvolutionAssembler
 
 
     protected:
-        const KernelFunction& kernel_;
+        KernelFunction kernel_;
         const int tQuadOrder_;
         const int aQuadOrder_;
 };

dune/fufem/assemblers/localassemblers/dunefunctionsl2functionalassembler.hh

@@ -34,7 +34,7 @@ class DuneFunctionsL2FunctionalAssembler :
         using GlobalCoordinate = typename Grid::template Codim<0>::Geometry::GlobalCoordinate;
         using LocalCoordinate = typename Grid::template Codim<0>::Geometry::LocalCoordinate;
         using Function = F;
-        using LocalFunction = std::decay_t<decltype(localFunction(std::declval<Function>()))>;
+        using LocalFunction = std::decay_t<decltype(localFunction(std::declval<Function&>()))>;
 
         static const int dim = Grid::dimension;
         static const int dimworld = Grid::dimensionworld;
@@ -82,9 +82,9 @@ class DuneFunctionsL2FunctionalAssembler :
          * \param f The L2 function
          */
         template<class FF, disableCopyMove<DuneFunctionsL2FunctionalAssembler, FF> = 0 >
-        DuneFunctionsL2FunctionalAssembler(F&& f) :
+        DuneFunctionsL2FunctionalAssembler(FF&& f) :
             f_(Dune::wrap_or_move<const Function>(std::forward<FF>(f))),
-            localF_(localFunction(f_)),
+            localF_(localFunction(*f_)),
             functionQuadKey_(dim, 0)
         {}
 

dune/fufem/assemblers/localassemblers/dunefunctionsmassassembler.hh

+#pragma once
+
+#warning This file is deprecated!  Please use Dune::Fufem::MassAssembler from the file massassembler.hh!
+
+#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
+
+namespace Dune::Fufem
+{
+  using DuneFunctionsLocalMassAssembler = Dune::Fufem::MassAssembler;
+} // namespace Dune::Fufem

dune/fufem/assemblers/localassemblers/gradientassembler.hh

@@ -8,8 +8,6 @@
 #include <dune/geometry/quadraturerules.hh>
 #include <dune/geometry/referenceelements.hh>
 
-#include <dune/localfunctions/common/virtualinterface.hh>
-
 #include <dune/istl/matrix.hh>
 
 #include <dune/fufem/functions/localbasisderivativefunction.hh>
@@ -57,9 +55,9 @@ class GradientAssembler
             // get transposed inverse of Jacobian of transformation
             const FMdimdim& invJacobian = geometry.jacobianInverseTransposed(pos);
 
-            typedef typename Dune::LocalFiniteElementFunctionBase<TrialLocalFE>::type FunctionBaseClass;
+            typedef typename TrialLocalFE::Traits::LocalBasisType::Traits FunctionTraits;
 
-            LocalBasisDerivativeFunction<AnsatzLocalFE, FunctionBaseClass> derivative(aFE.localBasis());
+            LocalBasisDerivativeFunction<AnsatzLocalFE, FunctionTraits> derivative(aFE.localBasis());
 
             std::vector<JacobianType> referenceGradients(tFE.localBasis().size());
 

dune/fufem/assemblers/localassemblers/interiorpenaltydgassembler.hh

@@ -45,8 +45,8 @@ class InteriorPenaltyDGAssembler : public LocalAssembler<GridType, TrialLocalFE,
     InteriorPenaltyDGAssembler() :
       sigma0_(0),
       dirichlet_(true),
-      sigma1_(0),
-      dgType_((double) DGType::SIPG) {}
+      dgType_((double) DGType::SIPG),
+      sigma1_(0) {}
 
     InteriorPenaltyDGAssembler(double penalty, bool dirichlet=true, DGType dgType=DGType::SIPG, double gradientPenalty=0) :
       sigma0_(penalty),
@@ -216,18 +216,10 @@ class InteriorPenaltyDGAssembler : public LocalAssembler<GridType, TrialLocalFE,
           // Depending on which elements i and j belong to, perform different operations.
           // In particular, i (resp. j) < insizeSize corresponds to functions from the inner element
           // while the others correspond to the outer element
-          if (i < insideSize){
-            if(j< insideSize)
-              Dune::MatrixVector::addToDiagonal(localMatrix[i][j], mc[0][0][i][j]);
+          if (i < insideSize)
+            localMatrix[i][j] += (j< insideSize) ?  mc[0][0][i][j] : mc[0][1][i][j-insideSize];
           else
-              Dune::MatrixVector::addToDiagonal(localMatrix[i][j], mc[0][1][i][j-insideSize]);
-          }
-          else {
-            if(j< insideSize)
-              Dune::MatrixVector::addToDiagonal(localMatrix[i][j], mc[1][0][i-insideSize][j]);
-            else
-              Dune::MatrixVector::addToDiagonal(localMatrix[i][j], mc[1][1][i-insideSize][j-insideSize]);
-          }
+            localMatrix[i][j] += (j< insideSize) ? mc[1][0][i-insideSize][j] : mc[1][1][i-insideSize][j-insideSize];
         }
       }
     }
@@ -362,16 +354,16 @@ class InteriorPenaltyDGAssembler : public LocalAssembler<GridType, TrialLocalFE,
      */
     double sigma0_; // Penalizes discontinuities
 
-    /** Penalty term \sigma_1
-     *
-     * This term penalizes jumps in the derivatives across the edges.
-     */
 
     /** Flag that checks if dirichlet_ terms must be assembled for the edges or not.*/
     bool dirichlet_;
 
     const double dgType_; // -1 leads to the symmetric interior penalty DG variant. {-1,0,1} are possible values.
 
+    /** Penalty term \sigma_1
+     *
+     * This term penalizes jumps in the derivatives across the edges.
+     */
     double sigma1_; // Penalizes jumps of the normal derivative along the non-boundary edges
 };
 #endif

dune/fufem/assemblers/localassemblers/laplaceassembler.hh

-#ifndef LAPLACE_ASSEMBLER_HH
-#define LAPLACE_ASSEMBLER_HH
+#ifndef DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_LAPLACEASSEMBLER_HH
+#define DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_LAPLACEASSEMBLER_HH
 
 
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
+#include <dune/common/version.hh>
+#if !DUNE_VERSION_GTE(DUNE_GEOMETRY, 2, 9)
+#include <dune/common/transpose.hh>
+#endif
 
 #include <dune/istl/matrix.hh>
 
 #include <dune/matrix-vector/addtodiagonal.hh>
 
+#include <dune/typetree/visitor.hh>
+#include <dune/typetree/traversal.hh>
+
 #include "dune/fufem/quadraturerules/quadraturerulecache.hh"
 
 #include "dune/fufem/assemblers/localoperatorassembler.hh"
 
-/** \brief Local assembler for the Laplace problem */
+
+namespace Dune::Fufem
+{
+
+  namespace Impl
+  {
+
+    /** \brief Computes the inner product of the jacobians of the local basis functions
+     *
+     * It computes at each leaf node the inner product (hence, stiffness matrix) of the
+     * jacobians of the local basis functions. It determines the local index and updates the
+     * corresponding entry of the local matrix.
+     */
+    template <class LV, class M>
+    class LeafNodeLaplaceAssembler
+    {
+    public:
+
+      using LocalView = LV;
+      using Matrix = M;
+      using field_type = typename M::field_type;
+
+      /** \brief Constructor with all necessary context information
+       *
+       * \param [in] lv (localView) element information and localBasis
+       * \param [out] m (matrix) resulting local stiffness matrix, wrapped in the ISTLBackend
+
+       */
+      LeafNodeLaplaceAssembler(const LV& lv, M& m)
+      : localView_(lv)
+      , matrix_(m)
+      {}
+
+      template<class Node, class TreePath>
+      void operator()(Node& node, [[maybe_unused]] TreePath treePath)
+      {
+        const auto& element = localView_.element();
+        const auto& geometry = element.geometry();
+        const auto& finiteElement = node.finiteElement();
+        const auto& localBasis = finiteElement.localBasis();
+        using JacobianType = typename std::decay_t<decltype(localBasis)>::Traits::JacobianType;
+        constexpr int gridDim = LocalView::GridView::dimension;
+
+        auto feQuadKey = QuadratureRuleKey(finiteElement);
+        // we have a product of jacobians
+        auto quadKey = feQuadKey.derivative().square();
+        const auto& quad = QuadratureRuleCache<double, gridDim>::rule(quadKey);
+
+        for( const auto& qp : quad )
+        {
+          const auto& quadPos = qp.position();
+          const auto integrationElement = geometry.integrationElement(quadPos);
+
+          std::vector<JacobianType> referenceJacobians;
+          localBasis.evaluateJacobian(quadPos, referenceJacobians);
+
+#if DUNE_VERSION_GTE(DUNE_GEOMETRY, 2, 9)
+          const auto& jacobianInverse = geometry.jacobianInverse(quadPos);
+#else
+          const auto& jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
+#endif
+
+          // transform jacobians
+#if DUNE_VERSION_GTE(DUNE_GEOMETRY, 2, 9)
+          std::vector<decltype(referenceJacobians[0] * jacobianInverse)> jacobians(referenceJacobians.size());
+#else
+          using Jacobian = Dune::FieldMatrix<double,JacobianType::rows,jacobianInverseTransposed.rows>;
+          std::vector<Jacobian> jacobians(referenceJacobians.size());
+#endif
+          for (size_t i=0; i<jacobians.size(); i++)
+#if DUNE_VERSION_GTE(DUNE_GEOMETRY, 2, 9)
+            jacobians[i] = referenceJacobians[i] * jacobianInverse;
+#else
+            for (size_t j=0; j<JacobianType::rows; j++)
+              jacobianInverseTransposed.mv(referenceJacobians[i][j], jacobians[i][j]);
+#endif
+
+          // compute matrix entries = inner products of jacobians
+          auto z = qp.weight() * integrationElement;
+          for (std::size_t i=0; i<localBasis.size(); ++i)
+          {
+            auto zi = z*jacobians[i];
+
+            auto rowIndex = node.localIndex(i);
+
+            // start off-diagonal, treat the diagonal extra
+            for (std::size_t j=i+1; j<localBasis.size(); ++j)
+            {
+              auto colIndex = node.localIndex(j);
+
+              // we need the Frobenius inner product here, since the matrix may be fully occupied (Raviart-Thomas, etc)
+              // I could not find an implementation anywhere in DUNE, therefore this ugly loop:
+              for ( std::size_t ii=0; ii<zi.N(); ii++ )
+              {
+                for ( std::size_t jj=0; jj<zi.M(); jj++ )
+                {
+                  // stiffness matrix is symmetric
+                  auto zij = zi[ii][jj] * jacobians[j][ii][jj];
+                  matrix_[rowIndex][colIndex] += zij;
+                  matrix_[colIndex][rowIndex] += zij;
+                }
+              }
+            }
+
+            // z * values[i]^2
+            for ( std::size_t ii=0; ii<zi.N(); ii++ )
+              for ( std::size_t jj=0; jj<zi.M(); jj++ )
+                matrix_[rowIndex][rowIndex] += zi[ii][jj] * jacobians[i][ii][jj];
+          }
+        }
+      }
+
+    private:
+      const LocalView& localView_;
+      Matrix& matrix_;
+    };
+
+  } // namespace Impl
+
+  /** \brief Local Laplace (stiffness) assembler for dune-functions basis
+   *
+   * We assume the stiffness matrix to be symmetric and hence ansatz and trial space to be equal!
+   * The implementation allows an arbitrary dune-functions compatible basis, as long as there
+   * is an ISTLBackend available for the resulting matrix.
+   */
+  class LaplaceAssembler
+  {
+  public:
+
+    template<class Element, class LocalMatrix, class LocalView>
+    void operator()(const Element& element, LocalMatrix& localMatrix, const LocalView& trialLocalView, const LocalView& ansatzLocalView) const
+    {
+      // the Laplace matrix operator assumes trial == ansatz
+      if (&trialLocalView != &ansatzLocalView)
+        DUNE_THROW(Dune::Exception,"The Laplace matrix operator assumes equal ansatz and trial functions");
+
+      // matrix was already resized before
+      localMatrix = 0.;
+
+      // create a tree visitor and compute the inner products of the local functions at the leaf nodes
+      Impl::LeafNodeLaplaceAssembler leadNodeLaplaceAssembler(trialLocalView,localMatrix);
+      TypeTree::forEachLeafNode(trialLocalView.tree(),leadNodeLaplaceAssembler);
+    }
+
+  };
+
+} // namespace Dune::Fufem
+
+
+/** \brief Local assembler for the Laplace problem
+
+  * \deprecated This assembler uses the old dune-fufem function space bases.  Please migrate
+                to the assembler based on dune-functions bases (above).
+*/
 template <class GridType, class TrialLocalFE, class AnsatzLocalFE, class T=Dune::FieldMatrix<double,1,1> >
 class LaplaceAssembler : public LocalOperatorAssembler <GridType, TrialLocalFE, AnsatzLocalFE, T >
 {
@@ -54,7 +214,7 @@ class LaplaceAssembler : public LocalOperatorAssembler <GridType, TrialLocalFE,
             quadOrder_(-1)
         {}
 
-        DUNE_DEPRECATED_MSG("Quadrature order is now selected automatically. you don't need to specify it anymore.")
+        [[deprecated("Quadrature order is now selected automatically. you don't need to specify it anymore.")]]
         LaplaceAssembler(int quadOrder) :
             quadOrder_(quadOrder)
         {}
@@ -237,5 +397,5 @@ class LaplaceAssembler : public LocalOperatorAssembler <GridType, TrialLocalFE,
 };
 
 
-#endif
+#endif  // DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_LAPLACEASSEMBLER_HH
 

dune/fufem/assemblers/localassemblers/lumpedmassassembler.hh

@@ -33,7 +33,7 @@ class LumpedMassAssembler : public LocalOperatorAssembler < GridType, TrialLocal
             quadOrder_(-1)
         {}
 
-        DUNE_DEPRECATED_MSG("Quadrature order is now selected automatically. you don't need to specify it anymore.")
+        [[deprecated("Quadrature order is now selected automatically. you don't need to specify it anymore.")]]
         LumpedMassAssembler(int quadOrder):
             quadOrder_(quadOrder)
         {}

dune/fufem/assemblers/localassemblers/massassembler.hh

-#ifndef MASS_ASSEMBLER_HH
-#define MASS_ASSEMBLER_HH
+#ifndef DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_MASSASSEMBLER_HH
+#define DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_MASSASSEMBLER_HH
 
+#include <vector>
 
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
 
 #include <dune/istl/matrix.hh>
 
-#include <dune/matrix-vector/addtodiagonal.hh>
+#include <dune/typetree/visitor.hh>
+#include <dune/typetree/traversal.hh>
 
-#include "dune/fufem/quadraturerules/quadraturerulecache.hh"
+#include <dune/matrix-vector/addtodiagonal.hh>
 
 #include "dune/fufem/assemblers/localoperatorassembler.hh"
+#include "dune/fufem/quadraturerules/quadraturerulecache.hh"
 
 
-//** \brief Local mass assembler **//
-template <class GridType, class TrialLocalFE, class AnsatzLocalFE, class T=Dune::FieldMatrix<double,1,1> >
-class MassAssembler : public LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE, T >
+namespace Dune::Fufem
 {
-    private:
-        typedef LocalOperatorAssembler < GridType, TrialLocalFE, AnsatzLocalFE ,T > Base;
-        static const int dim = GridType::dimension;
-        const int quadOrder_;
-
-
-    public:
-        typedef typename Base::Element Element;
-        typedef typename Element::Geometry Geometry;
-        typedef typename Base::BoolMatrix BoolMatrix;
-        typedef typename Base::LocalMatrix LocalMatrix;
-
-        MassAssembler():
-            quadOrder_(-1)
-        {}
-
-        DUNE_DEPRECATED_MSG("Quadrature order is now selected automatically. you don't need to specify it anymore.")
-        MassAssembler(int quadOrder):
-            quadOrder_(quadOrder)
-        {}
 
-        void indices(const Element& element, BoolMatrix& isNonZero, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
+  namespace Impl
   {
-            isNonZero = true;
-        }
 
-        template <class BoundaryIterator>
-        void indices(const BoundaryIterator& it, BoolMatrix& isNonZero, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
+    /** \brief Computes the inner product of the local basis functions
+     *
+     * It computes at each leaf node the inner product (hence, mass matrix) of the
+     * local basis function. It determines the local index and updates the
+     * corresponding entry of the local matrix.
+     */
+    template <class LV, class M>
+    class LeafNodeMassAssembler
     {
-            isNonZero = true;
-        }
+    public:
 
-        void assemble(const Element& element, LocalMatrix& localMatrix, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
-        {
-            typedef typename Dune::template FieldVector<double,dim> FVdim;
-            typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::RangeType RangeType;
+      using LocalView = LV;
+      using Matrix = M;
+      using field_type = typename M::field_type;
 
-            // Make sure we got suitable shape functions
-            assert(tFE.type() == element.type());
-            assert(aFE.type() == element.type());
+      /** \brief Constructor with all necessary context information
+       *
+       * \param [in] lv (localView) element information and localBasis
+       * \param [out] m (matrix) resulting local mass matrix, wrapped in the ISTLBackend
 
-            // check if ansatz local fe = test local fe
-//            if (not Base::isSameFE(tFE, aFE))
-//                DUNE_THROW(Dune::NotImplemented, "MassAssembler is only implemented for ansatz space=test space!");
+       */
+      LeafNodeMassAssembler(const LV& lv, M& m)
+      : localView_(lv)
+      , matrix_(m)
+      {}
 
-            int rows = localMatrix.N();
-            int cols = localMatrix.M();
+      template<class Node, class TreePath>
+      void operator()(Node& node, TreePath treePath)
+      {
+        const auto& element = localView_.element();
+        const auto& geometry = element.geometry();
+        const auto& finiteElement = node.finiteElement();
+        const auto& localBasis = finiteElement.localBasis();
+        using RangeType = typename std::decay_t<decltype(localBasis)>::Traits::RangeType;
+        constexpr int gridDim = LocalView::GridView::dimension;
 
-            // get geometry and store it
-            const Geometry geometry = element.geometry();
+        std::vector<RangeType> values(localBasis.size());
 
-            localMatrix = 0.0;
+        auto feQuadKey = QuadratureRuleKey(finiteElement);
+        // we have a product of FE functions
+        auto quadKey = feQuadKey.square();
+        const auto& quad = QuadratureRuleCache<double, gridDim>::rule(quadKey);
 
-            // get quadrature rule
-            QuadratureRuleKey quadKey = QuadratureRuleKey(tFE)
-                .product(QuadratureRuleKey(aFE));
-            if (quadOrder_>=0)
-                quadKey.setOrder(quadOrder_);
-            const Dune::template QuadratureRule<double, dim>& quad = QuadratureRuleCache<double, dim>::rule(quadKey);
+        for( const auto& qp : quad )
+        {
+          const auto& quadPos = qp.position();
+          const auto integrationElement = geometry.integrationElement(quadPos);
 
-            // store values of shape functions
-            std::vector<RangeType> values(tFE.localBasis().size());
+          localBasis.evaluateFunction(quadPos, values);
 
-            // loop over quadrature points
-            for (size_t pt=0; pt < quad.size(); ++pt)
+          // compute matrix entries = inner products
+          auto z = qp.weight() * integrationElement;
+          for (std::size_t i=0; i<localBasis.size(); ++i)
           {
-                // get quadrature point
-                const FVdim& quadPos = quad[pt].position();
-
-                // get integration factor
-                const double integrationElement = geometry.integrationElement(quadPos);
+            auto zi = values[i]*z;
 
-                // evaluate basis functions
-                tFE.localBasis().evaluateFunction(quadPos, values);
+            auto rowIndex = node.localIndex(i);
 
-                // compute matrix entries
-                double z = quad[pt].weight() * integrationElement;
-                for(int i=0; i<rows; ++i)
+            // start off-diagonal, treat the diagonal extra
+            for (std::size_t j=i+1; j<localBasis.size(); ++j)
             {
-                    double zi = values[i]*z;
+              auto zij = values[j] * zi;
 
-                    for (int j=i+1; j<cols; ++j)
-                    {
-                        double zij = values[j] * zi;
-                        Dune::MatrixVector::addToDiagonal(localMatrix[i][j],zij);
-                        Dune::MatrixVector::addToDiagonal(localMatrix[j][i],zij);
+              auto colIndex = node.localIndex(j);
+
+              // mass matrix is symmetric
+              matrix_[rowIndex][colIndex] += zij;
+              matrix_[colIndex][rowIndex] += zij;
             }
-                    Dune::MatrixVector::addToDiagonal(localMatrix[i][i], values[i] * zi);
+
+            // z * values[i]^2
+            matrix_[rowIndex][rowIndex] += values[i] * zi;
           }
         }
       }
 
+    private:
+      const LocalView& localView_;
+      Matrix& matrix_;
+    };
+
+  } // namespace Impl
 
-        /** \brief Assemble the local mass matrix for a given boundary face
+  /** \brief Local mass assembler for dune-functions basis
+   *
+   * We assume the mass matrix to be symmetric and hence ansatz and trial space to be equal!
+   * The implementation allows an arbitrary dune-functions compatible basis, as long as there
+   * is an ISTLBackend available for the resulting matrix.
    */
-        template <class BoundaryIterator>
-        void assemble(const BoundaryIterator& it, LocalMatrix& localMatrix, const TrialLocalFE& tFE, const AnsatzLocalFE& aFE) const
+  class MassAssembler
   {
-            typedef typename BoundaryIterator::Intersection::Geometry Geometry;
-            typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::RangeType RangeType;
-
-            // Make sure we got suitable shape functions
-            assert(tFE.type() == it->inside().type());
-            assert(aFE.type() == it->inside().type());
-
-            // check if ansatz local fe = test local fe
-            if (not Base::isSameFE(tFE, aFE))
-                DUNE_THROW(Dune::NotImplemented, "MassAssembler is only implemented for ansatz space=test space!");
-
-            int rows = localMatrix.N();
-            int cols = localMatrix.M();
-
-            // get geometry and store it
-            const Geometry intersectionGeometry = it->geometry();
-
-            localMatrix = 0.0;
-
-            // get quadrature rule
-            QuadratureRuleKey tFEquad(it->type(), tFE.localBasis().order());
-            QuadratureRuleKey quadKey = tFEquad.square();
-
-            const Dune::template QuadratureRule<double, dim-1>& quad = QuadratureRuleCache<double, dim-1>::rule(quadKey);
-
-            // store values of shape functions
-            std::vector<RangeType> values(tFE.localBasis().size());
+  public:
 
-            // loop over quadrature points
-            for (size_t pt=0; pt < quad.size(); ++pt)
+    template<class Element, class LocalMatrix, class LocalView>
+    void operator()(const Element& element, LocalMatrix& localMatrix, const LocalView& trialLocalView, const LocalView& ansatzLocalView) const
     {
-                // get quadrature point
-                const Dune::FieldVector<double,dim-1>& quadPos = quad[pt].position();
-
-                // get integration factor
-                const double integrationElement = intersectionGeometry.integrationElement(quadPos);
-
-                // get values of shape functions
-                tFE.localBasis().evaluateFunction(it->geometryInInside().global(quadPos), values);
+      // the mass matrix operator assumes trial == ansatz
+      if (&trialLocalView != &ansatzLocalView)
+        DUNE_THROW(Dune::Exception,"The mass matrix operator assumes equal ansatz and trial functions");
 
-                // compute matrix entries
-                double z = quad[pt].weight() * integrationElement;
-                for (int i=0; i<rows; ++i)
-                {
-                    for (int j=i+1; j<cols; ++j)
-                    {
-                        double zij = values[i] * values[j] * z;
-                        Dune::MatrixVector::addToDiagonal(localMatrix[i][j],zij);
-                        Dune::MatrixVector::addToDiagonal(localMatrix[j][i],zij);
-                    }
-                    Dune::MatrixVector::addToDiagonal(localMatrix[i][i], values[i] * values[i] * z);
-                }
-            }
+      // matrix was already resized before
+      localMatrix = 0.;
 
+      // create a tree visitor and compute the inner products of the local functions at the leaf nodes
+      Impl::LeafNodeMassAssembler leadNodeMassAssembler(trialLocalView,localMatrix);
+      TypeTree::forEachLeafNode(trialLocalView.tree(),leadNodeMassAssembler);
     }
 
   };
 
+} // namespace Dune::Fufem
 
-#endif
+#endif  // DUNE_FUFEM_ASSEMBLERS_LOCALASSEMBLERS_MASSASSEMBLER_HH
 

dune/fufem/assemblers/localassemblers/secondorderoperatorassembler.hh

@@ -9,7 +9,7 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/common/hybridutilities.hh>
-#include <dune/common/std/apply.hh>
+#include <dune/common/concept.hh>
 
 #include <dune/istl/matrix.hh>
 
@@ -24,8 +24,6 @@
 #include "dune/fufem/assemblers/localoperatorassembler.hh"
 
 
-#include "dune/fufem/concept.hh"
-
 namespace Dune {
 namespace Fufem {
 namespace Concept {
@@ -39,13 +37,13 @@ namespace Concept {
 
         template<class F>
         auto require(F&& f) -> decltype(
-            Dune::Std::apply(f, args)
+            std::apply(f, args)
         );
     };
 
     template<class ResultMat, class Scalar, class Coord, class... ReturnTypes>
     struct SecondOrderOperatorAssemblerContraction
-        : Refines<SecondOrderOperatorContraction<Coord, ReturnTypes...>>
+        : Dune::Concept::Refines<SecondOrderOperatorContraction<Coord, ReturnTypes...>>
     {
         using BaseConcept = SecondOrderOperatorContraction<Coord, ReturnTypes...>;
 
@@ -55,7 +53,7 @@ namespace Concept {
 
         template<class F>
         auto require(F&& f) -> decltype(
-            Dune::MatrixVector::addProduct(r, s, Dune::Std::apply(f, args))
+            Dune::MatrixVector::addProduct(r, s, std::apply(f, args))
         );
     };
 
@@ -164,7 +162,7 @@ public:
                 MatrixEntry, double, WorldCoordinate,
                 typename std::result_of<FunctionTypes(WorldCoordinate)>::type...>;
         static_assert(
-            Dune::Fufem::Concept::models<ContractionConcept, Contraction>(),
+            Dune::models<ContractionConcept, Contraction>(),
             "Contraction type does not support: Dune::MatrixVector::addProduct(MatrixEntry, double, contraction(WorldCoordinate, WorldCoordinate, ...))");
     }
     SecondOrderOperatorAssembler(const Contraction& contraction,
@@ -256,7 +254,7 @@ public:
                     };
                     Dune::MatrixVector::addProduct(
                         localMatrix[i][j], z,
-                        Dune::Std::apply(boundContraction, returnValues));
+                        std::apply(boundContraction, returnValues));
                 }
             }
         }

dune/fufem/assemblers/localassemblers/subgridl2functionalassembler.hh

@@ -58,7 +58,7 @@ class SubgridL2FunctionalAssembler :
           * \param grid the subgrid (!) type
           * \param order the quadrature order (DEFAULT 2)
           */
-        DUNE_DEPRECATED_MSG("Please select quadrature for SubgridL2FunctionalAssembler using a QuadratureRuleKey.")
+        [[deprecated("Please select quadrature for SubgridL2FunctionalAssembler using a QuadratureRuleKey.")]]
         SubgridL2FunctionalAssembler(const Function& f, const GridType& grid, int order=2) :
             L2FunctionalAssembler<GridType,TrialLocalFE, T>(f, order),
             quadOrder_(order),

dune/fufem/assemblers/localassemblers/vvlaplaceassembler.hh

@@ -35,13 +35,13 @@ class VVLaplaceAssembler : public LocalOperatorAssembler <GridType, TrialLocalFE
       : one_(one)
     {}
 
-    void indices(const Element& element DUNE_UNUSED, BoolMatrix& isNonZero, const TrialLocalFE& tFE DUNE_UNUSED, const AnsatzLocalFE& aFE DUNE_UNUSED) const
+    void indices([[maybe_unused]] const Element& element, BoolMatrix& isNonZero, [[maybe_unused]] const TrialLocalFE& tFE, [[maybe_unused]] const AnsatzLocalFE& aFE) const
     {
       isNonZero = true;
     }
 
     template <class BoundaryIterator>
-    void indices(const BoundaryIterator& it DUNE_UNUSED, BoolMatrix& isNonZero, const TrialLocalFE& tFE DUNE_UNUSED, const AnsatzLocalFE& aFE DUNE_UNUSED) const
+    void indices([[maybe_unused]] const BoundaryIterator& it, BoolMatrix& isNonZero, [[maybe_unused]] const TrialLocalFE& tFE, [[maybe_unused]] const AnsatzLocalFE& aFE) const
     {
       isNonZero = true;
     }

dune/fufem/assemblers/localassemblers/vvmassassembler.hh

@@ -32,13 +32,13 @@ class VVMassAssembler : public LocalOperatorAssembler < GridType, TrialLocalFE,
       : one_(one)
     {}
 
-    void indices(const Element& element DUNE_UNUSED, BoolMatrix& isNonZero, const TrialLocalFE& tFE DUNE_UNUSED, const AnsatzLocalFE& aFE DUNE_UNUSED) const
+    void indices([[maybe_unused]] const Element& element, BoolMatrix& isNonZero, [[maybe_unused]] const TrialLocalFE& tFE, [[maybe_unused]] const AnsatzLocalFE& aFE) const
     {
       isNonZero = true;
     }
 
     template <class BoundaryIterator>
-    void indices(const BoundaryIterator& it DUNE_UNUSED, BoolMatrix& isNonZero, const TrialLocalFE& tFE DUNE_UNUSED, const AnsatzLocalFE& aFE DUNE_UNUSED) const
+    void indices([[maybe_unused]] const BoundaryIterator& it, BoolMatrix& isNonZero, [[maybe_unused]] const TrialLocalFE& tFE, [[maybe_unused]] const AnsatzLocalFE& aFE) const
     {
       isNonZero = true;
     }

dune/fufem/assemblers/localfunctionalassembler.hh

@@ -9,7 +9,6 @@
  *
  *  \tparam GridType The grid we are assembling on
  *  \tparam TrialLocalFE the local finite element of the trial space
- *  \tparam AnsatzLocalFE the local finite element of the ansatz space
  *  \tparam T the block type
  */
 template <class GridType, class TrialLocalFE, typename T>

dune/fufem/assemblers/transferoperatorassembler.hh

@@ -29,7 +29,7 @@ class MultilevelBasis
         MultilevelBasis(const GridType& grid) :
             grid(grid)
         {
-            const auto& globalIdSet = grid.globalIdSet();
+            const auto& localIdSet = grid.localIdSet();
             const auto& leafIndexSet = grid.leafIndexSet();
 
             int maxLevel = grid.maxLevel();
@@ -43,7 +43,7 @@ class MultilevelBasis
 
             size_[maxLevel] = grid.size(dim);
 
-            // build extended level indices as map (globalId -> index)
+            // build extended level indices as map (localId -> index)
             // 1. enumerate vertices on each level
             // 2. enumerate leaf vertices on each finer level to emulate copies of them
             // 3. on maxlevel take the leaf indices
@@ -62,16 +62,16 @@ class MultilevelBasis
                 const auto& indexSet = grid.levelIndexSet(level);
                 size_[level] = indexSet.size(dim);
                 for (const auto& it : vertices(grid.levelGridView(level)))
-                    idToIndex[level][globalIdSet.id(it)] = indexSet.index(it);
+                    idToIndex[level][localIdSet.id(it)] = indexSet.index(it);
             }
 
             for (const auto& it : vertices(grid.leafGridView()))
             {
-                idToIndex[maxLevel][globalIdSet.id(it)] = leafIndexSet.index(it);
+                idToIndex[maxLevel][localIdSet.id(it)] = leafIndexSet.index(it);
 
                 for(int level=it.level()+1; level<maxLevel; ++level)
                 {
-                    idToIndex[level][globalIdSet.id(it)] = size_[level];
+                    idToIndex[level][localIdSet.id(it)] = size_[level];
                     ++size_[level];
                 }
             }
@@ -90,7 +90,7 @@ class MultilevelBasis
         int index(const Element& e, const int i, const int level) const
         {
             const Dune::LocalKey& localKey = getLocalFiniteElement(e).localCoefficients().localKey(i);
-            const IdType id = grid.globalIdSet().subId(e, localKey.subEntity(), localKey.codim());
+            const IdType id = grid.localIdSet().subId(e, localKey.subEntity(), localKey.codim());
             return idToIndex[level].find(id)->second;
         }
 
@@ -99,7 +99,7 @@ class MultilevelBasis
         const GridType& grid;
         std::vector<std::shared_ptr<LevelBasis> > levelBasis_;
 
-        typedef typename GridType::Traits::GlobalIdSet::IdType IdType;
+        typedef typename GridType::Traits::LocalIdSet::IdType IdType;
         std::vector< std::map<IdType,int> > idToIndex;
         std::vector<int> size_;
 };
@@ -205,7 +205,7 @@ class TransferOperatorAssembler {
 
             MultilevelBasis<GridType> multiLevelBasis(grid);
 #ifdef FE_VERBOSE
-            std::cout << "FE:" << "globalId -> index maps build in " << timer.elapsed() << " seconds." << std::endl;
+            std::cout << "FE:" << "localId -> index maps build in " << timer.elapsed() << " seconds." << std::endl;
 #endif
 
 

dune/fufem/concept.hh

@@ -3,6 +3,8 @@
 #ifndef DUNE_FUFEM_CONCEPT_HH
 #define DUNE_FUFEM_CONCEPT_HH
 
+#warning This file is deprecated and will be removed after 2.10!  Please use dune/common/concept.hh instead!
+
 #include <type_traits>
 #include <utility>
 

dune/fufem/dgindexset.hh

@@ -21,7 +21,7 @@ public:
 
     void setup(const GridView& gridView) {
 
-        mapper_.update();
+        mapper_.update(gridView);
 
         const int dim = GridView::dimension;
 

dune/fufem/dgpqkindexset.hh

@@ -21,7 +21,7 @@ public:
 
         typedef typename GridView::template Codim<0>::Iterator ElementIterator;
         
-        mapper_.update();
+        mapper_.update(gridview);
 
         const int dim = GridView::dimension;
         

dune/fufem/dunedataio.hh

@@ -16,7 +16,7 @@ namespace serialization {
  * mark the use of the operator & that works both ways, << or >> , depending on the type of Archive.
  */
 template <class Archive, class T, int n>
-void serialize(Archive& ar, Dune::FieldVector<T,n>& fvec, const unsigned int version DUNE_UNUSED)
+void serialize(Archive& ar, Dune::FieldVector<T,n>& fvec, [[maybe_unused]] const unsigned int version)
 {
     for (int i=0; i<n; ++i)
         ar & fvec[i];
@@ -26,7 +26,7 @@ void serialize(Archive& ar, Dune::FieldVector<T,n>& fvec, const unsigned int ver
  * two separate save() and load() functions later to be "unified" via boost::serialization::split_free (see below)
  */
 template <class Archive, class BlockType>
-void save(Archive& ar, const Dune::BlockVector<BlockType>& vec, const unsigned int version DUNE_UNUSED)
+void save(Archive& ar, const Dune::BlockVector<BlockType>& vec, [[maybe_unused]] const unsigned int version)
 {
     size_t size = vec.size();
     ar << size;
@@ -35,7 +35,7 @@ void save(Archive& ar, const Dune::BlockVector<BlockType>& vec, const unsigned i
 }
 
 template <class Archive, class BlockType>
-void load(Archive& ar, Dune::BlockVector<BlockType>& vec, const unsigned int version DUNE_UNUSED)
+void load(Archive& ar, Dune::BlockVector<BlockType>& vec, [[maybe_unused]] const unsigned int version)
 {
     size_t size;
     ar >> size;