.gitignore

@@ -4,3 +4,6 @@
 # Default cmake build directory
 /build-cmake
 
+# ignore Python files
+*.pyc
+

.gitlab-ci.yml

@@ -6,6 +6,24 @@ before_script:
   - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/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
+  script: duneci-standard-test
+  variables:
+    DUNECI_BRANCH: releases/2.8
+
+dune:2.8 gcc C++20:
+  image: registry.dune-project.org/docker/ci/dune:2.8-debian-11-gcc-9-20
+  script: duneci-standard-test
+  variables:
+    DUNECI_BRANCH: releases/2.8
+
+dune:2.8  clang  C++17:
+  image: registry.dune-project.org/docker/ci/dune:2.8-debian-10-clang-7-libcpp-17
+  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
   script: duneci-standard-test

CHANGELOG.md

-# Master (will become release 2.8)
+# Master (will become 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
+
+# 2.8 Release
 
 - Introduce class `LocalDensity` for material-specific implementations
 - Introduce class `LocalIntegralEnergy` to work with the densities
 - Local energies and FE assemblers use now `dune-functions` power bases instead of scalar `dune-fufem` bases; the key element is the `LocalView` which contains the information for each element
 - Introduce class `LocalHyperDualStiffness` and `hyperdual` to calculate gradient and hessian using hyper-dual numbers
