From bb45a8de10a2c1d565c59fdbd4a2cd405fc0a58b Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Tue, 6 Feb 2024 03:55:22 +0100
Subject: [PATCH] Make MixedLocalGeodesicFEADOLCStiffness accept a
ProductManifold
This is what the non-mixed ADOL-C assembler does, and I want to
merge the two.
---
.../assemblers/mixedlocalgfeadolcstiffness.hh | 49 ++++++++++---------
src/cosserat-continuum.cc | 10 ++--
src/film-on-substrate.cc | 7 ++-
src/simofoxshell.cc | 6 +--
test/cosseratcontinuumtest.cc | 3 +-
test/geodesicfeassemblerwrappertest.cc | 3 +-
6 files changed, 39 insertions(+), 39 deletions(-)
diff --git a/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh b/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
index 03e481f9..02d287d0 100644
--- a/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
@@ -17,36 +17,23 @@
/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
*/
-template<class Basis, class TargetSpace0, class TargetSpace1>
+template<class Basis, class TargetSpace>
class MixedLocalGFEADOLCStiffness
- : public LocalGeodesicFEStiffness<Basis,Dune::GFE::ProductManifold<TargetSpace0,TargetSpace1> >
+ : public LocalGeodesicFEStiffness<Basis,TargetSpace>
{
- using TargetSpace = Dune::GFE::ProductManifold<TargetSpace0,TargetSpace1>;
-
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
- typedef typename TargetSpace0::ctype RT;
- typedef typename GridView::template Codim<0>::Entity Entity;
+ typedef typename TargetSpace::ctype RT;
// The 'active' target spaces, i.e., the number type is replaced by adouble
using ATargetSpace = typename TargetSpace::template rebind<adouble>::other;
- typedef typename TargetSpace0::template rebind<adouble>::other ATargetSpace0;
- typedef typename TargetSpace1::template rebind<adouble>::other ATargetSpace1;
// some other sizes
constexpr static int gridDim = GridView::dimension;
public:
- //! Dimension of a tangent space
- constexpr static int blocksize0 = TargetSpace0::TangentVector::dimension;
- constexpr static int blocksize1 = TargetSpace1::TangentVector::dimension;
-
- //! Dimension of the embedding space
- constexpr static int embeddedBlocksize0 = TargetSpace0::EmbeddedTangentVector::dimension;
- constexpr static int embeddedBlocksize1 = TargetSpace1::EmbeddedTangentVector::dimension;
-
MixedLocalGFEADOLCStiffness(const Dune::GFE::LocalEnergy<Basis,ATargetSpace>* energy, bool adolcScalarMode = false)
: localEnergy_(energy),
adolcScalarMode_(adolcScalarMode)
@@ -89,26 +76,32 @@ public:
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeGradient& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const override;
- const Dune::GFE::LocalEnergy<Basis, Dune::GFE::ProductManifold<ATargetSpace0, ATargetSpace1> >* localEnergy_;
+ const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* localEnergy_;
const bool adolcScalarMode_;
};
-template <class Basis, class TargetSpace0, class TargetSpace1>
-typename MixedLocalGFEADOLCStiffness<Basis, TargetSpace0, TargetSpace1>::RT
-MixedLocalGFEADOLCStiffness<Basis, TargetSpace0, TargetSpace1>::
+template <class Basis, class TargetSpace>
+typename MixedLocalGFEADOLCStiffness<Basis, TargetSpace>::RT
+MixedLocalGFEADOLCStiffness<Basis, TargetSpace>::
energy(const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& localConfiguration) const
{
using namespace Dune::Indices;
int rank = Dune::MPIHelper::getCommunication().rank();
- double pureEnergy;
+
+ using TargetSpace0 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_0])>;
+ using TargetSpace1 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_1])>;
+
+ using ATargetSpace0 = typename TargetSpace0::template rebind<adouble>::other;
+ using ATargetSpace1 = typename TargetSpace1::template rebind<adouble>::other;
Dune::TupleVector<std::vector<ATargetSpace0>, std::vector<ATargetSpace1> > localAConfiguration;
localAConfiguration[_0].resize(localConfiguration[_0].size());
localAConfiguration[_1].resize(localConfiguration[_1].size());
+ double pureEnergy;
trace_on(rank);
adouble energy = 0;
@@ -159,8 +152,8 @@ energy(const typename Basis::LocalView& localView,
// To compute the Hessian we need to compute the gradient anyway, so we may
// as well return it. This saves assembly time.
// ///////////////////////////////////////////////////////////
-template <class Basis, class TargetSpace0, class TargetSpace1>
-void MixedLocalGFEADOLCStiffness<Basis, TargetSpace0, TargetSpace1>::
+template <class Basis, class TargetSpace>
+void MixedLocalGFEADOLCStiffness<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,
@@ -178,6 +171,16 @@ assembleGradientAndHessian(const typename Basis::LocalView& localView,
// into the correct coordinates.
/////////////////////////////////////////////////////////////////
+ //! Dimension of the embedding space
+ using TargetSpace0 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_0])>;
+ using TargetSpace1 = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_1])>;
+
+ constexpr static int blocksize0 = TargetSpace0::TangentVector::dimension;
+ constexpr static int blocksize1 = TargetSpace1::TangentVector::dimension;
+
+ constexpr static int embeddedBlocksize0 = TargetSpace0::EmbeddedTangentVector::dimension;
+ constexpr static int embeddedBlocksize1 = TargetSpace1::EmbeddedTangentVector::dimension;
+
// Compute the actual gradient
size_t nDofs0 = localConfiguration[_0].size();
size_t nDofs1 = localConfiguration[_1].size();
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index 9eae2869..6220cd31 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -97,10 +97,10 @@ const int rotationOrder = GFE_ORDER;
#if !MIXED_SPACE
static_assert(displacementOrder==rotationOrder, "displacement and rotation order do not match!");
+#endif
// Image space of the geodesic fe functions
using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
-#endif
int main (int argc, char *argv[]) try
@@ -497,9 +497,9 @@ int main (int argc, char *argv[]) try
&neumannBoundary,
neumannFunction,
volumeLoad);
- MixedLocalGFEADOLCStiffness<CompositeBasis,
- RealTuple<double,3>,
- Rotation<double,3> > localGFEADOLCStiffness(&localCosseratEnergy, adolcScalarMode);
+
+ MixedLocalGFEADOLCStiffness<CompositeBasis,TargetSpace> localGFEADOLCStiffness(&localCosseratEnergy,
+ adolcScalarMode);
MixedGFEAssembler<CompositeBasis,
RealTuple<double,3>,
Rotation<double,3> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
@@ -588,7 +588,7 @@ int main (int argc, char *argv[]) try
}
#endif
} else { //dim != dimworld
- using StiffnessType = MixedLocalGFEADOLCStiffness<CompositeBasis, RealTuple<double,3>, Rotation<double,3> >;
+ using StiffnessType = MixedLocalGFEADOLCStiffness<CompositeBasis, TargetSpace>;
std::shared_ptr<StiffnessType> localGFEStiffness;
#if HAVE_DUNE_CURVEDGEOMETRY && WORLD_DIM == 3 && GRID_DIM == 2
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 811b7a71..6faf4932 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -518,14 +518,14 @@ int main (int argc, char *argv[]) try
fThickness,
fLame);
+ using RBM = GFE::ProductManifold<RealTuple<double, dim>,Rotation<double,dim> >;
+
GFE::SumEnergy<CompositeBasis, RealTuple<ValueType,targetDim>, Rotation<ValueType,targetDim> > sumEnergy;
sumEnergy.addLocalEnergy(neumannEnergy);
sumEnergy.addLocalEnergy(elasticEnergy);
sumEnergy.addLocalEnergy(surfaceCosseratEnergy);
- MixedLocalGFEADOLCStiffness<CompositeBasis,
- RealTuple<double,dim>,
- Rotation<double,dim> > localGFEADOLCStiffness(&sumEnergy);
+ MixedLocalGFEADOLCStiffness<CompositeBasis,RBM> localGFEADOLCStiffness(&sumEnergy);
MixedGFEAssembler<CompositeBasis,
RealTuple<double,dim>,
Rotation<double,dim> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
@@ -574,7 +574,6 @@ int main (int argc, char *argv[]) try
//The MixedRiemannianTrustRegionSolver can treat the Deformation and Orientation Space as separate ones
//The RiemannianTrustRegionSolver can only treat the Deformation and Rotation together in a ProductManifold
//Therefore, x and the dirichletDofs are converted to a ProductManifold structure, as well as the Hessian and Gradient that are returned by the assembler
- using RBM = GFE::ProductManifold<RealTuple<double, dim>,Rotation<double,dim> >;
std::vector<RBM> xRBM(compositeBasis.size({0}));
BitSetVector<RBM::TangentVector::dimension> dirichletDofsRBM(compositeBasis.size({0}), false);
for (int i = 0; i < compositeBasis.size({0}); i++) {
diff --git a/src/simofoxshell.cc b/src/simofoxshell.cc
index 29cab86f..431732b7 100644
--- a/src/simofoxshell.cc
+++ b/src/simofoxshell.cc
@@ -320,9 +320,10 @@ int main(int argc, char *argv[]) try
neumannFunction,
nullptr, x0);
+ using TargetSpace = Dune::GFE::ProductManifold<RealTuple<double,3>,UnitVector<double,3> >;
+
MixedLocalGFEADOLCStiffness<decltype(compositeBasis),
- RealTuple<double,3>,
- UnitVector<double,3> > localGFEADOLCStiffness(&simoFoxEnergyADOLCLocalStiffness);
+ TargetSpace> localGFEADOLCStiffness(&simoFoxEnergyADOLCLocalStiffness);
MixedGFEAssembler<decltype(compositeBasis),
RealTuple<double,3>, UnitVector<double,3> > assembler(compositeBasis, &localGFEADOLCStiffness);
@@ -388,7 +389,6 @@ int main(int argc, char *argv[]) try
x = solver.getSol();
} else {
#if !MIXED_SPACE
- using TargetSpace = Dune::GFE::ProductManifold<RealTuple<double,3>,UnitVector<double,3> >;
std::vector<TargetSpace> xTargetSpace(compositeBasis.size({0}));
BitSetVector<TargetSpace::TangentVector::dimension> dirichletDofsTargetSpace(compositeBasis.size({0}), false);
for (std::size_t i = 0; i < compositeBasis.size({0}); i++) {
diff --git a/test/cosseratcontinuumtest.cc b/test/cosseratcontinuumtest.cc
index 2680242d..cd1fd2f4 100644
--- a/test/cosseratcontinuumtest.cc
+++ b/test/cosseratcontinuumtest.cc
@@ -208,8 +208,7 @@ int main (int argc, char *argv[])
sumEnergy.addLocalEnergy(neumannEnergy);
MixedLocalGFEADOLCStiffness<CompositeBasis,
- RealTuple<double,dim>,
- Rotation<double,dim> > localGFEADOLCStiffness(&sumEnergy);
+ GFE::ProductManifold<RealTuple<double,dim>,Rotation<double,dim> > > localGFEADOLCStiffness(&sumEnergy);
MixedGFEAssembler<CompositeBasis,
RealTuple<double,dim>,
Rotation<double,dim> > mixedAssembler(compositeBasis, &localGFEADOLCStiffness);
diff --git a/test/geodesicfeassemblerwrappertest.cc b/test/geodesicfeassemblerwrappertest.cc
index 79d0940d..a3f1ad32 100644
--- a/test/geodesicfeassemblerwrappertest.cc
+++ b/test/geodesicfeassemblerwrappertest.cc
@@ -147,8 +147,7 @@ int main (int argc, char *argv[])
neumannFunction,
nullptr);
MixedLocalGFEADOLCStiffness<CompositeBasis,
- RealTuple<double,dim>,
- Rotation<double,dim> > mixedLocalGFEADOLCStiffness(&cosseratEnergy);
+ GFE::ProductManifold<RealTuple<double,dim>,Rotation<double,dim> > > mixedLocalGFEADOLCStiffness(&cosseratEnergy);
MixedGFEAssembler<CompositeBasis,
RealTuple<double,dim>,
Rotation<double,dim> > mixedAssembler(compositeBasis, &mixedLocalGFEADOLCStiffness);
--
GitLab