From dcec0e99850702e641bd25e2999d5ee204e83240 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Thu, 1 Feb 2024 17:01:49 +0100
Subject: [PATCH] Make the local Hessian a local variable of the global
assembler
This will later simplify merging the mixed and non-mixed versions
of the assembler code.
---
dune/gfe/assemblers/geodesicfeassembler.hh | 5 ++--
.../localgeodesicfeadolcstiffness.hh | 16 +++++-----
.../assemblers/localgeodesicfefdstiffness.hh | 14 +++++----
.../assemblers/localgeodesicfestiffness.hh | 7 ++---
dune/gfe/assemblers/mixedgfeassembler.hh | 20 +++++++++----
.../mixedlocalgeodesicfestiffness.hh | 18 +++++++-----
.../assemblers/mixedlocalgfeadolcstiffness.hh | 29 +++++++++++--------
test/adolctest.cc | 8 +++--
8 files changed, 69 insertions(+), 48 deletions(-)
diff --git a/dune/gfe/assemblers/geodesicfeassembler.hh b/dune/gfe/assemblers/geodesicfeassembler.hh
index f75264cf..467db9e3 100644
--- a/dune/gfe/assemblers/geodesicfeassembler.hh
+++ b/dune/gfe/assemblers/geodesicfeassembler.hh
@@ -168,9 +168,10 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
localSolution[i] = sol[localView.index(i)];
std::vector<Dune::FieldVector<double,blocksize> > localGradient(numOfBaseFct);
+ Dune::Matrix<Dune::FieldMatrix<double,blocksize,blocksize> > localHessian(numOfBaseFct,numOfBaseFct);
// setup local matrix and gradient
- localStiffness_->assembleGradientAndHessian(localView, localSolution, localGradient);
+ localStiffness_->assembleGradientAndHessian(localView, localSolution, localGradient, localHessian);
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
@@ -180,7 +181,7 @@ assembleGradientAndHessian(const std::vector<TargetSpace>& sol,
for (int j=0; j<numOfBaseFct; j++ ) {
auto col = localView.index(j);
- hessian[row][col] += localStiffness_->A_[i][j];
+ hessian[row][col] += localHessian[i][j];
}
}
diff --git a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
index 49937c2c..74565e23 100644
--- a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
@@ -47,7 +47,7 @@ public:
/** \brief Compute the energy at the current configuration */
virtual RT energy (const typename Basis::LocalView& localView,
- const std::vector<TargetSpace>& localSolution) const;
+ const std::vector<TargetSpace>& localSolution) const override;
/** \brief Assemble the element gradient of the energy functional
@@ -55,7 +55,7 @@ public:
*/
virtual void assembleGradient(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& solution,
- std::vector<typename TargetSpace::TangentVector>& gradient) const;
+ std::vector<typename TargetSpace::TangentVector>& gradient) const override;
/** \brief Assemble the local stiffness matrix at the current position
@@ -63,7 +63,8 @@ public:
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
- std::vector<typename TargetSpace::TangentVector>& localGradient);
+ std::vector<typename TargetSpace::TangentVector>& localGradient,
+ Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> >& localHessian) const override;
const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* localEnergy_;
const bool adolcScalarMode_;
@@ -166,7 +167,8 @@ template <class Basis, class TargetSpace>
void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
- std::vector<typename TargetSpace::TangentVector>& localGradient)
+ std::vector<typename TargetSpace::TangentVector>& localGradient,
+ Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> >& localHessian) const
{
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localSolution);
@@ -302,7 +304,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
typedef typename TargetSpace::TangentVector TangentVector;
- this->A_.setSize(nDofs,nDofs);
+ localHessian.setSize(nDofs,nDofs);
for (size_t col=0; col<nDofs; col++) {
@@ -316,7 +318,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
embeddedHessian[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize; subRow++)
- this->A_[row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+ localHessian[row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
}
}
@@ -343,7 +345,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
TangentVector tmp2;
orthonormalFrame[row].mv(tmp1,tmp2);
- this->A_[row][row][subRow] += tmp2;
+ localHessian[row][row][subRow] += tmp2;
}
}
diff --git a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
index fd2a36c4..442fc88b 100644
--- a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
@@ -55,7 +55,8 @@ public:
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
- std::vector<typename TargetSpace::TangentVector>& localGradient) override;
+ std::vector<typename TargetSpace::TangentVector>& localGradient,
+ typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const override;
const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* localEnergy_;
@@ -131,15 +132,16 @@ 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<typename TargetSpace::TangentVector>& localGradient,
+ typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const
{
// Number of degrees of freedom for this element
size_t nDofs = localSolution.size();
// Clear assemble data
- this->A_.setSize(nDofs, nDofs);
+ localHessian.setSize(nDofs, nDofs);
- this->A_ = 0;
+ localHessian = 0;
#ifdef MULTIPRECISION
const field_type eps = 1e-10;
@@ -245,9 +247,9 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
field_type foo = 0.5 * (forwardValue - 2*centerValue + backwardValue) / (eps*eps);
#ifdef MULTIPRECISION
- this->A_[i][j][i2][j2] = this->A_[j][i][j2][i2] = foo.template convert_to<double>();
+ localHessian[i][j][i2][j2] = localHessian[j][i][j2][i2] = foo.template convert_to<double>();
#else
- this->A_[i][j][i2][j2] = this->A_[j][i][j2][i2] = foo;
+ localHessian[i][j][i2][j2] = localHessian[j][i][j2][i2] = foo;
#endif
}
}
diff --git a/dune/gfe/assemblers/localgeodesicfestiffness.hh b/dune/gfe/assemblers/localgeodesicfestiffness.hh
index 24dab769..389caf17 100644
--- a/dune/gfe/assemblers/localgeodesicfestiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfestiffness.hh
@@ -26,7 +26,8 @@ public:
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localSolution,
- std::vector<typename TargetSpace::TangentVector>& localGradient) = 0;
+ std::vector<typename TargetSpace::TangentVector>& localGradient,
+ Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> >& localHessian) const = 0;
/** \brief Compute the energy at the current configuration */
virtual RT energy (const typename Basis::LocalView& localView,
@@ -37,10 +38,6 @@ public:
virtual void assembleGradient(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& solution,
std::vector<typename TargetSpace::TangentVector>& gradient) const = 0;
-
- // assembled data
- Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> > A_;
-
};
#endif
diff --git a/dune/gfe/assemblers/mixedgfeassembler.hh b/dune/gfe/assemblers/mixedgfeassembler.hh
index ebfe80f2..d7f111df 100644
--- a/dune/gfe/assemblers/mixedgfeassembler.hh
+++ b/dune/gfe/assemblers/mixedgfeassembler.hh
@@ -194,10 +194,20 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
std::vector<Dune::FieldVector<double,blocksize0> > localGradient0(nDofs0);
std::vector<Dune::FieldVector<double,blocksize1> > localGradient1(nDofs1);
+ using Row0 = Dune::MultiTypeBlockVector<Dune::Matrix<Dune::FieldMatrix<double, blocksize0, blocksize0> >,
+ Dune::Matrix<Dune::FieldMatrix<double, blocksize0, blocksize1> > >;
+ using Row1 = Dune::MultiTypeBlockVector<Dune::Matrix<Dune::FieldMatrix<double, blocksize1, blocksize0> >,
+ Dune::Matrix<Dune::FieldMatrix<double, blocksize1, blocksize1> > >;
+
+ using HessianType = Dune::MultiTypeBlockMatrix<Row0, Row1>;
+
+ HessianType localHessian;
+
// setup local matrix and gradient
localStiffness_->assembleGradientAndHessian(localView,
localConfiguration0, localConfiguration1,
- localGradient0, localGradient1);
+ localGradient0, localGradient1,
+ localHessian);
// Add element matrix to global stiffness matrix
for (int i=0; i<nDofs0+nDofs1; i++)
@@ -227,16 +237,16 @@ assembleGradientAndHessian(const std::vector<TargetSpace0>& configuration0,
auto col = localView.index(localIndexCol);
if (row[0]==0 and col[0]==0)
- hessian[_0][_0][row[1]][col[1]] += localStiffness_->A_[_0][_0][i][j];
+ hessian[_0][_0][row[1]][col[1]] += localHessian[_0][_0][i][j];
if (row[0]==0 and col[0]==1)
- hessian[_0][_1][row[1]][col[1]] += localStiffness_->A_[_0][_1][i][j-nDofs0];
+ hessian[_0][_1][row[1]][col[1]] += localHessian[_0][_1][i][j-nDofs0];
if (row[0]==1 and col[0]==0)
- hessian[_1][_0][row[1]][col[1]] += localStiffness_->A_[_1][_0][i-nDofs0][j];
+ hessian[_1][_0][row[1]][col[1]] += localHessian[_1][_0][i-nDofs0][j];
if (row[0]==1 and col[0]==1)
- hessian[_1][_1][row[1]][col[1]] += localStiffness_->A_[_1][_1][i-nDofs0][j-nDofs0];
+ hessian[_1][_1][row[1]][col[1]] += localHessian[_1][_1][i-nDofs0][j-nDofs0];
}
// Add local gradient to global gradient
diff --git a/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh b/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
index b74d94fb..7175a144 100644
--- a/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
+++ b/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
@@ -26,13 +26,22 @@ public:
constexpr static int blocksize0 = DeformationTargetSpace::TangentVector::dimension;
constexpr static int blocksize1 = OrientationTargetSpace::TangentVector::dimension;
+ // Type of the local Hessian
+ using Row0 = Dune::MultiTypeBlockVector<Dune::Matrix<Dune::FieldMatrix<RT, blocksize0, blocksize0> >,
+ Dune::Matrix<Dune::FieldMatrix<RT, blocksize0, blocksize1> > >;
+ using Row1 = Dune::MultiTypeBlockVector<Dune::Matrix<Dune::FieldMatrix<RT, blocksize1, blocksize0> >,
+ Dune::Matrix<Dune::FieldMatrix<RT, blocksize1, blocksize1> > >;
+
+ using HessianType = Dune::MultiTypeBlockMatrix<Row0, Row1>;
+
/** \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,
std::vector<typename DeformationTargetSpace::TangentVector>& localDeformationGradient,
- std::vector<typename OrientationTargetSpace::TangentVector>& localOrientationGradient)
+ std::vector<typename OrientationTargetSpace::TangentVector>& localOrientationGradient,
+ HessianType& hessian)
{
DUNE_THROW(Dune::NotImplemented, "!");
}
@@ -42,13 +51,6 @@ public:
const std::vector<DeformationTargetSpace>& localDeformationConfiguration,
const std::vector<OrientationTargetSpace>& localOrientationConfiguration) const = 0;
- // assembled tangent matrix
- using Row0 = Dune::MultiTypeBlockVector<Dune::Matrix<Dune::FieldMatrix<RT, blocksize0, blocksize0> >,
- Dune::Matrix<Dune::FieldMatrix<RT, blocksize0, blocksize1> > >;
- using Row1 = Dune::MultiTypeBlockVector<Dune::Matrix<Dune::FieldMatrix<RT, blocksize1, blocksize0> >,
- Dune::Matrix<Dune::FieldMatrix<RT, blocksize1, blocksize1> > >;
-
- Dune::MultiTypeBlockMatrix<Row0, Row1> A_;
};
#endif
diff --git a/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh b/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
index 603cb811..f30ff1df 100644
--- a/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
@@ -34,6 +34,9 @@ class MixedLocalGFEADOLCStiffness
// some other sizes
constexpr static int gridDim = GridView::dimension;
+ using Base = MixedLocalGeodesicFEStiffness<Basis,Dune::GFE::ProductManifold<TargetSpace0,TargetSpace1> >;
+ using HessianType = typename Base::HessianType;
+
public:
//! Dimension of a tangent space
@@ -63,7 +66,8 @@ public:
const std::vector<TargetSpace0>& localConfiguration0,
const std::vector<TargetSpace1>& localConfiguration1,
std::vector<typename TargetSpace0::TangentVector>& localGradient0,
- std::vector<typename TargetSpace1::TangentVector>& localGradient1) override;
+ std::vector<typename TargetSpace1::TangentVector>& localGradient1,
+ HessianType& localHessian) override;
const MixedLocalGeodesicFEStiffness<Basis, Dune::GFE::ProductManifold<ATargetSpace0, ATargetSpace1> >* localEnergy_;
const bool adolcScalarMode_;
@@ -140,7 +144,8 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
const std::vector<TargetSpace0>& localConfiguration0,
const std::vector<TargetSpace1>& localConfiguration1,
std::vector<typename TargetSpace0::TangentVector>& localGradient0,
- std::vector<typename TargetSpace1::TangentVector>& localGradient1)
+ std::vector<typename TargetSpace1::TangentVector>& localGradient1,
+ HessianType& localHessian)
{
int rank = Dune::MPIHelper::getCommunication().rank();
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
@@ -353,7 +358,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
using namespace Dune::Indices;
- this->A_[_0][_0].setSize(nDofs0,nDofs0);
+ localHessian[_0][_0].setSize(nDofs0,nDofs0);
for (size_t col=0; col<nDofs0; col++) {
@@ -367,14 +372,14 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
embeddedHessian00[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize0; subRow++)
- this->A_[_0][_0][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+ localHessian[_0][_0][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
}
}
}
- this->A_[_0][_1].setSize(nDofs0,nDofs1);
+ localHessian[_0][_1].setSize(nDofs0,nDofs1);
for (size_t col=0; col<nDofs1; col++) {
@@ -388,14 +393,14 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
embeddedHessian01[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize0; subRow++)
- this->A_[_0][_1][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+ localHessian[_0][_1][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
}
}
}
- this->A_[_1][_0].setSize(nDofs1,nDofs0);
+ localHessian[_1][_0].setSize(nDofs1,nDofs0);
for (size_t col=0; col<nDofs0; col++) {
@@ -409,14 +414,14 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
embeddedHessian10[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize1; subRow++)
- this->A_[_1][_0][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+ localHessian[_1][_0][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
}
}
}
- this->A_[_1][_1].setSize(nDofs1,nDofs1);
+ localHessian[_1][_1].setSize(nDofs1,nDofs1);
for (size_t col=0; col<nDofs1; col++) {
@@ -430,7 +435,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
embeddedHessian11[row][col].mv(z,semiEmbeddedProduct);
for (int subRow=0; subRow<blocksize1; subRow++)
- this->A_[_1][_1][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
+ localHessian[_1][_1][row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
}
}
@@ -462,7 +467,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
typename TargetSpace0::TangentVector tmp2;
orthonormalFrame0[row].mv(tmp1,tmp2);
- this->A_[_0][_0][row][row][subRow] += tmp2;
+ localHessian[_0][_0][row][row][subRow] += tmp2;
}
}
@@ -477,7 +482,7 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
typename TargetSpace1::TangentVector tmp2;
orthonormalFrame1[row].mv(tmp1,tmp2);
- this->A_[_1][_1][row][row][subRow] += tmp2;
+ localHessian[_1][_1][row][row][subRow] += tmp2;
}
}
diff --git a/test/adolctest.cc b/test/adolctest.cc
index d53f84b7..8b59ac4e 100644
--- a/test/adolctest.cc
+++ b/test/adolctest.cc
@@ -580,14 +580,16 @@ int main (int argc, char *argv[]) try
localGFEADOLCStiffness.assembleGradientAndHessian(localView,
localSolution,
- localRiemannianADGradient);
+ localRiemannianADGradient,
+ localRiemannianADHessian);
localGFEFDStiffness.assembleGradientAndHessian(localView,
localSolution,
- localRiemannianFDGradient);
+ localRiemannianFDGradient,
+ localRiemannianFDHessian);
// compare
- compareMatrices(localGFEADOLCStiffness.A_, "Riemannian AD", localGFEFDStiffness.A_, "Riemannian FD");
+ compareMatrices(localRiemannianADHessian, "Riemannian AD", localRiemannianFDHessian, "Riemannian FD");
}
--
GitLab