+- The class `FEAssembler` now stores a `const` reference to the basis, rather than a value.
+  This allows to use the assembler with bases that are not copyable
+  (as, e.g., the `RaviartThomasBasis`, see [dune-functions issue 58](https://gitlab.dune-project.org/staging/dune-functions/-/issues/58)).
 
 ## Deprecations and removals
 

CMakeLists.txt

@@ -27,7 +27,7 @@ add_subdirectory("dune")
 add_subdirectory("src")
 add_subdirectory("cmake/modules")
 
-if(HAVE_AMIRAMESH AND HAVE_IPOPT AND HAVE_PSURFACE AND HAVE_UG)
+if(HAVE_AMIRAMESH AND HAVE_IPOPT AND HAVE_PSURFACE AND dune-uggrid_FOUND)
   set(programs viscoelast nonlinelast)
 
   foreach(_program ${programs})
@@ -40,7 +40,7 @@ if(HAVE_AMIRAMESH AND HAVE_IPOPT AND HAVE_PSURFACE AND HAVE_UG)
   endforeach()
 endif()
 
-if(HAVE_AMIRAMESH AND HAVE_IPOPT AND HAVE_PSURFACE AND HAVE_UG AND ADOLC_FOUND)
+if(HAVE_AMIRAMESH AND HAVE_IPOPT AND HAVE_PSURFACE AND dune-uggrid_FOUND AND ADOLC_FOUND)
   add_dune_adolc_flags(nonlinelast)
 endif()
 

dune.module

@@ -4,8 +4,8 @@
 
 #Name of the module
 Module:dune-elasticity
-Version: 2.8-git
-Maintainer: oliver.sander@tu-dresden.de, youett@mi.fu-berlin.de
+Version: 2.9-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
 Suggests: dune-uggrid dune-parmg

dune/elasticity/assemblers/feassembler.hh

@@ -13,6 +13,8 @@
 #include <dune/functions/functionspacebases/concepts.hh>
 #include <dune/functions/backends/istlvectorbackend.hh>
 
+#include <dune/grid/common/partitionset.hh>
+
 #include "localfestiffness.hh"
 
 
@@ -31,7 +33,7 @@ class FEAssembler {
 
 public:
 
-    const Basis basis_;
+    const Basis& basis_;
 
     /** \brief Partition type on which to assemble
      *
@@ -52,6 +54,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
@@ -69,6 +76,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>::
@@ -117,18 +165,18 @@ assembleGradientAndHessian(const VectorType& configuration,
     Matrix<field_type> localHessian;
     localHessian.setSize(size, size);
 
-    for (int i=0; i<size; i++)
+    for (size_t i=0; i<size; i++)
       localConfiguration[i] = configurationBackend[ localView.index(i) ];
 
     // setup local matrix and gradient
     localStiffness_->assembleGradientAndHessian(localView, localConfiguration, localGradient, localHessian);
 
-    for(int i=0; i<size; i++)
+    for (size_t i=0; i<size; i++)
     {
       auto row = localView.index(i);
       gradientBackend[row] += localGradient[i];
 
-      for (int j=0; j<size; j++ )
+      for (size_t j=0; j<size; j++ )
       {
         auto col = localView.index(j);
         // fufem ISTLBackend uses operator() for entry access
@@ -194,8 +242,8 @@ FEAssembler
     typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
 
 public:
-    const Basis
-      DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::FEAssembler with dune-functions bases.") basis_;
+    [[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
      *

dune/elasticity/assemblers/localadolcstiffness.hh

@@ -35,10 +35,20 @@ public:
     , useHessian2_(useHessian2)
     {}
 
+    virtual ~LocalADOLCStiffness() {}
+
     /** \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
+    *
+    *     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.
@@ -59,7 +69,7 @@ double LocalADOLCStiffness<LocalView>::
 energy(const LocalView& localView,
        const std::vector<double>& localConfiguration) const
 {
-    int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
+    int rank = Dune::MPIHelper::getCommunication().rank();
     double pureEnergy;
 
     std::vector<adouble> localAConfiguration(localConfiguration.size());
@@ -74,18 +84,43 @@ energy(const LocalView& localView,
     try {
         energy = localEnergy_->energy(localView,localAConfiguration);
     } catch (Dune::Exception &e) {
-        trace_off(rank);
-        throw e;
+        trace_off();
+        throw;
     }
 
     energy >>= pureEnergy;
 
-    trace_off(rank);
+    trace_off();
 
     return pureEnergy;
 }
 
 
+
+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,
@@ -94,25 +129,20 @@ assembleGradientAndHessian(const LocalView& localView,
                            Dune::Matrix<double>& localHessian
                           )
 {
-    int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
-    // Tape energy computation.  We may not have to do this every time, but it's comparatively cheap.
-    energy(localView, localConfiguration);
+    int rank = Dune::MPIHelper::getCommunication().rank();
 
     /////////////////////////////////////////////////////////////////
     // 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
@@ -180,7 +210,7 @@ public:
 
     LocalADOLCStiffness(const Dune::LocalEnergy<GridView, LocalFiniteElement, AVectorType>* energy)
     : localEnergy_(energy)
-    {} DUNE_DEPRECATED_MSG("dune-elasticity with dune-fufem bases is now deprecated. Use Elasticity::LocalADOLCStiffness with LocalView concept!")
+    {} [[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,
@@ -208,7 +238,7 @@ energy(const Entity& element,
        const LocalFiniteElement& localFiniteElement,
        const VectorType& localSolution) const
 {
-    int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
+    int rank = Dune::MPIHelper::getCommunication().rank();
     double pureEnergy;
 
     std::vector<Dune::FieldVector<adouble,blocksize> > localASolution(localSolution.size());
@@ -225,7 +255,7 @@ energy(const Entity& element,
 
     energy >>= pureEnergy;
 
-    trace_off(rank);
+    trace_off();
 
     return pureEnergy;
 }
@@ -244,7 +274,7 @@ assembleGradientAndHessian(const Entity& element,
                 const VectorType& localSolution,
                 VectorType& localGradient)
 {
-    int rank = Dune::MPIHelper::getCollectiveCommunication().rank();
+    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);
 

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,
@@ -55,8 +62,8 @@ public:
                                  VectorType& localGradient) = 0;
 
     // assembled data
-    Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> >
-      DUNE_DEPRECATED_MSG("Use dune-functions powerBases with LocalView concept. See Dune::Elasticity::LocalFEStiffness") A_;
+    [[deprecated("Use dune-functions powerBases with LocalView concept. See Dune::Elasticity::LocalFEStiffness")]]
+    Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> > A_;
 
 };
 

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/common/hyperdual.hh

@@ -58,63 +58,63 @@ public:
 
 	//basic manipulation
 	hyperdual operator+ () const;
-	hyperdual operator+ (const hyperdual rhs) const;
-	friend hyperdual operator+ (const double lhs, const hyperdual rhs);
+	hyperdual operator+ (const hyperdual& rhs) const;
+	friend hyperdual operator+ (const double lhs, const hyperdual& rhs);
 	hyperdual operator- () const;
-	hyperdual operator- (const hyperdual rhs) const;
-	friend hyperdual operator- (const double lhs, const hyperdual rhs);
-	hyperdual operator* (const hyperdual rhs)const;
-	friend hyperdual operator* (const double lhs, const hyperdual rhs);
-	friend hyperdual operator/ (const hyperdual lhs, const hyperdual rhs);
-	friend hyperdual operator/ (const double lhs, const hyperdual rhs);
-	friend hyperdual operator/ (const hyperdual lhs, const double rhs);
-	hyperdual& operator+= (hyperdual rhs);
-	hyperdual& operator-= (hyperdual rhs);
-	hyperdual& operator*= (hyperdual rhs);
+	hyperdual operator- (const hyperdual& rhs) const;
+	friend hyperdual operator- (const double lhs, const hyperdual& rhs);
+	hyperdual operator* (const hyperdual& rhs)const;
+	friend hyperdual operator* (const double lhs, const hyperdual& rhs);
+	friend hyperdual operator/ (const hyperdual& lhs, const hyperdual& rhs);
+	friend hyperdual operator/ (const double lhs, const hyperdual& rhs);
+	friend hyperdual operator/ (const hyperdual& lhs, const double rhs);
+	hyperdual& operator+= (const hyperdual& rhs);
+	hyperdual& operator-= (const hyperdual& rhs);
+	hyperdual& operator*= (const hyperdual& rhs);
 	hyperdual& operator*= (double rhs);
-	hyperdual& operator/= (hyperdual rhs);
+	hyperdual& operator/= (const hyperdual& rhs);
 	hyperdual& operator/= (double rhs);
 
 	//math.h functions
-	friend hyperdual pow (hyperdual x, double a);
-	friend hyperdual pow (hyperdual x, hyperdual a);
-	friend hyperdual exp(hyperdual x);
-	friend hyperdual log(hyperdual x);
-	friend hyperdual sin(hyperdual x);
-	friend hyperdual cos(hyperdual x);
-	friend hyperdual tan(hyperdual x);
-	friend hyperdual asin(hyperdual x);
-	friend hyperdual acos(hyperdual x);
-	friend hyperdual atan(hyperdual x);
-	friend hyperdual sqrt(hyperdual x);
-	friend hyperdual fabs(hyperdual x);
-	friend hyperdual abs(hyperdual x);  
-	friend hyperdual max(hyperdual x1, hyperdual x2);
-	friend hyperdual max(hyperdual x1, double x2);
-	friend hyperdual max(double x1, hyperdual x2);
-	friend hyperdual min(hyperdual x1, hyperdual x2);
-	friend hyperdual min(hyperdual x1, double x2);
-	friend hyperdual min(double x1, hyperdual x2);
+	friend hyperdual pow (const hyperdual& x, double a);
+	friend hyperdual pow (const hyperdual& x, const hyperdual& a);
+	friend hyperdual exp(const hyperdual& x);
+	friend hyperdual log(const hyperdual& x);
+	friend hyperdual sin(const hyperdual& x);
+	friend hyperdual cos(const hyperdual& x);
+	friend hyperdual tan(const hyperdual& x);
+	friend hyperdual asin(const hyperdual& x);
+	friend hyperdual acos(const hyperdual& x);
+	friend hyperdual atan(const hyperdual& x);
+	friend hyperdual sqrt(const hyperdual& x);
+	friend hyperdual fabs(const hyperdual& x);
+	friend hyperdual abs(const hyperdual& x);
+	friend hyperdual max(const hyperdual& x1, const hyperdual& x2);
+	friend hyperdual max(const hyperdual& x1, double x2);
+	friend hyperdual max(double x1, const hyperdual& x2);
+	friend hyperdual min(const hyperdual& x1, const hyperdual& x2);
+	friend hyperdual min(const hyperdual& x1, double x2);
+	friend hyperdual min(double x1, const hyperdual& x2);
 
 	//comparisons
-	friend bool operator> (hyperdual lhs, hyperdual rhs);
-	friend bool operator> (double lhs, hyperdual rhs);
-	friend bool operator> (hyperdual lhs, double rhs);
-	friend bool operator>= (hyperdual lhs, hyperdual rhs);
-	friend bool operator>= (double lhs, hyperdual rhs);
-	friend bool operator>= (hyperdual lhs, double rhs);
-	friend bool operator< (hyperdual lhs, hyperdual rhs);
-	friend bool operator< (double lhs, hyperdual rhs);
-	friend bool operator< (hyperdual lhs, double rhs);
-	friend bool operator<= (hyperdual lhs, hyperdual rhs);
-	friend bool operator<= (double lhs, hyperdual rhs);
-	friend bool operator<= (hyperdual lhs, double rhs);
-	friend bool operator== (hyperdual lhs, hyperdual rhs);
-	friend bool operator== (double lhs, hyperdual rhs);
-	friend bool operator== (hyperdual lhs, double rhs);
-	friend bool operator!= (hyperdual lhs, hyperdual rhs);
-	friend bool operator!= (double lhs, hyperdual rhs);
-	friend bool operator!= (hyperdual lhs, double rhs);
+	friend bool operator> (const hyperdual& lhs, const hyperdual& rhs);
+	friend bool operator> (double lhs, const hyperdual& rhs);
+	friend bool operator> (const hyperdual& lhs, double rhs);
+	friend bool operator>= (const hyperdual& lhs, const hyperdual& rhs);
+	friend bool operator>= (double lhs, const hyperdual& rhs);
+	friend bool operator>= (const hyperdual& lhs, double rhs);
+	friend bool operator< (const hyperdual& lhs, const hyperdual& rhs);
+	friend bool operator< (double lhs, const hyperdual& rhs);
+	friend bool operator< (const hyperdual& lhs, double rhs);
+	friend bool operator<= (const hyperdual& lhs, const hyperdual& rhs);
+	friend bool operator<= (double lhs, const hyperdual& rhs);
+	friend bool operator<= (const hyperdual& lhs, double rhs);
+	friend bool operator== (const hyperdual& lhs, const hyperdual& rhs);
+	friend bool operator== (double lhs, const hyperdual& rhs);
+	friend bool operator== (const hyperdual& lhs, double rhs);
+	friend bool operator!= (const hyperdual& lhs, const hyperdual& rhs);
+	friend bool operator!= (double lhs, const hyperdual& rhs);
+	friend bool operator!= (const hyperdual& lhs, double rhs);
 };
 
 
@@ -192,7 +192,7 @@ hyperdual hyperdual::operator+ () const
 {
 	return *this;
 }
-hyperdual hyperdual::operator+ (const hyperdual rhs) const
+hyperdual hyperdual::operator+ (const hyperdual& rhs) const
 {
 	hyperdual temp;
 	temp.f0 = f0 + rhs.f0;
@@ -201,7 +201,7 @@ hyperdual hyperdual::operator+ (const hyperdual rhs) const
 	temp.f12 = f12 + rhs.f12;
 	return temp;
 }
-hyperdual operator+ (const double lhs, const hyperdual rhs)
+hyperdual operator+ (const double lhs, const hyperdual& rhs)
 {
 	hyperdual temp;
 	temp.f0 = lhs + rhs.f0;
@@ -219,7 +219,7 @@ hyperdual hyperdual::operator- () const
 	temp.f12 = -f12;
 	return temp;
 }
-hyperdual hyperdual::operator- (const hyperdual rhs) const
+hyperdual hyperdual::operator- (const hyperdual& rhs) const
 {
 	hyperdual temp;
 	temp.f0 = f0 - rhs.f0;
@@ -228,7 +228,7 @@ hyperdual hyperdual::operator- (const hyperdual rhs) const
 	temp.f12 = f12 - rhs.f12;
 	return temp;
 }
-hyperdual operator- (const double lhs, const hyperdual rhs)
+hyperdual operator- (const double lhs, const hyperdual& rhs)
 {
 	hyperdual temp;
 	temp.f0 = lhs - rhs.f0;
@@ -237,7 +237,7 @@ hyperdual operator- (const double lhs, const hyperdual rhs)
 	temp.f12 = -rhs.f12;
 	return temp;
 }
-hyperdual hyperdual::operator* (const hyperdual rhs) const
+hyperdual hyperdual::operator* (const hyperdual& rhs) const
 {
 	hyperdual temp;
 	temp.f0 = f0*rhs.f0;
@@ -246,7 +246,7 @@ hyperdual hyperdual::operator* (const hyperdual rhs) const
 	temp.f12 = f0*rhs.f12 + f1*rhs.f2 + f2*rhs.f1 + f12*rhs.f0;
 	return temp;
 }
-hyperdual operator* (const double lhs, const hyperdual rhs)
+hyperdual operator* (const double lhs, const hyperdual& rhs)
 {
 	hyperdual temp;
 	temp.f0 = lhs*rhs.f0;
@@ -255,7 +255,7 @@ hyperdual operator* (const double lhs, const hyperdual rhs)
 	temp.f12 = lhs*rhs.f12;
 	return temp;
 }
-hyperdual operator/ (const hyperdual lhs, const hyperdual rhs)
+hyperdual operator/ (const hyperdual& lhs, const hyperdual& rhs)
 {
 	hyperdual inv;
 	double prod;
@@ -267,13 +267,13 @@ hyperdual operator/ (const hyperdual lhs, const hyperdual rhs)
 	inv.f12 = 2.0*rhs.f1*rhs.f2/(prod*rhs.f0) - rhs.f12/prod;
 	return lhs*inv;
 }
-hyperdual operator/ (const double lhs, const hyperdual rhs)
+hyperdual operator/ (const double lhs, const hyperdual& rhs)
 {
 	hyperdual temp,inv;
 	temp = hyperdual(lhs)/rhs;
 	return temp;
 }
-hyperdual operator/ (const hyperdual lhs, const double rhs)
+hyperdual operator/ (const hyperdual& lhs, const double rhs)
 {
 	hyperdual temp;
 	double inv;
@@ -284,7 +284,7 @@ hyperdual operator/ (const hyperdual lhs, const double rhs)
 	temp.f12 = inv*lhs.f12;
 	return temp;
 }
-hyperdual& hyperdual::operator+= (hyperdual rhs)
+hyperdual& hyperdual::operator+= (const hyperdual& rhs)
 {
 	f0 += rhs.f0;
 	f1 += rhs.f1;
@@ -292,7 +292,7 @@ hyperdual& hyperdual::operator+= (hyperdual rhs)
 	f12 += rhs.f12;
 	return *this;
 }
-hyperdual& hyperdual::operator-= (hyperdual rhs)
+hyperdual& hyperdual::operator-= (const hyperdual& rhs)
 {
 	f0 -= rhs.f0;
 	f1 -= rhs.f1;
@@ -300,7 +300,7 @@ hyperdual& hyperdual::operator-= (hyperdual rhs)
 	f12 -= rhs.f12;
 	return *this;
 }
-hyperdual& hyperdual::operator*= (hyperdual rhs)
+hyperdual& hyperdual::operator*= (const hyperdual& rhs)
 {
 	double tf0,tf1,tf2,tf12;
 	tf0=f0;
@@ -321,7 +321,7 @@ hyperdual& hyperdual::operator*= (double rhs)
 	f12 *= rhs;
 	return *this;
 }
-hyperdual& hyperdual::operator/= (hyperdual rhs) 
+hyperdual& hyperdual::operator/= (const hyperdual& rhs)
 {
     *this = *this / rhs;
     return *this;
@@ -334,7 +334,7 @@ hyperdual& hyperdual::operator/= (double rhs)
 	f12 /= rhs;
 	return *this;
 }
-hyperdual pow (hyperdual x, double a)
+hyperdual pow (const hyperdual& x, double a)
 {
 	hyperdual temp;
 	double deriv,xval,tol;
@@ -356,11 +356,11 @@ hyperdual pow (hyperdual x, double a)
 
 	return temp;
 }
-hyperdual pow (hyperdual x, hyperdual a)
+hyperdual pow (const hyperdual& x, const hyperdual& a)
 {
 	return exp(a*log(x));
 }
-hyperdual exp(hyperdual x)
+hyperdual exp(const hyperdual& x)
 {
 	hyperdual temp;
 	double deriv;
@@ -371,7 +371,7 @@ hyperdual exp(hyperdual x)
 	temp.f12 = deriv*(x.f12 + x.f1*x.f2);
 	return temp;
 }
-hyperdual log(hyperdual x)
+hyperdual log(const hyperdual& x)
 {
 	hyperdual temp;
 	double deriv1,deriv2;
@@ -383,7 +383,7 @@ hyperdual log(hyperdual x)
 	temp.f12 = x.f12/x.f0 - (deriv1*deriv2);
 	return temp;
 }
-hyperdual sin(hyperdual x)
+hyperdual sin(const hyperdual& x)
 {
 	hyperdual temp;
 	double funval,deriv;
@@ -395,7 +395,7 @@ hyperdual sin(hyperdual x)
 	temp.f12 = deriv*x.f12 - funval*x.f1*x.f2;
 	return temp;
 }
-hyperdual cos(hyperdual x)
+hyperdual cos(const hyperdual& x)
 {
 	hyperdual temp;
 	double funval,deriv;
@@ -407,7 +407,7 @@ hyperdual cos(hyperdual x)
 	temp.f12 = deriv*x.f12 - funval*x.f1*x.f2;
 	return temp;
 }
-hyperdual tan(hyperdual x)
+hyperdual tan(const hyperdual& x)
 {
 	hyperdual temp;
 	double funval,deriv;
@@ -419,7 +419,7 @@ hyperdual tan(hyperdual x)
 	temp.f12 = deriv*x.f12 + x.f1*x.f2*(2*funval*deriv);
 	return temp;
 }
-hyperdual asin(hyperdual x)
+hyperdual asin(const hyperdual& x)
 {
 	hyperdual temp;
 	double funval,deriv1,deriv;
@@ -432,7 +432,7 @@ hyperdual asin(hyperdual x)
 	temp.f12 = deriv*x.f12 + x.f1*x.f2*(x.f0*std::pow(deriv1,-1.5));
 	return temp;
 }
-hyperdual acos(hyperdual x)
+hyperdual acos(const hyperdual& x)
 {
 	hyperdual temp;
 	double funval,deriv1,deriv;
@@ -445,7 +445,7 @@ hyperdual acos(hyperdual x)
 	temp.f12 = deriv*x.f12 + x.f1*x.f2*(-x.f0*std::pow(deriv1,-1.5));
 	return temp;
 }
-hyperdual atan(hyperdual x)
+hyperdual atan(const hyperdual& x)
 {
 	hyperdual temp;
 	double funval,deriv1,deriv;
@@ -458,7 +458,7 @@ hyperdual atan(hyperdual x)
 	temp.f12 = deriv*x.f12 + x.f1*x.f2*(-2*x.f0/(deriv1*deriv1));
 	return temp;
 }
-hyperdual sqrt(hyperdual x)
+hyperdual sqrt(const hyperdual& x)
 {
 	hyperdual temp;
 	temp.f0 = std::sqrt(x.f0);
@@ -467,7 +467,7 @@ hyperdual sqrt(hyperdual x)
 	temp.f12 = 0.5*(x.f12 - 2.*temp.f1*temp.f2)/temp.f0;
 	return temp;
 }
-hyperdual fabs(hyperdual x)
+hyperdual fabs(const hyperdual& x)
 {
 	hyperdual temp;
 	if (x < 0.0)
@@ -477,11 +477,11 @@ hyperdual fabs(hyperdual x)
 	return temp;
 }
 
-hyperdual abs(hyperdual x)
+hyperdual abs(const hyperdual& x)
 {
 	return fabs(x);
 }
-hyperdual max(hyperdual x1, hyperdual x2)
+hyperdual max(const hyperdual& x1, const hyperdual& x2)
 {
 	hyperdual temp;
 	if (x1>x2)
@@ -490,7 +490,7 @@ hyperdual max(hyperdual x1, hyperdual x2)
 		temp = x2;
 	return temp;
 }
-hyperdual max(hyperdual x1, double x2)
+hyperdual max(const hyperdual& x1, double x2)
 {
 	hyperdual temp;
 	if (x1>x2)
@@ -499,7 +499,7 @@ hyperdual max(hyperdual x1, double x2)
 		temp = x2;
 	return temp;
 }
-hyperdual max(double x1, hyperdual x2)
+hyperdual max(double x1, const hyperdual& x2)
 {
 	hyperdual temp;
 	if (x1>x2)
@@ -508,7 +508,7 @@ hyperdual max(double x1, hyperdual x2)
 		temp = x2;
 	return temp;
 }
-hyperdual min(hyperdual x1, hyperdual x2)
+hyperdual min(const hyperdual& x1, const hyperdual& x2)
 {
 	hyperdual temp;
 	if (x1<x2)
@@ -517,7 +517,7 @@ hyperdual min(hyperdual x1, hyperdual x2)
 		temp = x2;
 	return temp;
 }
-hyperdual min(hyperdual x1, double x2)
+hyperdual min(const hyperdual& x1, double x2)
 {
 	hyperdual temp;
 	if (x1<x2)
@@ -526,7 +526,7 @@ hyperdual min(hyperdual x1, double x2)
 		temp = x2;
 	return temp;
 }
-hyperdual min(double x1, hyperdual x2)
+hyperdual min(double x1, const hyperdual& x2)
 {
 	hyperdual temp;
 	if (x1<x2)
@@ -536,79 +536,103 @@ hyperdual min(double x1, hyperdual x2)
 	return temp;
 }
 
-bool operator> (hyperdual lhs, hyperdual rhs)
+bool operator> (const hyperdual& lhs, const hyperdual& rhs)
 {
 	return (lhs.f0 > rhs.f0);
 }
-bool operator> (double lhs, hyperdual rhs)
+bool operator> (double lhs, const hyperdual& rhs)
 {
 	return (lhs > rhs.f0);
 }
-bool operator> (hyperdual lhs, double rhs)
+bool operator> (const hyperdual& lhs, double rhs)
 {
 	return (lhs.f0 > rhs);
 }
-bool operator>= (hyperdual lhs, hyperdual rhs)
+bool operator>= (const hyperdual& lhs, const hyperdual& rhs)
 {
 	return (lhs.f0 >= rhs.f0);
 }
-bool operator>= (double lhs, hyperdual rhs)
+bool operator>= (double lhs, const hyperdual& rhs)
 {
 	return (lhs >= rhs.f0);
 }
-bool operator>= (hyperdual lhs, double rhs)
+bool operator>= (const hyperdual& lhs, double rhs)
 {
 	return (lhs.f0 >= rhs);
 }
-bool operator< (hyperdual lhs, hyperdual rhs)
+bool operator< (const hyperdual& lhs, const hyperdual& rhs)
 {
 	return (lhs.f0 < rhs.f0);
 }
-bool operator< (double lhs, hyperdual rhs)
+bool operator< (double lhs, const hyperdual& rhs)
 {
 	return (lhs < rhs.f0);
 }
-bool operator< (hyperdual lhs, double rhs)
+bool operator< (const hyperdual& lhs, double rhs)
 {
 	return (lhs.f0 < rhs);
 }
-bool operator<= (hyperdual lhs, hyperdual rhs)
+bool operator<= (const hyperdual& lhs, const hyperdual& rhs)
 {
 	return (lhs.f0 <= rhs.f0);
 }
-bool operator<= (double lhs, hyperdual rhs)
+bool operator<= (double lhs, const hyperdual& rhs)
 {
 	return (lhs <= rhs.f0);
 }
-bool operator<= (hyperdual lhs, double rhs)
+bool operator<= (const hyperdual& lhs, double rhs)
 {
 	return (lhs.f0 <= rhs);
 }
-bool operator== (hyperdual lhs, hyperdual rhs)
+bool operator== (const hyperdual& lhs, const hyperdual& rhs)
 {
 	return (lhs.f0 == rhs.f0);
 }
-bool operator== (double lhs, hyperdual rhs)
+bool operator== (double lhs, const hyperdual& rhs)
 {
 	return (lhs == rhs.f0);
 }
-bool operator== (hyperdual lhs, double rhs)
+bool operator== (const hyperdual& lhs, double rhs)
 {
 	return (lhs.f0 == rhs);
 }
-bool operator!= (hyperdual lhs, hyperdual rhs)
+bool operator!= (const hyperdual& lhs, const hyperdual& rhs)
 {
 	return (lhs.f0 != rhs.f0);
 }
-bool operator!= (double lhs, hyperdual rhs)
+bool operator!= (double lhs, const hyperdual& rhs)
 {
 	return (lhs != rhs.f0);
 }
-bool operator!= (hyperdual lhs, double rhs)
+bool operator!= (const hyperdual& lhs, double rhs)
 {
 	return (lhs.f0 != rhs);
 }
 
+bool isnormal( const hyperdual& x )
+{
+  using std::isnormal;
+  return isnormal(x.real());
+}
+
+bool isfinite( const hyperdual& x )
+{
+  using std::isfinite;
+  return isfinite(x.real());
+}
+
+bool isnan( const hyperdual& x )
+{
+  using std::isnan;
+  return isnan(x.real());
+}
+
+bool isinf( const hyperdual& x )
+{
+  using std::isinf;
+  return isinf(x.real());
+}
+
 
 namespace std
 {

dune/elasticity/common/trustregionsolver.cc

@@ -62,7 +62,6 @@ setup(const typename BasisType::GridView::Grid& grid,
     baseTolerance_            = baseTolerance;
     damping_                  = damping;
 
-    int numLevels = grid_->maxLevel()+1;
     const auto dim = VectorType::value_type::dimension;
 
 #if HAVE_DUNE_PARMG
@@ -90,7 +89,7 @@ setup(const typename BasisType::GridView::Grid& grid,
     std::conditional_t< // do we have a dune-functions basis?
       Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
         BasisType,
-        DuneFunctionsBasis<BasisType> > basis(grid.leafGridView());
+        DuneFunctionsBasis<BasisType> > basis(assembler_->basis_);
     // ////////////////////////////////
     //   Create a multigrid solver
     // ////////////////////////////////
@@ -289,10 +288,14 @@ setup(const typename BasisType::GridView::Grid& grid,
       isP1Basis = std::is_same<Basis,Dune::Functions::LagrangeBasis<typename Basis::GridView, 1> >::value;
     }
 
+    int numLevels = grid_->maxLevel()+1;
+
     using TransferOperatorType = typename TruncatedCompressedMGTransfer<CorrectionType>::TransferOperatorType;
     std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<CorrectionType>>> transferOperators(isP1Basis ? numLevels-1 : numLevels);
 
     // Here we set up the restriction of the leaf grid space into the leaf grid P1/Q1 space
+    // TODO: IIRC that move to P1Basis was only implemented because I didn't have general
+    // prolongation operators between arbitratry FE spaces.
     if (not isP1Basis)
     {
       TransferOperatorType pkToP1TransferMatrix;
@@ -325,13 +328,8 @@ setup(const typename BasisType::GridView::Grid& grid,
     for (int i=0; i<numLevels-1; i++) {
 
         // Construct the local multigrid transfer matrix
-        TruncatedCompressedMGTransfer<CorrectionType>* newTransferOp = new TruncatedCompressedMGTransfer<CorrectionType>;
-        newTransferOp->setup(*grid_,i,i+1);
-
-        auto op = std::make_shared<TruncatedCompressedMGTransfer<CorrectionType>>();
-        transferOperators[i] = op;
-        std::shared_ptr<TransferOperatorType> transferOperatorMatrix = std::make_shared<TransferOperatorType>(newTransferOp->getMatrix());
-        op->setMatrix(transferOperatorMatrix);
+        transferOperators[i] = std::make_shared<TruncatedCompressedMGTransfer<CorrectionType> >();
+        transferOperators[i]->setup(*grid_,i,i+1);
     }
 
     mmgStep->setTransferOperators(transferOperators);
@@ -361,14 +359,11 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
         mgStep = dynamic_cast<MonotoneMGStep<MatrixType,CorrectionType>*>(&iterationStep);
 
     }
-
-    BasisType basis(grid_->leafGridView());
 #else
-    BasisType basis(grid_->levelGridView(grid_->maxLevel()));
     BasisType coarseBasis(grid_->levelGridView(0));
     std::vector<BoxConstraint<typename VectorType::field_type, blocksize>> coarseTrustRegionObstacles(coarseBasis.size());
 #endif
-    MaxNormTrustRegion<blocksize> trustRegion(basis.size(), initialTrustRegionRadius_);
+    MaxNormTrustRegion<blocksize> trustRegion(assembler_->basis_.size(), initialTrustRegionRadius_);
 
     std::vector<BoxConstraint<field_type,blocksize> > trustRegionObstacles;
 
@@ -556,16 +551,12 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
 
             if (correctionInfinityNorm < this->tolerance_) {
                 if (this->verbosity_ == NumProc::FULL and rank==0)
-                    std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
-
-                if (this->verbosity_ != NumProc::QUIET and rank==0)
-                    std::cout << i+1 << " trust-region steps were taken." << std::endl;
-                break;
-            }
-
-            if (trustRegion.radius() < this->tolerance_) {
-              if (this->verbosity_ == NumProc::FULL and rank==0)
-                    std::cout << "TRUST REGION RADIUS IS TOO SMALL, SMALLER THAN TOLERANCE, STOPPING NOW" << std::endl;
+                {
+                    if (correctionInfinityNorm < trustRegion.radius())
+                        std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
+                    else
+                        std::cout << "TRUST-REGION UNDERFLOW!" << std::endl;
+                }
 
                 if (this->verbosity_ != NumProc::QUIET and rank==0)
                     std::cout << i+1 << " trust-region steps were taken." << std::endl;

dune/elasticity/common/trustregionsolver.hh

@@ -26,6 +26,11 @@
 #include <dune/elasticity/assemblers/feassembler.hh>
 
 #if HAVE_DUNE_PARMG
+
+// enable a temporary hack to let parmg work with powerBases
+// which is needed for the current implementation of elasticty
+#define DUNE_PARMG_POWER_BASIS_HACK
+
 #include <dune/parmg/iterationstep/lambdastep.hh>
 #include <dune/parmg/iterationstep/multigrid.hh>
 #include <dune/parmg/parallel/dofmap.hh>
@@ -68,14 +73,7 @@ public:
     TrustRegionSolver()
         : IterativeSolver<VectorType, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL),
           hessianMatrix_(std::shared_ptr<MatrixType>(NULL)), h1SemiNorm_(NULL)
-    {
-#if HAVE_DUNE_PARMG
-      // TODO currently, we cannot use PARMG together with dune-functions bases
-      static_assert( not Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() ,
-        "Currently, dune-parmg and dune-functions bases cannot be used together."
-      );
-#endif
-    }
+    {}
 
     /** \brief Set up the solver using a monotone multigrid method as the inner solver */
     void setup(const GridType& grid,

dune/elasticity/materials/exphenckydensity.hh

@@ -85,7 +85,7 @@ public:
   }
 
   /** \brief Lame constants */
-  field_type mu_, lambda_, kappa_, k_, khat_, alpha_, beta_;
+  double mu_, lambda_, kappa_, k_, khat_, alpha_, beta_;
 };
 
 } // namespace Dune::Elasticity

dune/elasticity/materials/exphenckyenergy.hh

@@ -7,7 +7,7 @@
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::ExpHenckyDensity")
+class [[deprecated("Use Elasticity::LocalIntegralEnergy with Elasticity::ExpHenckyDensity")]]
 ExpHenckyEnergy
   : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
 {

dune/elasticity/materials/henckyenergy.hh

@@ -7,7 +7,7 @@
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::HenckyDensity")
+class [[deprecated("Use Elasticity::LocalIntegralEnergy with Elasticity::HenckyDensity")]]
 HenckyEnergy
   : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
 {

dune/elasticity/materials/localintegralenergy.hh

@@ -2,6 +2,10 @@
 #define DUNE_ELASTICITY_MATERIALS_LOCALINTEGRALENERGY_HH
 
 #include <dune/common/fmatrix.hh>
+#include <dune/common/version.hh>
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 8)
+#include <dune/common/transpose.hh>
+#endif
 
 #include <dune/geometry/quadraturerules.hh>
 
@@ -53,8 +57,9 @@ energy(const LocalView& localView,
   field_type energy = 0;
 
   // store gradients of shape functions and base functions
-  std::vector<FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
-  std::vector<FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
+  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;
@@ -65,23 +70,30 @@ energy(const LocalView& localView,
   {
     const DT integrationElement = element.geometry().integrationElement(qp.position());
 
-    const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(qp.position());
+    const auto geometryJacobianIT = element.geometry().jacobianInverseTransposed(qp.position());
 
     // Global position
     auto x = element.geometry().global(qp.position());
 
     // Get gradients of shape functions
-    localFiniteElement.localBasis().evaluateJacobian(qp.position(), referenceGradients);
+    localFiniteElement.localBasis().evaluateJacobian(qp.position(), jacobians);
 
     // compute gradients of base functions
-    for (size_t i=0; i<gradients.size(); ++i)
-      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+    for (size_t i=0; i<jacobians.size(); ++i)
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 8)
+      jacobians[i] = jacobians[i] * transpose(geometryJacobianIT);
+#else
+    {
+      auto referenceJacobian = jacobians[i];
+      geometryJacobianIT.mv(referenceJacobian[0], jacobians[i][0]);
+    }
+#endif
 
     // Deformation gradient
     FieldMatrix<field_type,gridDim,gridDim> deformationGradient(0);
-    for (size_t i=0; i<gradients.size(); i++)
+    for (size_t i=0; i<jacobians.size(); i++)
       for (int j=0; j<gridDim; j++)
-        deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], gradients[i]);
+        deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], jacobians[i][0]);
     // Integrate the energy density
     energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
   }
