Skip to content
Snippets Groups Projects
Commit dcec0e99 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

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.
parent 8a568780
No related branches found
No related tags found
1 merge request!143Various cleanup patches
Hide whitespace changes
Inline Side-by-side
Showing
with 69 additions and 48 deletions
dune/gfe/assemblers/geodesicfeassembler.hh
+ 3
− 2
View file @ dcec0e99
......@@ -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];
}
}
......
dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
+ 9
− 7
View file @ dcec0e99
......@@ -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;
}
}
......
dune/gfe/assemblers/localgeodesicfefdstiffness.hh
+ 8
− 6
View file @ dcec0e99
......@@ -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
}
}
......
dune/gfe/assemblers/localgeodesicfestiffness.hh
+ 2
− 5
View file @ dcec0e99
......@@ -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
dune/gfe/assemblers/mixedgfeassembler.hh
+ 15
− 5
View file @ dcec0e99
......@@ -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
......
dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
+ 10
− 8
View file @ dcec0e99
......@@ -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
dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
+ 17
− 12
View file @ dcec0e99
......@@ -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;
}
}
......
test/adolctest.cc
+ 5
− 3
View file @ dcec0e99
......@@ -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");
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment