WIP: Add test computation for bending isometries using reducedcubichermitetrianglebasis

dune/gfe/energies/discretekirchhoffbendingenergy.hh

@@ -23,7 +23,7 @@
 namespace Dune::GFE
 {
 
-  template<class Basis, class LocalDiscreteKirchhoffFunction, class TargetSpace>
+  template<class Basis, class LocalDiscreteKirchhoffFunction, class LocalForce, class TargetSpace>
   class DiscreteKirchhoffBendingEnergy
     : public Dune::GFE::LocalEnergy<Basis,TargetSpace>
   {
@@ -41,12 +41,14 @@ namespace Dune::GFE
     typedef typename Dune::LagrangeSimplexLocalFiniteElement<DT, double, gridDim, 2> P2LocalFiniteElement;
 
   public:
-    DiscreteKirchhoffBendingEnergy(LocalDiscreteKirchhoffFunction &localFunction)
-        : localFunction_(localFunction)
+    DiscreteKirchhoffBendingEnergy(LocalDiscreteKirchhoffFunction &localFunction,
+                                   LocalForce &localForce)
+        : localFunction_(localFunction),
+          localForce_(localForce)
     {}
 
     /** \brief Assemble the energy for a single element */
-    virtual RT energy (const typename Basis::LocalView& localView,          // actually not needed?!
+    virtual RT energy (const typename Basis::LocalView& localView,         
                const std::vector<TargetSpace>& localSolution) const        // actually not needed?!
     {
       RT energy = 0;
@@ -65,14 +67,18 @@ namespace Dune::GFE
                                                             : (localFiniteElement.localBasis().order() * gridDim - 1) * 2;
 #endif
 
+      const auto &quad = QuadratureRules<double, gridDim>::rule(lagrangeLFE_.type(), quadOrder);
+
       const auto element = localView.element();
       localFunction_.bind(element);
+      localForce_.bind(element);
 
       auto geometry = element.geometry();
       auto gridView = localFunction_.gridView();
 
-
-      //! \brief Compute energy contribution from the harmonic energy:
+      /**
+       * @brief Compute energy contribution from the harmonic energy:
+       */
       RT harmonicEnergy = 0;
 
       /**
@@ -205,8 +211,7 @@ namespace Dune::GFE
        * The entries of discrete Hessian are stored in a Fieldmatrix<double,6,2> 
        * to compute the squared Frobenius norm via the method .frobenius_norm2()
        */
-      const auto &quad = QuadratureRules<double, gridDim>::rule(lagrangeLFE_.type(), quadOrder);
-
+      
       for (auto&& quadPoint : quad)
       {
         // Get Jacobians of the P2-Basis functions on current quadrature point
@@ -236,7 +241,7 @@ namespace Dune::GFE
         /**
          * @brief Compute squared Frobenius norm of the discrete Hessian
          */
-        auto quadResult = discreteHessian.frobenius_norm2();
+        // auto quadResult = discreteHessian.frobenius_norm2();
         // std::cout << "quadResult: " << weight*quadResult << std::endl;
 
         harmonicEnergy += weight * discreteHessian.frobenius_norm2();
@@ -245,20 +250,40 @@ namespace Dune::GFE
       }
 
 
-      // Compute contribution of the harmonic energy:
-      RT forcingTermEnergy = 0;
       /**
-       * @brief TODO: Add missing/optional force term contributation in energy
+       * @brief Compute contribution of the force Term
        * 
+       * Integrate the scalar product of the force 'f' 
+       * with the local deformation function 'localFunction_'
+       * over the current element.
        */
+      RT forceTermEnergy = 0;
 
-        return 0.5 * harmonicEnergy - forcingTermEnergy;
+      for (auto&& quadPoint : quad)
+      {
+          auto deformationValue = localFunction_.evaluate(quadPoint.position());
+          auto forceValue = localForce_(quadPoint.position());
+
+          // printvector(std::cout, deformationValue.globalCoordinates(), "deformationValue", "--");
+          // printvector(std::cout, forceValue, "forceValue", "--");
+
+          const auto jacobian = geometry.jacobianInverseTransposed(quadPoint.position());
+          auto weight = quadPoint.weight() * element.geometry().integrationElement(quadPoint.position());
+
+          // std::cout << "weight:" << weight << std::endl;
+          forceTermEnergy += deformationValue.globalCoordinates() * forceValue * weight;
+          // std::cout << "forceTermEnergy:" << forceTermEnergy << std::endl;
+      }
+      
+      // std::cout << "energy output:" << 0.5 * harmonicEnergy - forceTermEnergy << std::endl;
+        return 0.5 * harmonicEnergy - forceTermEnergy;
     }
 
   private:
-
     LocalDiscreteKirchhoffFunction& localFunction_;
 
+    LocalForce& localForce_;
+
     P2LocalFiniteElement lagrangeLFE_;
   };