@@ -124,8 +136,8 @@ public:
                      const std::vector<Dune::FieldVector<field_type,gridDim> >& localConfiguration) const;
 
 protected:
-  const std::shared_ptr<Elasticity::LocalDensity<gridDim,field_type,DT>>
-    DUNE_DEPRECATED_MSG("Use dune-functions powerBases with LocalView concept. See Elasticity::LocalIntegralEnergy") localDensity_ = nullptr;
+  [[deprecated("Use dune-functions powerBases with LocalView concept. See Elasticity::LocalIntegralEnergy")]]
+  const std::shared_ptr<Elasticity::LocalDensity<gridDim,field_type,DT>> localDensity_ = nullptr;
 
 };
 

dune/elasticity/materials/mooneyrivlindensity.hh

@@ -2,7 +2,6 @@
 #define DUNE_ELASTICITY_MATERIALS_MOONEYRIVLINDENSITY_HH
 
 #include <dune/common/fmatrix.hh>
-#include <dune/common/fmatrixev.hh>
 
 #include <dune/elasticity/materials/localdensity.hh>
 
@@ -25,7 +24,9 @@ public:
       mooneyrivlin_a = parameters.get<double>("mooneyrivlin_a");
       mooneyrivlin_b = parameters.get<double>("mooneyrivlin_b");
       mooneyrivlin_c = parameters.get<double>("mooneyrivlin_c");
