MixedLocalGeodesicFEStiffness: Use a tuple for the coefficients

dune/gfe/assemblers/mixedgfeassembler.hh

@@ -169,8 +169,9 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
     // nDofs1 are the degrees of freedom for *one* subspacebasis of the power basis of the rotational part
     // this is why the indices (_0,0) and (_1,0) are used: _0 takes the whole displacement part and _1 the whole rotational part; and 0 the first subspacebasis respectively
     // Extract local solution
-    std::vector<TargetSpace0> localConfiguration0(nDofs0);
-    std::vector<TargetSpace1> localConfiguration1(nDofs1);
+    Dune::TupleVector<std::vector<TargetSpace0>, std::vector<TargetSpace1> > localConfiguration;
+    localConfiguration[_0].resize(nDofs0);
+    localConfiguration[_1].resize(nDofs1);
 
     for (int i=0; i<nDofs0+nDofs1; i++)
     {
@@ -186,9 +187,9 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
       //CompositeBasis number is contained in multiIndex[0], the Subspacebasis is contained in multiIndex[2]
       //multiIndex[1] contains the actual index
       if (multiIndex[0] == 0)
-        localConfiguration0[i] = configuration0[multiIndex[1]];
+        localConfiguration[_0][i] = configuration0[multiIndex[1]];
       else if (multiIndex[0] == 1)
-        localConfiguration1[i-nDofs0] = configuration1[multiIndex[1]];
+        localConfiguration[_1][i-nDofs0] = configuration1[multiIndex[1]];
     }
 
     std::vector<Dune::FieldVector<double,blocksize0> > localGradient0(nDofs0);
@@ -205,7 +206,7 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
 
     // setup local matrix and gradient
     localStiffness_->assembleGradientAndHessian(localView,
-                                                localConfiguration0, localConfiguration1,
+                                                localConfiguration,
                                                 localGradient0, localGradient1,
                                                 localHessian);
 
@@ -282,7 +283,7 @@ computeEnergy(const std::vector<TargetSpace0>& configuration0,
     localView.bind(element);
 
     // Number of degrees of freedom on this element
-    using namespace Dune::TypeTree::Indices;
+    using namespace Dune::Indices;
     const int nDofs0 = localView.tree().child(_0,0).finiteElement().size();
     const int nDofs1 = localView.tree().child(_1,0).finiteElement().size();
 

dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh

@@ -14,6 +14,7 @@ namespace Dune::GFE
   {
     template<class TargetSpace>
     class MixedLocalStiffnessTypes
+      : public LocalFirstOrderModelTypes<TargetSpace>
     {
       // Number type
       typedef typename TargetSpace::ctype RT;
@@ -61,14 +62,10 @@ public:
   /** \brief Assemble the local stiffness matrix at the current position
    */
   virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
-                                          const std::vector<DeformationTargetSpace>& localDisplacementConfiguration,
-                                          const std::vector<OrientationTargetSpace>& localOrientationConfiguration,
+                                          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>::MixedHessian& localHessian)
-  {
-    DUNE_THROW(Dune::NotImplemented, "!");
-  }
+                                          typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::MixedHessian& localHessian) const = 0;
 };
 
 #endif

dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh

@@ -79,11 +79,10 @@ public:
      This uses the automatic differentiation toolbox ADOL_C.
    */
   virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
-                                          const std::vector<TargetSpace0>& localConfiguration0,
-                                          const std::vector<TargetSpace1>& localConfiguration1,
+                                          const typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
                                           std::vector<typename TargetSpace0::TangentVector>& localGradient0,
                                           std::vector<typename TargetSpace1::TangentVector>& localGradient1,
-                                          HessianType& localHessian) override;
+                                          HessianType& localHessian) const override;
 
   const Dune::GFE::LocalEnergy<Basis, Dune::GFE::ProductManifold<ATargetSpace0, ATargetSpace1> >* localEnergy_;
   const bool adolcScalarMode_;
