From 25239378ad947b219eed7bb88eaabfc9a5cf2a15 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Sun, 14 Jan 2024 06:29:12 +0100
Subject: [PATCH] Replace RigidBodyMotion by ProductManifold
The ProductManifold class generalizes RigidBodyMotion, and can do
everything that the RigidBodyMotion class can. Therefore there is
no point in keeping RigidBodyMotion any longer. Having two
implementations for the same thing will just confuse people.
---
CHANGELOG.md | 3 +
.../gfe/assemblers/cosseratenergystiffness.hh | 32 +-
dune/gfe/assemblers/cosseratrodenergy.hh | 56 +--
.../assemblers/geodesicfeassemblerwrapper.hh | 10 +-
dune/gfe/assemblers/localintegralenergy.hh | 5 +-
.../nonplanarcosseratshellenergy.hh | 33 +-
dune/gfe/assemblers/surfacecosseratenergy.hh | 14 +-
dune/gfe/cosseratvtkreader.hh | 10 +-
dune/gfe/cosseratvtkwriter.hh | 34 +-
dune/gfe/localgeodesicfefunction.hh | 32 +-
dune/gfe/localprojectedfefunction.hh | 25 +-
dune/gfe/rodfactory.hh | 32 +-
dune/gfe/spaces/CMakeLists.txt | 1 -
dune/gfe/spaces/rigidbodymotion.hh | 431 ------------------
src/compute-disc-error.cc | 20 +-
src/cosserat-continuum.cc | 37 +-
src/film-on-substrate.cc | 22 +-
src/harmonicmaps.cc | 3 +-
src/rod3d.cc | 34 +-
src/rodobstacle.cc | 9 +-
test/CMakeLists.txt | 3 +-
test/adolctest-scalar-and-vector-mode.cc | 4 +-
test/adolctest.cc | 12 +-
test/cosseratenergytest.cc | 17 +-
test/cosseratrodenergytest.cc | 34 +-
test/cosseratrodtest.cc | 25 +-
test/geodesicfeassemblerwrappertest.cc | 12 +-
test/localgeodesicfefunctiontest.cc | 3 -
test/localprojectedfefunctiontest.cc | 9 -
test/nonplanarcosseratenergytest.cc | 15 +-
test/rigidbodymotiontest.cc | 129 ------
test/targetspacetest.cc | 2 -
test/valuefactory.hh | 15 +-
33 files changed, 302 insertions(+), 821 deletions(-)
delete mode 100644 dune/gfe/spaces/rigidbodymotion.hh
delete mode 100644 test/rigidbodymotiontest.cc
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 52d877d9..b51284fe 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,8 @@
# Master
+- The `RigidBodyMotion` class has been removed. Please use
+ `ProductManifold<RealTuple,Rotation>` from now on.
+
- Building the module requires CMake version 3.16 now, to be in line
with the current core modules.
diff --git a/dune/gfe/assemblers/cosseratenergystiffness.hh b/dune/gfe/assemblers/cosseratenergystiffness.hh
index 3caffc35..08831326 100644
--- a/dune/gfe/assemblers/cosseratenergystiffness.hh
+++ b/dune/gfe/assemblers/cosseratenergystiffness.hh
@@ -23,7 +23,9 @@
#include <dune/gfe/tensor3.hh>
#include <dune/gfe/cosseratstrain.hh>
#include <dune/gfe/spaces/orthogonalmatrix.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
#define DONT_USE_CURL
@@ -69,7 +71,7 @@ public:
template<class Basis, int dim, class field_type=double>
class CosseratEnergyLocalStiffness
- : public Dune::GFE::LocalEnergy<Basis,RigidBodyMotion<field_type,dim> >,
+ : public Dune::GFE::LocalEnergy<Basis,Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > >,
public MixedLocalGeodesicFEStiffness<Basis,
RealTuple<field_type,dim>,
Rotation<field_type,dim> >
@@ -77,7 +79,7 @@ class CosseratEnergyLocalStiffness
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
- typedef RigidBodyMotion<field_type,dim> TargetSpace;
+ typedef Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > TargetSpace;
typedef typename TargetSpace::ctype RT;
typedef typename GridView::template Codim<0>::Entity Entity;
@@ -342,12 +344,14 @@ public:
};
template <class Basis, int dim, class field_type>
-[[deprecated("Use an std::vector<RealTuple<field_type,dim>> and an std::vector<Rotation<field_type,dim>> together with the MixedGFEAssembler and the GFEAssemblerWrapper instead of std::vector<RigidBodyMotion<field_type,dim>>.")]]
+[[deprecated("Use an std::vector<RealTuple<field_type,dim>> and an std::vector<Rotation<field_type,dim>> together with the MixedGFEAssembler and the GFEAssemblerWrapper instead of std::vector<Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> >>.")]]
typename CosseratEnergyLocalStiffness<Basis,dim,field_type>::RT
CosseratEnergyLocalStiffness<Basis,dim,field_type>::
energy(const typename Basis::LocalView& localView,
- const std::vector<RigidBodyMotion<field_type,dim> >& localSolution) const
+ const std::vector<Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > >& localSolution) const
{
+ using namespace Dune::Indices;
+
RT energy = 0;
auto element = localView.element();
@@ -378,7 +382,7 @@ energy(const typename Basis::LocalView& localView,
DT weight = quad[pt].weight() * integrationElement;
// The value of the local function
- RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
+ Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > value = localGeodesicFEFunction.evaluate(quadPos);
// The derivative of the local function defined on the reference element
typename LocalGFEFunctionType::DerivativeType referenceDerivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
@@ -396,7 +400,7 @@ energy(const typename Basis::LocalView& localView,
//
Dune::FieldMatrix<field_type,dim,dim> R;
- value.q.matrix(R);
+ value[_1].matrix(R);
Dune::GFE::CosseratStrain<field_type,dim,gridDim> U(derivative,R);
@@ -405,7 +409,13 @@ energy(const typename Basis::LocalView& localView,
// Note: we need it in matrix coordinates
//////////////////////////////////////////////////////////
- Tensor3<field_type,3,3,gridDim> DR = value.quaternionTangentToMatrixTangent(derivative);
+ // Extract the orientation derivative
+ Dune::FieldMatrix<field_type,4,gridDim> orientationDerivative;
+ for (size_t i=0; i<4; ++i)
+ for (size_t j=0; j<gridDim; ++j)
+ orientationDerivative[i][j] = derivative[3+i][j];
+
+ Tensor3<field_type,3,3,gridDim> DR = value[_1].quaternionTangentToMatrixTangent(orientationDerivative);
// Add the local energy density
if (gridDim==2) {
@@ -440,7 +450,7 @@ energy(const typename Basis::LocalView& localView,
// Only translational dofs are affected by the volume load
for (size_t i=0; i<volumeLoadDensity.size(); i++)
- energy += (volumeLoadDensity[i] * value.r[i]) * quad[pt].weight() * integrationElement;
+ energy += (volumeLoadDensity[i] * value[_0].globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
}
@@ -467,14 +477,14 @@ energy(const typename Basis::LocalView& localView,
const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
// The value of the local function
- RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
+ Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > value = localGeodesicFEFunction.evaluate(quadPos);
// Value of the Neumann data at the current position
auto neumannValue = neumannFunction_(it.geometry().global(quad[pt].position()));
// Only translational dofs are affected by the Neumann force
for (size_t i=0; i<neumannValue.size(); i++)
- energy += (neumannValue[i] * value.r[i]) * quad[pt].weight() * integrationElement;
+ energy += (neumannValue[i] * value[_0].globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
}
diff --git a/dune/gfe/assemblers/cosseratrodenergy.hh b/dune/gfe/assemblers/cosseratrodenergy.hh
index b2c6dced..5a336135 100644
--- a/dune/gfe/assemblers/cosseratrodenergy.hh
+++ b/dune/gfe/assemblers/cosseratrodenergy.hh
@@ -17,15 +17,17 @@
#include <dune/fufem/boundarypatch.hh>
#include <dune/gfe/assemblers/localenergy.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
namespace Dune::GFE {
template<class Basis, class LocalInterpolationRule, class RT>
class CosseratRodEnergy
- : public LocalEnergy<Basis, RigidBodyMotion<RT,3> >
+ : public LocalEnergy<Basis, ProductManifold<RealTuple<RT,3>,Rotation<RT,3> > >
{
- typedef RigidBodyMotion<RT,3> TargetSpace;
+ using TargetSpace = ProductManifold<RealTuple<RT,3>,Rotation<RT,3> >;
// grid types
using GridView = typename Basis::GridView;
@@ -50,7 +52,7 @@ namespace Dune::GFE {
The referenceConfiguration is not a variable, and we don't
want to use `adouble` for it.
*/
- std::vector<RigidBodyMotion<double,3> > referenceConfiguration_;
+ std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > > referenceConfiguration_;
public:
@@ -99,13 +101,13 @@ namespace Dune::GFE {
/** \brief Set the stress-free configuration
*/
- void setReferenceConfiguration(const std::vector<RigidBodyMotion<double,3> >& referenceConfiguration) {
+ void setReferenceConfiguration(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& referenceConfiguration) {
referenceConfiguration_ = referenceConfiguration;
}
/** \brief Compute local element energy */
virtual RT energy (const typename Basis::LocalView& localView,
- const std::vector<RigidBodyMotion<RT,3> >& localSolution) const override;
+ const std::vector<TargetSpace>& localSolution) const override;
/** \brief Get the rod strain at one point in the rod
*
@@ -121,16 +123,16 @@ namespace Dune::GFE {
* \tparam Number This is a member template because the method has to work for double and adouble
*/
template<class Number>
- auto getStress(const std::vector<RigidBodyMotion<Number,3> >& localSolution,
+ auto getStress(const std::vector<ProductManifold<RealTuple<Number,3>,Rotation<Number,3> > >& localSolution,
const Entity& element,
const FieldVector<double,1>& pos) const;
/** \brief Get average strain for each element */
- void getStrain(const std::vector<RigidBodyMotion<double,3> >& sol,
+ void getStrain(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
BlockVector<FieldVector<double, blocksize> >& strain) const;
/** \brief Get average stress for each element */
- void getStress(const std::vector<RigidBodyMotion<double,3> >& sol,
+ void getStress(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
BlockVector<FieldVector<double, blocksize> >& stress) const;
/** \brief Return resultant force across boundary in canonical coordinates
@@ -138,14 +140,14 @@ namespace Dune::GFE {
\note Linear run-time in the size of the grid */
template <class PatchGridView>
auto getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
- const std::vector<RigidBodyMotion<double,3> >& sol) const;
+ const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol) const;
protected:
- std::vector<RigidBodyMotion<double,3> > getLocalReferenceConfiguration(const typename Basis::LocalView& localView) const
+ std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > > getLocalReferenceConfiguration(const typename Basis::LocalView& localView) const
{
unsigned int numOfBaseFct = localView.size();
- std::vector<RigidBodyMotion<double,3> > localReferenceConfiguration(numOfBaseFct);
+ std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > > localReferenceConfiguration(numOfBaseFct);
for (size_t i=0; i<numOfBaseFct; i++)
localReferenceConfiguration[i] = referenceConfiguration_[localView.index(i)];
@@ -170,7 +172,7 @@ namespace Dune::GFE {
template<class Basis, class LocalInterpolationRule, class RT>
RT CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
energy(const typename Basis::LocalView& localView,
- const std::vector<RigidBodyMotion<RT,3> >& localCoefficients) const
+ const std::vector<TargetSpace>& localCoefficients) const
{
const auto& localFiniteElement = localView.tree().finiteElement();
LocalInterpolationRule localConfiguration(localFiniteElement, localCoefficients);
@@ -179,8 +181,8 @@ namespace Dune::GFE {
RT energy = 0;
- std::vector<RigidBodyMotion<double,3> > localReferenceCoefficients = getLocalReferenceConfiguration(localView);
- using InactiveLocalInterpolationRule = typename LocalInterpolationRule::template rebind<RigidBodyMotion<double,3> >::other;
+ 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(localFiniteElement, localReferenceCoefficients);
// ///////////////////////////////////////////////////////////////////////////////
@@ -268,15 +270,17 @@ namespace Dune::GFE {
// /////////////////////////////////////////////
// Strain defined on each element
+ using namespace Dune::Indices;
+
// Part I: the shearing and stretching strain
FieldVector<Number, 6> strain(0);
- strain[0] = r_s * value.q.director(0); // shear strain
- strain[1] = r_s * value.q.director(1); // shear strain
- strain[2] = r_s * value.q.director(2); // stretching strain
+ strain[0] = r_s * value[_1].director(0); // shear strain
+ strain[1] = r_s * value[_1].director(1); // shear strain
+ strain[2] = r_s * value[_1].director(2); // stretching strain
// Part II: the Darboux vector
- FieldVector<Number,3> u = darboux<Number>(value.q, q_s);
+ FieldVector<Number,3> u = darboux<Number>(value[_1], q_s);
strain[3] = u[0];
strain[4] = u[1];
strain[5] = u[2];
@@ -287,7 +291,7 @@ namespace Dune::GFE {
template<class Basis, class LocalInterpolationRule, class RT>
template <class Number>
auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
- getStress(const std::vector<RigidBodyMotion<Number,3> >& localSolution,
+ getStress(const std::vector<ProductManifold<RealTuple<Number,3>,Rotation<Number,3> > >& localSolution,
const Entity& element,
const FieldVector<double, 1>& pos) const
{
@@ -309,7 +313,7 @@ namespace Dune::GFE {
template<class Basis, class LocalInterpolationRule, class RT>
void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
- getStrain(const std::vector<RigidBodyMotion<double,3> >& sol,
+ getStrain(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
BlockVector<FieldVector<double, blocksize> >& strain) const
{
const typename GridView::Traits::IndexSet& indexSet = this->basis_.gridView().indexSet();
@@ -330,7 +334,7 @@ namespace Dune::GFE {
Dune::LagrangeSimplexLocalFiniteElement<double, double, 1, 1> localFiniteElement;
int numOfBaseFct = localFiniteElement.localCoefficients().size();
- std::vector<RigidBodyMotion<double,3> > localSolution(2);
+ std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > > localSolution(2);
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[indexSet.subIndex(element,i,1)];
@@ -364,7 +368,7 @@ namespace Dune::GFE {
template<class Basis, class LocalInterpolationRule, class RT>
void CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
- getStress(const std::vector<RigidBodyMotion<double,3> >& sol,
+ getStress(const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol,
BlockVector<FieldVector<double, blocksize> >& stress) const
{
// Get the strain
@@ -389,7 +393,7 @@ namespace Dune::GFE {
template <class PatchGridView>
auto CosseratRodEnergy<Basis, LocalInterpolationRule, RT>::
getResultantForce(const BoundaryPatch<PatchGridView>& boundary,
- const std::vector<RigidBodyMotion<double,3> >& sol) const
+ const std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& sol) const
{
const typename GridView::Traits::IndexSet& indexSet = this->basis_.gridView().indexSet();
@@ -408,11 +412,11 @@ namespace Dune::GFE {
double pos = facet.geometryInInside().corner(0);
- std::vector<RigidBodyMotion<double,3> > localSolution(2);
+ std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > > localSolution(2);
localSolution[0] = sol[indexSet.subIndex(*facet.inside(),0,1)];
localSolution[1] = sol[indexSet.subIndex(*facet.inside(),1,1)];
- std::vector<RigidBodyMotion<double,3> > localRefConf(2);
+ std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > > localRefConf(2);
localRefConf[0] = referenceConfiguration_[indexSet.subIndex(*facet.inside(),0,1)];
localRefConf[1] = referenceConfiguration_[indexSet.subIndex(*facet.inside(),1,1)];
diff --git a/dune/gfe/assemblers/geodesicfeassemblerwrapper.hh b/dune/gfe/assemblers/geodesicfeassemblerwrapper.hh
index e6a3fdda..12d2b789 100644
--- a/dune/gfe/assemblers/geodesicfeassemblerwrapper.hh
+++ b/dune/gfe/assemblers/geodesicfeassemblerwrapper.hh
@@ -4,7 +4,6 @@
#include <dune/gfe/assemblers/mixedgfeassembler.hh>
#include <dune/common/tuplevector.hh>
#include <dune/gfe/spaces/productmanifold.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
namespace Dune::GFE {
@@ -87,13 +86,8 @@ splitVector(const std::vector<TargetSpace>& sol) const {
solutionSplit[_0].resize(n);
solutionSplit[_1].resize(n);
for (std::size_t i = 0; i < n; i++) {
- if constexpr(std::is_base_of<TargetSpace,RigidBodyMotion<double, 3> >::value) {
- solutionSplit[_0][i] = sol[i].r; // Deformation part
- solutionSplit[_1][i] = sol[i].q; // Rotational part
- } else if constexpr(std::is_base_of<TargetSpace,ProductManifold<RealTuple<double,3>,UnitVector<double,3> > >::value) {
- solutionSplit[_0][i] = sol[i][_0]; // Deformation part
- solutionSplit[_1][i] = sol[i][_1]; // Rotational part
- }
+ solutionSplit[_0][i] = sol[i][_0]; // Deformation part
+ solutionSplit[_1][i] = sol[i][_1]; // Rotational part
}
return solutionSplit;
}
diff --git a/dune/gfe/assemblers/localintegralenergy.hh b/dune/gfe/assemblers/localintegralenergy.hh
index 40fb241b..23c8dd8d 100644
--- a/dune/gfe/assemblers/localintegralenergy.hh
+++ b/dune/gfe/assemblers/localintegralenergy.hh
@@ -10,7 +10,7 @@
#include <dune/gfe/cosseratstrain.hh>
#include <dune/gfe/densities/localdensity.hh>
#include <dune/gfe/spaces/realtuple.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/rotation.hh>
#ifdef PROJECTED_INTERPOLATION
#include <dune/gfe/localprojectedfefunction.hh>
#else
@@ -65,9 +65,6 @@ namespace Dune::GFE {
using TargetSpaceDeformation = typename std::tuple_element<0, std::tuple<TargetSpaces...> >::type;
using TargetSpaceRotation = typename std::tuple_element<1, std::tuple<TargetSpaces...> >::type;
- static_assert( (std::is_same<TargetSpaceDeformation, RigidBodyMotion<RT,gridDim> >::value)
- or (std::is_same<TargetSpaceDeformation, RealTuple<RT,gridDim> >::value), "The first TargetSpace of LocalGeodesicIntegralEnergy needs to be RigidBodyMotion or RealTuple!" );
-
const std::vector<TargetSpaceDeformation>& localDeformationConfiguration = std::get<0>(std::forward_as_tuple(localSolutions ...));
const std::vector<TargetSpaceRotation>& localOrientationConfiguration = std::get<1>(std::forward_as_tuple(localSolutions ...));
diff --git a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
index fd1a9eb0..98d941b5 100644
--- a/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
+++ b/dune/gfe/assemblers/nonplanarcosseratshellenergy.hh
@@ -21,8 +21,9 @@
#include <dune/gfe/tensor3.hh>
#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
-#include <dune/gfe/spaces/unitvector.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
#if HAVE_DUNE_CURVEDGEOMETRY
#include <dune/curvedgeometry/curvedgeometry.hh>
@@ -38,7 +39,7 @@
*/
template<class Basis, int dim, class field_type, class StressFreeStateGridFunction>
class NonplanarCosseratShellEnergy
- : public Dune::GFE::LocalEnergy<Basis,RigidBodyMotion<field_type,dim> >,
+ : public Dune::GFE::LocalEnergy<Basis,Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > >,
public MixedLocalGeodesicFEStiffness<Basis,
RealTuple<field_type,dim>,
Rotation<field_type,dim> >
@@ -46,7 +47,7 @@ class NonplanarCosseratShellEnergy
// grid types
typedef typename Basis::GridView GridView;
typedef typename GridView::ctype DT;
- typedef RigidBodyMotion<field_type,dim> TargetSpace;
+ typedef Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > TargetSpace;
typedef typename TargetSpace::ctype RT;
typedef typename GridView::template Codim<0>::Entity Entity;
@@ -191,12 +192,14 @@ public:
};
template <class Basis, int dim, class field_type, class StressFreeStateGridFunction>
-[[deprecated("Use an std::vector<RealTuple<field_type,dim>> and an std::vector<Rotation<field_type,dim>> together with the MixedGFEAssembler and the GFEAssemblerWrapper instead of std::vector<RigidBodyMotion<field_type,dim>>.")]]
+[[deprecated("Use an std::vector<RealTuple<field_type,dim>> and an std::vector<Rotation<field_type,dim>> together with the MixedGFEAssembler and the GFEAssemblerWrapper instead of std::vector<Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> >>.")]]
typename NonplanarCosseratShellEnergy<Basis, dim, field_type, StressFreeStateGridFunction>::RT
NonplanarCosseratShellEnergy<Basis,dim,field_type, StressFreeStateGridFunction>::
energy(const typename Basis::LocalView& localView,
- const std::vector<RigidBodyMotion<field_type,dim> >& localSolution) const
+ const std::vector<Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > >& localSolution) const
{
+ using namespace Dune::Indices;
+
// The element geometry
auto element = localView.element();
@@ -240,7 +243,7 @@ energy(const typename Basis::LocalView& localView,
const DT integrationElement = geometry.integrationElement(quadPos);
// The value of the local function
- RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
+ Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > value = localGeodesicFEFunction.evaluate(quadPos);
// The derivative of the local function w.r.t. the coordinate system of the tangent space
auto derivative = localGeodesicFEFunction.evaluateDerivative(quadPos,value);
@@ -251,11 +254,17 @@ energy(const typename Basis::LocalView& localView,
//////////////////////////////////////////////////////////
Dune::FieldMatrix<field_type,dim,dim> R;
- value.q.matrix(R);
+ value[_1].matrix(R);
auto RT = Dune::GFE::transpose(R);
+ // Extract the orientation derivative
+ Dune::FieldMatrix<field_type,4,gridDim> orientationDerivative;
+ for (size_t i=0; i<4; ++i)
+ for (size_t j=0; j<gridDim; ++j)
+ orientationDerivative[i][j] = derivative[3+i][j];
+
//Derivative of the rotation w.r.t. the coordinate system of the tangent space
- Tensor3<field_type,3,3,gridDim> DR = value.quaternionTangentToMatrixTangent(derivative);
+ Tensor3<field_type,3,3,gridDim> DR = value[_1].quaternionTangentToMatrixTangent(orientationDerivative);
//////////////////////////////////////////////////////////
// Fundamental forms and curvature
@@ -389,7 +398,7 @@ energy(const typename Basis::LocalView& localView,
// Only translational dofs are affected by the volume load
for (size_t i=0; i<volumeLoadDensity.size(); i++)
- energy += (volumeLoadDensity[i] * value.r[i]) * quad[pt].weight() * integrationElement;
+ energy += (volumeLoadDensity[i] * value[_0].globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
}
@@ -415,14 +424,14 @@ energy(const typename Basis::LocalView& localView,
const DT integrationElement = it.geometry().integrationElement(quad[pt].position());
// The value of the local function
- RigidBodyMotion<field_type,dim> value = localGeodesicFEFunction.evaluate(quadPos);
+ Dune::GFE::ProductManifold<RealTuple<field_type,dim>,Rotation<field_type,dim> > value = localGeodesicFEFunction.evaluate(quadPos);
// Value of the Neumann data at the current position
Dune::FieldVector<double,3> neumannValue = neumannFunction_(it.geometry().global(quad[pt].position()));
// Only translational dofs are affected by the Neumann force
for (size_t i=0; i<neumannValue.size(); i++)
- energy += (neumannValue[i] * value.r[i]) * quad[pt].weight() * integrationElement;
+ energy += (neumannValue[i] * value[_0].globalCoordinates()[i]) * quad[pt].weight() * integrationElement;
}
}
diff --git a/dune/gfe/assemblers/surfacecosseratenergy.hh b/dune/gfe/assemblers/surfacecosseratenergy.hh
index b5d6c16b..e70d8415 100644
--- a/dune/gfe/assemblers/surfacecosseratenergy.hh
+++ b/dune/gfe/assemblers/surfacecosseratenergy.hh
@@ -13,7 +13,9 @@
#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/assemblers/mixedlocalgeodesicfestiffness.hh>
#include <dune/gfe/tensor3.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
#include <dune/curvedgeometry/curvedgeometry.hh>
#include <dune/localfunctions/lagrange/lfecache.hh>
@@ -35,7 +37,7 @@ namespace Dune::GFE {
using Entity = typename GridView::template Codim<0>::Entity ;
using RBM0 = RealTuple<RT,GridView::dimensionworld> ;
using RBM1 = Rotation<RT,GridView::dimensionworld> ;
- using RBM = RigidBodyMotion<RT,GridView::dimensionworld> ;
+ using RBM = ProductManifold<RBM0,RBM1> ;
constexpr static int dimWorld = GridView::dimensionworld;
constexpr static int gridDim = GridView::dimension;
@@ -60,7 +62,7 @@ namespace Dune::GFE {
//
// So, DR[i][j][k] contains \partial R_ij / \partial k
Tensor3<RT, dimWorld, dimWorld, 4> dd_dq;
- value.q.getFirstDerivativesOfDirectors(dd_dq);
+ value[Indices::_1].getFirstDerivativesOfDirectors(dd_dq);
derivative_rotation = RT(0);
for (int i=0; i<dimWorld; i++)
@@ -144,8 +146,8 @@ namespace Dune::GFE {
std::vector<RBM> localSolutionRBM(localSolution0.size());
for (int i = 0; i < localSolution0.size(); i++) {
for (int j = 0; j < dimWorld; j++)
- localSolutionRBM[i].r[j] = localSolution0[i][j];
- localSolutionRBM[i].q = localSolution1[i];
+ localSolutionRBM[i][_0][j] = localSolution0[i][j];
+ localSolutionRBM[i][_1] = localSolution1[i];
}
typedef LocalGeodesicFEFunction<gridDim, DT, decltype(deformationLocalFiniteElement), RBM> LocalGFEFunctionType;
LocalGFEFunctionType localGeodesicFEFunction(deformationLocalFiniteElement,localSolutionRBM);
@@ -232,7 +234,7 @@ namespace Dune::GFE {
//////////////////////////////////////////////////////////
Dune::FieldMatrix<RT,dimWorld,dimWorld> R;
- value.q.matrix(R);
+ value[_1].matrix(R);
auto rt = Dune::GFE::transpose(R);
Tensor3<RT,dimWorld,dimWorld,boundaryDim> derivative_rotation;
diff --git a/dune/gfe/cosseratvtkreader.hh b/dune/gfe/cosseratvtkreader.hh
index 7e0e6109..659634ff 100644
--- a/dune/gfe/cosseratvtkreader.hh
+++ b/dune/gfe/cosseratvtkreader.hh
@@ -1,7 +1,9 @@
#ifndef DUNE_GFE_COSSERATVTKREADER_HH
#define DUNE_GFE_COSSERATVTKREADER_HH
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
namespace Dune
{
@@ -13,7 +15,7 @@ namespace Dune
{
public:
- static void read(std::vector<RigidBodyMotion<double,3> >& configuration,
+ static void read(std::vector<ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& configuration,
const std::string& filename)
{
VTKFile vtkFile;
@@ -23,14 +25,14 @@ namespace Dune
for (size_t i=0; i<configuration.size(); i++)
{
- configuration[i].r = vtkFile.points_[i];
+ configuration[i][Indices::_0].globalCoordinates() = vtkFile.points_[i];
FieldMatrix<double,3,3> R;
for (int j=0; j<3; j++)
for (int k=0; k<3; k++)
R[j][k] = vtkFile.directors_[k][i][j];
- configuration[i].q.set(R);
+ configuration[i][Indices::_1].set(R);
}
}
diff --git a/dune/gfe/cosseratvtkwriter.hh b/dune/gfe/cosseratvtkwriter.hh
index feb41c41..1bd51347 100644
--- a/dune/gfe/cosseratvtkwriter.hh
+++ b/dune/gfe/cosseratvtkwriter.hh
@@ -11,7 +11,9 @@
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/gfe/vtkfile.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
/** \brief Write the configuration of a Cosserat material in VTK format */
@@ -22,8 +24,8 @@ class CosseratVTKWriter
static const int dim = GridType::dimension;
template <typename Basis1, typename Basis2>
- static void downsample(const Basis1& basis1, const std::vector<RigidBodyMotion<double,3> >& v1,
- const Basis2& basis2, std::vector<RigidBodyMotion<double,3> >& v2)
+ static void downsample(const Basis1& basis1, const std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& v1,
+ const Basis2& basis2, std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& v2)
{
// Embed v1 into R^7
std::vector<Dune::FieldVector<double,7> > v1Embedded(v1.size());
@@ -35,10 +37,10 @@ class CosseratVTKWriter
std::vector<Dune::FieldVector<double,7> > v2Embedded;
Dune::Functions::interpolate(basis2, v2Embedded, function);
- // Copy back from R^7 into RigidBodyMotions
+ // Copy back from R^7 into ProductManifold
v2.resize(v2Embedded.size());
for (size_t i=0; i<v2.size(); i++)
- v2[i] = RigidBodyMotion<double,3>(v2Embedded[i]);
+ v2[i] = Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >(v2Embedded[i]);
}
template <typename Basis1, typename Basis2>
@@ -118,11 +120,11 @@ public:
const std::string& filename)
{
using namespace Dune::TypeTree::Indices;
- std::vector<RigidBodyMotion<double,3> > xRBM(basis.size());
+ std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > > xRBM(basis.size());
for (std::size_t i = 0; i < basis.size(); i++) {
for (int j = 0; j < 3; j ++) // Displacement part
- xRBM[i].r[j] = configuration[_0][i][j];
- xRBM[i].q = configuration[_1][i]; // Rotation part
+ xRBM[i][_0].globalCoordinates()[j] = configuration[_0][i][j];
+ xRBM[i][_1] = configuration[_1][i]; // Rotation part
}
write(basis,xRBM,filename);
}
@@ -134,10 +136,10 @@ public:
const std::string& filename)
{
using namespace Dune::TypeTree::Indices;
- std::vector<RigidBodyMotion<double,3> > xRBM(basis.size());
+ std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > > xRBM(basis.size());
for (std::size_t i = 0; i < basis.size(); i++) {
for (int j = 0; j < 3; j ++) // Displacement part
- xRBM[i].r[j] = configuration[i][j];
+ xRBM[i][_0].globalCoordinates()[j] = configuration[i][j];
}
write(basis,xRBM,filename);
}
@@ -146,9 +148,11 @@ public:
*/
template <typename Basis>
static void write(const Basis& basis,
- const std::vector<RigidBodyMotion<double,3> >& configuration,
+ const std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& configuration,
const std::string& filename)
{
+ using namespace Dune::Indices;
+
assert(basis.size() == configuration.size());
auto gridView = basis.gridView();
@@ -173,7 +177,7 @@ public:
blockedInterleaved()
));
- std::vector<RigidBodyMotion<double,3> > downsampledConfig;
+ std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > > downsampledConfig;
downsample(basis, configuration, blockedP2Basis, downsampledConfig);
@@ -198,7 +202,7 @@ public:
// Enter vertex coordinates
std::vector<Dune::FieldVector<double, 3> > points(configuration.size());
for (size_t i=0; i<configuration.size(); i++)
- points[i] = configuration[i].r;
+ points[i] = configuration[i][_0].globalCoordinates();
vtkFile.points_ = points;
@@ -363,7 +367,7 @@ public:
// Z coordinate for better visualization of wrinkles
std::vector<double> zCoord(points.size());
for (size_t i=0; i<configuration.size(); i++)
- zCoord[i] = configuration[i].r[2];
+ zCoord[i] = configuration[i][_0].globalCoordinates()[2];
vtkFile.zCoord_ = zCoord;
@@ -372,7 +376,7 @@ public:
{
vtkFile.directors_[i].resize(configuration.size());
for (size_t j=0; j<configuration.size(); j++)
- vtkFile.directors_[i][j] = configuration[j].q.director(i);
+ vtkFile.directors_[i][j] = configuration[j][_1].director(i);
}
// Actually write the VTK file to disk
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 76e0d864..6892e936 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -10,7 +10,9 @@
#include <dune/gfe/averagedistanceassembler.hh>
#include <dune/gfe/targetspacertrsolver.hh>
#include <dune/gfe/localquickanddirtyfefunction.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
#include <dune/gfe/tensor3.hh>
#include <dune/gfe/tensorssd.hh>
@@ -551,18 +553,17 @@ evaluateFDDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>
/** \brief A function defined by simplicial geodesic interpolation
- from the reference element to a RigidBodyMotion.
+ from the reference element to a ProductManifold<RealTuple,Rotation>.
This is a specialization for speeding up the code.
- We use that a RigidBodyMotion is a product manifold.
\tparam dim Dimension of the reference element
\tparam ctype Type used for coordinates on the reference element
*/
template <int dim, class ctype, class LocalFiniteElement, class field_type>
-class LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,RigidBodyMotion<field_type,3> >
+class LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement, Dune::GFE::ProductManifold<RealTuple<field_type,3>, Rotation<field_type,3> > >
{
- typedef RigidBodyMotion<field_type,3> TargetSpace;
+ using TargetSpace = Dune::GFE::ProductManifold<RealTuple<field_type,3>, Rotation<field_type,3> >;
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
static const int embeddedDim = EmbeddedTangentVector::dimension;
@@ -583,14 +584,15 @@ public:
: localFiniteElement_(localFiniteElement),
translationCoefficients_(coefficients.size())
{
+ using namespace Dune::Indices;
assert(localFiniteElement.localBasis().size() == coefficients.size());
for (size_t i=0; i<coefficients.size(); i++)
- translationCoefficients_[i] = coefficients[i].r;
+ translationCoefficients_[i] = coefficients[i][_0].globalCoordinates();
std::vector<Rotation<field_type,3> > orientationCoefficients(coefficients.size());
for (size_t i=0; i<coefficients.size(); i++)
- orientationCoefficients[i] = coefficients[i].q;
+ orientationCoefficients[i] = coefficients[i][_1];
orientationFEFunction_ = std::unique_ptr<LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3> > > (new LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3> >(localFiniteElement,orientationCoefficients));
@@ -618,17 +620,19 @@ public:
/** \brief Evaluate the function */
TargetSpace evaluate(const Dune::FieldVector<ctype, dim>& local) const
{
+ using namespace Dune::Indices;
+
TargetSpace result;
// Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
std::vector<Dune::FieldVector<ctype,1> > w;
localFiniteElement_.localBasis().evaluateFunction(local,w);
- result.r = 0;
+ result[_0] = Dune::FieldVector<field_type,3>(0.0);
for (size_t i=0; i<w.size(); i++)
- result.r.axpy(w[i][0], translationCoefficients_[i]);
+ result[_0].globalCoordinates().axpy(w[i][0], translationCoefficients_[i]);
- result.q = orientationFEFunction_->evaluate(local);
+ result[_1] = orientationFEFunction_->evaluate(local);
return result;
}
@@ -661,6 +665,8 @@ public:
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
const TargetSpace& q) const
{
+ using namespace Dune::Indices;
+
DerivativeType result(0);
// get translation part
@@ -672,7 +678,7 @@ public:
result[j].axpy(translationCoefficients_[i][j], sfDer[i][0]);
// get orientation part
- Dune::FieldMatrix<field_type,4,dim> qResult = orientationFEFunction_->evaluateDerivative(local,q.q);
+ Dune::FieldMatrix<field_type,4,dim> qResult = orientationFEFunction_->evaluateDerivative(local,q[_1]);
for (int i=0; i<4; i++)
for (int j=0; j<dim; j++)
result[3+i][j] = qResult[i][j];
@@ -777,9 +783,7 @@ private:
*/
const LocalFiniteElement& localFiniteElement_;
- // The two factors of a RigidBodyMotion
- //LocalGeodesicFEFunction<dim,ctype,RealTuple<3> > translationFEFunction_;
-
+ // Coefficients for the two factors of the product manifold
std::vector<Dune::FieldVector<field_type,3> > translationCoefficients_;
std::unique_ptr<LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3> > > orientationFEFunction_;
diff --git a/dune/gfe/localprojectedfefunction.hh b/dune/gfe/localprojectedfefunction.hh
index 46acb95a..5f1a0c49 100644
--- a/dune/gfe/localprojectedfefunction.hh
+++ b/dune/gfe/localprojectedfefunction.hh
@@ -11,7 +11,7 @@
#include <dune/gfe/linearalgebra.hh>
#include <dune/gfe/polardecomposition.hh>
#include <dune/gfe/spaces/realtuple.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
#include <dune/gfe/spaces/rotation.hh>
namespace Dune {
@@ -414,10 +414,10 @@ namespace Dune {
* \tparam LocalFiniteElement A Lagrangian finite element whose shape functions define the interpolation weights
*/
template <int dim, class ctype, class LocalFiniteElement, class field_type>
- class LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,RigidBodyMotion<field_type,3> >
+ class LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> > >
{
public:
- typedef RigidBodyMotion<field_type,3> TargetSpace;
+ using TargetSpace = ProductManifold<RealTuple<field_type,3>,Rotation<field_type,3> >;
private:
typedef typename TargetSpace::ctype RT;
@@ -441,13 +441,14 @@ namespace Dune {
translationCoefficients_(coefficients.size())
{
assert(localFiniteElement.localBasis().size() == coefficients.size());
+ using namespace Dune::Indices;
for (size_t i=0; i<coefficients.size(); i++)
- translationCoefficients_[i] = coefficients[i].r;
+ translationCoefficients_[i] = coefficients[i][_0].globalCoordinates();
std::vector<Rotation<field_type,3> > orientationCoefficients(coefficients.size());
for (size_t i=0; i<coefficients.size(); i++)
- orientationCoefficients[i] = coefficients[i].q;
+ orientationCoefficients[i] = coefficients[i][_1];
orientationFunction_ = std::make_unique<LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,Rotation<field_type,3> > > (localFiniteElement,orientationCoefficients);
}
@@ -474,17 +475,19 @@ namespace Dune {
/** \brief Evaluate the function */
TargetSpace evaluate(const Dune::FieldVector<ctype, dim>& local) const
{
- RigidBodyMotion<field_type,3> result;
+ using namespace Dune::Indices;
+
+ TargetSpace result;
// Evaluate the weighting factors---these are the Lagrangian shape function values at 'local'
std::vector<Dune::FieldVector<ctype,1> > w;
localFiniteElement_.localBasis().evaluateFunction(local,w);
- result.r = 0;
+ result[_0] = FieldVector<field_type,3>(0.0);
for (size_t i=0; i<w.size(); i++)
- result.r.axpy(w[i][0], translationCoefficients_[i]);
+ result[_0].globalCoordinates().axpy(w[i][0], translationCoefficients_[i]);
- result.q = orientationFunction_->evaluate(local);
+ result[_1] = orientationFunction_->evaluate(local);
return result;
}
@@ -506,6 +509,8 @@ namespace Dune {
DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
const TargetSpace& q) const
{
+ using namespace Dune::Indices;
+
DerivativeType result(0);
// get translation part
@@ -517,7 +522,7 @@ namespace Dune {
result[j].axpy(translationCoefficients_[i][j], sfDer[i][0]);
// get orientation part
- Dune::FieldMatrix<field_type,4,dim> qResult = orientationFunction_->evaluateDerivative(local,q.q);
+ Dune::FieldMatrix<field_type,4,dim> qResult = orientationFunction_->evaluateDerivative(local,q[_1]);
for (int i=0; i<4; i++)
for (int j=0; j<dim; j++)
diff --git a/dune/gfe/rodfactory.hh b/dune/gfe/rodfactory.hh
index 3f27c7d5..bea34db0 100644
--- a/dune/gfe/rodfactory.hh
+++ b/dune/gfe/rodfactory.hh
@@ -7,11 +7,13 @@
#include <dune/matrix-vector/crossproduct.hh>
-#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/rodassembler.hh>
#include <dune/gfe/riemanniantrsolver.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
/** \brief A factory class that implements various ways to create rod configurations
*/
@@ -33,7 +35,7 @@ public:
\param[in] n The number of vertices
*/
template <int dim>
- void create(std::vector<RigidBodyMotion<double,dim> >& rod,
+ void create(std::vector<Dune::GFE::ProductManifold<RealTuple<double,dim>,Rotation<double,dim> > >& rod,
const Dune::FieldVector<double,3>& beginning, const Dune::FieldVector<double,3>& end)
{
// Compute the correct orientation
@@ -54,7 +56,9 @@ public:
orientation = Rotation<double,3>(axis, angle);
// Set the values
- create(rod, RigidBodyMotion<double,dim>(beginning,orientation), RigidBodyMotion<double,dim>(end,orientation));
+ create(rod,
+ Dune::GFE::ProductManifold<RealTuple<double,dim>,Rotation<double,dim> > (beginning,orientation),
+ Dune::GFE::ProductManifold<RealTuple<double,dim>,Rotation<double,dim> > (end,orientation));
}
@@ -63,9 +67,9 @@ public:
\param[out] rod The new rod
*/
template <int spaceDim>
- void create(std::vector<RigidBodyMotion<double,spaceDim> >& rod,
- const RigidBodyMotion<double,spaceDim>& beginning,
- const RigidBodyMotion<double,spaceDim>& end)
+ void create(std::vector<Dune::GFE::ProductManifold<RealTuple<double,spaceDim>,Rotation<double,spaceDim> > >& rod,
+ const Dune::GFE::ProductManifold<RealTuple<double,spaceDim>,Rotation<double,spaceDim> > & beginning,
+ const Dune::GFE::ProductManifold<RealTuple<double,spaceDim>,Rotation<double,spaceDim> > & end)
{
static const int dim = GridView::dimension; // de facto: 1
@@ -106,8 +110,8 @@ public:
\param[out] rod The new rod
*/
template <int spaceDim>
- void create(std::vector<RigidBodyMotion<double,spaceDim> >& rod,
- const RigidBodyMotion<double,spaceDim>& value)
+ void create(std::vector<Dune::GFE::ProductManifold<RealTuple<double,spaceDim>,Rotation<double,spaceDim> > >& rod,
+ const Dune::GFE::ProductManifold<RealTuple<double,spaceDim>,Rotation<double,spaceDim> > & value)
{
rod.resize(gridView_.size(1));
std::fill(rod.begin(), rod.end(), value);
@@ -119,7 +123,7 @@ public:
\param[in,out] rod The new rod
*/
template <int spaceDim>
- void create(std::vector<RigidBodyMotion<double,spaceDim> >& rod)
+ void create(std::vector<Dune::GFE::ProductManifold<RealTuple<double,spaceDim>,Rotation<double,spaceDim> > >& rod)
{
static const int dim = GridView::dimension; // de facto: 1
assert(gridView_.size(dim)==rod.size());
@@ -133,7 +137,7 @@ public:
double min = std::numeric_limits<double>::max();
double max = -std::numeric_limits<double>::max();
- RigidBodyMotion<double,spaceDim> beginning, end;
+ Dune::GFE::ProductManifold<RealTuple<double,spaceDim>,Rotation<double,spaceDim> > beginning, end;
for (; vIt != vEndIt; ++vIt) {
if (vIt->geometry().corner(0)[0] < min) {
@@ -174,10 +178,10 @@ public:
\param[out] rod The new rod
*/
template <int spaceDim>
- void create(std::vector<RigidBodyMotion<double,spaceDim> >& rod,
+ void create(std::vector<Dune::GFE::ProductManifold<RealTuple<double,spaceDim>,Rotation<double,spaceDim> > >& rod,
double radius, double E, double nu,
- const RigidBodyMotion<double,spaceDim>& beginning,
- const RigidBodyMotion<double,spaceDim>& end)
+ const Dune::GFE::ProductManifold<RealTuple<double,spaceDim>,Rotation<double,spaceDim> > & beginning,
+ const Dune::GFE::ProductManifold<RealTuple<double,spaceDim>,Rotation<double,spaceDim> > & end)
{
// Make Dirichlet bitfields for the rods as well
@@ -207,7 +211,7 @@ public:
rod.back() = end;
// Trust--Region solver
- RiemannianTrustRegionSolver<RodP1Basis, RigidBodyMotion<double,spaceDim> > rodSolver;
+ RiemannianTrustRegionSolver<RodP1Basis, Dune::GFE::ProductManifold<RealTuple<double,spaceDim>,Rotation<double,spaceDim> > > rodSolver;
rodSolver.setup(gridView_.grid(), &assembler, rod,
rodDirichletNodes,
1e-10, 100, // TR tolerance and iterations
diff --git a/dune/gfe/spaces/CMakeLists.txt b/dune/gfe/spaces/CMakeLists.txt
index f158d605..bb85104d 100644
--- a/dune/gfe/spaces/CMakeLists.txt
+++ b/dune/gfe/spaces/CMakeLists.txt
@@ -5,7 +5,6 @@ install(FILES
productmanifold.hh
quaternion.hh
realtuple.hh
- rigidbodymotion.hh
rotation.hh
unitvector.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/gfe/spaces)
diff --git a/dune/gfe/spaces/rigidbodymotion.hh b/dune/gfe/spaces/rigidbodymotion.hh
deleted file mode 100644
index e284c18f..00000000
--- a/dune/gfe/spaces/rigidbodymotion.hh
+++ /dev/null
@@ -1,431 +0,0 @@
-#ifndef DUNE_GFE_SPACES_RIGIDBODYMOTION_HH
-#define DUNE_GFE_SPACES_RIGIDBODYMOTION_HH
-
-#include <dune/common/fvector.hh>
-
-#include <dune/gfe/spaces/realtuple.hh>
-#include <dune/gfe/spaces/rotation.hh>
-
-/** \brief A rigid-body motion in R^N, i.e., a member of SE(N) */
-template <class T, int N>
-struct RigidBodyMotion
-{
-public:
-
- /** \brief Dimension of manifold */
- static const int dim = N + Rotation<T,N>::dim;
-
- /** \brief Dimension of the embedding space */
- static const int embeddedDim = N + Rotation<T,N>::embeddedDim;
-
- /** \brief Type of an infinitesimal rigid body motion */
- typedef Dune::FieldVector<T, dim> TangentVector;
-
- /** \brief Type of an infinitesimal rigid body motion */
- typedef Dune::FieldVector<T, embeddedDim> EmbeddedTangentVector;
-
- /** \brief The type used for coordinates */
- typedef T ctype;
- typedef T field_type;
-
- /** \brief The type used for global coordinates */
- typedef Dune::FieldVector<T,embeddedDim> CoordinateType;
-
- /** \brief The global convexity radius of the rigid body motions */
- static constexpr double convexityRadius = Rotation<T,N>::convexityRadius;
-
- /** \brief Default constructor */
- RigidBodyMotion()
- {}
-
- /** \brief Constructor from a translation and a rotation */
- RigidBodyMotion(const Dune::FieldVector<ctype, N>& translation,
- const Rotation<ctype,N>& rotation)
- : r(translation), q(rotation)
- {}
-
- RigidBodyMotion(const CoordinateType& globalCoordinates)
- {
- for (int i=0; i<N; i++)
- r[i] = globalCoordinates[i];
-
- for (int i=N; i<embeddedDim; i++)
- q[i-N] = globalCoordinates[i];
-
- // Turn this into a unit quaternion if it isn't already
- q.normalize();
- }
-
- /** \brief Assignment from RigidBodyMotion with different type -- used for automatic differentiation with ADOL-C */
- template <class T2>
- RigidBodyMotion& operator <<= (const RigidBodyMotion<T2,N>& other) {
- for (int i=0; i<N; i++)
- r[i] <<= other.r[i];
- q <<= other.q;
- return *this;
- }
-
- /** \brief Rebind the RigidBodyMotion to another coordinate type */
- template<class U>
- struct rebind
- {
- typedef RigidBodyMotion<U,N> other;
- };
-
- /** \brief The exponential map from a given point $p \in SE(d)$.
-
- Why the template parameter? Well, it should work with both TangentVector and EmbeddedTangentVector.
- In general these differ and we could just have two exp methods. However in 2d they do _not_ differ,
- and then the compiler complains about having two methods with the same signature.
- */
- template <class TVector>
- static RigidBodyMotion<ctype,N> exp(const RigidBodyMotion<ctype,N>& p, const TVector& v) {
-
- RigidBodyMotion<ctype,N> result;
-
- // Add translational correction
- for (int i=0; i<N; i++)
- result.r[i] = p.r[i] + v[i];
-
- // Add rotational correction
- typedef typename std::conditional<std::is_same<TVector,TangentVector>::value,
- typename Rotation<ctype,N>::TangentVector,
- typename Rotation<ctype,N>::EmbeddedTangentVector>::type RotationTangentVector;
- RotationTangentVector qCorr;
- for (int i=0; i<RotationTangentVector::dimension; i++)
- qCorr[i] = v[N+i];
-
- result.q = Rotation<ctype,N>::exp(p.q, qCorr);
- return result;
- }
-
- /** \brief Compute difference vector from a to b on the tangent space of a */
- static EmbeddedTangentVector log(const RigidBodyMotion<ctype,N>& a,
- const RigidBodyMotion<ctype,N>& b)
- {
- EmbeddedTangentVector result;
-
- // Usual linear difference
- for (int i=0; i<N; i++)
- result[i] = b.r[i] - a.r[i];
-
- // Subtract orientations on the tangent space of 'a'
- typename Rotation<ctype,N>::EmbeddedTangentVector v = Rotation<ctype,N>::log(a.q, b.q);
-
- // Compute difference on T_a SO(3)
- for (int i=0; i<Rotation<ctype,N>::EmbeddedTangentVector::dimension; i++)
- result[i+N] = v[i];
-
- return result;
- }
-
- /** \brief Compute geodesic distance from a to b */
- static T distance(const RigidBodyMotion<ctype,N>& a, const RigidBodyMotion<ctype,N>& b) {
-
- T euclideanDistanceSquared = (a.r - b.r).two_norm2();
-
- T rotationDistance = Rotation<ctype,N>::distance(a.q, b.q);
-
- return std::sqrt(euclideanDistanceSquared + rotationDistance*rotationDistance);
- }
-
- /** \brief Compute difference vector from a to b on the tangent space of a
-
- * \warning The method is buggy! See https://gitlab.mn.tu-dresden.de/osander/dune-gfe/-/merge_requests/2
- * \deprecated Use the log method instead!
- */
- [[deprecated("Use RigidBodyMotion::log instead of RigidBodyMotion::difference!")]]
- static TangentVector difference(const RigidBodyMotion<ctype,N>& a,
- const RigidBodyMotion<ctype,N>& b) {
-
- TangentVector result;
-
- // Usual linear difference
- for (int i=0; i<N; i++)
- result[i] = a.r[i] - b.r[i];
-
- // Subtract orientations on the tangent space of 'a'
- typename Rotation<ctype,N>::TangentVector v = Rotation<ctype,N>::difference(a.q, b.q).axial();
-
- // Compute difference on T_a SO(3)
- for (int i=0; i<Rotation<ctype,N>::TangentVector::dimension; i++)
- result[i+N] = v[i];
-
- return result;
- }
-
- static EmbeddedTangentVector derivativeOfDistanceSquaredWRTSecondArgument(const RigidBodyMotion<ctype,N>& a,
- const RigidBodyMotion<ctype,N>& b) {
-
- // linear part
- Dune::FieldVector<ctype,N> linearDerivative = a.r;
- linearDerivative -= b.r;
- linearDerivative *= -2;
-
- // rotation part
- typename Rotation<ctype,N>::EmbeddedTangentVector rotationDerivative
- = Rotation<ctype,N>::derivativeOfDistanceSquaredWRTSecondArgument(a.q, b.q);
-
- return concat(linearDerivative, rotationDerivative);
- }
-
- /** \brief Compute the Hessian of the squared distance function keeping the first argument fixed */
- static Dune::SymmetricMatrix<T,embeddedDim> secondDerivativeOfDistanceSquaredWRTSecondArgument(const RigidBodyMotion<ctype,N> & p, const RigidBodyMotion<ctype,N> & q)
- {
- Dune::SymmetricMatrix<T,embeddedDim> result;
-
- // The linear part
- Dune::SymmetricMatrix<T,N> linearPart = RealTuple<T,N>::secondDerivativeOfDistanceSquaredWRTSecondArgument(p.r,q.r);
- for (int i=0; i<N; i++)
- for (int j=0; j<=i; j++)
- result(i,j) = linearPart(i,j);
-
- // The rotation part
- Dune::SymmetricMatrix<T,Rotation<T,N>::embeddedDim> rotationPart
- = Rotation<ctype,N>::secondDerivativeOfDistanceSquaredWRTSecondArgument(p.q,q.q);
- for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
- {
- for (int j=0; j<N; j++)
- result(N+i,j) = 0;
- for (int j=0; j<=i; j++)
- result(N+i,N+j) = rotationPart(i,j);
- }
-
- return result;
- }
-
- /** \brief Compute the mixed second derivate \partial d^2 / \partial da db
-
- Unlike the distance itself the squared distance is differentiable at zero
- */
- static Dune::FieldMatrix<T,embeddedDim,embeddedDim> secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(const RigidBodyMotion<ctype,N> & p, const RigidBodyMotion<ctype,N> & q)
- {
- Dune::FieldMatrix<T,embeddedDim,embeddedDim> result(0);
-
- // The linear part
- Dune::FieldMatrix<T,N,N> linearPart = RealTuple<T,N>::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(p.r,q.r);
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- result[i][j] = linearPart[i][j];
-
- // The rotation part
- Dune::FieldMatrix<T,Rotation<T,N>::embeddedDim,Rotation<T,N>::embeddedDim> rotationPart
- = Rotation<ctype,N>::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(p.q,q.q);
- for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
- for (int j=0; j<Rotation<T,N>::embeddedDim; j++)
- result[N+i][N+j] = rotationPart[i][j];
-
- return result;
- }
-
- /** \brief Compute the third derivative \partial d^3 / \partial dq^3
-
- Unlike the distance itself the squared distance is differentiable at zero
- */
- static Tensor3<T,embeddedDim,embeddedDim,embeddedDim> thirdDerivativeOfDistanceSquaredWRTSecondArgument(const RigidBodyMotion<ctype,N> & p, const RigidBodyMotion<ctype,N> & q)
- {
- Tensor3<T,embeddedDim,embeddedDim,embeddedDim> result(0);
-
- // The linear part
- Tensor3<T,N,N,N> linearPart = RealTuple<T,N>::thirdDerivativeOfDistanceSquaredWRTSecondArgument(p.r,q.r);
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- for (int k=0; k<N; k++)
- result[i][j][k] = linearPart[i][j][k];
-
- // The rotation part
- Tensor3<T,Rotation<T,N>::embeddedDim,Rotation<T,N>::embeddedDim,Rotation<T,N>::embeddedDim> rotationPart
- = Rotation<ctype,N>::thirdDerivativeOfDistanceSquaredWRTSecondArgument(p.q,q.q);
-
- for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
- for (int j=0; j<Rotation<T,N>::embeddedDim; j++)
- for (int k=0; k<Rotation<T,N>::embeddedDim; k++)
- result[N+i][N+j][N+k] = rotationPart[i][j][k];
-
- return result;
- }
-
- /** \brief Compute the mixed third derivative \partial d^3 / \partial da db^2
-
- Unlike the distance itself the squared distance is differentiable at zero
- */
- static Tensor3<T,embeddedDim,embeddedDim,embeddedDim> thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(const RigidBodyMotion<ctype,N> & p, const RigidBodyMotion<ctype,N> & q)
- {
- Tensor3<T,embeddedDim,embeddedDim,embeddedDim> result(0);
-
- // The linear part
- Tensor3<T,N,N,N> linearPart = RealTuple<T,N>::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(p.r,q.r);
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- for (int k=0; k<N; k++)
- result[i][j][k] = linearPart[i][j][k];
-
- // The rotation part
- Tensor3<T,Rotation<T,N>::embeddedDim,Rotation<T,N>::embeddedDim,Rotation<T,N>::embeddedDim> rotationPart = Rotation<ctype,N>::thirdDerivativeOfDistanceSquaredWRTFirst1AndSecond2Argument(p.q,q.q);
- for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
- for (int j=0; j<Rotation<T,N>::embeddedDim; j++)
- for (int k=0; k<Rotation<T,N>::embeddedDim; k++)
- result[N+i][N+j][N+k] = rotationPart[i][j][k];
-
- return result;
- }
-
-
-
- /** \brief Project tangent vector of R^n onto the tangent space */
- EmbeddedTangentVector projectOntoTangentSpace(const EmbeddedTangentVector& v) const {
- EmbeddedTangentVector result;
-
- // translation part
- for (int i=0; i<N; i++)
- result[i] = v[i];
-
- // rotation part
- typename Rotation<T,N>::EmbeddedTangentVector rotV;
- for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
- rotV[i] = v[i+N];
-
- rotV = q.projectOntoTangentSpace(rotV);
-
- for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
- result[i+N] = rotV[i];
-
- return result;
- }
-
- /** \brief Project tangent vector of R^n onto the normal space space */
- EmbeddedTangentVector projectOntoNormalSpace(const EmbeddedTangentVector& v) const {
-
- EmbeddedTangentVector result;
-
- // translation part
- for (int i=0; i<N; i++)
- result[i] = 0;
-
- // rotation part
- T sp = 0;
- for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
- sp += v[i+N] * q[i];
-
- for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
- result[i+N] = sp * q[i];
-
- return result;
- }
-
- /** \brief Transform tangent vectors from quaternion representation to matrix representation
- *
- * This class represents rotations as unit quaternions, and therefore variations of rotations
- * (i.e., tangent vector) are represented in quaternion space, too. However, some applications
- * require the tangent vectors as matrices. To obtain matrix coordinates we use the
- * chain rule, which means that we have to multiply the given derivative with
- * the derivative of the embedding of the unit quaternion into the space of 3x3 matrices.
- * This second derivative is almost given by the method getFirstDerivativesOfDirectors.
- * However, since the directors of a given unit quaternion are the _columns_ of the
- * corresponding orthogonal matrix, we need to invert the i and j indices
- *
- * As a typical GFE assembler will require this for several tangent vectors at once,
- * the implementation here is a vector one: It allows to treat several tangent vectors
- * together, which have to come as the columns of the `derivative` parameter.
- *
- * So, if I am not mistaken, result[i][j][k] contains \partial R_ij / \partial k
- *
- * \param derivative Tangent vector in quaternion coordinates
- * \returns DR Tangent vector in matrix coordinates
- */
- template <int blocksize>
- Tensor3<T,3,3,blocksize> quaternionTangentToMatrixTangent(const Dune::FieldMatrix<T,7,blocksize>& derivative) const
- {
- Tensor3<T,3,3,blocksize> result = T(0);
-
- Tensor3<T,3 , 3, 4> dd_dq;
- q.getFirstDerivativesOfDirectors(dd_dq);
-
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++)
- for (int k=0; k<blocksize; k++)
- for (int l=0; l<4; l++)
- result[i][j][k] += dd_dq[j][i][l] * derivative[l+3][k];
-
- return result;
- }
-
- /** \brief The Weingarten map */
- EmbeddedTangentVector weingarten(const EmbeddedTangentVector& z, const EmbeddedTangentVector& v) const {
-
- EmbeddedTangentVector result;
-
- // translation part: nothing, the space is flat
- for (int i=0; i<N; i++)
- result[i] = 0;
-
- // rotation part
- T sp = 0;
- for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
- sp += v[i+N] * q[i];
-
- for (int i=0; i<Rotation<T,N>::embeddedDim; i++)
- result[i+N] = -sp * z[i+N];
-
- return result;
- }
-
- /** \brief Compute an orthonormal basis of the tangent space of SE(3).
-
- This basis may not be globally continuous.
- */
- Dune::FieldMatrix<T,dim,embeddedDim> orthonormalFrame() const {
- Dune::FieldMatrix<T,dim,embeddedDim> result(0);
-
- // Get the R^d part
- for (int i=0; i<N; i++)
- result[i][i] = 1;
-
- Dune::FieldMatrix<T,Rotation<T,N>::dim,Rotation<T,N>::embeddedDim> SO3Part = q.orthonormalFrame();
-
- for (int i=0; i<Rotation<T,N>::dim; i++)
- for (int j=0; j<Rotation<T,N>::embeddedDim; j++)
- result[N+i][N+j] = SO3Part[i][j];
-
- return result;
- }
-
- /** \brief The global coordinates, if you really want them */
- CoordinateType globalCoordinates() const {
- return concat(r, q.globalCoordinates());
- }
-
-
-
- // Translational part
- Dune::FieldVector<ctype, N> r;
-
- // Rotational part
- Rotation<ctype,N> q;
-
-private:
-
- /** \brief Concatenate two FieldVectors */
- template <int NN, int M>
- static Dune::FieldVector<ctype,NN+M> concat(const Dune::FieldVector<ctype,NN>& a,
- const Dune::FieldVector<ctype,M>& b)
- {
- Dune::FieldVector<ctype,NN+M> result;
- for (int i=0; i<NN; i++)
- result[i] = a[i];
- for (int i=0; i<M; i++)
- result[i+NN] = b[i];
- return result;
- }
-
-};
-
-//! Send configuration to output stream
-template <int N, class ctype>
-std::ostream& operator<< (std::ostream& s, const RigidBodyMotion<ctype,N>& c)
-{
- s << "(" << c.r << ") (" << c.q << ")";
- return s;
-}
-
-#endif
diff --git a/src/compute-disc-error.cc b/src/compute-disc-error.cc
index f5b5b277..db771c37 100644
--- a/src/compute-disc-error.cc
+++ b/src/compute-disc-error.cc
@@ -222,7 +222,7 @@ void measureDiscreteEOC(const GridView gridView,
auto localReferenceSolution = localFunction(referenceSolution);
auto localNumericalSolution = localFunction(numericalSolution);
- if (std::is_same<TargetSpace,RigidBodyMotion<double,3> >::value)
+ if (std::is_same<TargetSpace,GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >::value)
{
double deformationL2ErrorSquared = 0;
double orientationL2ErrorSquared = 0;
@@ -852,9 +852,9 @@ int main (int argc, char *argv[]) try
parameterSet);
} else if (targetSpace=="RigidBodyMotion")
{
- measureEOC<GridType,RigidBodyMotion<double,2> >(grid,
- referenceGrid,
- parameterSet);
+ measureEOC<GridType,GFE::ProductManifold<RealTuple<double,2>,Rotation<double,2> > >(grid,
+ referenceGrid,
+ parameterSet);
#endif
} else
DUNE_THROW(NotImplemented, "Target space '" << targetSpace << "' is not implemented");
@@ -878,9 +878,9 @@ int main (int argc, char *argv[]) try
parameterSet);
} else if (targetSpace=="RigidBodyMotion")
{
- measureEOC<GridType,RigidBodyMotion<double,3> >(grid,
- referenceGrid,
- parameterSet);
+ measureEOC<GridType,GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >(grid,
+ referenceGrid,
+ parameterSet);
} else
DUNE_THROW(NotImplemented, "Target space '" << targetSpace << "' is not implemented");
break;
@@ -904,9 +904,9 @@ int main (int argc, char *argv[]) try
parameterSet);
} else if (targetSpace=="RigidBodyMotion")
{
- measureEOC<GridType,RigidBodyMotion<double,4> >(grid,
- referenceGrid,
- parameterSet);
+ measureEOC<GridType,GFE::ProductManifold<RealTuple<double,4>,Rotation<double,4> > >(grid,
+ referenceGrid,
+ parameterSet);
#endif
} else
DUNE_THROW(NotImplemented, "Target space '" << targetSpace << "' is not implemented");
diff --git a/src/cosserat-continuum.cc b/src/cosserat-continuum.cc
index c81a9192..8906f527 100644
--- a/src/cosserat-continuum.cc
+++ b/src/cosserat-continuum.cc
@@ -70,7 +70,9 @@
#include <dune/gfe/assemblers/geodesicfeassemblerwrapper.hh>
#include <dune/gfe/riemannianpnsolver.hh>
#include <dune/gfe/riemanniantrsolver.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
#endif
#if HAVE_DUNE_VTK
@@ -80,6 +82,8 @@
#include <dune/gmsh4/gmsh4reader.hh>
#include <dune/gmsh4/gridcreators/lagrangegridcreator.hh>
+using namespace Dune;
+
// grid dimension
const int dim = GRID_DIM;
const int dimworld = WORLD_DIM;
@@ -94,10 +98,9 @@ const int rotationOrder = GFE_ORDER;
static_assert(displacementOrder==rotationOrder, "displacement and rotation order do not match!");
// Image space of the geodesic fe functions
-using TargetSpace = RigidBodyMotion<double, 3>;
+using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
#endif
-using namespace Dune;
int main (int argc, char *argv[]) try
{
@@ -381,8 +384,8 @@ int main (int argc, char *argv[]) try
for (size_t i=0; i<x.size(); i++) {
auto vTargetSpace = TargetSpace(v[i]);
- x[_0][i] = vTargetSpace.r;
- x[_1][i] = vTargetSpace.q;
+ x[_0][i] = vTargetSpace[_0];
+ x[_1][i] = vTargetSpace[_1];
}
}
#endif
@@ -523,16 +526,16 @@ int main (int argc, char *argv[]) try
x = solver.getSol();
#else
//The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
- //The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a RigidBodyMotion
- //Therefore, x and the dirichletDofs are converted to a RigidBodyMotion structure, as well as the Hessian and Gradient that are returned by the assembler
+ //The RiemannianTrustRegionSolver can only treat the Displacement 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
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++) {
for (int j = 0; j < 3; j ++) { // Displacement part
- xTargetSpace[i].r[j] = x[_0][i][j];
+ xTargetSpace[i][_0].globalCoordinates()[j] = x[_0][i][j];
dirichletDofsTargetSpace[i][j] = deformationDirichletDofs[i][j];
}
- xTargetSpace[i].q = x[_1][i]; // Rotation part
+ xTargetSpace[i][_1] = x[_1][i]; // Rotation part
for (int j = 3; j < TargetSpace::TangentVector::dimension; j ++)
dirichletDofsTargetSpace[i][j] = orientationDirichletDofs[i][j-3];
}
@@ -575,8 +578,8 @@ int main (int argc, char *argv[]) try
xTargetSpace = solver.getSol();
}
for (std::size_t i = 0; i < xTargetSpace.size(); i++) {
- x[_0][i] = xTargetSpace[i].r;
- x[_1][i] = xTargetSpace[i].q;
+ x[_0][i] = xTargetSpace[i][_0];
+ x[_1][i] = xTargetSpace[i][_1];
}
#endif
} else { //dim != dimworld
@@ -623,16 +626,16 @@ int main (int argc, char *argv[]) try
x = solver.getSol();
#else
//The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
- //The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a RigidBodyMotion
- //Therefore, x and the dirichletDofs are converted to a RigidBodyMotion structure, as well as the Hessian and Gradient that are returned by the assembler
+ //The RiemannianTrustRegionSolver can only treat the Displacement 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
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++) {
for (int j = 0; j < 3; j ++) { // Displacement part
- xTargetSpace[i].r[j] = x[_0][i][j];
+ xTargetSpace[i][_0].globalCoordinates()[j] = x[_0][i][j];
dirichletDofsTargetSpace[i][j] = deformationDirichletDofs[i][j];
}
- xTargetSpace[i].q = x[_1][i]; // Rotation part
+ xTargetSpace[i][_1] = x[_1][i]; // Rotation part
for (int j = 3; j < TargetSpace::TangentVector::dimension; j ++)
dirichletDofsTargetSpace[i][j] = orientationDirichletDofs[i][j-3];
}
@@ -676,8 +679,8 @@ int main (int argc, char *argv[]) try
}
for (std::size_t i = 0; i < xTargetSpace.size(); i++) {
- x[_0][i] = xTargetSpace[i].r;
- x[_1][i] = xTargetSpace[i].q;
+ x[_0][i] = xTargetSpace[i][_0];
+ x[_1][i] = xTargetSpace[i][_1];
}
#endif
}
diff --git a/src/film-on-substrate.cc b/src/film-on-substrate.cc
index 016275d8..33fb40eb 100644
--- a/src/film-on-substrate.cc
+++ b/src/film-on-substrate.cc
@@ -69,7 +69,7 @@
#include <dune/gfe/assemblers/geodesicfeassemblerwrapper.hh>
#include <dune/gfe/riemannianpnsolver.hh>
#include <dune/gfe/riemanniantrsolver.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
#endif
#include <dune/istl/multitypeblockvector.hh>
@@ -564,18 +564,18 @@ int main (int argc, char *argv[]) try
x[_1][i].set(dOV[i]);
#if !MIXED_SPACE
- //The MixedRiemannianTrustRegionSolver can treat the Displacement and Orientation Space as separate ones
- //The RiemannianTrustRegionSolver can only treat the Displacement and Rotation together in a RigidBodyMotion
- //Therefore, x and the dirichletDofs are converted to a RigidBodyMotion structure, as well as the Hessian and Gradient that are returned by the assembler
- typedef RigidBodyMotion<double, dim> RBM;
+ //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++) {
for (int j = 0; j < dim; j ++) { // Displacement part
- xRBM[i].r[j] = x[_0][i][j];
+ xRBM[i][_0].globalCoordinates()[j] = x[_0][i][j];
dirichletDofsRBM[i][j] = dirichletDofs[_0][i][j];
}
- xRBM[i].q = x[_1][i]; // Rotation part
+ xRBM[i][_1] = x[_1][i]; // Rotation part
for (int j = dim; j < RBM::TangentVector::dimension; j ++)
dirichletDofsRBM[i][j] = dirichletDofs[_1][i][j-dim];
}
@@ -632,8 +632,8 @@ int main (int argc, char *argv[]) try
solver.solve();
xRBM = solver.getSol();
for (int i = 0; i < xRBM.size(); i++) {
- x[_0][i] = xRBM[i].r;
- x[_1][i] = xRBM[i].q;
+ x[_0][i] = xRBM[i][_0];
+ x[_1][i] = xRBM[i][_1];
}
#endif
} else { //parameterSet.get<std::string>("solvertype") == "proximalNewton"
@@ -655,8 +655,8 @@ int main (int argc, char *argv[]) try
solver.solve();
xRBM = solver.getSol();
for (int i = 0; i < xRBM.size(); i++) {
- x[_0][i] = xRBM[i].r;
- x[_1][i] = xRBM[i].q;
+ x[_0][i] = xRBM[i][_0];
+ x[_1][i] = xRBM[i][_1];
}
#endif
}
diff --git a/src/harmonicmaps.cc b/src/harmonicmaps.cc
index 4033c842..940e55da 100644
--- a/src/harmonicmaps.cc
+++ b/src/harmonicmaps.cc
@@ -47,7 +47,7 @@
#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/embeddedglobalgfefunction.hh>
#include <dune/gfe/spaces/realtuple.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
#include <dune/gfe/spaces/rotation.hh>
#include <dune/gfe/spaces/unitvector.hh>
@@ -58,7 +58,6 @@ const int dimworld = 2;
// Image space of the geodesic fe functions
// typedef Rotation<double,2> TargetSpace;
// typedef Rotation<double,3> TargetSpace;
-// typedef RigidBodyMotion<double,3> TargetSpace;
// typedef UnitVector<double,2> TargetSpace;
typedef UnitVector<double,3> TargetSpace;
// typedef UnitVector<double,4> TargetSpace;
diff --git a/src/rod3d.cc b/src/rod3d.cc
index d934c906..e0e25e63 100644
--- a/src/rod3d.cc
+++ b/src/rod3d.cc
@@ -39,16 +39,20 @@
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/riemanniantrsolver.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
-typedef RigidBodyMotion<double,3> TargetSpace;
+using namespace Dune;
+using namespace Dune::Indices;
+
+using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
const int blocksize = TargetSpace::TangentVector::dimension;
// Approximation order of the finite element space
constexpr int order = 2;
-using namespace Dune;
int main (int argc, char *argv[]) try
{
@@ -112,15 +116,13 @@ int main (int argc, char *argv[]) try
Functions::interpolate(scalarBasis, referenceConfigurationX, identity);
- using Configuration = std::vector<RigidBodyMotion<double,3> >;
+ using Configuration = std::vector<TargetSpace>;
Configuration referenceConfiguration(scalarBasis.size());
for (std::size_t i=0; i<referenceConfiguration.size(); i++)
{
- referenceConfiguration[i].r[0] = 0;
- referenceConfiguration[i].r[1] = 0;
- referenceConfiguration[i].r[2] = referenceConfigurationX[i];
- referenceConfiguration[i].q = Rotation<double,3>::identity();
+ referenceConfiguration[i][_0] = {0, 0, referenceConfigurationX[i]};
+ referenceConfiguration[i][_1] = Rotation<double,3>::identity();
}
// Select the reference configuration as initial iterate
@@ -157,18 +159,18 @@ int main (int argc, char *argv[]) try
}
// Set Dirichlet values
- x[rightBoundaryDof].r = parameterSet.get<FieldVector<double,3> >("dirichletValue");
+ x[rightBoundaryDof][_0] = parameterSet.get<FieldVector<double,3> >("dirichletValue");
auto axis = parameterSet.get<FieldVector<double,3> >("dirichletAxis");
double angle = parameterSet.get<double>("dirichletAngle");
- x[rightBoundaryDof].q = Rotation<double,3>(axis, M_PI*angle/180);
+ x[rightBoundaryDof][_1] = Rotation<double,3>(axis, M_PI*angle/180);
// backup for error measurement later
std::cout << "Right boundary orientation:" << std::endl;
- std::cout << "director 0: " << x[rightBoundaryDof].q.director(0) << std::endl;
- std::cout << "director 1: " << x[rightBoundaryDof].q.director(1) << std::endl;
- std::cout << "director 2: " << x[rightBoundaryDof].q.director(2) << std::endl;
+ std::cout << "director 0: " << x[rightBoundaryDof][_1].director(0) << std::endl;
+ std::cout << "director 1: " << x[rightBoundaryDof][_1].director(1) << std::endl;
+ std::cout << "director 2: " << x[rightBoundaryDof][_1].director(2) << std::endl;
//////////////////////////////////////////////
// Create the energy and assembler
@@ -207,7 +209,7 @@ int main (int argc, char *argv[]) try
// Create a solver for the rod problem
/////////////////////////////////////////////
- RiemannianTrustRegionSolver<ScalarBasis,RigidBodyMotion<double,3> > rodSolver;
+ RiemannianTrustRegionSolver<ScalarBasis,TargetSpace> rodSolver;
rodSolver.setup(grid,
&rodAssembler,
@@ -253,7 +255,7 @@ int main (int argc, char *argv[]) try
// The rod displacement field
BlockVector<FieldVector<double,3> > displacement(worldBasis.size());
for (std::size_t i=0; i<x.size(); i++)
- displacement[i] = x[i].r;
+ displacement[i] = x[i][_0].globalCoordinates();
std::vector<double> xEmbedding;
Functions::interpolate(scalarBasis, xEmbedding, [](FieldVector<double,1> x){
@@ -277,7 +279,7 @@ int main (int argc, char *argv[]) try
{
director[i].resize(worldBasis.size());
for (std::size_t j=0; j<x.size(); j++)
- director[i][j] = x[j].q.director(i);
+ director[i][j] = x[j][_1].director(i);
directorFunction[i] = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,3> >(worldBasis, std::move(director[i]));
vtkWriter.addPointData(*directorFunction[i], "director " + std::to_string(i), 3);
diff --git a/src/rodobstacle.cc b/src/rodobstacle.cc
index 6267e6d2..c1b58b70 100644
--- a/src/rodobstacle.cc
+++ b/src/rodobstacle.cc
@@ -19,6 +19,9 @@
#include <dune/gfe/rodwriter.hh>
#include <dune/gfe/rodassembler.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
// 3 (x, y, theta) for a planar rod
@@ -71,7 +74,7 @@ int main (int argc, char *argv[]) try
// Some types that I need
typedef BCRSMatrix<FieldMatrix<double, blocksize, blocksize> > MatrixType;
typedef BlockVector<FieldVector<double, blocksize> > CorrectionType;
- typedef std::vector<RigidBodyMotion<double,2> > SolutionType;
+ typedef std::vector<GFE::ProductManifold<RealTuple<double,2>,Rotation<doublee,2> > > SolutionType;
// parse data file
ParameterTree parameterSet;
@@ -293,7 +296,7 @@ int main (int argc, char *argv[]) try
SolutionType newIterate = x;
for (int j=0; j<newIterate.size(); j++)
- newIterate[j] = RigidBodyMotion<double,2>::exp(newIterate[j], corr[j]);
+ newIterate[j] = exp(newIterate[j], corr[j]);
/** \todo Don't always recompute oldEnergy */
double oldEnergy = rodAssembler.computeEnergy(x);
@@ -304,7 +307,7 @@ int main (int argc, char *argv[]) try
// Add correction to the current solution
for (int j=0; j<x.size(); j++)
- x[j] = RigidBodyMotion<double,2>::exp(x[j], corr[j]);
+ x[j] = exp(x[j], corr[j]);
// Subtract correction from the current obstacle
for (int k=0; k<corr.size(); k++)
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 296f756e..3ce0039e 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -12,7 +12,6 @@ set(TESTS
localprojectedfefunctiontest
orthogonalmatrixtest
polardecompositiontest
- rigidbodymotiontest
rotationtest
svdtest
targetspacetest
@@ -62,4 +61,4 @@ file(COPY stressPlotData/stressPlotTestRotation DESTINATION ${CMAKE_CURRENT_BINA
file(COPY stressPlotData/stressPlotTestInitialDeformation DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/stressPlotData)
# Copy an .adolcrc file into the build dir to make sure the tapes stay in memory and do not get written to disk because this causes errors sometimes
-file(COPY .adolcrc DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
\ No newline at end of file
+file(COPY .adolcrc DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
diff --git a/test/adolctest-scalar-and-vector-mode.cc b/test/adolctest-scalar-and-vector-mode.cc
index 61d8f1e5..a06d3b98 100644
--- a/test/adolctest-scalar-and-vector-mode.cc
+++ b/test/adolctest-scalar-and-vector-mode.cc
@@ -55,7 +55,7 @@ using MatrixRow1 = MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,dimCR,dim>
using MatrixTypeMixed = MultiTypeBlockMatrix<MatrixRow0,MatrixRow1>;
//Types for the Non-mixed space
-using RBM = RigidBodyMotion<double, dim>;
+using RBM = GFE::ProductManifold<RealTuple<double,dim>,Rotation<double,dim> >;
using RBMVector = std::vector<RBM>;
const static int blocksize = RBM::TangentVector::dimension;
using CorrectionType = BlockVector<FieldVector<double, blocksize> >;
@@ -170,7 +170,7 @@ int main (int argc, char *argv[])
initialDeformation[i][j] = std::cos(identity[i][j]);
initialDeformation[i][2] = identity[i][0]*identity[i][0];
x[_0][i] = initialDeformation[i];
- xRBM[i].r = initialDeformation[i];
+ xRBM[i][_0] = initialDeformation[i];
}
//////////////////////////////////////////////////////////////////////////////
diff --git a/test/adolctest.cc b/test/adolctest.cc
index b43c9062..d53f84b7 100644
--- a/test/adolctest.cc
+++ b/test/adolctest.cc
@@ -42,15 +42,17 @@ typedef double FDType;
#include <dune/gfe/assemblers/cosseratenergystiffness.hh>
#include <dune/gfe/assemblers/localgeodesicfeadolcstiffness.hh>
#include <dune/gfe/assemblers/localgeodesicfefdstiffness.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
+
+using namespace Dune;
// grid dimension
const int dim = 2;
// Image space of the geodesic fe functions
-typedef RigidBodyMotion<double,3> TargetSpace;
-
-using namespace Dune;
+using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
/** \brief Assembles energy gradient and Hessian with ADOL-C
*/
@@ -477,7 +479,7 @@ int main (int argc, char *argv[]) try
Functions::interpolate(powerBasis, v, identity);
for (size_t i=0; i<x.size(); i++)
- x[i].r = v[i];
+ x[i][Indices::_0] = v[i];
#endif
// ////////////////////////////////////////////////////////////
diff --git a/test/cosseratenergytest.cc b/test/cosseratenergytest.cc
index f9348a2f..9589bb56 100644
--- a/test/cosseratenergytest.cc
+++ b/test/cosseratenergytest.cc
@@ -8,15 +8,18 @@
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/gfe/assemblers/cosseratenergystiffness.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
#include "multiindex.hh"
#include "valuefactory.hh"
-typedef RigidBodyMotion<double,3> TargetSpace;
-
using namespace Dune;
+using namespace Dune::Indices;
+
+using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
// ////////////////////////////////////////////////////////
// Make a test grid consisting of a single simplex
@@ -97,10 +100,10 @@ void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficien
for (size_t j=0; j<coefficients.size(); j++) {
FieldVector<double,3> tmp;
- matrix.mv(coefficients[j].r, tmp);
- rotatedCoefficients[j].r = tmp;
+ matrix.mv(coefficients[j][_0].globalCoordinates(), tmp);
+ rotatedCoefficients[j][_0].globalCoordinates() = tmp;
- rotatedCoefficients[j].q = testRotations[i].mult(coefficients[j].q);
+ rotatedCoefficients[j][_1] = testRotations[i].mult(coefficients[j][_1]);
}
double energy = assembler.energy(localView,
@@ -143,7 +146,7 @@ void testFrameInvariance()
bool identicalPoints = false;
for (int j=0; j<domainDim+1; j++)
for (int k=0; k<domainDim+1; k++)
- if (j!=k and (testPoints[index[j]].r - testPoints[index[k]].r).two_norm() < 1e-5)
+ if (j!=k and (testPoints[index[j]][_0].globalCoordinates() - testPoints[index[k]][_0].globalCoordinates()).two_norm() < 1e-5)
identicalPoints = true;
if (identicalPoints)
diff --git a/test/cosseratrodenergytest.cc b/test/cosseratrodenergytest.cc
index 95cbc183..7a875cac 100644
--- a/test/cosseratrodenergytest.cc
+++ b/test/cosseratrodenergytest.cc
@@ -6,16 +6,20 @@
#include <dune/gfe/assemblers/cosseratrodenergy.hh>
#include <dune/gfe/localgeodesicfefunction.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
using namespace Dune;
+using namespace Dune::Indices;
int main (int argc, char *argv[]) try
{
- // Some types that I need
- typedef std::vector<RigidBodyMotion<double,3> > SolutionType;
+ // Type used for algebraic rod configurations
+ using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
+ using SolutionType = std::vector<TargetSpace>;
// Problem settings
const int numRodBaseElements = 100;
@@ -40,15 +44,13 @@ int main (int argc, char *argv[]) try
for (size_t i=0; i<x.size(); i++)
{
double s = double(i)/(x.size()-1);
- x[i].r[0] = 0.1*std::cos(2*M_PI*s);
- x[i].r[1] = 0.1*std::sin(2*M_PI*s);
- x[i].r[2] = s;
- x[i].q = Rotation<double,3>::identity();
+ x[i][_0] = {0.1*std::cos(2*M_PI*s), 0.1*std::sin(2*M_PI*s), s};
+ x[i][_1] = Rotation<double,3>::identity();
//x[i].q = Quaternion<double>(zAxis, (double(i)*M_PI)/(2*(x.size()-1)) );
}
FieldVector<double,3> zAxis(0); zAxis[2]=1;
- x.back().q = Rotation<double,3>(zAxis, M_PI/4);
+ x.back()[_1] = Rotation<double,3>(zAxis, M_PI/4);
// /////////////////////////////////////////////////////////////////////
// Create a second, rotated copy of the configuration
@@ -63,31 +65,29 @@ int main (int argc, char *argv[]) try
for (size_t i=0; i<rotatedX.size(); i++)
{
- rotatedX[i].r = rotation.rotate(x[i].r);
- rotatedX[i].r += displacement;
+ rotatedX[i][_0] = rotation.rotate(x[i][_0].globalCoordinates());
+ rotatedX[i][_0].globalCoordinates() += displacement;
- rotatedX[i].q = rotation.mult(x[i].q);
+ rotatedX[i][_1] = rotation.mult(x[i][_1]);
}
using GeodesicInterpolationRule = LocalGeodesicFEFunction<1, double,
FEBasis::LocalView::Tree::FiniteElement,
- RigidBodyMotion<double,3> >;
+ TargetSpace>;
GFE::CosseratRodEnergy<FEBasis,
GeodesicInterpolationRule,
double> localRodEnergy(gridView,
1,1,1,1e6,0.3);
- std::vector<RigidBodyMotion<double,3> > referenceConfiguration(gridView.size(1));
+ SolutionType referenceConfiguration(gridView.size(1));
for (const auto& vertex : vertices(gridView))
{
auto idx = gridView.indexSet().index(vertex);
- referenceConfiguration[idx].r[0] = 0;
- referenceConfiguration[idx].r[1] = 0;
- referenceConfiguration[idx].r[2] = vertex.geometry().corner(0)[0];
- referenceConfiguration[idx].q = Rotation<double,3>::identity();
+ referenceConfiguration[idx][_0] = {0.0, 0.0, vertex.geometry().corner(0)[0]};
+ referenceConfiguration[idx][_1] = Rotation<double,3>::identity();
}
localRodEnergy.setReferenceConfiguration(referenceConfiguration);
diff --git a/test/cosseratrodtest.cc b/test/cosseratrodtest.cc
index 59fc5a58..2bb10581 100644
--- a/test/cosseratrodtest.cc
+++ b/test/cosseratrodtest.cc
@@ -22,17 +22,20 @@
#include <dune/gfe/localgeodesicfefunction.hh>
#include <dune/gfe/localprojectedfefunction.hh>
#include <dune/gfe/riemanniantrsolver.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
-typedef RigidBodyMotion<double,3> TargetSpace;
+using namespace Dune;
+using namespace Dune::Indices;
+
+using TargetSpace = GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> >;
const int blocksize = TargetSpace::TangentVector::dimension;
// Approximation order of the finite element space
constexpr int order = 2;
-using namespace Dune;
-
int main (int argc, char *argv[]) try
{
MPIHelper::instance(argc, argv);
@@ -70,15 +73,13 @@ int main (int argc, char *argv[]) try
Functions::interpolate(scalarBasis, referenceConfigurationX, identity);
- using Configuration = std::vector<RigidBodyMotion<double,3> >;
+ using Configuration = std::vector<TargetSpace>;
Configuration referenceConfiguration(scalarBasis.size());
for (std::size_t i=0; i<referenceConfiguration.size(); i++)
{
- referenceConfiguration[i].r[0] = 0;
- referenceConfiguration[i].r[1] = 0;
- referenceConfiguration[i].r[2] = referenceConfigurationX[i];
- referenceConfiguration[i].q = Rotation<double,3>::identity();
+ referenceConfiguration[i][_0] = {0.0, 0.0, referenceConfigurationX[i]};
+ referenceConfiguration[i][_1] = Rotation<double,3>::identity();
}
// Select the reference configuration as initial iterate
@@ -114,12 +115,12 @@ int main (int argc, char *argv[]) try
}
// Set Dirichlet values
- x[rightBoundaryDof].r = {1,0,0};
+ x[rightBoundaryDof][_0] = {1,0,0};
FieldVector<double,3> axis = {1,0,0};
double angle = 0;
- x[rightBoundaryDof].q = Rotation<double,3>(axis, M_PI*angle/180);
+ x[rightBoundaryDof][_1] = Rotation<double,3>(axis, M_PI*angle/180);
//////////////////////////////////////////////
// Create the energy and assembler
@@ -160,7 +161,7 @@ int main (int argc, char *argv[]) try
// Create a solver for the rod problem
/////////////////////////////////////////////
- RiemannianTrustRegionSolver<ScalarBasis,RigidBodyMotion<double,3> > solver;
+ RiemannianTrustRegionSolver<ScalarBasis,TargetSpace> solver;
solver.setup(grid,
&rodAssembler,
diff --git a/test/geodesicfeassemblerwrappertest.cc b/test/geodesicfeassemblerwrappertest.cc
index ae08c341..79d0940d 100644
--- a/test/geodesicfeassemblerwrappertest.cc
+++ b/test/geodesicfeassemblerwrappertest.cc
@@ -30,6 +30,9 @@
#include <dune/gfe/assemblers/mixedgfeassembler.hh>
#include <dune/gfe/assemblers/mixedlocalgfeadolcstiffness.hh>
#include <dune/gfe/riemanniantrsolver.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
#include <dune/gfe/assemblers/geodesicfeassemblerwrapper.hh>
@@ -64,7 +67,7 @@ using MatrixRow1 = MultiTypeBlockVector<BCRSMatrix<FieldMatrix<double,dimCR,dim>
using MatrixType = MultiTypeBlockMatrix<MatrixRow0,MatrixRow1>;
//Types for the Non-mixed space
-using RBM = RigidBodyMotion<double, dim>;
+using RBM = GFE::ProductManifold<RealTuple<double,dim>,Rotation<double, dim> >;
const static int blocksize = RBM::TangentVector::dimension;
using CorrectionTypeWrapped = BlockVector<FieldVector<double, blocksize> >;
using MatrixTypeWrapped = BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >;
@@ -150,7 +153,6 @@ int main (int argc, char *argv[])
RealTuple<double,dim>,
Rotation<double,dim> > mixedAssembler(compositeBasis, &mixedLocalGFEADOLCStiffness);
- using RBM = RigidBodyMotion<double, dim>;
using DeformationFEBasis = Functions::LagrangeBasis<GridView,displacementOrder>;
DeformationFEBasis deformationFEBasis(gridView);
using GFEAssemblerWrapper = GFE::GeodesicFEAssemblerWrapper<CompositeBasis, DeformationFEBasis, RBM, RealTuple<double, dim>, Rotation<double,dim> >;
@@ -179,10 +181,8 @@ int main (int argc, char *argv[])
for (int j = 0; j < gridDim; j++)
initialDeformation[i][j] = identity[i][j];
x[_0][i] = initialDeformation[i];
- for (int j = 0; j < dim; j ++) { // Displacement part
- xRBM[i].r[j] = x[_0][i][j];
- }
- xRBM[i].q = x[_1][i]; // Rotation part
+ xRBM[i][_0] = x[_0][i];
+ xRBM[i][_1] = x[_1][i]; // Rotation part
}
//////////////////////////////////////////////////////////////////////////////
diff --git a/test/localgeodesicfefunctiontest.cc b/test/localgeodesicfefunctiontest.cc
index 1ff9526d..400493ea 100644
--- a/test/localgeodesicfefunctiontest.cc
+++ b/test/localgeodesicfefunctiontest.cc
@@ -303,7 +303,6 @@ int main()
test<UnitVector<double,2>,1>(GeometryTypes::simplex(1));
test<UnitVector<double,3>,1>(GeometryTypes::simplex(1));
test<Rotation<double,3>,1>(GeometryTypes::simplex(1));
- test<RigidBodyMotion<double,3>,1>(GeometryTypes::simplex(1));
typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > CrazyManifold;
test<CrazyManifold,1>(GeometryTypes::simplex(1));
@@ -315,7 +314,6 @@ int main()
test<UnitVector<double,2>,2>(GeometryTypes::simplex(2));
test<UnitVector<double,3>,2>(GeometryTypes::simplex(2));
test<Rotation<double,3>,2>(GeometryTypes::simplex(2));
- test<RigidBodyMotion<double,3>,2>(GeometryTypes::simplex(2));
typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > CrazyManifold;
test<CrazyManifold,2>(GeometryTypes::simplex(2));
@@ -327,7 +325,6 @@ int main()
test<UnitVector<double,2>,2>(GeometryTypes::cube(2));
test<UnitVector<double,3>,2>(GeometryTypes::cube(2));
test<Rotation<double,3>,2>(GeometryTypes::cube(2));
- test<RigidBodyMotion<double,3>,2>(GeometryTypes::cube(2));
typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > CrazyManifold;
test<CrazyManifold,2>(GeometryTypes::cube(2));
diff --git a/test/localprojectedfefunctiontest.cc b/test/localprojectedfefunctiontest.cc
index 58373aba..31d3ff6c 100644
--- a/test/localprojectedfefunctiontest.cc
+++ b/test/localprojectedfefunctiontest.cc
@@ -98,12 +98,6 @@ void testDerivativeTangentiality(const Rotation<double,vectorDim-1>& x,
const FieldMatrix<double,vectorDim,domainDim>& derivative)
{}
-// the columns of the derivative must be tangential to the manifold
-template <int domainDim, int vectorDim>
-void testDerivativeTangentiality(const RigidBodyMotion<double,3>& x,
- const FieldMatrix<double,vectorDim,domainDim>& derivative)
-{}
-
// the columns of the derivative must be tangential to the manifold
template <int domainDim, int vectorDim,typename ... TargetSpaces>
void testDerivativeTangentiality(const Dune::GFE::ProductManifold<TargetSpaces...>& x,
@@ -270,7 +264,6 @@ int main()
test<UnitVector<double,2>,1>(GeometryTypes::line);
test<UnitVector<double,3>,1>(GeometryTypes::line);
test<Rotation<double,3>,1>(GeometryTypes::line);
- test<RigidBodyMotion<double,3>,1>(GeometryTypes::line);
typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > CrazyManifold;
test<CrazyManifold, 1>(GeometryTypes::line);
@@ -283,7 +276,6 @@ int main()
test<RealTuple<double,3>,2>(GeometryTypes::triangle);
test<UnitVector<double,3>,2>(GeometryTypes::triangle);
test<Rotation<double,3>,2>(GeometryTypes::triangle);
- test<RigidBodyMotion<double,3>,2>(GeometryTypes::triangle);
typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > CrazyManifold;
test<CrazyManifold, 2>(GeometryTypes::triangle);
@@ -295,7 +287,6 @@ int main()
test<UnitVector<double,2>,2>(GeometryTypes::quadrilateral);
test<UnitVector<double,3>,2>(GeometryTypes::quadrilateral);
test<Rotation<double,3>,2>(GeometryTypes::quadrilateral);
- test<RigidBodyMotion<double,3>,2>(GeometryTypes::quadrilateral);
typedef Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > CrazyManifold;
test<CrazyManifold, 2>(GeometryTypes::quadrilateral);
diff --git a/test/nonplanarcosseratenergytest.cc b/test/nonplanarcosseratenergytest.cc
index ddb2ee83..5753614b 100644
--- a/test/nonplanarcosseratenergytest.cc
+++ b/test/nonplanarcosseratenergytest.cc
@@ -13,7 +13,9 @@
#include <dune/gfe/cosseratvtkwriter.hh>
#include <dune/gfe/assemblers/nonplanarcosseratshellenergy.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/productmanifold.hh>
+#include <dune/gfe/spaces/realtuple.hh>
+#include <dune/gfe/spaces/rotation.hh>
#include "multiindex.hh"
#include "valuefactory.hh"
@@ -24,7 +26,7 @@ static const int dim = 2;
static const int dimworld = 3;
using GridType = FoamGrid<dim,dimworld>;
-using TargetSpace = RigidBodyMotion<double,dimworld>;
+using TargetSpace = GFE::ProductManifold<RealTuple<double,dimworld>,Rotation<double,dimworld> >;
//////////////////////////////////////////////////////////
// Make a test grid consisting of a single triangle
@@ -102,16 +104,15 @@ double calculateEnergy(const int numLevels, const F1 referenceConfigurationFunct
BlockVector<FieldVector<double,3> > helperVector2(feBasis.size());
Dune::Functions::interpolate(deformationPowerBasis, helperVector2, configurationFunction);
for (std::size_t i = 0; i < feBasis.size(); i++) {
- for (int j = 0; j < dimworld; j++)
- sol[i].r[j] = helperVector2[i][j];
+ sol[i][_0].globalCoordinates() = helperVector2[i];
FieldVector<double,4> idRotation = {0, 0, 0, 1}; //set rotation = Id everywhere
Rotation<double,dimworld> rotation(idRotation);
FieldMatrix<double,dimworld,dimworld> rotationMatrix(0);
rotation.matrix(rotationMatrix);
- sol[i].q.set(rotationMatrix);
- solTuple[_0][i] = sol[i].r;
- solTuple[_1][i] = sol[i].q;
+ sol[i][_1].set(rotationMatrix);
+ solTuple[_0][i] = sol[i][_0];
+ solTuple[_1][i] = sol[i][_1];
}
CosseratVTKWriter<GridType>::write<FEBasis>(feBasis, solTuple, "configuration_l" + std::to_string(numLevels));
diff --git a/test/rigidbodymotiontest.cc b/test/rigidbodymotiontest.cc
deleted file mode 100644
index af824ad2..00000000
--- a/test/rigidbodymotiontest.cc
+++ /dev/null
@@ -1,129 +0,0 @@
-#include "config.h"
-
-#include <dune/common/parallel/mpihelper.hh>
-
-#include <dune/geometry/type.hh>
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
-
-#include <dune/gfe/localgeodesicfefunction.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
-
-#include "multiindex.hh"
-#include "valuefactory.hh"
-
-const double eps = 1e-4;
-
-typedef RigidBodyMotion<double,3> TargetSpace;
-
-using namespace Dune;
-
-
-template <int domainDim,class LocalFiniteElement>
-Tensor3<double,3,3,domainDim> evaluateDerivativeFD(const LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace>& f,
- const Dune::FieldVector<double,domainDim>& local)
-{
- const double stepSize = std::sqrt(std::numeric_limits<double>::epsilon());
- Tensor3<double,3,3,domainDim> result(0);
-
- for (int i=0; i<domainDim; i++) {
-
- Dune::FieldVector<double, domainDim> forward = local;
- Dune::FieldVector<double, domainDim> backward = local;
-
- forward[i] += stepSize;
- backward[i] -= stepSize;
-
- TargetSpace forwardValue = f.evaluate(forward);
- TargetSpace backwardValue = f.evaluate(backward);
-
- FieldMatrix<double,3,3> forwardMatrix, backwardMatrix;
- forwardValue.q.matrix(forwardMatrix);
- backwardValue.q.matrix(backwardMatrix);
-
- FieldMatrix<double,3,3> fdDer = (forwardMatrix - backwardMatrix) / (2*stepSize);
-
- for (int j=0; j<3; j++)
- for (int k=0; k<3; k++)
- result[j][k][i] = fdDer[j][k];
-
- }
-
- return result;
-
-}
-
-template <int domainDim>
-void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners)
-{
- // Make local fe function to be tested
- LagrangeLocalFiniteElementCache<double,double,domainDim,1> feCache;
- typedef typename LagrangeLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement;
-
- LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement, TargetSpace> f(feCache.get(GeometryTypes::simplex(domainDim)), corners);
-
- // A quadrature rule as a set of test points
- int quadOrder = 3;
-
- const auto& quad = Dune::QuadratureRules<double, domainDim>::rule(GeometryTypes::simplex(domainDim), quadOrder);
-
- for (size_t pt=0; pt<quad.size(); pt++) {
-
- const Dune::FieldVector<double,domainDim>& quadPos = quad[pt].position();
-
- // evaluate actual derivative
- Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::dimension, domainDim> derivative = f.evaluateDerivative(quadPos);
-
- Tensor3<double,3,3,domainDim> DR = f.evaluate(quadPos).quaternionTangentToMatrixTangent(derivative);
-
- // evaluate fd approximation of derivative
- Tensor3<double,3,3,domainDim> DR_fd = evaluateDerivativeFD(f,quadPos);
-
- double maxDiff = 0;
- for (int i=0; i<3; i++)
- for (int j=0; j<3; j++)
- for (int k=0; k<domainDim; k++)
- maxDiff = std::max(maxDiff, std::abs(DR[i][j][k] - DR_fd[i][j][k]));
-
- if ( maxDiff > 100*eps ) {
- std::cout << className(corners[0]) << ": Analytical gradient does not match fd approximation." << std::endl;
- std::cout << "Analytical:\n " << DR << std::endl;
- std::cout << "FD :\n " << DR_fd << std::endl;
- assert(false);
- }
-
- }
-}
-
-int main(int argc, char** argv)
-{
- MPIHelper::instance(argc, argv);
-
- const int domainDim = 2;
-
- ////////////////////////////////////////////////////////////////////////////
- // Create a local assembler object
- ////////////////////////////////////////////////////////////////////////////
- std::cout << " --- Testing derivative of rotation matrix, domain dimension: " << domainDim << " ---" << std::endl;
-
- std::vector<Rotation<double,3> > testPoints;
-
- ValueFactory<Rotation<double,3> >::get(testPoints);
-
- int nTestPoints = testPoints.size();
-
- // Set up elements of SO(3)
- std::vector<TargetSpace> corners(domainDim+1);
-
- ::MultiIndex index(domainDim+1, nTestPoints);
- int numIndices = index.cycle();
-
- for (int i=0; i<numIndices; i++, ++index)
- {
- for (int j=0; j<domainDim+1; j++)
- corners[j].q = testPoints[index[j]];
-
- testDerivativeOfRotationMatrix<domainDim>(corners);
- }
-}
diff --git a/test/targetspacetest.cc b/test/targetspacetest.cc
index 5c2d855c..a0eb6b03 100644
--- a/test/targetspacetest.cc
+++ b/test/targetspacetest.cc
@@ -412,8 +412,6 @@ int main() try
test<Dune::GFE::ProductManifold<RealTuple<double,1>,Rotation<double,3>,UnitVector<double,2> > >();
test<Dune::GFE::ProductManifold<Rotation<double,3>,UnitVector<double,5> > >();
- // Test the RigidBodyMotion class
- test<RigidBodyMotion<double,3> >();
//
// test<HyperbolicHalfspacePoint<double,2> >();
diff --git a/test/valuefactory.hh b/test/valuefactory.hh
index 303f5e09..ed853b9e 100644
--- a/test/valuefactory.hh
+++ b/test/valuefactory.hh
@@ -6,7 +6,7 @@
#include <dune/gfe/spaces/hyperbolichalfspacepoint.hh>
#include <dune/gfe/spaces/orthogonalmatrix.hh>
#include <dune/gfe/spaces/productmanifold.hh>
-#include <dune/gfe/spaces/rigidbodymotion.hh>
+#include <dune/gfe/spaces/realtuple.hh>
#include <dune/gfe/spaces/rotation.hh>
#include <dune/gfe/spaces/unitvector.hh>
@@ -174,13 +174,15 @@ public:
/** \brief A class that creates sets of values of various types, to be used in unit tests
*
- * This is the specialization for RigidBodyMotion<3>
+ * This is the specialization for a ProductManifold<RealTuple,Rotation>
*/
template <>
-class ValueFactory<RigidBodyMotion<double,3> >
+class ValueFactory<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >
{
public:
- static void get(std::vector<RigidBodyMotion<double,3> >& values) {
+ static void get(std::vector<Dune::GFE::ProductManifold<RealTuple<double,3>,Rotation<double,3> > >& values) {
+
+ using namespace Dune::Indices;
std::vector<RealTuple<double,3> > rValues;
ValueFactory<RealTuple<double,3> >::get(rValues);
@@ -194,7 +196,10 @@ public:
// Set up elements of SE(3)
for (int i=0; i<nTestPoints; i++)
- values[i] = RigidBodyMotion<double,3>(rValues[i].globalCoordinates(),qValues[i]);
+ {
+ values[i][_0] = rValues[i];
+ values[i][_1] = qValues[i];
+ }
}
--
GitLab