-    } else if (mooneyrivlin_energy == "log" or mooneyrivlin_energy == "square") {
+    }
+    else if (mooneyrivlin_energy == "log" or mooneyrivlin_energy == "square")
+    {
       mooneyrivlin_10 = parameters.get<double>("mooneyrivlin_10");
       mooneyrivlin_01 = parameters.get<double>("mooneyrivlin_01");
       mooneyrivlin_20 = parameters.get<double>("mooneyrivlin_20");
@@ -35,8 +36,10 @@ public:
       mooneyrivlin_03 = parameters.get<double>("mooneyrivlin_03");
       mooneyrivlin_21 = parameters.get<double>("mooneyrivlin_21");
       mooneyrivlin_12 = parameters.get<double>("mooneyrivlin_12");
-      mooneyrivlin_k = parameters.get<double>("mooneyrivlin_k");
-    } else {
+      mooneyrivlin_k  = parameters.get<double>("mooneyrivlin_k");
+    }
+    else
+    {
       DUNE_THROW(Exception, "Error: Selected mooneyrivlin implementation (" << mooneyrivlin_energy << ") not available!");
     }
   }
@@ -58,78 +61,96 @@ public:
         for (int k=0; k<dim; k++)
           C[i][j] += gradient[k][i] * gradient[k][j];
 
