.gitlab-ci.yml

 ---
 # Install external dependencies
 before_script:
-  - patch -d /duneci/modules/dune-common -p1 < dune-common-densematrix-adolc-workaround.diff
-  - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-matrix-vector.git
-  - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-fufem.git
+  - duneci-install-module https://gitlab.dune-project.org/extensions/dune-vtk.git
+  - duneci-install-module https://gitlab.dune-project.org/fufem/dune-matrix-vector.git
+  - duneci-install-module https://gitlab.dune-project.org/fufem/dune-fufem.git
   - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-solvers.git
 
-dune:2.8 gcc C++17:
-  image: registry.dune-project.org/docker/ci/dune:2.8-debian-10-gcc-8-17
+dune:2.9 gcc C++20:
+  image: registry.dune-project.org/docker/ci/dune:2.9-debian-11-gcc-10-20
   script: duneci-standard-test
   variables:
-    DUNECI_BRANCH: releases/2.8
+    DUNECI_BRANCH: releases/2.9
 
-dune:2.8 gcc C++20:
-  image: registry.dune-project.org/docker/ci/dune:2.8-debian-11-gcc-9-20
+dune:2.9  clang  C++20:
+  image: registry.dune-project.org/docker/ci/dune:2.9-debian-11-clang-11-20
   script: duneci-standard-test
   variables:
-    DUNECI_BRANCH: releases/2.8
+    DUNECI_BRANCH: releases/2.9
 
-dune:2.8  clang  C++17:
-  image: registry.dune-project.org/docker/ci/dune:2.8-debian-10-clang-7-libcpp-17
+dune:git  clang-11  C++20:
+  image: registry.dune-project.org/docker/ci/dune:git-debian-11-clang-11-20
   script: duneci-standard-test
-  variables:
-    DUNECI_BRANCH: releases/2.8
 
-dune:git  clang  C++17:
-  image: registry.dune-project.org/docker/ci/dune:git-debian-10-clang-7-libcpp-17
+dune:git  gcc-10  C++20:
+  image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
   script: duneci-standard-test
 
-dune:git  gcc-8  C++17:
-  image: registry.dune-project.org/docker/ci/dune:git-debian-10-gcc-8-17
-  script: duneci-standard-test
+# Check for spelling mistakes in text
+code-spelling-check:
+  stage: .pre
+  # Avoid the global 'before_script'
+  before_script: ""
+  image: registry.dune-project.org/docker/ci/debian:11
+  script:
+  - codespell
+    --ignore-words-list afe,tre

CHANGELOG.md

-# Master (will become release 2.9)
+# Master (will be release 2.11)
+
+- The `LocalDensity` abstract base class and all its implementations are now
+  in a dedicated `densities` subdirectory.  This matches the structure
+  used in the `dune-gfe` module.
+
+- All classes and methods that were not already in the `Dune::Elasticity`
+  namespace have been moved there.
+
+- The file `directionalspringrobinfunction.hh` has been removed.
+  It was not used anywhere.  It used obsolete `dune-fufem` interfaces
+  and would have to be rewritten from scratch anyway.
+
+- The file `stresswriter.hh` has been removed.  It was based on the
+  AmiraMesh file format and the AmiraMesh writing support of _dune-grid_.
+  However, that support has been removed a long time ago, and therefore
+  the `StressWriter` class would need to be rewritten from scratch anyway.
+
+# Release 2.10
+
+- The class `Dune::FEAssembler` has been removed -- please use `Dune::Elasticity::FEAssembler`
+  from now on.  The removed code still used the old `dune-fufem` function space basis interface.
+
+- The class `LocalADOLCStiffness` (which used `dune-fufem`-style bases) has been removed.
+  Please use the class of the same name from the `Dune::Elasticity` namespace from now on.
+
+- The files `exphenckyenergy.hh`, `henckyenergy.hh`, `mooneyrivlinenergy.hh`,
+  `neohookeenergy.hh`, and `stvenantkirchhoffenergy.hh` have been removed.
+  Use the corresponding density classes together with `LocalIntegralEnergy` from now on.
+
+- The old implementations of `LocalEnergy` and `LocalIntegralEnergy` (the ones *not*
+  in the `Dune::Elasticity` namespace) have been removed.
+
+- The old implementation of `LocalFEStiffness` (the one *not* in the `Dune::Elasticity` namespace)
+  has been removed.
+
+- The module requires CMake 3.16 now.  This is what the Dune core modules currently require.
+
+# Release 2.9
 
 ## Deprecations and removals
 
-- The support of `Amiramesh` is dropped in dune-grid and thus the program `linear-elasticity` is removed from CMakeLists.txt
+- The support of `Amiramesh` is dropped in dune-grid and thus the program `linear-elasticity` is modified to compile with `Gmsh` but it will not run.
 
 # 2.8 Release
 

CMakeLists.txt

