diff --git a/dune/gfe/assemblers/CMakeLists.txt b/dune/gfe/assemblers/CMakeLists.txt
index a890d6ab41268334e9380c202ed4ceb2ff70783b..e0dfa540445ee4ad122b247b682505a90a70846a 100644
--- a/dune/gfe/assemblers/CMakeLists.txt
+++ b/dune/gfe/assemblers/CMakeLists.txt
@@ -14,7 +14,6 @@ install(FILES
localgeodesicfestiffness.hh
localintegralenergy.hh
mixedgfeassembler.hh
- mixedlocalgeodesicfestiffness.hh
mixedlocalgfeadolcstiffness.hh
nonplanarcosseratshellenergy.hh
simofoxenergy.hh
diff --git a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
index fb59c5cd9437a4a758ea331daf3578031e08c929..93b6723a66aac4fde8402cac012b7023b58e97cc 100644
--- a/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh
@@ -63,17 +63,27 @@ public:
This uses the automatic differentiation toolbox ADOL_C.
*/
virtual void assembleGradient(const typename Basis::LocalView& localView,
- const std::vector<TargetSpace>& solution,
- std::vector<typename TargetSpace::TangentVector>& gradient) const override;
+ const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
+ typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Gradient& gradient) const override;
/** \brief Assemble the local stiffness matrix at the current position
This uses the automatic differentiation toolbox ADOL_C.
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
- const std::vector<TargetSpace>& localSolution,
- std::vector<typename TargetSpace::TangentVector>& localGradient,
- Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> >& localHessian) const override;
+ const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
+ typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Gradient& 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,
+ typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const override
+ {
+ DUNE_THROW(Dune::NotImplemented, "!");
+ }
const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* localEnergy_;
const bool adolcScalarMode_;
@@ -130,8 +140,8 @@ energy(const typename Basis::LocalView& localView,
template <class Basis, class TargetSpace>
void LocalGeodesicFEADOLCStiffness<Basis, TargetSpace>::
assembleGradient(const typename Basis::LocalView& localView,
- const std::vector<TargetSpace>& localSolution,
- std::vector<typename TargetSpace::TangentVector>& localGradient) const
+ const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& localSolution,
+ typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Gradient& localGradient) const
{
// Tape energy computation. We may not have to do this every time, but it's comparatively cheap.
energy(localView, localSolution);
@@ -175,9 +185,9 @@ assembleGradient(const typename Basis::LocalView& localView,
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,
- Dune::Matrix<Dune::FieldMatrix<RT,blocksize,blocksize> >& localHessian) const
+ const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& localSolution,
+ typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Gradient& 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.
energy(localView, localSolution);
diff --git a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
index cf400921d0c5636a164ea64f85533cabf78ba8c9..6f4fe03af2ea7882a6bfbf8e85f3b1b8d56b2b63 100644
--- a/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfefdstiffness.hh
@@ -45,7 +45,7 @@ public:
RT energy (const typename Basis::LocalView& localView,
const typename Dune::GFE::Impl::LocalEnergyTypes<TargetSpace>::CompositeCoefficients& coefficients) const override
{
- DUNE_THROW(Dune::NotImplemented, "!");
+ return localEnergy_->energy(localView,coefficients);
}
/** \brief Assemble the element gradient of the energy functional
@@ -64,6 +64,17 @@ public:
std::vector<typename TargetSpace::TangentVector>& localGradient,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const override;
+ /** \brief Assemble the local tangent matrix and gradient at the current position
+
+ This implementation uses finite-difference approximations
+ */
+ 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,
+ typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const override
+ {
+ DUNE_THROW(Dune::NotImplemented, "!");
+ }
const Dune::GFE::LocalEnergy<Basis, ATargetSpace>* localEnergy_;
diff --git a/dune/gfe/assemblers/localgeodesicfestiffness.hh b/dune/gfe/assemblers/localgeodesicfestiffness.hh
index d43f18ecafc61c9dcae6a09b595c02137a1f5074..c32bf68382af1849dc7dd25122df93cef20c51b2 100644
--- a/dune/gfe/assemblers/localgeodesicfestiffness.hh
+++ b/dune/gfe/assemblers/localgeodesicfestiffness.hh
@@ -3,6 +3,7 @@
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrix.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
#include <dune/gfe/assemblers/localfirstordermodel.hh>
@@ -10,6 +11,10 @@ namespace Dune::GFE
{
namespace Impl
{
+ /** \brief A class exporting container types for sets local Hesse matrices
+ *
+ * This generic template handles TargetSpaces that are not product manifolds.
+ */
template<class TargetSpace>
class LocalStiffnessTypes
: public LocalFirstOrderModelTypes<TargetSpace>
@@ -24,6 +29,46 @@ namespace Dune::GFE
// Type of the local Hessian
using Hessian = Matrix<FieldMatrix<RT, blocksize, blocksize> >;
+
+ using Row = MultiTypeBlockVector<Matrix<FieldMatrix<RT, blocksize, blocksize> > >;
+ using CompositeHessian = MultiTypeBlockMatrix<Row>;
+ };
+
+ /** \brief A class exporting container types for sets local Hesse matrices
+ *
+ * This is the specialization for product manifolds.
+ */
+ template<class ... Factors>
+ class LocalStiffnessTypes<ProductManifold<Factors...> >
+ : public LocalFirstOrderModelTypes<ProductManifold<Factors...> >
+ {
+ using TargetSpace = ProductManifold<Factors...>;
+
+ // Number type
+ typedef typename ProductManifold<Factors...>::ctype RT;
+
+ using DeformationTargetSpace = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_0])>;
+ using OrientationTargetSpace = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_1])>;
+
+ // Dimension of the product tangent space
+ constexpr static auto blocksize = TargetSpace::TangentVector::dimension;
+
+ // Dimensions of the individual factor tangent spaces
+ constexpr static auto blocksize0 = DeformationTargetSpace::TangentVector::dimension;
+ constexpr static auto blocksize1 = OrientationTargetSpace::TangentVector::dimension;
+
+ public:
+
+ // Type of the local Hessian
+ using Hessian = Matrix<FieldMatrix<RT, blocksize, blocksize> >;
+
+ // Type of the local Hessian
+ using Row0 = MultiTypeBlockVector<Matrix<FieldMatrix<RT, blocksize0, blocksize0> >,
+ Matrix<FieldMatrix<RT, blocksize0, blocksize1> > >;
+ using Row1 = MultiTypeBlockVector<Matrix<FieldMatrix<RT, blocksize1, blocksize0> >,
+ Matrix<FieldMatrix<RT, blocksize1, blocksize1> > >;
+
+ using CompositeHessian = MultiTypeBlockMatrix<Row0, Row1>;
};
}
}
@@ -41,6 +86,13 @@ public:
const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Coefficients& coefficients,
typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Gradient& 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,
+ typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const = 0;
};
#endif
diff --git a/dune/gfe/assemblers/mixedgfeassembler.hh b/dune/gfe/assemblers/mixedgfeassembler.hh
index f4182de4bf2113a77da8b48460feabe23dc37bf1..6b134d41193dfc001ed4311fa3daefd5279163a0 100644
--- a/dune/gfe/assemblers/mixedgfeassembler.hh
+++ b/dune/gfe/assemblers/mixedgfeassembler.hh
@@ -7,7 +7,7 @@
#include <dune/istl/matrix.hh>
#include <dune/istl/multitypeblockmatrix.hh>
-#include <dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh>
+#include <dune/gfe/assemblers/localgeodesicfestiffness.hh>
/** \brief A global FE assembler for problems involving functions that map into non-Euclidean spaces
@@ -36,13 +36,13 @@ public:
Dune::MultiTypeBlockVector<MatrixBlock10,MatrixBlock11> > MatrixType;
const Basis basis_;
- MixedLocalGeodesicFEStiffness<Basis, Dune::GFE::ProductManifold<TargetSpace0,TargetSpace1> >* localStiffness_;
+ LocalGeodesicFEStiffness<Basis, Dune::GFE::ProductManifold<TargetSpace0,TargetSpace1> >* localStiffness_;
public:
/** \brief Constructor for a given grid */
MixedGFEAssembler(const Basis& basis,
- MixedLocalGeodesicFEStiffness<Basis, Dune::GFE::ProductManifold<TargetSpace0, TargetSpace1> >* localStiffness)
+ LocalGeodesicFEStiffness<Basis, Dune::GFE::ProductManifold<TargetSpace0, TargetSpace1> >* localStiffness)
: basis_(basis),
localStiffness_(localStiffness)
{}
diff --git a/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh b/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
deleted file mode 100644
index ae5620df418df9eecfa7ba56c93a6256598c7877..0000000000000000000000000000000000000000
--- a/dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh
+++ /dev/null
@@ -1,59 +0,0 @@
-#ifndef DUNE_GFE_MIXEDLOCALGEODESICFESTIFFNESS_HH
-#define DUNE_GFE_MIXEDLOCALGEODESICFESTIFFNESS_HH
-
-#include <dune/common/fmatrix.hh>
-
-#include <dune/istl/matrix.hh>
-#include <dune/istl/multitypeblockmatrix.hh>
-
-#include <dune/gfe/assemblers/localfirstordermodel.hh>
-
-namespace Dune::GFE
-{
- namespace Impl
- {
- template<class TargetSpace>
- class MixedLocalStiffnessTypes
- : public LocalFirstOrderModelTypes<TargetSpace>
- {
- // Number type
- typedef typename TargetSpace::ctype RT;
-
- using DeformationTargetSpace = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_0])>;
- using OrientationTargetSpace = std::decay_t<decltype(std::declval<TargetSpace>()[Dune::Indices::_1])>;
-
- //! Dimension of a tangent space
- constexpr static int blocksize0 = DeformationTargetSpace::TangentVector::dimension;
- constexpr static int blocksize1 = OrientationTargetSpace::TangentVector::dimension;
-
- public:
-
- // Type of the local Hessian
- using Row0 = MultiTypeBlockVector<Matrix<FieldMatrix<RT, blocksize0, blocksize0> >,
- Matrix<FieldMatrix<RT, blocksize0, blocksize1> > >;
- using Row1 = MultiTypeBlockVector<Matrix<FieldMatrix<RT, blocksize1, blocksize0> >,
- Matrix<FieldMatrix<RT, blocksize1, blocksize1> > >;
-
- using MixedHessian = Dune::MultiTypeBlockMatrix<Row0, Row1>;
- };
- }
-}
-
-/** \brief Abstract base class for second-order energy approximations on one grid element
- *
- * \tparam TargetSpace The space we map into. MUST be a ProductManifold with two factors
- */
-template<class Basis, class TargetSpace>
-class MixedLocalGeodesicFEStiffness
- : public Dune::GFE::LocalFirstOrderModel<Basis,TargetSpace>
-{
-public:
- /** \brief Assemble the local stiffness matrix at the current position
- */
- virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
- const typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
- typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeGradient& localGradient,
- typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::MixedHessian& localHessian) const = 0;
-};
-
-#endif
diff --git a/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh b/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
index d62186934244861875b5d72509af1b63adf1fd88..03e481f9d2545db4f088006994a8cd55bae7adf8 100644
--- a/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
+++ b/dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh
@@ -12,14 +12,14 @@
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrix.hh>
-#include <dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh>
+#include <dune/gfe/assemblers/localgeodesicfestiffness.hh>
#include <dune/gfe/spaces/productmanifold.hh>
/** \brief Assembles energy gradient and Hessian with ADOL-C (automatic differentiation)
*/
template<class Basis, class TargetSpace0, class TargetSpace1>
class MixedLocalGFEADOLCStiffness
- : public MixedLocalGeodesicFEStiffness<Basis,Dune::GFE::ProductManifold<TargetSpace0,TargetSpace1> >
+ : public LocalGeodesicFEStiffness<Basis,Dune::GFE::ProductManifold<TargetSpace0,TargetSpace1> >
{
using TargetSpace = Dune::GFE::ProductManifold<TargetSpace0,TargetSpace1>;
@@ -37,8 +37,6 @@ class MixedLocalGFEADOLCStiffness
// some other sizes
constexpr static int gridDim = GridView::dimension;
- using HessianType = typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::MixedHessian;
-
public:
//! Dimension of a tangent space
@@ -74,14 +72,22 @@ public:
DUNE_THROW(Dune::NotImplemented, "!");
}
+ 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,
+ typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::Hessian& localHessian) const override
+ {
+ DUNE_THROW(Dune::NotImplemented, "!");
+ }
+
/** \brief Assemble the local stiffness matrix at the current position
This uses the automatic differentiation toolbox ADOL_C.
*/
virtual void assembleGradientAndHessian(const typename Basis::LocalView& localView,
- const typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
- typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeGradient& localGradient,
- HessianType& localHessian) const override;
+ const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
+ 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 bool adolcScalarMode_;
@@ -156,9 +162,9 @@ 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 typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
- typename Dune::GFE::Impl::MixedLocalStiffnessTypes<TargetSpace>::CompositeGradient& localGradient,
- HessianType& localHessian) const
+ const typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeCoefficients& localConfiguration,
+ typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeGradient& localGradient,
+ typename Dune::GFE::Impl::LocalStiffnessTypes<TargetSpace>::CompositeHessian& localHessian) const
{
int rank = Dune::MPIHelper::getCommunication().rank();