-    //////////////////////////////////////////////////////////
-    //  Eigenvalues of C = F^TF
-    //////////////////////////////////////////////////////////
-
-
-    // eigenvalues of C, i.e., squared singular values \sigma_i of F
-    Dune::FieldVector<field_type, dim> sigmaSquared;
-    FMatrixHelp::eigenValues(C, sigmaSquared);
-
-    field_type normFSquared = gradient.frobenius_norm2();
+    /* The Mooney-Rivlin-Density is given as a function of the eigenvalues of the right Cauchy-Green-Deformation tensor C = F^TF
+       or the right Cauchy-Green-Deformation tensor B = FF^T.
+       C = F^TF and B = FF^T have the same eigenvalues - we can use either one of them.
+
+       There are three Mooney-Rivlin-Variants:
+       ciarlet: W(F) = mooneyrivlin_a*(normF)^2 + mooneyrivlin_b*(normFinv)^2*detF^2 + mooneyrivlin_c*(detF)^2 -
+                       (2*mooneyrivlin_a + 4*mooneyrivlin_b + 2*mooneyrivlin_c)*ln(detF),
+                       where the last term is chosen s.t. W( t*I ) is minimal for t=1
+       log:     W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * ln(det(F))^2
+       square:  W(F) = \sum_{i,j=0}^{i+j<=n} mooneyrivlin_ij * (I1 - 3)^i * (I2 - 3)^j + mooneyrivlin_k * 0.5 * (det(F) - 1)^2
+
+       For log and square: I1 and I2 are the first two invariants of C (or B), multiplied with a factor depending on det(F):
+       I1 = (det(F)^(-2/dim)) * [ first invariant of C ]
+          = (det(F)^(-2/dim)) * (sum of all eigenvalues of C)
+          = (det(F)^(-2/dim)) * trace(C)
+          = (det(F)^(-2/dim)) * trace(F^TF)
+          = (det(F)^(-2/dim)) * (frobenius_norm(F))^2
+       I2 = (det(F)^(-4/dim)) * [ second invariant of C ]
+          = (det(F)^(-4/dim)) * 0.5 * [ trace(C)^2 - tr(C^2) ]
+          = (det(F)^(-4/dim)) * [C_11C_22 + C_11C_33 + C_22C_33 - C_12C_21 - C_13C_31 - C_23C_32]
+          = (det(F)^(-4/dim)) * [C_11C_22 + C_11C_33 + C_22C_33 - C_12^2 - C_13^2 - C_23^2]
+    */
+
+    field_type frobeniusNormFsquared = gradient.frobenius_norm2();
     field_type detF = gradient.determinant();
 
     if (detF < 0)
