From 3359c543f84dff0f3ecf8656d8742f91da9ea787 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Tue, 6 Feb 2024 01:56:34 +0100
Subject: [PATCH] Store local gradients in a TupleVector when using a product
 space

This makes the interface more consistent.
---
 dune/gfe/assemblers/localfirstordermodel.hh        |  2 ++
 dune/gfe/assemblers/mixedgfeassembler.hh           | 12 +++++++-----
 .../assemblers/mixedlocalgeodesicfestiffness.hh    | 14 +-------------
 dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh | 10 ++++------
 4 files changed, 14 insertions(+), 24 deletions(-)

diff --git a/dune/gfe/assemblers/localfirstordermodel.hh b/dune/gfe/assemblers/localfirstordermodel.hh
index e9bcb5319..6e2f41328 100644
--- a/dune/gfe/assemblers/localfirstordermodel.hh
+++ b/dune/gfe/assemblers/localfirstordermodel.hh
@@ -16,6 +16,7 @@ namespace Dune::GFE
       : public LocalEnergyTypes<TargetSpace>
     {
       using Gradient = std::vector<typename TargetSpace::TangentVector>;
+      using CompositeGradient = TupleVector<std::vector<typename TargetSpace::TangentVector> >;
     };
 
     /** \brief A class exporting container types for sets of tangent vectors -- specialization for product manifolds
@@ -25,6 +26,7 @@ namespace Dune::GFE
       : public LocalEnergyTypes<ProductManifold<Factors...> >
     {
       using Gradient = std::vector<typename ProductManifold<Factors...>::TangentVector>;
+      using CompositeGradient = TupleVector<std::vector<typename Factors::TangentVector>... >;
     };
 
   }  // namespace Impl
diff --git a/dune/gfe/assemblers/mixedgfeassembler.hh b/dune/gfe/assemblers/mixedgfeassembler.hh
index d3e8a8f8f..f4182de4b 100644
--- a/dune/gfe/assemblers/mixedgfeassembler.hh
+++ b/dune/gfe/assemblers/mixedgfeassembler.hh
@@ -192,8 +192,10 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
         localConfiguration[_1][i-nDofs0] = configuration1[multiIndex[1]];
     }
 
-    std::vector<Dune::FieldVector<double,blocksize0> > localGradient0(nDofs0);
-    std::vector<Dune::FieldVector<double,blocksize1> > localGradient1(nDofs1);
+    Dune::TupleVector<std::vector<Dune::FieldVector<double,blocksize0> >,
+        std::vector<Dune::FieldVector<double,blocksize1> > > localGradient;
+    localGradient[_0].resize(nDofs0);
+    localGradient[_1].resize(nDofs1);
 
     using Row0 = Dune::MultiTypeBlockVector<Dune::Matrix<Dune::FieldMatrix<double, blocksize0, blocksize0> >,
         Dune::Matrix<Dune::FieldMatrix<double, blocksize0, blocksize1> > >;
@@ -207,7 +209,7 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
     // setup local matrix and gradient
     localStiffness_->assembleGradientAndHessian(localView,
                                                 localConfiguration,
-                                                localGradient0, localGradient1,
+                                                localGradient,
                                                 localHessian);
 
     // Add element matrix to global stiffness matrix
@@ -252,9 +254,9 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
 
       // Add local gradient to global gradient
       if (row[0] == 0)
-        gradient0[row[1]] += localGradient0[i];
+        gradient0[row[1]] += localGradient[_0][i];
       else
-        gradient1[row[1]] += localGradient1[i-nDofs0];
+        gradient1[row[1]] += localGradient[_1][i-nDofs0];
     }
 
   }
diff --git a/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh b/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
index cd640f226..ae5620df4 100644
--- a/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
+++ b/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
@@ -47,24 +47,12 @@ template<class Basis, class TargetSpace>
 class MixedLocalGeodesicFEStiffness
   : public Dune::GFE::LocalFirstOrderModel<Basis,TargetSpace>
 {
-  using DeformationTargetSpace = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_0])>;
-  using OrientationTargetSpace = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_1])>;
-
-  // Number type
-  typedef typename TargetSpace::ctype RT;
-
 public:
-
-  //! Dimension of a tangent space
-  constexpr static int blocksize0 = DeformationTargetSpace::TangentVector::dimension;
-  constexpr static int blocksize1 = OrientationTargetSpace::TangentVector::dimension;
-
   /** \brief Assemble the local stiffness matrix at the current position
    */
   virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
                                           const typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
-                                          std::vector<typename DeformationTargetSpace::TangentVector>& localDeformationGradient,
-                                          std::vector<typename OrientationTargetSpace::TangentVector>& localOrientationGradient,
+                                          typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeGradient& localGradient,
                                           typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::MixedHessian& localHessian) const = 0;
 };
 
diff --git a/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh b/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
index bdcedd6b7..d62186934 100644
--- a/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
@@ -80,8 +80,7 @@ public:
    */
   virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
                                           const typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
-                                          std::vector<typename TargetSpace0::TangentVector>& localGradient0,
-                                          std::vector<typename TargetSpace1::TangentVector>& localGradient1,
+                                          typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeGradient& localGradient,
                                           HessianType& localHessian) const override;
 
   const Dune::GFE::LocalEnergy<Basis, Dune::GFE::ProductManifold<ATargetSpace0, ATargetSpace1> >* localEnergy_;
@@ -158,8 +157,7 @@ template <class Basis, class TargetSpace0, class TargetSpace1>
 void MixedLocalGFEADOLCStiffness<Basis, TargetSpace0, TargetSpace1>::
 assembleGradientAndHessian(const typename Basis::LocalView& localView,
                            const typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
-                           std::vector<typename TargetSpace0::TangentVector>& localGradient0,
-                           std::vector<typename TargetSpace1::TangentVector>& localGradient1,
+                           typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeGradient& localGradient,
                            HessianType& localHessian) const
 {
   int rank = Dune::MPIHelper::getCommunication().rank();
@@ -203,7 +201,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
       localEmbeddedGradient0[i][j] = g[idx++];
 
     // Express gradient in local coordinate system
-    localConfiguration[_0][i].orthonormalFrame().mv(localEmbeddedGradient0[i],localGradient0[i]);
+    localConfiguration[_0][i].orthonormalFrame().mv(localEmbeddedGradient0[i],localGradient[_0][i]);
   }
 
   for (size_t i=0; i<localConfiguration[_1].size(); i++) {
@@ -211,7 +209,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
       localEmbeddedGradient1[i][j] = g[idx++];
 
     // Express gradient in local coordinate system
-    localConfiguration[_1][i].orthonormalFrame().mv(localEmbeddedGradient1[i],localGradient1[i]);
+    localConfiguration[_1][i].orthonormalFrame().mv(localEmbeddedGradient1[i],localGradient[_1][i]);
   }
 
   /////////////////////////////////////////////////////////////////
-- 
GitLab