diff --git a/CHANGELOG.md b/CHANGELOG.md
index ccc354c82fb15ad0b4b2a2f6f74c07a5138882bc..f452d9055c0fd7990f1adf7c8824f81c6d5ba3b1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,11 @@
- The class `LocalIntegralEnergy` now does not construct the local interpolation
rule itself anymore. Rather, it expects to be given one in its constructor.
This will allow interpolation rules with state.
+
+- The same holds for `CosseratRodEnergy`. Unfortunately, that move
+ implies that `CosseratRodEnergy` gains an extra template argument
+ which specifies the GFE function used to implement the reference
+ configuration.
- The classes `LocalGeodesicFEFunction` and `LocalProjectedFEFunction`
now take a `dune-functions` basis as their first arguments, instead
diff --git a/dune/gfe/assemblers/cosseratrodenergy.hh b/dune/gfe/assemblers/cosseratrodenergy.hh
index d8dde8cd048f7246a1a0f9900f3228ece1d0d08d..d47fa74a9ec5ed5cc6d042dbb6159e6468f30362 100644
--- a/dune/gfe/assemblers/cosseratrodenergy.hh
+++ b/dune/gfe/assemblers/cosseratrodenergy.hh
@@ -2,6 +2,7 @@
#define DUNE_GFE_COSSERATRODENERGY_HH
#include <array>
+#include <memory>
#include <dune/common/fmatrix.hh>
#include <dune/common/version.hh>
@@ -21,7 +22,7 @@
namespace Dune::GFE
{
- template<class Basis, class LocalInterpolationRule, class RT>
+ template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
class CosseratRodEnergy
: public LocalEnergy<Basis, ProductManifold<RealTuple<RT,3>,Rotation<RT,3> > >
{
@@ -43,6 +44,14 @@ namespace Dune::GFE
public:
+ /** \brief The current configuration whose energy is being computed
+ */
+ const std::shared_ptr<LocalInterpolationRule> localGFEFunction_;
+
+ /** \brief The GFE function that implements the stress-free reference configuration
+ */
+ const std::shared_ptr<ReferenceInterpolationRule> localReferenceConfiguration_;
+
/** \brief The stress-free configuration
The number type cannot be RT, because RT become `adouble` when
@@ -68,22 +77,68 @@ namespace Dune::GFE
GridView gridView_;
//! Constructor
- CosseratRodEnergy(const GridView& gridView,
+ CosseratRodEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
+ std::shared_ptr<ReferenceInterpolationRule> localReferenceConfiguration,
+ const GridView& gridView,
const std::array<double,3>& K, const std::array<double,3>& A)
- : K_(K),
+ : localGFEFunction_(localGFEFunction)
+ , localReferenceConfiguration_(localReferenceConfiguration)
+ , K_(K),
A_(A),
gridView_(gridView)
{}
+ //! Constructor
+ CosseratRodEnergy(LocalInterpolationRule&& localGFEFunction,
+ ReferenceInterpolationRule&& localReferenceConfiguration,
+ const GridView& gridView,
+ const std::array<double,3>& K, const std::array<double,3>& A)
+ : localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction)))
+ , localReferenceConfiguration_(std::make_shared<ReferenceInterpolationRule>(std::move(localReferenceConfiguration)))
+ , K_(K),
+ A_(A),
+ gridView_(gridView)
+ {}
+
+ /** \brief Constructor setting shape constants and material parameters
+ \param A The rod section area
+ \param J1, J2 The geometric moments (Flächenträgheitsmomente)
+ \param E Young's modulus
+ \param nu Poisson number
+ */
+ CosseratRodEnergy(std::shared_ptr<LocalInterpolationRule> localGFEFunction,
+ std::shared_ptr<ReferenceInterpolationRule> localReferenceConfiguration,
+ const GridView& gridView,
+ double A, double J1, double J2, double E, double nu)
+ : localGFEFunction_(localGFEFunction)
+ , localReferenceConfiguration_(localReferenceConfiguration)
+ , gridView_(gridView)
+ {
+ // shear modulus
+ double G = E/(2+2*nu);
+
+ K_[0] = E * J1;
+ K_[1] = E * J2;
+ K_[2] = G * (J1 + J2);
+
+ A_[0] = G * A;
+ A_[1] = G * A;
+ A_[2] = E * A;
+ }
+
/** \brief Constructor setting shape constants and material parameters
\param A The rod section area
\param J1, J2 The geometric moments (Flächenträgheitsmomente)
\param E Young's modulus
\param nu Poisson number
*/
- CosseratRodEnergy(const GridView& gridView,
+ CosseratRodEnergy(LocalInterpolationRule&& localGFEFunction,
+ ReferenceInterpolationRule&& localReferenceConfiguration,
+ const GridView& gridView,
double A, double J1, double J2, double E, double nu)
- : gridView_(gridView)
+ : localGFEFunction_(std::make_shared<LocalInterpolationRule>(std::move(localGFEFunction)))
+ , localReferenceConfiguration_(std::make_shared<ReferenceInterpolationRule>(std::move(localReferenceConfiguration)))
+ , gridView_(gridView)
{
// shear modulus
double G = E/(2+2*nu);
@@ -173,23 +228,20 @@ namespace Dune::GFE
};
- template<class Basis, class LocalInterpolationRule, class RT>
- RT CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+ template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
+ RT CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
energy(const typename Basis::LocalView& localView,
const std::vector<TargetSpace>& localCoefficients) const
{
const auto& localFiniteElement = localView.tree().finiteElement();
- LocalInterpolationRule localConfiguration;
- localConfiguration.bind(localFiniteElement, localCoefficients);
+ localGFEFunction_->bind(localFiniteElement, localCoefficients);
const auto& element = localView.element();
RT energy = 0;
std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > > localReferenceCoefficients = getLocalReferenceConfiguration(localView);
- using InactiveLocalInterpolationRule = typename LocalInterpolationRule::template rebind<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >::other;
- InactiveLocalInterpolationRule localReferenceConfiguration;
- localReferenceConfiguration.bind(localFiniteElement, localReferenceCoefficients);
+ localReferenceConfiguration_->bind(localFiniteElement, localReferenceCoefficients);
// ///////////////////////////////////////////////////////////////////////////////
// The following two loops are a reduced integration scheme. We integrate
@@ -208,10 +260,10 @@ namespace Dune::GFE
double weight = shearingQuad[pt].weight() * integrationElement;
- auto strain = getStrain(localConfiguration, element, quadPos);
+ auto strain = getStrain(*localGFEFunction_, element, quadPos);
// The reference strain
- auto referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
+ auto referenceStrain = getStrain(*localReferenceConfiguration_, element, quadPos);
for (int i=0; i<3; i++)
energy += weight * 0.5 * A_[i] * (strain[i] - referenceStrain[i]) * (strain[i] - referenceStrain[i]);
@@ -228,10 +280,10 @@ namespace Dune::GFE
double weight = bendingQuad[pt].weight() * element.geometry().integrationElement(quadPos);
- auto strain = getStrain(localConfiguration, element, quadPos);
+ auto strain = getStrain(*localGFEFunction_, element, quadPos);
// The reference strain
- auto referenceStrain = getStrain(localReferenceConfiguration, element, quadPos);
+ auto referenceStrain = getStrain(*localReferenceConfiguration_, element, quadPos);
// Part II: the bending and twisting energy
for (int i=0; i<3; i++)
@@ -243,9 +295,9 @@ namespace Dune::GFE
}
- template<class Basis, class LocalInterpolationRule, class RT>
+ template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
template <class ReboundLocalInterpolationRule>
- auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+ auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStrain(const ReboundLocalInterpolationRule& localInterpolation,
const Entity& element,
const FieldVector<double,1>& pos) const
@@ -288,9 +340,9 @@ namespace Dune::GFE
return strain;
}
- template<class Basis, class LocalInterpolationRule, class RT>
+ template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
template <class Number>
- auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+ auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStress(const std::vector<ProductManifold<RealTuple<Number,3>,Rotation<Number,3> > >& localSolution,
const Entity& element,
const FieldVector<double, 1>& pos) const
@@ -311,8 +363,8 @@ namespace Dune::GFE
return stress;
}
- template<class Basis, class LocalInterpolationRule, class RT>
- void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+ template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
+ void CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStrain(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
BlockVector<FieldVector<double, blocksize> >& strain) const
{
@@ -366,8 +418,8 @@ namespace Dune::GFE
}
}
- template<class Basis, class LocalInterpolationRule, class RT>
- void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+ template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
+ void CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getStress(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
BlockVector<FieldVector<double, blocksize> >& stress) const
{
@@ -389,9 +441,9 @@ namespace Dune::GFE
}
}
- template<class Basis, class LocalInterpolationRule, class RT>
+ template<class Basis, class LocalInterpolationRule, class ReferenceInterpolationRule, class RT>
template <class PatchGridView>
- auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
+ auto CosseratRodEnergy<Basis, LocalInterpolationRule, ReferenceInterpolationRule, RT>::
getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol) const
{
diff --git a/src/cosserat-rod.cc b/src/cosserat-rod.cc
index 3df2a0ccc2c2bf1e9e539d2347fe983086dc56ff..646dc9427ed9506aa082705283a6b754922e5100 100644
--- a/src/cosserat-rod.cc
+++ b/src/cosserat-rod.cc
@@ -219,24 +219,36 @@ int main (int argc, char *argv[]) try
// Create the energy and assembler
//////////////////////////////////////////////
+ // GFE function describing the reference configuration
+ using ReferenceGeodesicInterpolationRule = GFE::LocalGeodesicFEFunction<ScalarBasis, TargetSpace>;
+ using ReferenceProjectedInterpolationRule = GFE::LocalProjectedFEFunction<ScalarBasis, TargetSpace>;
+ ReferenceGeodesicInterpolationRule referenceGeodesicFEFunction;
+ ReferenceProjectedInterpolationRule referenceProjectedFEFunction;
+
using ATargetSpace = TargetSpace::rebind<adouble>::other;
using GeodesicInterpolationRule = GFE::LocalGeodesicFEFunction<ScalarBasis, ATargetSpace>;
using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<ScalarBasis, ATargetSpace>;
+ GeodesicInterpolationRule localGeodesicFEFunction;
+ ProjectedInterpolationRule localProjectedFEFunction;
// Assembler using ADOL-C
std::shared_ptr<GFE::LocalEnergy<ScalarBasis,ATargetSpace> > localRodEnergy;
if (parameterSet["interpolationMethod"] == "geodesic")
{
- auto energy = std::make_shared<GFE::CosseratRodEnergy<ScalarBasis, GeodesicInterpolationRule, adouble> >(gridView,
- A, J1, J2, E, nu);
+ auto energy = std::make_shared<GFE::CosseratRodEnergy<ScalarBasis, GeodesicInterpolationRule, ReferenceGeodesicInterpolationRule, adouble> >(std::move(localGeodesicFEFunction),
+ std::move(referenceGeodesicFEFunction),
+ gridView,
+ A, J1, J2, E, nu);
energy->setReferenceConfiguration(referenceConfiguration);
localRodEnergy = energy;
}
else if (parameterSet["interpolationMethod"] == "projected")
{
- auto energy = std::make_shared<GFE::CosseratRodEnergy<ScalarBasis, ProjectedInterpolationRule, adouble> >(gridView,
- A, J1, J2, E, nu);
+ auto energy = std::make_shared<GFE::CosseratRodEnergy<ScalarBasis, ProjectedInterpolationRule, ReferenceProjectedInterpolationRule, adouble> >(std::move(localProjectedFEFunction),
+ std::move(referenceProjectedFEFunction),
+ gridView,
+ A, J1, J2, E, nu);
energy->setReferenceConfiguration(referenceConfiguration);
localRodEnergy = energy;
}
diff --git a/test/cosseratrodenergytest.cc b/test/cosseratrodenergytest.cc
index 72f08f5f9e2ed6fda12f75e59189b267085e4df6..0ef37b2aed8875a28246f6010a7d117b793ea2ad 100644
--- a/test/cosseratrodenergytest.cc
+++ b/test/cosseratrodenergytest.cc
@@ -73,10 +73,15 @@ int main (int argc, char *argv[]) try
using GeodesicInterpolationRule = GFE::LocalGeodesicFEFunction<FEBasis,
TargetSpace>;
+ GeodesicInterpolationRule localGFEFunction;
+ GeodesicInterpolationRule referenceGFEFunction;
GFE::CosseratRodEnergy<FEBasis,
GeodesicInterpolationRule,
- double> localRodEnergy(gridView,
+ GeodesicInterpolationRule,
+ double> localRodEnergy(std::move(localGFEFunction),
+ std::move(referenceGFEFunction),
+ gridView,
1,1,1,1e6,0.3);
SolutionType referenceConfiguration(gridView.size(1));
diff --git a/test/cosseratrodtest.cc b/test/cosseratrodtest.cc
index 51b42cd5ccc97273e438b8ae77c6fad7c6087cb5..630b8e536fb0fc1400eef24dd478c72c86fc1a36 100644
--- a/test/cosseratrodtest.cc
+++ b/test/cosseratrodtest.cc
@@ -134,6 +134,13 @@ int main (int argc, char *argv[]) try
using ATargetSpace = TargetSpace::rebind<adouble>::other;
using GeodesicInterpolationRule = GFE::LocalGeodesicFEFunction<ScalarBasis, ATargetSpace>;
using ProjectedInterpolationRule = GFE::LocalProjectedFEFunction<ScalarBasis, ATargetSpace>;
+ GeodesicInterpolationRule localGeodesicFEFunction;
+ ProjectedInterpolationRule localProjectedFEFunction;
+
+ using ReferenceGeodesicInterpolationRule = GFE::LocalGeodesicFEFunction<ScalarBasis, TargetSpace>;
+ using ReferenceProjectedInterpolationRule = GFE::LocalProjectedFEFunction<ScalarBasis, TargetSpace>;
+ ReferenceGeodesicInterpolationRule localReferenceGeodesicFEFunction;
+ ReferenceProjectedInterpolationRule localReferenceProjectedFEFunction;
// Assembler using ADOL-C
std::shared_ptr<GFE::LocalEnergy<ScalarBasis,ATargetSpace> > localRodEnergy;
@@ -142,15 +149,17 @@ int main (int argc, char *argv[]) try
if (interpolationMethod == "geodesic")
{
- auto energy = std::make_shared<GFE::CosseratRodEnergy<ScalarBasis, GeodesicInterpolationRule, adouble> >(gridView,
- A, J1, J2, E, nu);
+ auto energy = std::make_shared<GFE::CosseratRodEnergy<ScalarBasis, GeodesicInterpolationRule, ReferenceGeodesicInterpolationRule, adouble> >(std::move(localGeodesicFEFunction),
+ std::move(localReferenceGeodesicFEFunction), gridView,
+ A, J1, J2, E, nu);
energy->setReferenceConfiguration(referenceConfiguration);
localRodEnergy = energy;
}
else if (interpolationMethod == "projected")
{
- auto energy = std::make_shared<GFE::CosseratRodEnergy<ScalarBasis, ProjectedInterpolationRule, adouble> >(gridView,
- A, J1, J2, E, nu);
+ auto energy = std::make_shared<GFE::CosseratRodEnergy<ScalarBasis, ProjectedInterpolationRule, ReferenceProjectedInterpolationRule, adouble> >(std::move(localProjectedFEFunction),
+ std::move(localReferenceProjectedFEFunction), gridView,
+ A, J1, J2, E, nu);
energy->setReferenceConfiguration(referenceConfiguration);
localRodEnergy = energy;
}