-# We require version CMake version 3.1 to prevent issues
-# with dune_enable_all_packages and older CMake versions.
-cmake_minimum_required(VERSION 3.1)
+cmake_minimum_required(VERSION 3.16)
 project(dune-elasticity CXX)
 
 if(NOT (dune-common_DIR OR dune-common_ROOT OR

README

-Getting started
-===============
-
-You have to create the configure-script on your own!
-
-First, you'll need the followings programs installed:
-
-  automake >= 1.5
-
-  autoconf >= 2.50
-
-  libtool
-
-Then run
-
-  ./autogen.sh DUNEDIR
-
-where DUNEDIR is the (absolute or relative) path of the dune/-directory.
-
-The directory is needed because autogen.sh needs access to the tests
-stored in DUNEDIR/m4
-
-Now call
-
-  ./configure --with-dune=DIR
-
-where DIR is the directory _above_ dune/. If you want to include third
-party packages check
-
-  ./configure --help
-
-for options on how to include Albert, Grape, UG, ...
-
-If configure checked your system without problems you can use
-
-  make
-
-to build all examples.
-

dune-common-densematrix-adolc-workaround.diff

-diff --git a/dune/common/densematrix.hh b/dune/common/densematrix.hh
-index b03bbb0b..917ecef9 100644
---- a/dune/common/densematrix.hh
-+++ b/dune/common/densematrix.hh
-@@ -897,7 +897,7 @@ namespace Dune
-         for (size_type k=i+1; k<A.rows(); k++)
-         {
-           auto abs = fvmeta::absreal(A[k][i]);
--          auto mask = abs > pivmax;
-+          bool mask = abs > pivmax;
-           pivmax = Simd::cond(mask, abs, pivmax);
-           imax   = Simd::cond(mask, simd_index_type(k), imax);
-         }

dune.module

@@ -4,8 +4,8 @@
 
 #Name of the module
 Module:dune-elasticity
-Version: 2.9-git
+Version: 2.11-git
 Maintainer: oliver.sander@tu-dresden.de, patrick.jaap@tu-dresden.de
 #depending on
-Depends:dune-common dune-grid dune-istl dune-solvers dune-fufem dune-geometry dune-functions
+Depends:dune-common dune-grid dune-istl dune-solvers dune-fufem dune-geometry dune-functions dune-vtk
 Suggests: dune-uggrid dune-parmg

dune/elasticity/CMakeLists.txt

 add_subdirectory(assemblers)
 add_subdirectory(common)
+add_subdirectory(densities)
 add_subdirectory(estimators)
 add_subdirectory(materials)
-add_subdirectory(utilities)

dune/elasticity/assemblers/CMakeLists.txt

 install(FILES
     feassembler.hh
+    functionalassembler.hh
+    geomexactstvenantkirchhofffunctionalassembler.hh
+    geomexactstvenantkirchhoffoperatorassembler.hh
+    localadolcintegralstiffness.hh
     localadolcstiffness.hh
+    localfunctionalassembler.hh
     localhyperdualstiffness.hh
+    localoperatorassembler.hh
+    operatorassembler.hh
     localenergy.hh
     localfestiffness.hh
+    localstiffnesssum.hh
     neohookefunctionalassembler.hh
     neohookeoperatorassembler.hh
     ogdenfunctionalassembler.hh

dune/elasticity/assemblers/feassembler.hh

@@ -2,8 +2,7 @@
 #define DUNE_ELASTICITY_ASSEMBLERS_FEASSEMBLER_HH
 
 #include <dune/common/fmatrix.hh>
-
-#include <dune/fufem/assemblers/istlbackend.hh>
+#include <dune/common/version.hh>
 
 #include <dune/istl/bvector.hh>
 #include <dune/istl/bcrsmatrix.hh>
@@ -15,6 +14,12 @@
 
 #include <dune/grid/common/partitionset.hh>
 
+#if DUNE_VERSION_GT(DUNE_COMMON, 2, 9)
+#include <dune/fufem/backends/istlmatrixbackend.hh>
+#else
+#include <dune/fufem/assemblers/istlbackend.hh>
+#endif
+
 #include "localfestiffness.hh"
 
 
@@ -54,6 +59,11 @@ public:
           localStiffness_(localStiffness)
     {}
 
+    /** \brief Assemble the functional gradient */
+    void assembleGradient(const VectorType& configuration,
+                          VectorType& gradient) const;
+
+
     /** \brief Assemble the tangent stiffness matrix and the functional gradient together
      *
      * This may be more efficient than computing them separately
@@ -71,6 +81,47 @@ public:
 }; // end class
 
 
+
+template <class Basis, class VectorType>
+void FEAssembler<Basis,VectorType>::
+assembleGradient(const VectorType& configuration,
+                 VectorType& gradient) const
+{
+  // create backends for multi index access
+  auto configurationBackend      = Dune::Functions::istlVectorBackend(configuration);
+  auto gradientBackend = Dune::Functions::istlVectorBackend(gradient);
+
+  gradientBackend.resize(basis_);
+
+  gradient = 0;
+
+  auto localView = basis_.localView();
+
+  for (const auto& element : elements(basis_.gridView(), partitionType))
+  {
+    localView.bind(element);
+    const auto size = localView.size();
+
+    // Extract local configuration
+    std::vector<field_type> localConfiguration(size);
+    std::vector<field_type> localGradient(size);
+
+    for (size_t i=0; i<size; i++)
+      localConfiguration[i] = configurationBackend[ localView.index(i) ];
+
+    // setup local matrix and gradient
+    localStiffness_->assembleGradient(localView, localConfiguration, localGradient);
+
+    for (size_t i=0; i<size; i++)
+    {
+      auto row = localView.index(i);
+      gradientBackend[row] += localGradient[i];
+    }
+  }
+}
+
+
+
 template <class Basis, class VectorType>
 template <class MatrixType>
 void FEAssembler<Basis,VectorType>::
@@ -172,188 +223,4 @@ computeEnergy(const VectorType& configuration) const
 
 } // namespace Dune::Elasticity
 
-
-
-namespace Dune
-{
-
-/** \brief A global FE assembler for variational problems (old fufem bases version)
- */
-template <class Basis, class VectorType >
-class
-FEAssembler
-{
-
-    typedef typename Basis::GridView GridView;
-
-    //! Dimension of the grid.
-    enum { gridDim = GridView::dimension };
-
-    //! Dimension of a tangent space
-    enum { blocksize = VectorType::value_type::dimension };
-
-    //!
-    typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
-
-public:
-    [[deprecated("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::FEAssembler with dune-functions bases.")]]
-    const Basis basis_;
-
-    /** \brief Partition type on which to assemble
-     *
-     * We want a fully nonoverlapping partition, and therefore need to skip ghost elements.
-     */
-    static constexpr auto partitionType = Dune::Partitions::interiorBorder;
-
-protected:
-
-    LocalFEStiffness<GridView,
-                             typename Basis::LocalFiniteElement,
-                             VectorType>* localStiffness_;
-
-public:
-
-    /** \brief Constructor for a given grid */
-    FEAssembler(const Basis& basis,
-                LocalFEStiffness<GridView,typename Basis::LocalFiniteElement, VectorType>* localStiffness)
-        : basis_(basis),
-          localStiffness_(localStiffness)
-    {}
-
-    /** \brief Assemble the tangent stiffness matrix and the functional gradient together
-     *
-     * This may be more efficient than computing them separately
-     */
-    virtual void assembleGradientAndHessian(const VectorType& sol,
-                                            Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
-                                            Dune::BCRSMatrix<MatrixBlock>& hessian,
-                                            bool computeOccupationPattern=true) const;
-
-    /** \brief Compute the energy of a deformation state */
-    virtual double computeEnergy(const VectorType& sol) const;
-
-    //protected:
-    void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const;
-
-}; // end class
-
-
-
-template <class Basis, class TargetSpace>
-void FEAssembler<Basis,TargetSpace>::
-getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const
-{
-    int n = basis_.size();
-
-    nb.resize(n, n);
-
-    for (const auto& element : elements(basis_.getGridView(), partitionType))
-    {
-        const typename Basis::LocalFiniteElement& lfe = basis_.getLocalFiniteElement(element);
-
-        for (size_t i=0; i<lfe.localBasis().size(); i++) {
-
-            for (size_t j=0; j<lfe.localBasis().size(); j++) {
-
-                int iIdx = basis_.index(element,i);
-                int jIdx = basis_.index(element,j);
-
-                nb.add(iIdx, jIdx);
-
-            }
-
-        }
-
-    }
-
-}
-
-template <class Basis, class VectorType>
-void FEAssembler<Basis,VectorType>::
-assembleGradientAndHessian(const VectorType& sol,
-                           Dune::BlockVector<Dune::FieldVector<double, blocksize> >& gradient,
-                           Dune::BCRSMatrix<MatrixBlock>& hessian,
-                           bool computeOccupationPattern) const
-{
-    if (computeOccupationPattern) {
-
-        Dune::MatrixIndexSet neighborsPerVertex;
-        getNeighborsPerVertex(neighborsPerVertex);
-        neighborsPerVertex.exportIdx(hessian);
-
-    }
-
-    hessian = 0;
-
-    gradient.resize(sol.size());
-    gradient = 0;
-
-    for (const auto& element : elements(basis_.getGridView(), partitionType))
-    {
-        const int numOfBaseFct = basis_.getLocalFiniteElement(element).localBasis().size();
-
-        // Extract local solution
-        VectorType localSolution(numOfBaseFct);
-        VectorType localPointLoads(numOfBaseFct);
-
-        for (int i=0; i<numOfBaseFct; i++)
-            localSolution[i]   = sol[basis_.index(element,i)];
-
-        VectorType localGradient(numOfBaseFct);
-
-        // setup local matrix and gradient
-        localStiffness_->assembleGradientAndHessian(element, basis_.getLocalFiniteElement(element), localSolution, localGradient);
-
-        // Add element matrix to global stiffness matrix
-        for(int i=0; i<numOfBaseFct; i++) {
-
-            int row = basis_.index(element,i);
-
-            for (int j=0; j<numOfBaseFct; j++ ) {
-
-                int col = basis_.index(element,j);
-                hessian[row][col] += localStiffness_->A_[i][j];
-
-            }
-        }
-
-        // Add local gradient to global gradient
-        for (int i=0; i<numOfBaseFct; i++)
-            gradient[basis_.index(element,i)] += localGradient[i];
-
-    }
-
-}
-
-
-template <class Basis, class VectorType>
-double FEAssembler<Basis, VectorType>::
-computeEnergy(const VectorType& sol) const
-{
-    double energy = 0;
-
-    if (sol.size()!=basis_.size())
-        DUNE_THROW(Dune::Exception, "Solution vector doesn't match the basis!");
-
-    // Loop over all elements
-    for (const auto& element : elements(basis_.getGridView(), partitionType))
-    {
-        // Number of degrees of freedom on this element
-        size_t nDofs = basis_.getLocalFiniteElement(element).localBasis().size();
-
-        VectorType localSolution(nDofs);
-
-        for (size_t i=0; i<nDofs; i++)
-            localSolution[i]   = sol[basis_.index(element,i)];
-
-        energy += localStiffness_->energy(element, basis_.getLocalFiniteElement(element), localSolution);
-
-    }
-
-    return energy;
-
-}
-
-} // namespace Dune
-
 #endif

dune/elasticity/assemblers/functionalassembler.hh

+#ifndef FUNCTIONAL_ASSEMBLER_HH
+#define FUNCTIONAL_ASSEMBLER_HH
+
+#warning functionalassembler.hh is deprecated!
+
+#include <dune/istl/bvector.hh>
+
+#include "dune/fufem/functionspacebases/functionspacebasis.hh"
+
+//! Generic global assembler for functionals on a gridview
+template <class TestBasis>
+class FunctionalAssembler
+{
+    private:
+        typedef typename TestBasis::GridView GridView;
+
+    public:
+        //! create assembler for gridview
+        FunctionalAssembler(const TestBasis& tBasis) :
+            tBasis_(tBasis)
+        {}
+
+        /** \brief Assemble
+         *
+         * \param localAssembler local assembler
+         * \param[out] b target vector
+         * \param initializeVector If this is set the output vector is
+         * set to the correct length and initialized to zero before
+         * assembly. Otherwise the assembled values are just added to
+         * the vector.
+         */
+        template <class LocalFunctionalAssemblerType, class GlobalVectorType>
+        void assemble(LocalFunctionalAssemblerType& localAssembler, GlobalVectorType& b, bool initializeVector=true) const
+        {
+            typedef typename LocalFunctionalAssemblerType::LocalVector LocalVector;
+            typedef typename TestBasis::LinearCombination LinearCombination;
+
+            int rows = tBasis_.size();
+
+            if (initializeVector)
+            {
+                b.resize(rows);
+                b=0.0;
+            }
+
+            for (const auto& element : elements(tBasis_.getGridView()))
+            {
+                // get shape functions
+                const typename TestBasis::LocalFiniteElement& tFE = tBasis_.getLocalFiniteElement(element);
+
+                LocalVector localB(tFE.localBasis().size());
+                localAssembler.assemble(element, localB, tFE);
+
+                for (size_t i=0; i<tFE.localBasis().size(); ++i)
+                {
+                    int idx = tBasis_.index(element, i);
+                    const LinearCombination& constraints = tBasis_.constraints(idx);
+                    bool isConstrained = tBasis_.isConstrained(idx);
+                    if (isConstrained)
+                    {
+                        for (size_t w=0; w<constraints.size(); ++w)
+                            b[constraints[w].index].axpy(constraints[w].factor, localB[i]);
+                    }
+                    else
+                        b[idx] += localB[i];
+                }
+            }
+            return;
+        }
+
+    private:
+        const TestBasis& tBasis_;
+};
+
+#endif
+

dune/elasticity/assemblers/geomexactstvenantkirchhofffunctionalassembler.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef DUNE_ELASTICITY_ASSEMBLERS_GEOMETRICALLY_EXACT_ST_VENANT_KIRCHHOFF_FUNCTIONAL_ASSEMBLER_HH
+#define DUNE_ELASTICITY_ASSEMBLERS_GEOMETRICALLY_EXACT_ST_VENANT_KIRCHHOFF_FUNCTIONAL_ASSEMBLER_HH
+
+#include <array>
+#include <memory>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
+#include <dune/fufem/symmetrictensor.hh>
+#include <dune/fufem/functions/virtualgridfunction.hh>
+
+#include <dune/elasticity/assemblers/localfunctionalassembler.hh>
+
+namespace Dune::Elasticity
+{
+
+/** \brief Local assembler of the first derivative of the geometrically exact St.Venant material at a displacement u
+ *
+ *  \f[
+ *      J(u) = \int_{\Omega} E(u):C:E(u) dx
+ *  \f]
+ *
+ *  which results in the linearform
+ *
+ *    \f[
+ *      l_u(v) := \int_{\Omega} D\varphi(u)\sigma(u):Dv dx = \int_{\Omega}\sigma(u): E'_u (v) dx
+ *    \f]
+ *   with
+ *      - \f$D\varphi\f$: the deformation gradient
+ *      - \f$E\f$: the nonlinear strain tensor
+ *      - \f$E'_u(v)\f$: the linearized strain tensor at the point u in direction v
+ *      - \f$\sigma(u) = C:E(u)\f$: the second Piola-Kirchhoff stress tensor
+ *      - \f$C\f$: the Hooke tensor
+ *      - \f$Dv\f$: the gradient of the test function
+*/
+template <class GridType, class TestLocalFE>
+class GeomExactStVenantKirchhoffFunctionalAssembler :
+    public LocalFunctionalAssembler<GridType, TestLocalFE, Dune::FieldVector<typename GridType::ctype ,GridType::dimension> >
+
+{
+    static const int dim = GridType::dimension;
+    typedef typename GridType::ctype ctype;
+    typedef typename Dune::FieldVector<ctype,dim> T;
+    typedef LocalFunctionalAssembler<GridType, TestLocalFE, T> Base;
+
+public:
+    typedef typename Base::Element Element;
+    typedef typename Base::LocalVector LocalVector;
+    typedef VirtualGridFunction<GridType, T> GridFunction;
+
+     //! Constructor
+    GeomExactStVenantKirchhoffFunctionalAssembler(ctype E,ctype nu, std::shared_ptr<GridFunction> configuration)
+       : E_(E), nu_(nu),
+        configuration_(configuration)
+    {}
+
+
+    //! Constructor
+    GeomExactStVenantKirchhoffFunctionalAssembler(ctype E, ctype nu) :
+        E_(E), nu_(nu), configuration_(NULL)
+    {}
+
+    void assemble(const Element& element, LocalVector& localVector, const TestLocalFE& tFE) const
+    {
+
+        typedef typename Dune::FieldMatrix<ctype,dim,dim> FMdimdim;
+        typedef typename TestLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType;
+        typedef typename Element::Geometry::LocalCoordinate LocalCoordinate;
+        int rows = tFE.localBasis().size();
+        localVector = 0.0;
+
+        // get quadrature rule
+        const int order = (element.type().isSimplex())
+            ? (tFE.localBasis().order())
+            : (tFE.localBasis().order()*dim);
+
+        // get quadrature rule
+        const Dune::template QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(element.type(), order, IsRefinedLocalFiniteElement<TestLocalFE>::value(tFE) );
+
+        // store gradients of shape functions and base functions
+        std::vector<JacobianType> referenceGradients(tFE.localBasis().size());
+        std::vector<T> gradients(tFE.localBasis().size());
+
+        // store geometry mapping of the entity
+        const typename Element::Geometry geometry = element.geometry();
+
+        // loop over quadrature points
+        for (size_t pt=0; pt < quad.size(); ++pt) {
+
+            // get quadrature point
+            const LocalCoordinate& quadPos = quad[pt].position();
+
+            // get transposed inverse of Jacobian of transformation
+            const FMdimdim& invJacobian = geometry.jacobianInverseTransposed(quadPos);
+
+            // get integration factor
+            const ctype integrationElement = geometry.integrationElement(quadPos);
+
+            // get gradients of shape functions
+            tFE.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+            // compute gradients of base functions
+            for (size_t i=0; i<gradients.size(); ++i) {
+
+                // transform gradients
+                gradients[i] = 0.0;
+                invJacobian.umv(referenceGradients[i][0], gradients[i]);
+            }
+
+            // evaluate the displacement gradients of the configuration at the quadrature point
+            typename GridFunction::DerivativeType localConfGrad;
+            if (configuration_->isDefinedOn(element))
+                configuration_->evaluateDerivativeLocal(element, quadPos, localConfGrad);
+            else
+                configuration_->evaluateDerivative(geometry.global(quadPos),localConfGrad);
+
+
+            // /////////////////////////////////////////////
+            //   Compute stress of configuration at quad point
+            // /////////////////////////////////////////////
+
+            /* Compute the nonlinear strain tensor from the displacement gradient*/
+            SymmetricTensor<dim> strain;
+            computeStrain(localConfGrad,strain);
+
+            // and the stress
+            SymmetricTensor<dim> stress = hookeTimesStrain(strain);
+
+            // make deformation gradient out of the discplacement
+            for (int i=0;i<dim;i++)
+                localConfGrad[i][i] += 1;
+
+            //compute linearized strain for all shapefunctions
+            std::vector<std::array<SymmetricTensor<dim>,dim> > linStrain(rows);
+
+            for (int i=0; i<rows; i++)
+                for (int k=0; k<dim; k++) {
+
+                    // The deformation gradient of the shape function
+                    FMdimdim deformationGradient(0);
+                    deformationGradient[k] = gradients[i];
+
+                    // The linearized strain is the symmetric product of defGrad and (Id+localConf)
+                    linStrain[i][k]=symmetricProduct(deformationGradient,localConfGrad);
+                }
+
+            ctype z = quad[pt].weight()*integrationElement;
+
+            // add vector entries
+            for(int row=0; row<rows; row++)
+                for (int rcomp=0;rcomp<dim;rcomp++)
+                    localVector[row][rcomp] +=(stress*linStrain[row][rcomp])*z;
+
+        }
+        return;
+    }
+
+    /** \brief Set new configuration. In Newton iterations this needs to be assembled more than one time.*/
+    void setConfiguration(std::shared_ptr<GridFunction> newConfig) {
+        configuration_ = newConfig;
+
+        if (!configuration_)
+            DUNE_THROW(Dune::Exception,"In [GeomExactStVenantKirchhoffFunctionalAssembler]: GridFunction cast failed!");
+    }
+
+protected:
+    /** \brief Young's modulus */
+    ctype E_;
+
+    /** \brief The Poisson ratio */
+    ctype nu_;
+
+    /** \brief The configuration at which the functional is evaluated.*/
+    std::shared_ptr<GridFunction> configuration_;
+
+    /** \brief Compute nonlinear strain tensor from the deformation gradient. */
+    void computeStrain(const Dune::FieldMatrix<ctype, dim, dim>& grad, SymmetricTensor<dim>& strain) const {
+        strain = 0;
+        for (int i=0; i<dim ; ++i) {
+            strain(i,i) +=grad[i][i];
+
+            for (int k=0;k<dim;++k)
+                strain(i,i) += 0.5*grad[k][i]*grad[k][i];
+
+            for (int j=i+1; j<dim; ++j) {
+                strain(i,j) += 0.5*(grad[i][j] + grad[j][i]);
+                for (int k=0;k<dim;++k)
+                    strain(i,j) += 0.5*grad[k][i]*grad[k][j];
+            }
+        }
+    }
+
+    /** \brief Compute the stress tensor from the nonlinear strain tensor. */
+    SymmetricTensor<dim> hookeTimesStrain(const SymmetricTensor<dim>& strain) const {
+
+        SymmetricTensor<dim> stress = strain;
+        stress *= E_/(1+nu_);
+
+        ctype f = E_*nu_ / ((1+nu_)*(1-2*nu_)) * strain.trace();
+
+        stress.addToDiag(f);
+
+        return stress;
+    }
+
+    /** \brief Compute the symmetric product of two matrices, i.e.
+     *      0.5*(A^T * B + B^T * A )
+     */
+    SymmetricTensor<dim> symmetricProduct(const Dune::FieldMatrix<ctype, dim, dim>& a, const Dune::FieldMatrix<ctype, dim, dim>& b) const {
+
+        SymmetricTensor<dim> result(0.0);
+        for (int i=0;i<dim; i++)
+            for (int j=i;j<dim;j++)
+                for (int k=0;k<dim;k++)
+                    result(i,j)+= a[k][i]*b[k][j]+b[k][i]*a[k][j];
+        result *= 0.5;
+
+        return result;
+    }
+};
+
+} // namespace Dune::Elasticity
+
+#endif

dune/elasticity/assemblers/geomexactstvenantkirchhoffoperatorassembler.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef DUNE_ELASTICITY_ASSEMBLERS_GEOMETRICALLY_EXACT_ST_VENANT_KIRCHHOFF_OPERATOR_ASSEMBLER_HH
+#define DUNE_ELASTICITY_ASSEMBLERS_GEOMETRICALLY_EXACT_ST_VENANT_KIRCHHOFF_OPERATOR_ASSEMBLER_HH
+
+#include <array>
+#include <memory>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/istl/matrix.hh>
+
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
+#include <dune/fufem/staticmatrixtools.hh>
+#include <dune/fufem/symmetrictensor.hh>
+#include <dune/fufem/functions/virtualgridfunction.hh>
+#include <dune/matrix-vector/addtodiagonal.hh>
+
+#include <dune/elasticity/assemblers/localoperatorassembler.hh>
+
+namespace Dune::Elasticity
+{
+
+/** \brief Assemble the second derivative of the geometrically exact St.-Venant-Kirchhoff material
+ *
+ *  \f[
+ *      J(u) = \int_{\Omega} E(u):C:E(u) dx
+ *  \f]
+ *
+ *  which results in the bilinearform
+ *
+ *    \f[
+ *      a_u(v,w) = \int_{\Omega} \sigma(u)E''_u(v,w) + (C:E'_u(v)):E'_u(w) dx
+ *    \f]
+ *   with
+ *      - \f$D\varphi\f$: the deformation gradient
+ *      - \f$E\f$: the nonlinear strain tensor
+ *      - \f$E'_u(v)\f$: the linearized strain tensor at the point u in direction v
+ *      - \f$\sigma(u) = C:E(u)\f$: the second Piola-Kirchhoff stress tensor
+ *      - \f$C\f$: the Hooke tensor
+ *      - \f$Dv\f$: the gradient of the test function
+ *
+ */
+template <class GridType, class TestLocalFE, class AnsatzLocalFE>
+class GeomExactStVenantKirchhoffOperatorAssembler
+    : public LocalOperatorAssembler < GridType, TestLocalFE, AnsatzLocalFE, Dune::FieldMatrix<typename GridType::ctype ,GridType::dimension,GridType::dimension> >
+{
+
+
+    static const int dim = GridType::dimension;
+    typedef typename GridType::ctype ctype;
+    using field_type = typename AnsatzLocalFE::Traits::LocalBasisType::Traits::RangeFieldType;
+
+    typedef VirtualGridFunction<GridType, Dune::FieldVector<ctype, dim> > GridFunction;
+
+    typedef Dune::FieldMatrix<field_type, dim, dim> T;
+    using Base = LocalOperatorAssembler < GridType, TestLocalFE, AnsatzLocalFE, T >;
+
+public:
+    using MatrixEntry = typename Base::MatrixEntry;
+    using Element = typename Base::Element;
+    using BoolMatrix = typename Base::BoolMatrix;
+    using LocalMatrix = typename Base::LocalMatrix;
+
+   GeomExactStVenantKirchhoffOperatorAssembler(ctype E, ctype nu, std::shared_ptr<GridFunction> configuration)
+        : E_(E), nu_(nu), configuration_(configuration)
+    {}
+
+    GeomExactStVenantKirchhoffOperatorAssembler(ctype E, ctype nu) :
+        E_(E), nu_(nu)
+    {}
+
+    void indices(const Element& element, BoolMatrix& isNonZero, const TestLocalFE& tFE, const AnsatzLocalFE& aFE) const
+    {
+        isNonZero = true;
+    }
+
+    void assemble(const Element& element, LocalMatrix& localMatrix, const TestLocalFE& tFE, const AnsatzLocalFE& aFE) const
+    {
+
+        int rows = tFE.localBasis().size();
+        int cols = aFE.localBasis().size();
+
+        localMatrix.setSize(rows,cols);
+        localMatrix = 0.0;
+
+        const int order = (element.type().isSimplex())
+            ? (tFE.localBasis().order())
+            : (tFE.localBasis().order()*dim);
+        const auto& quad = QuadratureRuleCache<ctype, dim>::rule(element.type(), order, IsRefinedLocalFiniteElement<TestLocalFE>::value(tFE) );
+
+        using ReferenceLfeJacobian = typename TestLocalFE::Traits::LocalBasisType::Traits::JacobianType;
+        std::vector<ReferenceLfeJacobian> referenceGradients(tFE.localBasis().size());
+
+        const auto& geometry = element.geometry();
+
+        for (size_t pt=0; pt < quad.size(); ++pt) {
+
+            const auto& quadPos = quad[pt].position();
+
+            const auto& invJacobian = geometry.jacobianInverseTransposed(quadPos);
+
+            tFE.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+            using GlobalLfeJacobian = Dune::FieldVector<ctype, dim>;
+            std::vector<GlobalLfeJacobian> gradients(tFE.localBasis().size());
+            for (size_t i=0; i < gradients.size(); ++i) {
+                // transform gradients
+                gradients[i] = 0.0;
+                invJacobian.umv(referenceGradients[i][0], gradients[i]);
+            }
+
+            using Jacobian = typename GridFunction::DerivativeType;
+            Jacobian localConfGrad;
+            if (configuration_->isDefinedOn(element))
+                configuration_->evaluateDerivativeLocal(element, quadPos, localConfGrad);
+            else
+                configuration_->evaluateDerivative(geometry.global(quadPos), localConfGrad);
+
+            SymmetricTensor<dim, field_type> strain;
+            computeStrain(localConfGrad, strain);
+            SymmetricTensor<dim, field_type> stress = hookeTimesStrain(strain);
+
+            Dune::MatrixVector::addToDiagonal(localConfGrad, 1.0);
+
+            std::vector<std::array<Jacobian, dim > > defJac(rows);
+            std::vector<std::array<SymmetricTensor<dim>,dim> > linStrain(rows);
+
+            for (int i=0; i<rows; i++) {
+
+                for (int k=0; k<dim; k++) {
+
+                    defJac[i][k] = 0;
+                    defJac[i][k][k] = gradients[i];
+
+                    // The linearized strain is the symmetric product of defGrad and (Id+localConf)
+                    linStrain[i][k]=symmetricProduct(defJac[i][k], localConfGrad);
+                }
+            }
+
+            const auto integrationElement = geometry.integrationElement(quadPos);
+            ctype z = quad[pt].weight() * integrationElement;
+            for (int row=0; row < rows; row++)
+                for (int rcomp=0; rcomp<dim; rcomp++) {
+
+                    auto linStress = hookeTimesStrain(linStrain[row][rcomp]);
+
+                    for (int col=0; col < cols; col++)
+                        for (int ccomp=0; ccomp < dim; ccomp++) {
+                            localMatrix[row][col][rcomp][ccomp] += (linStress*linStrain[col][ccomp])*z;
+                            localMatrix[row][col][rcomp][ccomp] += (stress*symmetricProduct(defJac[row][rcomp], defJac[col][ccomp]))*z;
+                        }
+                }
+        }
+    }
+
+    /** \brief Set new configuration to be assembled at. Needed e.g. during Newton iteration.*/
+    void setConfiguration(std::shared_ptr<GridFunction> newConfig) {
+        configuration_ = newConfig;
+        if (!configuration_)
+            DUNE_THROW(Dune::Exception,"In [GeomExactStVenantKirchhoffOperatorAssembler]: You need to provide a GridFunction describing the displacement!");
+    }
+
+protected:
+    /** \brief Young's modulus */
+    ctype E_;
+
+    /** \brief The Poisson ratio */
+    ctype nu_;
+
+    /** \brief The configuration at which the linearized operator is evaluated.*/
+    std::shared_ptr<GridFunction> configuration_;
+
+    /** \brief Compute nonlinear strain tensor from the deformation gradient. */
+    void computeStrain(const Dune::FieldMatrix<ctype, dim, dim>& grad, SymmetricTensor<dim>& strain) const {
+        strain = 0;
+        for (int i=0; i<dim ; ++i) {
+            strain(i,i) +=grad[i][i];
+
+            for (int k=0;k<dim;++k)
+                strain(i,i) += 0.5*grad[k][i]*grad[k][i];
+
+            for (int j=i+1; j<dim; ++j) {
+                strain(i,j) += 0.5*(grad[i][j] + grad[j][i]);
+                for (int k=0;k<dim;++k)
+                    strain(i,j) += 0.5*grad[k][i]*grad[k][j];
+            }
+        }
+    }
+
+    /** \brief Compute the stress tensor from the nonlinear strain tensor. */
+    SymmetricTensor<dim> hookeTimesStrain(const SymmetricTensor<dim>& strain) const {
+
+        SymmetricTensor<dim> stress = strain;
+        stress *= E_/(1+nu_);
+
+        ctype f = E_*nu_ / ((1+nu_)*(1-2*nu_)) * strain.trace();
+
+        stress.addToDiag(f);
+
+        return stress;
+    }
+
+    /** \brief Compute the symmetric product of two matrices, i.e.
+     *      0.5*(A^T * B + B^T * A )
+     */
+    SymmetricTensor<dim> symmetricProduct(const Dune::FieldMatrix<ctype, dim, dim>& a, const Dune::FieldMatrix<ctype, dim, dim>& b) const {
+
+        SymmetricTensor<dim> result(0.0);
+        for (int i=0;i<dim; i++)
+            for (int j=i;j<dim;j++)
+                for (int k=0;k<dim;k++)
+                    result(i,j)+= a[k][i]*b[k][j]+b[k][i]*a[k][j];
+        result *= 0.5;
+
+        return result;
+    }
+
+};
+
+} // namespace Dune::Elasticity
+
+#endif
+

dune/elasticity/assemblers/localadolcintegralstiffness.hh

+#ifndef DUNE_ELASTICITY_ASSEMBLERS_LOCALADOLCINTEGRALSTIFFNESS_HH
+#define DUNE_ELASTICITY_ASSEMBLERS_LOCALADOLCINTEGRALSTIFFNESS_HH
+
+#include <adolc/adouble.h>            // use of active doubles
+
+#include <dune/common/fmatrix.hh>
+
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/istl/matrix.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+#include <dune/elasticity/densities/localdensity.hh>
+
+namespace Dune::Elasticity {
+
+/** \brief This class assembles the energy gradient and Hessian for an integral over a given density function.
+
+           For the gradient:
+           This class takes the derivative of the density function (using ADOL-C), evaluates this at the shape functions,
+           multiplies it with the derivative of the shape functions, and then integrates this over the local element
+
+           For the Hessian:
+           This class calculates the Hessian of the density function (using ADOL-C), evaluates this at the shape functions,
+           multiplies it with the derivative of the shape functions (left and right), and then integrates this over the local element
+ */
+template<class LocalView>
+class LocalADOLCIntegralStiffness
+    : public LocalFEStiffness<LocalView,double>
+{
+    static constexpr auto gridDim = LocalView::GridView::dimension;
+    static constexpr auto derivativeSize = gridDim*gridDim; //derivativeSize is the size of an argument of the LocalDensity
+    static constexpr auto targetDim = gridDim;
+    using FiniteElement = typename LocalView::Tree::ChildType::FiniteElement;
+    using Jacobian = typename FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
+    using ctype = typename LocalView::GridView::Grid::ctype;
+
+  public:
+
+    /** \brief Constructor from a local density function
+     *  \param[in]  ld            Pointer to the local density function which we take the derivative of
+     *  \param[in]  vectorMode    Boolean indicator whether to use the scalar or vector mode of ADOL-C
+     *
+     *  \note       There are different basic drivers for computing the hessian by ADOL-C, namely "hessian" and "hessian2",
+     *              where hessian is the scalar mode of ADOL-C, hessian2 is the vector mode.
+     *              They yield the same results, but "hessian2" seems to be faster. The user can switch to "hessian" with
+     *              vectorMode=false. See ADOL-C spec for more info about the drivers 
+   **/
+  LocalADOLCIntegralStiffness(const std::shared_ptr<LocalDensity<gridDim,adouble,ctype>>& ld, bool vectorMode = true)
+    : localDensity_(ld)
+    , vectorMode_(vectorMode)
+    {}
+      
+    virtual ~LocalADOLCIntegralStiffness() {}
+
+    virtual ctype energy (const LocalView& localView,
+                          const std::vector<ctype>& localConfiguration) const override final;
+
+    /** \brief Assemble the local gradient at the current position
+     */
+    virtual void assembleGradient(const LocalView& localView,
+                                  const std::vector<ctype>& localConfiguration,
+                                  std::vector<ctype>& localGradient);
+
+    /** \brief Assemble the local tangent matrix at the current position
+     */
+    virtual void assembleGradientAndHessian(const LocalView& localView,
+                                            const std::vector<ctype>& localConfiguration,
+                                            std::vector<ctype>& localGradient,
+                                            Matrix<ctype>& localHessian);
+
+  protected:
+    const std::shared_ptr<LocalDensity<gridDim,adouble,ctype>> localDensity_ = nullptr;
+    const bool vectorMode_;
+
+  private:
+    /** \brief Trace the call to the density function, evaluated at the current quadrature point
+               with the localConfiguration to construct the deformation finite element function
+     *  \param[in]  localView                    The local view
+     *  \param[in]  localConfiguration           The local configuration
+     *  \param[in]  qp                           The quadrature point
+     *  \param[out] jacobians                    Jacobians of the shape functions evaluated at the current quadrature point 
+     *  \param[out] localDeformationGradientVec  Gradient of the deformation finite element function evaluated at the current quadrature point 
+     */
+    void evaluateActiveDensity(const LocalView& localView,
+                               const std::vector<ctype>& localConfiguration,
+                               const QuadraturePoint<ctype,gridDim>& qp,
+                               std::vector<Jacobian>& jacobians,
+                               std::vector<ctype>& localDeformationGradientVec)
+    {
+      int rank = MPIHelper::getCommunication().rank();
+
+      const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+      const auto& element = localView.element();
+
+      auto x = element.geometry().global(qp.position());
+
+      localFiniteElement.localBasis().evaluateJacobian(qp.position(), jacobians);
+
+      const auto geometryJacobianI = element.geometry().jacobianInverse(qp.position());
+      for (size_t i=0; i<jacobians.size(); ++i)
+        jacobians[i] = jacobians[i] * geometryJacobianI;
+
+      FieldMatrix<ctype,targetDim,gridDim> localDeformationGradient(0);
+      for (size_t j=0; j<jacobians.size(); j++)
+        for (size_t i=0; i<targetDim; i++)
+          localDeformationGradient[i].axpy(localConfiguration[ localView.tree().child(i).localIndex(j) ], jacobians[j][0]);
+
+      trace_on(rank);
+      localDeformationGradientVec.resize(derivativeSize);
+      FieldMatrix<adouble,gridDim,gridDim> localADeformationGradient(0);
+      for (size_t i=0; i<targetDim; i++)
+        for (size_t j=0; j<gridDim; j++) {
+          // localDeformationGradientVec then contains packets of size gridDim, for j = 0, ..., targetDim
+          localDeformationGradientVec[i*gridDim + j] = localDeformationGradient[i][j];
+          localADeformationGradient[i][j] <<= localDeformationGradientVec[i*gridDim + j];
+        }
+
+      //Evaluate the density at the localDeformationGradient at the current quadrature point - to then take the derivative of that
+      adouble density = (*localDensity_)(x, localADeformationGradient);
+      ctype pureDensity;
+      density >>= pureDensity;
+
+      trace_off();
+    }
+
+};
+
+template<class LocalView>
+typename LocalADOLCIntegralStiffness<LocalView>::ctype LocalADOLCIntegralStiffness<LocalView>::
+energy(const LocalView& localView,
+       const std::vector<ctype>& localConfiguration) const
+{
+  const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+  const auto& element = localView.element();
+
+  ctype energy = 0;
+
+  using FiniteElement = typename LocalView::Tree::ChildType::FiniteElement;
+  using Jacobian = typename FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
+  std::vector<Jacobian> jacobians(localFiniteElement.size());
+
+  int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+                                               : localFiniteElement.localBasis().order() * gridDim;
+
+  const auto& quad = QuadratureRules<ctype, gridDim>::rule(element.type(), quadOrder);
+
+  for (const auto& qp : quad)
+  {
+    const auto integrationElement = element.geometry().integrationElement(qp.position());
+
+    const auto geometryJacobianI = element.geometry().jacobianInverse(qp.position());
+
+    auto x = element.geometry().global(qp.position());
+
+    localFiniteElement.localBasis().evaluateJacobian(qp.position(), jacobians);
+
+    for (size_t i=0; i<jacobians.size(); ++i)
+      jacobians[i] = jacobians[i] * geometryJacobianI;
+
+    FieldMatrix<adouble,gridDim,gridDim> deformationGradient(0);
+    for (size_t i=0; i<jacobians.size(); i++)
+      for (size_t j=0; j<gridDim; j++)
+        deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], jacobians[i][0]);
+
+    // TODO: Actually, the energy calculation does not have to use adoubles.
+    // We calculate the energy using adoubles, just because the density is of the type LocalDensity<gridDim,adouble,ctype>.
+    adouble aEnergyAtThisQP = qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
+    ctype energyAtThisQP;
+    aEnergyAtThisQP >>= energyAtThisQP;
+    energy += energyAtThisQP;
+  }
+
+  return energy;
+}
+
+template<class LocalView>
+void LocalADOLCIntegralStiffness<LocalView>::
+assembleGradient(const LocalView& localView,
+                 const std::vector<ctype>& localConfiguration,
+                 std::vector<ctype>& localGradient)
+{
+  int rank = MPIHelper::getCommunication().rank();
+
+  size_t nDoubles = localConfiguration.size();
+
+  localGradient.resize(nDoubles);
+  std::fill(localGradient.begin(), localGradient.end(), 0.0);
+
+  // powerBasis: grab the finite element of the first child
+  const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+  const auto& element = localView.element();
+
+  // store gradients of shape functions and basis functions
+  std::vector<Jacobian> jacobians(localFiniteElement.size());
+
+  int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+                                               : localFiniteElement.localBasis().order() * gridDim;
+
+  const auto& quad = QuadratureRules<ctype, gridDim>::rule(element.type(), quadOrder);
+
+  for (const auto& qp : quad)
+  {
+    std::vector<ctype> localDeformationGradientVec;
+    evaluateActiveDensity(localView,localConfiguration, qp, jacobians, localDeformationGradientVec);
+
+    const auto integrationElement = element.geometry().integrationElement(qp.position());
+
+    std::vector<ctype> gradientDensity(localDeformationGradientVec.size());
+    gradient(rank,derivativeSize,localDeformationGradientVec.data(),gradientDensity.data());
+
+    std::vector<ctype> rawGradient(nDoubles);
+    std::fill(rawGradient.begin(), rawGradient.end(), 0.0);
+
+    //Chain rule: Multiply gradientDensity with the gradient of the shape functions from the right
+    for (size_t k = 0; k < jacobians.size(); k++)
+      for (size_t i = 0; i < targetDim; i++)
+        for (size_t j = 0; j < gridDim; j++)
+          rawGradient[i*jacobians.size() + k] += gradientDensity[i*gridDim + j] * jacobians[k][0][j];
+
+    for (size_t i = 0; i < nDoubles; i++)
+      localGradient[i] += qp.weight() * integrationElement * rawGradient[i];
+  }
+
+}
+
+
+template<class LocalView>
+void LocalADOLCIntegralStiffness<LocalView>::
+assembleGradientAndHessian(const LocalView& localView,
+                           const std::vector<ctype>& localConfiguration,
+                           std::vector<ctype>& localGradient,
+                           Matrix<ctype>& localHessian)
+{
+  int rank = MPIHelper::getCommunication().rank();
+
+  size_t nDoubles = localConfiguration.size();
+
+  localGradient.resize(nDoubles);
+  localHessian.setSize(nDoubles,nDoubles);
+
+  std::fill(localGradient.begin(), localGradient.end(), 0.0);
+
+  for (size_t i = 0; i < nDoubles; i++)
+    for (size_t j = 0; j < nDoubles; j++)
+      localHessian[i][j] = 0;
+
+  // powerBasis: grab the finite element of the first child
+  const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+  const auto& element = localView.element();
+
+  // store gradients of shape functions and basis functions
+  std::vector<Jacobian> jacobians(localFiniteElement.size());
+
+  int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+                                               : localFiniteElement.localBasis().order() * gridDim;
+
+  const auto& quad = QuadratureRules<ctype, gridDim>::rule(element.type(), quadOrder);
+
+  ctype* hessianDensity[derivativeSize];
+  // we need only the lower half of the hessian, thus allocate i+1 doubles only
+  for(size_t i=0; i<derivativeSize; i++)
+    hessianDensity[i] = (ctype*)malloc((i+1)*sizeof(ctype));
+
+  for (const auto& qp : quad)
+  {
+    std::vector<ctype> localDeformationGradientVec;
+    evaluateActiveDensity(localView,localConfiguration, qp, jacobians, localDeformationGradientVec);
+
+    const auto integrationElement = element.geometry().integrationElement(qp.position());
+
+    std::vector<ctype> gradientDensity(localDeformationGradientVec.size());
+    gradient(rank,derivativeSize,localDeformationGradientVec.data(),gradientDensity.data());
+
+    // ADOL-C error code
+    int rc;
+
+    if ( vectorMode_ )
+      rc = hessian2(rank,derivativeSize,localDeformationGradientVec.data(),hessianDensity);
+    else
+      rc = hessian(rank,derivativeSize,localDeformationGradientVec.data(),hessianDensity);
+
+    if (rc < 0)
+      DUNE_THROW(Exception, "ADOL-C has returned with error code " << rc << "!");
+
+
+    std::vector<ctype> rawGradient(nDoubles);
+    std::fill(rawGradient.begin(), rawGradient.end(), 0.0);
+    
+    std::vector<ctype> rawHessian(nDoubles*nDoubles);
+    std::fill(rawHessian.begin(), rawHessian.end(), 0.0);
+    
+    //Chain rule: Multiply gradientDensity and gradientHessian with the gradient of the shape functions
+    //The gradient gets multiplied from the right, the Hessian from the left and the right
+    for (size_t k = 0; k < jacobians.size(); k++) {
+      for (size_t i = 0; i < targetDim; i++) {
+        for (size_t j = 0; j < gridDim; j++) {
+          rawGradient[i*jacobians.size() + k] += gradientDensity[i*gridDim + j] * jacobians[k][0][j];
+
+          auto index1 = i*jacobians.size() + k;
+          for (size_t kk = 0; kk < jacobians.size(); kk++) {
+            for (size_t ii = 0; ii < gridDim; ii++) {
+              for (size_t jj = 0; jj < gridDim; jj++) {
+                auto index2 = ii*jacobians.size() + kk;
+                if (ii*gridDim + jj < i*gridDim + j) { // lower half and upper half
+                  rawHessian[index1*nDoubles + index2] += hessianDensity[i*gridDim + j][ii*gridDim + jj]*jacobians[k][0][j]*jacobians[kk][0][jj];
+                  rawHessian[index2*nDoubles + index1] += hessianDensity[i*gridDim + j][ii*gridDim + jj]*jacobians[k][0][j]*jacobians[kk][0][jj];
+                }
+                if (ii*gridDim + jj == i*gridDim + j) { // diagonal
+                  rawHessian[index1*nDoubles + index2] += hessianDensity[i*gridDim + j][ii*gridDim + jj]*jacobians[k][0][j]*jacobians[kk][0][jj];
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+
+    for (size_t i = 0; i < nDoubles; i++) {
+      localGradient[i] += qp.weight() * integrationElement * rawGradient[i];
+      for (size_t j = 0; j < nDoubles; j++) {
+        localHessian[i][j] += qp.weight() * integrationElement * rawHessian[i*nDoubles + j];
+      }
+    }
+
+  }
+  for(size_t i=0; i<derivativeSize; i++)
+    free(hessianDensity[i]);
+}
+
+
+} // namespace Dune::Elasticity
+
+#endif

dune/elasticity/assemblers/localadolcstiffness.hh

@@ -41,6 +41,14 @@ public:
     virtual double energy (const LocalView& localView,
                            const std::vector<double>& localConfiguration) const;
 
+    /** \brief Assemble the local gradient at the current position
+    *
+    *     This uses the automatic differentiation toolbox ADOL_C.
+    */
+    virtual void assembleGradient(const LocalView& localView,
+                                  const std::vector<double>& localConfiguration,
+                                  std::vector<double>& localGradient);
+
     /** \brief Assemble the local stiffness matrix at the current position
      *
      *     This uses the automatic differentiation toolbox ADOL_C.
@@ -77,7 +85,7 @@ energy(const LocalView& localView,
         energy = localEnergy_->energy(localView,localAConfiguration);
     } catch (Dune::Exception &e) {
         trace_off();
-        throw e;
+        throw;
     }
 
     energy >>= pureEnergy;
@@ -88,6 +96,31 @@ energy(const LocalView& localView,
 }
 
 
+
+template<class LocalView>
+void LocalADOLCStiffness<LocalView>::
+assembleGradient(const LocalView& localView,
+                 const std::vector<double>& localConfiguration,
+                 std::vector<double>& localGradient
+                )
+{
+  int rank = Dune::MPIHelper::getCommunication().rank();
+  // Tape energy computation.  We may not have to do this every time, but it's comparatively cheap.
+  energy(localView, localConfiguration);
+
+  /////////////////////////////////////////////////////////////////
+  // Compute the energy gradient
+  /////////////////////////////////////////////////////////////////
+
+  // Compute the actual gradient
+  size_t nDoubles = localConfiguration.size();
+
+  // Compute gradient
+  localGradient.resize(nDoubles);
+  gradient(rank,nDoubles,localConfiguration.data(),localGradient.data());
+}
+
+
 template<class LocalView>
 void LocalADOLCStiffness<LocalView>::
 assembleGradientAndHessian(const LocalView& localView,
@@ -97,24 +130,19 @@ assembleGradientAndHessian(const LocalView& localView,
                           )
 {
     int rank = Dune::MPIHelper::getCommunication().rank();
-    // Tape energy computation.  We may not have to do this every time, but it's comparatively cheap.
-    energy(localView, localConfiguration);
 
     /////////////////////////////////////////////////////////////////
     // Compute the energy gradient
     /////////////////////////////////////////////////////////////////
 
-    // Compute the actual gradient
-    size_t nDoubles = localConfiguration.size();
-
-    // Compute gradient
-    localGradient.resize(nDoubles);
-    gradient(rank,nDoubles,localConfiguration.data(),localGradient.data());
+    // this tapes the energy computation
+    assembleGradient(localView, localConfiguration, localGradient);
 
     /////////////////////////////////////////////////////////////////
     // Compute Hessian
     /////////////////////////////////////////////////////////////////
 
+    size_t nDoubles = localConfiguration.size();
     localHessian.setSize(nDoubles,nDoubles);
 
     // allocate raw doubles: ADOL-C requires "**double" arguments
@@ -149,194 +177,11 @@ assembleGradientAndHessian(const LocalView& localView,
       localHessian[i][i] = rawHessian[i][i];
     }
 
-    // don't forget the clean everything up
+    // don't forget to clean up everything
     for(size_t i=0; i<nDoubles; i++)
       free(rawHessian[i]);
 }
 
 } // namespace Dune::Elasticity
 
-
-/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
- */
-template<class GridView, class LocalFiniteElement, class VectorType>
-class
-LocalADOLCStiffness
-    : public LocalFEStiffness<GridView,LocalFiniteElement,VectorType>
-{
-    // grid types
-    typedef typename GridView::Grid::ctype DT;
-    typedef typename VectorType::value_type::field_type RT;
-    typedef typename GridView::template Codim<0>::Entity Entity;
-
-    // some other sizes
-    enum {gridDim=GridView::dimension};
-
-    // Hack
-    typedef std::vector<Dune::FieldVector<adouble,gridDim> > AVectorType;
-
-public:
-
-    //! Dimension of a tangent space
-    enum { blocksize = VectorType::value_type::dimension };
-
-    LocalADOLCStiffness(const Dune::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* energy)
-    : localEnergy_(energy)
-    {} [[deprecated("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::LocalADOLCStiffness with LocalView concept!")]]
-
-    /** \brief Compute the energy at the current configuration */
-    virtual RT energy (const Entity& e,
-               const LocalFiniteElement& localFiniteElement,
-               const VectorType& localConfiguration) const;
-
-    /** \brief Assemble the local stiffness matrix at the current position
-
-    This uses the automatic differentiation toolbox ADOL_C.
-    */
-    virtual void assembleGradientAndHessian(const Entity& e,
-                         const LocalFiniteElement& localFiniteElement,
-                         const VectorType& localConfiguration,
-                         VectorType& localGradient);
-
-    const Dune::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* localEnergy_;
-
-};
-
-
-template <class GridView, class LocalFiniteElement, class VectorType>
-typename LocalADOLCStiffness<GridView, LocalFiniteElement, VectorType>::RT
-LocalADOLCStiffness<GridView, LocalFiniteElement, VectorType>::
-energy(const Entity& element,
-       const LocalFiniteElement& localFiniteElement,
-       const VectorType& localSolution) const
-{
-    int rank = Dune::MPIHelper::getCommunication().rank();
-    double pureEnergy;
-
-    std::vector<Dune::FieldVector<adouble,blocksize> > localASolution(localSolution.size());
-
-    trace_on(rank);
-
-    adouble energy = 0;
-
-    for (size_t i=0; i<localSolution.size(); i++)
-      for (size_t j=0; j<localSolution[i].size(); j++)
-        localASolution[i][j] <<= localSolution[i][j];
-
-    energy = localEnergy_->energy(element,localFiniteElement,localASolution);
-
-    energy >>= pureEnergy;
-
-    trace_off();
-
-    return pureEnergy;
-}
-
-
-
-// ///////////////////////////////////////////////////////////
-//   Compute gradient and Hessian together
-//   To compute the Hessian we need to compute the gradient anyway, so we may
-//   as well return it.  This saves assembly time.
-// ///////////////////////////////////////////////////////////
-template <class GridType, class LocalFiniteElement, class VectorType>
-void LocalADOLCStiffness<GridType, LocalFiniteElement, VectorType>::
-assembleGradientAndHessian(const Entity& element,
-                const LocalFiniteElement& localFiniteElement,
-                const VectorType& localSolution,
-                VectorType& localGradient)
-{
-    int rank = Dune::MPIHelper::getCommunication().rank();
-    // Tape energy computation.  We may not have to do this every time, but it's comparatively cheap.
-    energy(element, localFiniteElement, localSolution);
-
-    /////////////////////////////////////////////////////////////////
-    // Compute the energy gradient
-    /////////////////////////////////////////////////////////////////
-
-    // Compute the actual gradient
-    size_t nDofs = localSolution.size();
-    size_t nDoubles = nDofs*blocksize;
-    std::vector<double> xp(nDoubles);
-    int idx=0;
-    for (size_t i=0; i<nDofs; i++)
-        for (size_t j=0; j<blocksize; j++)
-            xp[idx++] = localSolution[i][j];
-
-  // Compute gradient
-    std::vector<double> g(nDoubles);
-    gradient(rank,nDoubles,xp.data(),g.data());                  // gradient evaluation
-
-    // Copy into Dune type
-    localGradient.resize(localSolution.size());
-
-    idx=0;
-    for (size_t i=0; i<nDofs; i++)
-        for (size_t j=0; j<blocksize; j++)
-            localGradient[i][j] = g[idx++];
-
-    /////////////////////////////////////////////////////////////////
-    // Compute Hessian
-    /////////////////////////////////////////////////////////////////
-
-    this->A_.setSize(nDofs,nDofs);
-
-#ifndef ADOLC_VECTOR_MODE
-    std::vector<double> v(nDoubles);
-    std::vector<double> w(nDoubles);
-
-    std::fill(v.begin(), v.end(), 0.0);
-
-    for (size_t i=0; i<nDofs; i++)
-      for (size_t ii=0; ii<blocksize; ii++)
-      {
-        // Evaluate Hessian in the direction of each vector of the orthonormal frame
-        for (size_t k=0; k<blocksize; k++)
-          v[i*blocksize + k] = (k==ii);
-
-        int rc= 3;
-        MINDEC(rc, hess_vec(rank, nDoubles, xp.data(), v.data(), w.data()));
-        if (rc < 0)
-          DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
-
-        for (size_t j=0; j<nDoubles; j++)
-          this->A_[i][j/blocksize][ii][j%blocksize] = w[j];
-
-        // Make v the null vector again
-        std::fill(&v[i*blocksize], &v[(i+1)*blocksize], 0.0);
-      }
-#else
-    int n = nDoubles;
-    int nDirections = nDofs * blocksize;
-    double* tangent[nDoubles];
-    for(size_t i=0; i<nDoubles; i++)
-        tangent[i] = (double*)malloc(nDirections*sizeof(double));
-
-    double* rawHessian[nDoubles];
-    for(size_t i=0; i<nDoubles; i++)
-        rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
-
-    for (int j=0; j<nDirections; j++)
-    {
-      for (int i=0; i<n; i++)
-        tangent[i][j] = 0.0;
-
-      for (int i=0; i<embeddedBlocksize; i++)
-        tangent[(j/blocksize)*embeddedBlocksize+i][j] = orthonormalFrame[j/blocksize][j%blocksize][i];
-    }
-
-    hess_mat(rank,nDoubles,nDirections,xp.data(),tangent,rawHessian);
-
-    // Copy Hessian into Dune data type
-    for(size_t i=0; i<nDoubles; i++)
-      for (size_t j=0; j<nDirections; j++)
-        this->A_[j/blocksize][i/embeddedBlocksize][j%blocksize][i%embeddedBlocksize] = rawHessian[i][j];
-
-    for(size_t i=0; i<nDoubles; i++) {
-        free(rawHessian[i]);
-        free(tangent[i]);
-    }
-#endif
-}
-
 #endif

dune/elasticity/assemblers/localenergy.hh

@@ -24,34 +24,5 @@ public:
 
 } // namespace Dune::Elasticity
 
-namespace Dune {
-
-// WARNING: This interface is deprecated and will be removed!
-
-/** \brief Base class for energies defined by integrating over one grid element */
-template<class GridView, class LocalFiniteElement, class VectorType>
-class
-LocalEnergy
-{
-    typedef typename VectorType::value_type::field_type RT;
-    typedef typename GridView::template Codim<0>::Entity Element;
-
-public:
-
-  /** \brief Compute the energy
-   *
-   * \param element A grid element
-   * \param LocalFiniteElement A finite element on the reference element
-   * \param localConfiguration The coefficients of a FE function on the current element
-   */
-  /** \brief Compute the energy at the current configuration */
-  virtual RT energy (const Element& element,
-                     const LocalFiniteElement& localFiniteElement,
-                     const VectorType& localConfiguration) const = 0;
-
-};
-
-}  // namespace Dune
-
 #endif  // DUNE_ELASTICITY_ASSEMBLERS_LOCALENERGY_HH
 

dune/elasticity/assemblers/localfestiffness.hh

@@ -16,6 +16,13 @@ class LocalFEStiffness
 
 public:
 
+    /** \brief Assemble the local gradient at the current position
+    */
+    virtual void assembleGradient(const LocalView& localView,
+                                  const std::vector<RT>& localConfiguration,
+                                  std::vector<RT>& localGradient
+                                 ) = 0;
+
     /** \brief Assemble the local gradient and tangent matrix at the current position
     */
     virtual void assembleGradientAndHessian(const LocalView& localView,
@@ -28,37 +35,5 @@ public:
 
 } // namespace Dune::Elasticity
 
-
-template<class GridView, class LocalFiniteElement, class VectorType>
-class
-LocalFEStiffness
-: public Dune::LocalEnergy<GridView, LocalFiniteElement, VectorType>
-{
-    // grid types
-    typedef typename GridView::Grid::ctype DT;
-    typedef typename VectorType::value_type::field_type RT;
-    typedef typename GridView::template Codim<0>::Entity Entity;
-
-    // some other sizes
-    enum {gridDim=GridView::dimension};
-
-public:
-
-    //! Dimension of a tangent space
-    enum { blocksize = VectorType::value_type::dimension };
-
-    /** \brief Assemble the local gradient and tangent matrix at the current position
-    */
-    virtual void assembleGradientAndHessian(const Entity& e,
-                                 const LocalFiniteElement& localFiniteElement,
-                                 const VectorType& localConfiguration,
-                                 VectorType& localGradient) = 0;
-
-    // assembled data
-    [[deprecated("Use dune-functions powerBases with LocalView concept. See Dune::Elasticity::LocalFEStiffness")]]
-    Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> > A_;
-
-};
-
 #endif
 

dune/elasticity/assemblers/localfunctionalassembler.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef DUNE_ELASTICITY_ASSEMBLERS_LOCALFUNCTIONALASSEMBLER_HH
+#define DUNE_ELASTICITY_ASSEMBLERS_LOCALFUNCTIONALASSEMBLER_HH
+
+#warning localfunctionalassembler.hh is deprecated!
+
+#include <dune/istl/bvector.hh>
+
+/** \brief Abstract base class for local operator assemblers
+ *
+ *  \tparam GridType The grid we are assembling on
+ *  \tparam TestLocalFE the local finite element of the test space
+ *  \tparam T the block type
+ */
+template <class GridType, class TestLocalFE, typename T>
+class LocalFunctionalAssembler
+{
+    public:
+        using Element = typename GridType::template Codim<0>::Entity;
+        using LocalVector = Dune::BlockVector<T>;
+
+        /** \brief Assemble a local problem
+         *
+         *  \param element the grid element on which to operate
+         *  \param localVector will contain the assembled element vector
+         *  \param tFE the local finite element in the test space used on element
+         */
+        virtual void assemble(const Element& element, LocalVector& localVector,
+                              const TestLocalFE& tFE) const  = 0;
+
+        template<class TestLocalView>
+        void operator()(
+                const Element& element,
+                LocalVector& localVector,
+                const TestLocalView& testLocalView) const
+        {
+            assemble(element, localVector, testLocalView.tree().finiteElement());
+        }
+};
+
+#endif
+

dune/elasticity/assemblers/localhyperdualstiffness.hh

@@ -8,7 +8,7 @@
 
 namespace Dune::Elasticity {
 
-/** \brief Assembles energy gradient and Hessian using Hyperdual Numbers 
+/** \brief Assembles energy gradient and Hessian using Hyperdual Numbers
  */
 template<class LocalView>
 class LocalHyperDualStiffness
@@ -27,6 +27,15 @@ public:
     virtual double energy (const LocalView& localView,
                            const std::vector<double>& localConfiguration) const;
 
+
+    /** \brief Assemble the local gradient at the current position
+     *
+     *     This uses the automatic differentiation using hyperdual numbers
+     */
+    virtual void assembleGradient(const LocalView& localView,
+                                  const std::vector<double>& localConfiguration,
+                                  std::vector<double>& localGradient);
+
     /** \brief Assemble the local stiffness matrix at the current position
      *
      *     This uses the automatic differentiation using hyperdual numbers
@@ -36,7 +45,7 @@ public:
                                             std::vector<double>& localGradient,
                                             Dune::Matrix<double>& localHessian);
 
-    std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, hyperdual>> localEnergy_; 
+    std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, hyperdual>> localEnergy_;
 
 };
 
@@ -50,10 +59,10 @@ energy(const LocalView& localView,
 
     std::vector<hyperdual> localHyperDualConfiguration(localConfiguration.size());
     hyperdual energy;
-    
+
     for (size_t i=0; i<localConfiguration.size(); i++)
-        localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);        
-        
+        localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
+
     energy = localEnergy_->energy(localView,localHyperDualConfiguration);
 
     pureEnergy = energy.real();
@@ -62,6 +71,37 @@ energy(const LocalView& localView,
 }
 
 
+
+template<class LocalView>
+void LocalHyperDualStiffness<LocalView>::
+assembleGradient(const LocalView& localView,
+                 const std::vector<double>& localConfiguration,
+                 std::vector<double>& localGradient
+)
+{
+    size_t nDoubles = localConfiguration.size();
+
+    localGradient.resize(nDoubles);
+
+    // Create hyperdual vector localHyperDualConfiguration = double vector localConfiguration
+    std::vector<hyperdual> localHyperDualConfiguration(nDoubles);
+    for (size_t i=0; i<nDoubles; i++)
+        localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
+
+    // Compute gradient and hessian (symmetry is used) using hyperdual function evaluation
+    // localHyperDualConfiguration puts Ones in the eps1 and eps2 component to evaluate the different partial derivatives
+    for (size_t i=0; i<nDoubles; i++)
+    {
+        localHyperDualConfiguration[i].seteps1(1.0);
+
+        hyperdual temp = localEnergy_->energy(localView, localHyperDualConfiguration);
+        localGradient[i] = temp.eps1();
+
+        localHyperDualConfiguration[i].seteps1(0.0);
+    }
+}
+
+
 template<class LocalView>
 void LocalHyperDualStiffness<LocalView>::
 assembleGradientAndHessian(const LocalView& localView,
@@ -70,7 +110,7 @@ assembleGradientAndHessian(const LocalView& localView,
                            Dune::Matrix<double>& localHessian
                           )
 {
-    size_t nDoubles = localConfiguration.size(); 
+    size_t nDoubles = localConfiguration.size();
 
     localGradient.resize(nDoubles);
 
@@ -79,28 +119,28 @@ assembleGradientAndHessian(const LocalView& localView,
     std::vector<hyperdual> localHyperDualConfiguration(nDoubles);
     for (size_t i=0; i<nDoubles; i++)
        localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
-    
+
     // Compute gradient and hessian (symmetry is used) using hyperdual function evaluation
     // localHyperDualConfiguration puts Ones in the eps1 and eps2 component to evaluate the different partial derivatives
-    for (size_t i=0; i<nDoubles; i++) 
+    for (size_t i=0; i<nDoubles; i++)
     {
       localHyperDualConfiguration[i].seteps1(1.0);
       localHyperDualConfiguration[i].seteps2(1.0);
-    
+
       hyperdual temp = localEnergy_->energy(localView, localHyperDualConfiguration);
-      localGradient[i] = temp.eps1();                   
+      localGradient[i] = temp.eps1();
       localHessian[i][i] = temp.eps1eps2();
-      
+
       localHyperDualConfiguration[i].seteps2(0.0);
 
-      for (size_t j=i+1; j<nDoubles; j++) 
-      { 
+      for (size_t j=i+1; j<nDoubles; j++)
+      {
         localHyperDualConfiguration[j].seteps2(1.0);
         localHessian[i][j] = localEnergy_->energy(localView, localHyperDualConfiguration).eps1eps2();
-        localHessian[j][i] = localHessian[i][j];        // Use symmetry of hessian 
+        localHessian[j][i] = localHessian[i][j];        // Use symmetry of hessian
         localHyperDualConfiguration[j].seteps2(0.0);
       }
-        
+
       localHyperDualConfiguration[i].seteps1(0.0);
     }
 

dune/elasticity/assemblers/localoperatorassembler.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef LOCAL_OPERATOR_ASSEMBLER_HH
+#define LOCAL_OPERATOR_ASSEMBLER_HH
+
+#warning localoperatorassembler.hh is deprecated!
+
+#include <type_traits>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/matrix.hh>
+
+
+/** \brief Abstract base class for local operator assemblers
+ *
+ *  \tparam GridType The grid we are assembling on
+ *  \tparam TestLocalFE the local finite element of the test space
+ *  \tparam AnsatzLocalFE the local finite element of the ansatz space
+ *  \tparam T the block type
+ */
+template <class GridType, class TestLocalFE, class AnsatzLocalFE, typename T>
+class LocalOperatorAssembler
+{
+
+    public:
+
+        typedef typename GridType::template Codim<0>::Entity Element;
+        typedef Dune::Matrix< Dune::FieldMatrix<int,1,1> > BoolMatrix;
+        typedef typename Dune::Matrix<T> LocalMatrix;
+        typedef typename T::field_type field_type;
+        typedef T MatrixEntry;
+
+        enum {ndofs = T::cols};
+
+        virtual ~LocalOperatorAssembler() {}
+
+        /** \brief After return BoolMatrix isNonZero contains the bit pattern of coupling basis functions
+         *
+         *  \param element the grid element on which to operate
+         *  \param isNonZero will contain bit pattern of coupling basis functions after return
+         *  \param tFE the local finite element in the test space used on element
+         *  \param aFE the local finite element in the ansatz space used on element
+         */
+        virtual void indices(const Element& element, BoolMatrix& isNonZero,
+                const TestLocalFE& tFE,
+                const AnsatzLocalFE& aFE) const = 0;
+
+        /** \brief Assemble a local problem providing a local solution.
+         *
+         *  Since this is a linear assembler the solution is simply discarded.
+         */
+        virtual void assemble(const Element& element,
+                const Dune::BlockVector<Dune::FieldVector<field_type,ndofs> >& localSolution,
+                LocalMatrix& localMatrix,
+                const TestLocalFE& tFE,
+                const AnsatzLocalFE& aFE) const
+        {
+            assemble(element, localMatrix, tFE, aFE);
+        }
+
+        /** \brief Assemble a local problem
+         *
+         *  \param element the grid element on which to operate
+         *  \param localMatrix will contain the assembled element matrix
+         *  \param tFE the local finite element in the test space used on element
+         *  \param aFE the local finite element in the ansatz space used on element
+         */
+        virtual void assemble(const Element& element, LocalMatrix& localMatrix,
+                const TestLocalFE& tFE,
+                const AnsatzLocalFE& aFE) const = 0;
+
+
+
+        template<class TestLocalView, class AnsatzLocalView>
+        void operator()(
+                const Element& element,
+                LocalMatrix& localMatrix,
+                const TestLocalView& testLocalView,
+                const AnsatzLocalView& ansatzLocalView) const
+        {
+            assemble(element, localMatrix, testLocalView.tree().finiteElement(), ansatzLocalView.tree().finiteElement());
+        }
+
+
+
+        /** \brief Check if test and ansatz local finite element are the same
+         *
+         *  \param tFE a local finite element in the test space
+         *  \param aFE a local finite element in the ansatz space
+         */
+        static bool isSameFE(const TestLocalFE& tFE, const AnsatzLocalFE& aFE)
+        {
+            if (not std::is_same<TestLocalFE, AnsatzLocalFE>::value)
+                return false;
+            return &tFE == reinterpret_cast<const TestLocalFE*>(&aFE);
+        }
+
+};
+
+
+#endif
+

dune/elasticity/assemblers/localstiffnesssum.hh

+#ifndef DUNE_ELASTICITY_ASSEMBLERS_LOCALSTIFFNESSSUM_HH
+#define DUNE_ELASTICITY_ASSEMBLERS_LOCALSTIFFNESSSUM_HH
+
+#include <vector>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/matrix.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+
+namespace Dune::Elasticity {
+
+/** \brief Sums the assembled energy gradient and Hessian of the given stiffness classes
+ */
+template<class LocalView>
+class LocalStiffnessSum
+    : public LocalFEStiffness<LocalView,double>
+{
+  public:
+
+    /** \brief Constructor with a vector of local energies (compatibility for old interface)*/
+    LocalStiffnessSum(std::vector<std::shared_ptr<LocalFEStiffness<LocalView, double>>> stiffnessVector)
+    {
+      for (auto localStiffness : stiffnessVector)
+        addLocalStiffness(localStiffness);
+    }
+
+    /** \brief Default Constructor*/
+    LocalStiffnessSum()
+    {}
+
+    virtual ~LocalStiffnessSum() {}
+
+    void addLocalStiffness(std::shared_ptr<LocalFEStiffness<LocalView, double>> localStiffess) {
+      localStiffnesses_.push_back(localStiffess);
+    }
+
+    /** \brief Compute the energy at the current configuration */
+    virtual double energy (const LocalView& localView,
+                           const std::vector<double>& localConfiguration) const;
+
+    /** \brief Assemble the local gradient at the current position
+    */
+    virtual void assembleGradient(const LocalView& localView,
+                                  const std::vector<double>& localConfiguration,
+                                  std::vector<double>& localGradient);
+
+    /** \brief Assemble the local stiffness matrix at the current position
+     */
+    virtual void assembleGradientAndHessian(const LocalView& localView,
+                                            const std::vector<double>& localConfiguration,
+                                            std::vector<double>& localGradient,
+                                            Matrix<double>& localHessian);
+
+  private:
+      std::vector<std::shared_ptr<LocalFEStiffness<LocalView, double>>> localStiffnesses_;
+
+};
+
+template<class LocalView>
+double LocalStiffnessSum<LocalView>::
+energy(const LocalView& localView,
+       const std::vector<double>& localConfiguration) const
+{
+  double energy = 0;
+  for (const auto& localStiffness: localStiffnesses_ )
+    energy += localStiffness->energy(localView,localConfiguration);
+  return energy;
+}
+
+template<class LocalView>
+void LocalStiffnessSum<LocalView>::
+assembleGradient(const LocalView& localView,
+                 const std::vector<double>& localConfiguration,
+                 std::vector<double>& localGradient
+                )
+{
+  DUNE_THROW(Exception, "Error: Use assembleGradientAndHessian instead of assembleGradient!");
+}
+
+
+template<class LocalView>
+void LocalStiffnessSum<LocalView>::
+assembleGradientAndHessian(const LocalView& localView,
+                           const std::vector<double>& localConfiguration,
+                           std::vector<double>& localGradient,
+                           Matrix<double>& localHessian
+                          )
+{
+  size_t nDoubles = localConfiguration.size();
+
+  localGradient.resize(nDoubles);
+  localHessian.setSize(nDoubles,nDoubles);
+
+  for (size_t i = 0; i < nDoubles; i++) {
+    localGradient[i] = 0;
+    for (size_t j = 0; j < nDoubles; j++) {
+      localHessian[i][j] = 0;
+    }
+  }
+  for (const auto& localStiffness: localStiffnesses_ ){
+    std::vector<double> innerLocalGradient;
+    Matrix<double> innerLocalHessian;
+    localStiffness->assembleGradientAndHessian(localView, localConfiguration, innerLocalGradient, innerLocalHessian);
+    for (size_t i = 0; i < nDoubles; i++) {
+      localGradient[i] += innerLocalGradient[i];
+      for (size_t j = 0; j < nDoubles; j++) {
+        localHessian[i][j] += innerLocalHessian[i][j];
+      }
+    }
+  }
+}
+
+} // namespace Dune::Elasticity
+
+#endif