Move the density of the harmonic energy into a separate class

The rest of the HarmonicEnergy class is now independent of
the density.  It is functionally equivalent to the LocalIntegralEnergy
class, and shall be replaced by that eventually.  However,
LocalIntegralEnergy still has the deformation/rotation-split
hard-coded in a few places, and that needs more work.

dune/gfe/assemblers/harmonicenergy.hh

@@ -2,9 +2,11 @@
 #define DUNE_GFE_HARMONICENERGY_HH
 
 #include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
 #include <dune/geometry/quadraturerules.hh>
 
 #include <dune/gfe/assemblers/localenergy.hh>
+#include <dune/gfe/densities/harmonicdensity.hh>
 
 template<class Basis, class LocalInterpolationRule, class TargetSpace>
 class HarmonicEnergy
@@ -32,6 +34,8 @@ HarmonicEnergy<Basis, LocalInterpolationRule, TargetSpace>::
 energy(const typename Basis::LocalView& localView,
        const std::vector<TargetSpace>& localSolution) const
 {
+  Dune::GFE::HarmonicDensity<Dune::FieldVector<DT,gridDim>,TargetSpace> density;
+
   RT energy = 0;
 
   const auto& localFiniteElement = localView.tree().finiteElement();
@@ -51,27 +55,24 @@ energy(const typename Basis::LocalView& localView,
 
     const auto integrationElement = element.geometry().integrationElement(quadPos);
 
-    const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(quadPos);
+    const auto jacobianInverse = element.geometry().jacobianInverse(quadPos);
 
     auto weight = quad[pt].weight() * integrationElement;
 
     // The derivative of the local function defined on the reference element
-    auto referenceDerivative = localInterpolationRule.evaluateDerivative(quadPos);
-
-    // Compute the Frobenius norm squared of the derivative.  This is the correct term
-    // if both domain and target space use the metric inherited from an embedding.
-    for (size_t i=0; i<jacobianInverseTransposed.N(); i++)
-      for (int j=0; j<TargetSpace::embeddedDim; j++)
-      {
-        RT entry = 0;
-        for (size_t k=0; k<jacobianInverseTransposed.M(); k++)
-          entry += jacobianInverseTransposed[i][k] * referenceDerivative[j][k];
-        energy += weight * entry * entry;
-      }
 
+    // Function value at the point where we are evaluating the derivative
+    // (not needed by the energy, but needed to compute the derivative)
+    TargetSpace q = localInterpolationRule.evaluate(quadPos);
+
+    // Compute the derivative
+    auto referenceDerivative = localInterpolationRule.evaluateDerivative(quadPos, q);
+    auto derivative = referenceDerivative * jacobianInverse;
+
+    energy += weight * density(quadPos,q,derivative);
   }
 
-  return 0.5 * energy;
+  return energy;
 }
 
 #endif

dune/gfe/densities/CMakeLists.txt

 #install headers
 install(FILES
         bulkcosseratdensity.hh
+        harmonicdensity.hh
         localdensity.hh
         DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/gfe/densities)

dune/gfe/densities/bulkcosseratdensity.hh

@@ -110,6 +110,18 @@ namespace Dune::GFE {
                                   + b2_ * GFE::skew(wryness).frobenius_norm2() + b3_ * GFE::traceSquared(wryness));
     }
 
+    /** \brief Evaluate the density
+     *
+     * \todo Currently this method is only here to please the compiler.
+     * Eventually, it will replace the other operator()-implementation.
+     */
+    virtual field_type operator() (const Position& x,
+                                   const GFE::ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> >& value,
+                                   const FieldMatrix<field_type,7,gridDim>& derivative) const override
+    {
+      DUNE_THROW(NotImplemented, "!");
+    }
+
     /** \brief Evaluation with the current position, the deformation value and its derivative, the orientation value and its derivative
      *
      * \param x The current position

dune/gfe/densities/harmonicdensity.hh

+#ifndef DUNE_GFE_DENSITIES_HARMONICDENSITY_HH
+#define DUNE_GFE_DENSITIES_HARMONICDENSITY_HH
+
+#include <dune/common/fmatrix.hh>
+
+#include <dune/gfe/densities/localdensity.hh>
+
+namespace Dune::GFE
+{
+  template<class Position, class TargetSpace>
+  class HarmonicDensity final
+    : public LocalDensity<Position,TargetSpace>
+  {
+    using field_type = typename TargetSpace::field_type;
+
+    constexpr static auto dim = Position::size();
+    constexpr static auto embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
+
+  public:
+
+    /** \brief Evaluate the density
+     *
+     * \param x The current position
+     * \param value The deformation at the current position
+     * \param derivative The derivative of the deformation at the current position
+     */
+    virtual field_type operator() (const Position& x,
+                                   const TargetSpace& value,
+                                   const FieldMatrix<field_type,embeddedBlocksize,dim>& derivative) const override
+    {
+      return 0.5 * derivative.frobenius_norm2();
+    }
+  };
+
+}  // namespace Dune:GFE
+
+#endif

dune/gfe/densities/localdensity.hh

@@ -16,8 +16,20 @@ namespace Dune::GFE {
   class LocalDensity
   {
     using field_type = typename TargetSpace::field_type;
+    using DerivativeType = FieldMatrix<field_type,TargetSpace::EmbeddedTangentVector::dimension,Position::size()>;
+
   public:
 
+    /** \brief Evaluation with the current position, the deformation function, the deformation gradient, the rotation and the rotation gradient
+     *
+     * \param x The current position
+     * \param value The deformation at the current position
+     * \param derivative The derivative of the deformation at the current position
+     */
+    virtual field_type operator() (const Position& x,
+                                   const TargetSpace& value,
+                                   const DerivativeType& derivative) const = 0;
+
     /** \brief Evaluation with the current position, the deformation function, the deformation gradient, the rotation and the rotation gradient
      *
      * \param x The current position
@@ -25,12 +37,18 @@ namespace Dune::GFE {
      * \param deformationDerivative The derivative of the deformation at the current position
      * \param orientationValue The orientation at the current position
      * \param orientationDerivative The derivative of the orientation at the current position
+     *
+     * \deprecated This is still used by Cosserat models, but shell be removed
+     * of the method above eventually.
      */
     virtual field_type operator() (const Position& x,
                                    const RealTuple<field_type,3>& deformation,
                                    const FieldMatrix<field_type,3,3>& gradient,
                                    const Rotation<field_type,3>& rotation,
-                                   const FieldMatrix<field_type, 4, 3>& rotationGradient) const = 0;
+                                   const FieldMatrix<field_type, 4, 3>& rotationGradient) const
+    {
+      std::terminate();
+    }
 
   };