@@ -158,19 +157,15 @@ energy(const typename Basis::LocalView& localView,
 template <class Basis, class TargetSpace0, class TargetSpace1>
 void MixedLocalGFEADOLCStiffness<Basis, TargetSpace0, TargetSpace1>::
 assembleGradientAndHessian(const typename Basis::LocalView& localView,
-                           const std::vector<TargetSpace0>& localConfiguration0,
-                           const std::vector<TargetSpace1>& localConfiguration1,
+                           const typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
                            std::vector<typename TargetSpace0::TangentVector>& localGradient0,
                            std::vector<typename TargetSpace1::TangentVector>& localGradient1,
-                           HessianType& localHessian)
+                           HessianType& localHessian) const
 {
   int rank = Dune::MPIHelper::getCommunication().rank();
 
   // Tape energy computation.  We may not have to do this every time, but it's comparatively cheap.
   using namespace Dune::Indices;
-  typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients localConfiguration;
-  localConfiguration[_0] = localConfiguration0;
-  localConfiguration[_1] = localConfiguration1;
 
   energy(localView, localConfiguration);
 
@@ -180,43 +175,43 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
   /////////////////////////////////////////////////////////////////
 
   // Compute the actual gradient
-  size_t nDofs0 = localConfiguration0.size();
-  size_t nDofs1 = localConfiguration1.size();
+  size_t nDofs0 = localConfiguration[_0].size();
+  size_t nDofs1 = localConfiguration[_1].size();
   size_t nDoubles = nDofs0*embeddedBlocksize0 + nDofs1*embeddedBlocksize1;
 
   std::vector<double> xp(nDoubles);
   int idx=0;
-  for (size_t i=0; i<localConfiguration0.size(); i++)
+  for (size_t i=0; i<localConfiguration[_0].size(); i++)
     for (size_t j=0; j<embeddedBlocksize0; j++)
-      xp[idx++] = localConfiguration0[i].globalCoordinates()[j];
+      xp[idx++] = localConfiguration[_0][i].globalCoordinates()[j];
 
-  for (size_t i=0; i<localConfiguration1.size(); i++)
+  for (size_t i=0; i<localConfiguration[_1].size(); i++)
     for (size_t j=0; j<embeddedBlocksize1; j++)
-      xp[idx++] = localConfiguration1[i].globalCoordinates()[j];
+      xp[idx++] = localConfiguration[_1][i].globalCoordinates()[j];
 
   // Compute gradient
   std::vector<double> g(nDoubles);
   gradient(rank,nDoubles,xp.data(),g.data());                    // gradient evaluation
 
   // Copy into Dune type
-  std::vector<typename TargetSpace0::EmbeddedTangentVector> localEmbeddedGradient0(localConfiguration0.size());
-  std::vector<typename TargetSpace1::EmbeddedTangentVector> localEmbeddedGradient1(localConfiguration1.size());
+  std::vector<typename TargetSpace0::EmbeddedTangentVector> localEmbeddedGradient0(localConfiguration[_0].size());
+  std::vector<typename TargetSpace1::EmbeddedTangentVector> localEmbeddedGradient1(localConfiguration[_1].size());
 
   idx=0;
-  for (size_t i=0; i<localConfiguration0.size(); i++) {
+  for (size_t i=0; i<localConfiguration[_0].size(); i++) {
     for (size_t j=0; j<embeddedBlocksize0; j++)
       localEmbeddedGradient0[i][j] = g[idx++];
 
     // Express gradient in local coordinate system
-    localConfiguration0[i].orthonormalFrame().mv(localEmbeddedGradient0[i],localGradient0[i]);
+    localConfiguration[_0][i].orthonormalFrame().mv(localEmbeddedGradient0[i],localGradient0[i]);
   }
 
-  for (size_t i=0; i<localConfiguration1.size(); i++) {
+  for (size_t i=0; i<localConfiguration[_1].size(); i++) {
     for (size_t j=0; j<embeddedBlocksize1; j++)
       localEmbeddedGradient1[i][j] = g[idx++];
 
     // Express gradient in local coordinate system
-    localConfiguration1[i].orthonormalFrame().mv(localEmbeddedGradient1[i],localGradient1[i]);
+    localConfiguration[_1][i].orthonormalFrame().mv(localEmbeddedGradient1[i],localGradient1[i]);
   }
 
   /////////////////////////////////////////////////////////////////
@@ -249,12 +244,12 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
   std::vector<Dune::FieldMatrix<RT,blocksize0,embeddedBlocksize0> > orthonormalFrame0(nDofs0);
 
   for (size_t i=0; i<nDofs0; i++)
-    orthonormalFrame0[i] = localConfiguration0[i].orthonormalFrame();
+    orthonormalFrame0[i] = localConfiguration[_0][i].orthonormalFrame();
 
   std::vector<Dune::FieldMatrix<RT,blocksize1,embeddedBlocksize1> > orthonormalFrame1(nDofs1);
 
   for (size_t i=0; i<nDofs1; i++)
-    orthonormalFrame1[i] = localConfiguration1[i].orthonormalFrame();
+    orthonormalFrame1[i] = localConfiguration[_1][i].orthonormalFrame();
 
   Dune::Matrix<Dune::FieldMatrix<double,blocksize0, embeddedBlocksize0> > embeddedHessian00(nDofs0,nDofs0);
   Dune::Matrix<Dune::FieldMatrix<double,blocksize0, embeddedBlocksize1> > embeddedHessian01(nDofs0,nDofs1);
@@ -473,11 +468,11 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
   // Project embedded gradient onto normal space
   std::vector<typename TargetSpace0::EmbeddedTangentVector> projectedGradient0(nDofs0);
   for (size_t i=0; i<nDofs0; i++)
-    projectedGradient0[i] = localConfiguration0[i].projectOntoNormalSpace(localEmbeddedGradient0[i]);
+    projectedGradient0[i] = localConfiguration[_0][i].projectOntoNormalSpace(localEmbeddedGradient0[i]);
 
   std::vector<typename TargetSpace1::EmbeddedTangentVector> projectedGradient1(nDofs1);
   for (size_t i=0; i<nDofs1; i++)
-    projectedGradient1[i] = localConfiguration1[i].projectOntoNormalSpace(localEmbeddedGradient1[i]);
+    projectedGradient1[i] = localConfiguration[_1][i].projectOntoNormalSpace(localEmbeddedGradient1[i]);
 
   // The Weingarten map has only diagonal entries
   for (size_t row=0; row<nDofs0; row++) {
@@ -485,7 +480,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
     for (size_t subRow=0; subRow<blocksize0; subRow++) {
 
       typename TargetSpace0::EmbeddedTangentVector z = orthonormalFrame0[row][subRow];
-      typename TargetSpace0::EmbeddedTangentVector tmp1 = localConfiguration0[row].weingarten(z,projectedGradient0[row]);
+      typename TargetSpace0::EmbeddedTangentVector tmp1 = localConfiguration[_0][row].weingarten(z,projectedGradient0[row]);
 
       typename TargetSpace0::TangentVector tmp2;
       orthonormalFrame0[row].mv(tmp1,tmp2);
@@ -500,7 +495,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
     for (size_t subRow=0; subRow<blocksize1; subRow++) {
 
       typename TargetSpace1::EmbeddedTangentVector z = orthonormalFrame1[row][subRow];
-      typename TargetSpace1::EmbeddedTangentVector tmp1 = localConfiguration1[row].weingarten(z,projectedGradient1[row]);
+      typename TargetSpace1::EmbeddedTangentVector tmp1 = localConfiguration[_1][row].weingarten(z,projectedGradient1[row]);
 
       typename TargetSpace1::TangentVector tmp2;
       orthonormalFrame1[row].mv(tmp1,tmp2);