-      DUNE_THROW( Dune::Exception, "det(F) < 0, so it is not possible to calculate the MooneyRivlinEnergy.");
-
-    field_type normFinvSquared = 0;
-
-    field_type c2Tilde = 0;
-    for (int i = 0; i < dim; i++) {
-      normFinvSquared += 1/sigmaSquared[i];
-      // compute D, which is the sum of the squared eigenvalues
-      for (int j = i+1; j < dim; j++)
-        c2Tilde += sigmaSquared[j]*sigmaSquared[i];
-      }
-
-    field_type trCTildeMinus3 = normFSquared/pow(detF, 2.0/dim) - 3;
-    // \tilde{D} = \frac{1}{\det{F}^{4/3}}D -> divide by det(F)^(4/3)= det(C)^(2/3)
-    c2Tilde /= pow(detF, 4.0/dim);
-    field_type c2TildeMinus3 = c2Tilde - 3;
+      DUNE_THROW( Dune::MathError, "det(F) < 0, so it is not possible to calculate the MooneyRivlinEnergy.");
 
     // Add the local energy density
     field_type strainEnergy = 0;
 
     if (mooneyrivlin_energy == "ciarlet")
     {
+      //To calculate the squared frobeniusnorm of F^(-1), we actually invert F^(-1)
+      auto gradientInverse = gradient;
+      gradientInverse.invert();
+      field_type frobeinusNormFInverseSquared = gradientInverse.frobenius_norm2();
       using std::log;
-      return mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*log(detF);
-    } else {
-      strainEnergy = mooneyrivlin_10 * trCTildeMinus3 +
-                        mooneyrivlin_01 * c2TildeMinus3  +
-                        mooneyrivlin_20 * trCTildeMinus3 * trCTildeMinus3 +
-                        mooneyrivlin_02 * c2TildeMinus3  * c2TildeMinus3  +
-                        mooneyrivlin_11 * trCTildeMinus3 * c2TildeMinus3  +
-                        mooneyrivlin_30 * trCTildeMinus3 * trCTildeMinus3 * trCTildeMinus3 +
-                        mooneyrivlin_21 * trCTildeMinus3 * trCTildeMinus3 * c2TildeMinus3  +
-                        mooneyrivlin_12 * trCTildeMinus3 * c2TildeMinus3  * c2TildeMinus3  +
-                        mooneyrivlin_03 * c2TildeMinus3  * c2TildeMinus3  * c2TildeMinus3;
-      if (mooneyrivlin_energy == "log") {
+      return mooneyrivlin_a*frobeniusNormFsquared + mooneyrivlin_b*frobeinusNormFInverseSquared*detF*detF + mooneyrivlin_c*detF*detF
+             - (2.0*mooneyrivlin_a + 4.0*mooneyrivlin_b + 2.0*mooneyrivlin_c)*log(detF);
+    }
+    else
+    {
+      // mooneyrivlin_energy is "log" or "square"
+      field_type a = pow(detF, 2.0/dim);
+      field_type invariant1Minus3 = frobeniusNormFsquared/a - 3;
+      field_type secondInvariantOfC = 0;
+      for (int i=0; i<dim-1; i++)
+        for (int j=i+1; j<dim; j++)
+          secondInvariantOfC += C[i][i]*C[j][j] - C[i][j]*C[i][j];
+      field_type invariant2Minus3 = secondInvariantOfC/(a*a) - 3;
+      strainEnergy = mooneyrivlin_10 * invariant1Minus3 +
+                     mooneyrivlin_01 * invariant2Minus3 +
+                     mooneyrivlin_20 * invariant1Minus3 * invariant1Minus3 +
+                     mooneyrivlin_02 * invariant2Minus3 * invariant2Minus3 +
+                     mooneyrivlin_11 * invariant1Minus3 * invariant2Minus3 +
+                     mooneyrivlin_30 * invariant1Minus3 * invariant1Minus3 * invariant1Minus3 +
+                     mooneyrivlin_21 * invariant1Minus3 * invariant1Minus3 * invariant2Minus3 +
+                     mooneyrivlin_12 * invariant1Minus3 * invariant2Minus3 * invariant2Minus3 +
+                     mooneyrivlin_03 * invariant2Minus3 * invariant2Minus3 * invariant2Minus3;
+      if (mooneyrivlin_energy == "log")
+      {
         using std::log;
         field_type logDetF = log(detF);
         return strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF;
-      } else { //mooneyrivlin_energy is "square"
+      }
+      else
+      {
+        //mooneyrivlin_energy is "square"
         field_type detFMinus1 = detF - 1;
-        return strainEnergy + mooneyrivlin_k* detFMinus1 * detFMinus1;
+        return strainEnergy + 0.5 * mooneyrivlin_k* detFMinus1 * detFMinus1;
       }
     }
   }
 
   /** \brief Lame constants */
