From a7ed421501ae3c811db1289060b45d1100c62566 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Fri, 5 Apr 2024 06:39:50 +0200
Subject: [PATCH] Use plain std::vector<double> for element energy gradients
This is what dune-function mandates. We don't actually use
dune-functions here a lot, but it is still nice to stay in the spirit.
Also, I am thinking about using dune-functions bases for *some*
indexing work in the future.
Plus, I hope that getting rid of the blocking will improve run-times
of debug builds a little bit, in particular in the
LocalIntegralStiffness class.
---
dune/gfe/assemblers/geodesicfeassembler.hh | 10 +++---
dune/gfe/assemblers/localfirstordermodel.hh | 31 ++---------------
.../localgeodesicfeadolcstiffness.hh | 34 +++++++++++++------
.../assemblers/localgeodesicfefdstiffness.hh | 20 +++++------
.../assemblers/localgeodesicfestiffness.hh | 8 ++---
dune/gfe/assemblers/mixedgfeassembler.hh | 11 +++---
test/adolctest.cc | 4 +--
7 files changed, 52 insertions(+), 66 deletions(-)
diff --git a/dune/gfe/assemblers/geodesicfeassembler.hh b/dune/gfe/assemblers/geodesicfeassembler.hh
index 467db9e3..070c4164 100644
--- a/dune/gfe/assemblers/geodesicfeassembler.hh
+++ b/dune/gfe/assemblers/geodesicfeassembler.hh
@@ -167,7 +167,7 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[localView.index(i)];
- std::vector<Dune::FieldVector<double,blocksize> > localGradient(numOfBaseFct);
+ std::vector<double> localGradient(numOfBaseFct*blocksize);
Dune::Matrix<Dune::FieldMatrix<double,blocksize,blocksize> > localHessian(numOfBaseFct,numOfBaseFct);
// setup local matrix and gradient
@@ -188,7 +188,8 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
// Add local gradient to global gradient
for (int i=0; i<numOfBaseFct; i++)
- gradient[localView.index(i)] += localGradient[i];
+ for (int j=0; j<blocksize; j++)
+ gradient[localView.index(i)][j] += localGradient[i*blocksize + j];
}
@@ -223,13 +224,14 @@ assembleGradient(const std::vector<TargetSpace>& sol,
localSolution[i] = sol[localView.index(i)];
// Assemble local gradient
- std::vector<Dune::FieldVector<double,blocksize> > localGradient(nDofs);
+ std::vector<double> localGradient(nDofs*blocksize);
localStiffness_->assembleGradient(localView, localSolution, localGradient);
// Add to global gradient
for (size_t i=0; i<nDofs; i++)
- grad[localView.index(i)[0]] += localGradient[i];
+ for (size_t j=0; j<blocksize; j++)
+ grad[localView.index(i)[0]][j] += localGradient[i*blocksize+j];
}
}
diff --git a/dune/gfe/assemblers/localfirstordermodel.hh b/dune/gfe/assemblers/localfirstordermodel.hh
index 6e2f4132..7ab06435 100644
--- a/dune/gfe/assemblers/localfirstordermodel.hh
+++ b/dune/gfe/assemblers/localfirstordermodel.hh
@@ -5,33 +5,6 @@
namespace Dune::GFE
{
- namespace Impl
- {
- /** \brief A class exporting container types for sets of tangent vectors
- *
- * This generic template handles TargetSpaces that are not product manifolds.
- */
- template <class TargetSpace>
- struct LocalFirstOrderModelTypes
- : 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
- */
- template <class ... Factors>
- struct LocalFirstOrderModelTypes<ProductManifold<Factors...> >
- : public LocalEnergyTypes<ProductManifold<Factors...> >
- {
- using Gradient = std::vector<typename ProductManifold<Factors...>::TangentVector>;
- using CompositeGradient = TupleVector<std::vector<typename Factors::TangentVector>... >;
- };
-
- } // namespace Impl
-
-
/** \brief Base class for problems that have an energy and a first derivative
*/
template<class Basis, class TargetSpace>
@@ -42,8 +15,8 @@ namespace Dune::GFE
/** \brief Assemble the element gradient of the energy functional */
virtual void assembleGradient(const typename Basis::LocalView& localView,
- const typename Impl::LocalFirstOrderModelTypes<TargetSpace>::Coefficients& coefficients,
- typename Impl::LocalFirstOrderModelTypes<TargetSpace>::Gradient& gradient) const = 0;
+ const typename Impl::LocalEnergyTypes<TargetSpace>::Coefficients& coefficients,
+ std::vector<double>& gradient) const = 0;
};
diff --git a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
index 25eea316..9220b8d1 100644
--- a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
@@ -54,7 +54,7 @@ public:
*/
virtual void assembleGradient(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
- typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Gradient& gradient) const override;
+ std::vector<double>& gradient) const override;
/** \brief Assemble the local stiffness matrix at the current position
@@ -62,14 +62,14 @@ public:
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
- typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Gradient& localGradient,
+ std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const override;
/** \brief Assemble the local stiffness matrix at the current position
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
- typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeGradient& localGradient,
+ std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const override;
const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* localEnergy_;
@@ -204,7 +204,7 @@ template <class Basis, class TargetSpace>
void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
assembleGradient(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& localSolution,
- typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Gradient& localGradient) const
+ std::vector<double>& localGradient) const
{
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localSolution);
@@ -250,7 +250,7 @@ template <class Basis, class TargetSpace>
void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& localSolution,
- typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Gradient& localGradient,
+ std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const
{
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
@@ -287,8 +287,10 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
// Express gradient in local coordinate system
for (size_t i=0; i<nDofs; i++) {
- Dune::FieldMatrix<RT,blocksize,embeddedBlocksize> orthonormalFrame = localSolution[i].orthonormalFrame();
- orthonormalFrame.mv(localEmbeddedGradient[i],localGradient[i]);
+ typename TargetSpace::TangentVector tmp;
+ localSolution[i].orthonormalFrame().mv(localEmbeddedGradient[i],tmp);
+ for (size_t j=0; j<blocksize; j++)
+ localGradient[i*blocksize+j] = tmp[j];
}
/////////////////////////////////////////////////////////////////
@@ -445,7 +447,7 @@ template <class Basis, class TargetSpace>
void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
- typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeGradient& localGradient,
+ std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const
{
using namespace Dune::Indices;
@@ -456,7 +458,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
//static_assert(localConfiguration.size()==1);
return assembleGradientAndHessian(localView,
localConfiguration[_0],
- localGradient[_0],
+ localGradient,
localHessian[_0][_0]);
}
else
@@ -513,15 +515,25 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
localEmbeddedGradient0[i][j] = g[idx++];
// Express gradient in local coordinate system
- localConfiguration[_0][i].orthonormalFrame().mv(localEmbeddedGradient0[i],localGradient[_0][i]);
+ typename TargetSpace0::TangentVector tmp;
+ localConfiguration[_0][i].orthonormalFrame().mv(localEmbeddedGradient0[i],tmp);
+
+ for (size_t j=0; j<blocksize0; j++)
+ localGradient[i*blocksize0+j] = tmp[j];
}
+ auto offset = localConfiguration[_0].size()*blocksize0;
+
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
- localConfiguration[_1][i].orthonormalFrame().mv(localEmbeddedGradient1[i],localGradient[_1][i]);
+ typename TargetSpace1::TangentVector tmp;
+ localConfiguration[_1][i].orthonormalFrame().mv(localEmbeddedGradient1[i],tmp);
+
+ for (size_t j=0; j<blocksize1; j++)
+ localGradient[offset + i*blocksize1 + j] = tmp[j];
}
/////////////////////////////////////////////////////////////////
diff --git a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
index 6f4fe03a..b8c50b4c 100644
--- a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
@@ -53,7 +53,7 @@ public:
The default implementation in this class uses a finite difference approximation */
virtual void assembleGradient(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& solution,
- std::vector<typename TargetSpace::TangentVector>& gradient) const override;
+ std::vector<field_type>& gradient) const override;
/** \brief Assemble the local tangent matrix and gradient at the current position
@@ -61,7 +61,7 @@ public:
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
- std::vector<typename TargetSpace::TangentVector>& localGradient,
+ std::vector<field_type>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const override;
/** \brief Assemble the local tangent matrix and gradient at the current position
@@ -70,7 +70,7 @@ public:
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localSolution,
- typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeGradient& localGradient,
+ std::vector<field_type>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const override
{
DUNE_THROW(Dune::NotImplemented, "!");
@@ -84,7 +84,7 @@ template <class Basis, class TargetSpace, class field_type>
void LocalGeodesicFEFDStiffness<Basis, TargetSpace, field_type>::
assembleGradient(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
- std::vector<typename TargetSpace::TangentVector>& localGradient) const
+ std::vector<field_type>& localGradient) const
{
// ///////////////////////////////////////////////////////////
@@ -125,9 +125,9 @@ assembleGradient(const typename Basis::LocalView& localView,
field_type foo = (localEnergy_->energy(localView,forwardSolution) - localEnergy_->energy(localView, backwardSolution)) / (2*eps);
#ifdef MULTIPRECISION
- localGradient[i][j] = foo.template convert_to<double>();
+ localGradient[i*blocksize+j] = foo.template convert_to<double>();
#else
- localGradient[i][j] = foo;
+ localGradient[i*blocksize+j] = foo;
#endif
}
@@ -149,7 +149,7 @@ template <class Basis, class TargetSpace, class field_type>
void LocalGeodesicFEFDStiffness<Basis, TargetSpace, field_type>::
assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
- std::vector<typename TargetSpace::TangentVector>& localGradient,
+ std::vector<field_type>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const
{
// Number of degrees of freedom for this element
@@ -212,16 +212,16 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
// Compute gradient by finite-difference approximation
//////////////////////////////////////////////////////////////
- localGradient.resize(localSolution.size());
+ localGradient.resize(nDofs*blocksize);
for (size_t i=0; i<localSolution.size(); i++)
for (int j=0; j<blocksize; j++)
{
field_type foo = (forwardEnergy[i][j] - backwardEnergy[i][j]) / (2*eps);
#ifdef MULTIPRECISION
- localGradient[i][j] = foo.template convert_to<double>();
+ localGradient[i*blocksize+j] = foo.template convert_to<double>();
#else
- localGradient[i][j] = foo;
+ localGradient[i*blocksize+j] = foo;
#endif
}
diff --git a/dune/gfe/assemblers/localgeodesicfestiffness.hh b/dune/gfe/assemblers/localgeodesicfestiffness.hh
index c32bf683..b7848210 100644
--- a/dune/gfe/assemblers/localgeodesicfestiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfestiffness.hh
@@ -17,7 +17,7 @@ namespace Dune::GFE
*/
template<class TargetSpace>
class LocalStiffnessTypes
- : public LocalFirstOrderModelTypes<TargetSpace>
+ : public LocalEnergyTypes<TargetSpace>
{
// Number type
typedef typename TargetSpace::ctype RT;
@@ -40,7 +40,7 @@ namespace Dune::GFE
*/
template<class ... Factors>
class LocalStiffnessTypes<ProductManifold<Factors...> >
- : public LocalFirstOrderModelTypes<ProductManifold<Factors...> >
+ : public LocalEnergyTypes<ProductManifold<Factors...> >
{
using TargetSpace = ProductManifold<Factors...>;
@@ -84,14 +84,14 @@ public:
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
- typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Gradient& localGradient,
+ std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const = 0;
/** \brief Assemble the local gradient and stiffness matrix at the current position -- Composite version
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& coefficients,
- typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeGradient& localGradient,
+ std::vector<double>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const = 0;
};
diff --git a/dune/gfe/assemblers/mixedgfeassembler.hh b/dune/gfe/assemblers/mixedgfeassembler.hh
index 6b134d41..4df0b50a 100644
--- a/dune/gfe/assemblers/mixedgfeassembler.hh
+++ b/dune/gfe/assemblers/mixedgfeassembler.hh
@@ -192,10 +192,7 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
localConfiguration[_1][i-nDofs0] = configuration1[multiIndex[1]];
}
- Dune::TupleVector<std::vector<Dune::FieldVector<double,blocksize0> >,
- std::vector<Dune::FieldVector<double,blocksize1> > > localGradient;
- localGradient[_0].resize(nDofs0);
- localGradient[_1].resize(nDofs1);
+ std::vector<double> localGradient(nDofs0*blocksize0 + nDofs1*blocksize1);
using Row0 = Dune::MultiTypeBlockVector<Dune::Matrix<Dune::FieldMatrix<double, blocksize0, blocksize0> >,
Dune::Matrix<Dune::FieldMatrix<double, blocksize0, blocksize1> > >;
@@ -254,9 +251,11 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
// Add local gradient to global gradient
if (row[0] == 0)
- gradient0[row[1]] += localGradient[_0][i];
+ for (std::size_t j=0; j<blocksize0; j++)
+ gradient0[row[1]][j] += localGradient[i*blocksize0 + j];
else
- gradient1[row[1]] += localGradient[_1][i-nDofs0];
+ for (std::size_t j=0; j<blocksize1; j++)
+ gradient1[row[1]][j] += localGradient[nDofs0*blocksize0 + (i-nDofs0)*blocksize1 + j];
}
}
diff --git a/test/adolctest.cc b/test/adolctest.cc
index 8b59ac4e..0dd4ec71 100644
--- a/test/adolctest.cc
+++ b/test/adolctest.cc
@@ -572,8 +572,8 @@ int main (int argc, char *argv[]) try
compareMatrices(localADHessian, "AD scalar", localADVMHessian, "AD vector");
// Assemble Riemannian derivatives
- std::vector<Dune::FieldVector<double,blocksize> > localRiemannianADGradient(numOfBaseFct);
- std::vector<Dune::FieldVector<double,blocksize> > localRiemannianFDGradient(numOfBaseFct);
+ std::vector<double> localRiemannianADGradient(numOfBaseFct*blocksize);
+ std::vector<double> localRiemannianFDGradient(numOfBaseFct*blocksize);
Matrix<FieldMatrix<double,blocksize,blocksize> > localRiemannianADHessian;
Matrix<FieldMatrix<double,blocksize,blocksize> > localRiemannianFDHessian;
--
GitLab