-  field_type mooneyrivlin_a,
-             mooneyrivlin_b,
-             mooneyrivlin_c,
-             mooneyrivlin_10,
-             mooneyrivlin_01,
-             mooneyrivlin_20,
-             mooneyrivlin_02,
-             mooneyrivlin_11,
-             mooneyrivlin_30,
-             mooneyrivlin_21,
-             mooneyrivlin_12,
-             mooneyrivlin_03,
-             mooneyrivlin_k;
+  double mooneyrivlin_a,
+         mooneyrivlin_b,
+         mooneyrivlin_c,
+         mooneyrivlin_10,
+         mooneyrivlin_01,
+         mooneyrivlin_20,
+         mooneyrivlin_02,
+         mooneyrivlin_11,
+         mooneyrivlin_30,
+         mooneyrivlin_21,
+         mooneyrivlin_12,
+         mooneyrivlin_03,
+         mooneyrivlin_k;
   std::string mooneyrivlin_energy;
 };
 

dune/elasticity/materials/mooneyrivlinenergy.hh

@@ -7,7 +7,7 @@
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::MooneyRivlinDensity")
+class [[deprecated("Use Elasticity::LocalIntegralEnergy with Elasticity::MooneyRivlinDensity")]]
 MooneyRivlinEnergy
   : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
 {

dune/elasticity/materials/neohookedensity.hh

@@ -76,7 +76,7 @@ public:
   }
 
   /** \brief Lame constants */
-  field_type mu_, lambda_, kappa_;
+  double mu_, lambda_, kappa_;
 };
 
 } // namespace Dune::Elasticity

dune/elasticity/materials/neohookeenergy.hh

@@ -7,7 +7,7 @@
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::NeoHookeDensity")
+class [[deprecated("Use Elasticity::LocalIntegralEnergy with Elasticity::NeoHookeDensity")]]
 NeoHookeEnergy
   : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
 {