diff --git a/CHANGELOG.md b/CHANGELOG.md
index 52d877d9b7d4c41dcc873e5b62a2d1b3e5f79c4a..b51284fef3cdeb1516e22e0a254a3e4b8e81f1f1 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 3caffc35ad99fa4ef61f53d22aaa00ddd34880aa..08831326982b1b697f3fd8062b650e38e0fd06fb 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 b2c6dced2911729df81ac55800222bac65e008db..5a3361354d1429d7bdfa5c23aea02443c07682a6 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 e6a3fdda5f6629bd2b8c7a60944e24fa7534e42b..12d2b789530bdbeb07227514992cedb4354d9ee9 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 40fb241b0586f89a51817707b9661232e6c7c522..23c8dd8d2b34656fdb1dd3e07e3694d63f05a46e 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 fd1a9eb04bda19594ae45814abbed3f2c1da1020..98d941b502a7cd2e5e6eafa75aa68ed40793788a 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 b5d6c16b39a4a3e8c04b733ed3865bb2d88d7a6a..e70d84150e64900b765f8acc5a3198688b6b80a1 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 7e0e6109999849cf1fe6b2f1d5577efaa212b65b..659634ff549deab6a2e56441ad37e8f3918c2589 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 feb41c418e052da2c6c8209828a2183393fd40fd..1bd513471403dc999432e53feafb49e675024833 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 76e0d86417eee26bdf66c38d57a607954f65d781..6892e93604dfd22ed4872f30835899739ea1e9d3 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 46acb95ad1710e5fa1136c5dff9d7be7b2082b95..5f1a0c49d8ea8493529d82bf6b76fcb6d7045955 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 3f27c7d544db4b507a7776de58596695abc401a2..bea34db0b847c2aa602a59babb5e5c053d1662d5 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 f158d6058cc5376c7eafc6889e595d9595678227..bb85104de0a572b451e47267359b0047bc58db45 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 e284c18f904f380184bca94dbf64cbab861ca35a..0000000000000000000000000000000000000000
--- 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 f5b5b2776560481484d2d83b5d5025de0dcc3afb..db771c377dc1ff0cfdc8e348a519b8a27a9567fd 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 c81a9192ca8a3a29e7d5ff7cdf98ba5c43c4c294..8906f52794c4407784e626f6f220a299e0e1d8e8 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 016275d86a3feb39c810eec7f8d044295df7114a..33fb40eb7741f4765326e3687c739b38a9e08e8f 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 4033c84227f1845277bc2574f364ee9a953925bc..940e55dab9f448e16b486768e5314a55959f6960 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 d934c906df95da45c888d15e0d427fb8eb6985a2..e0e25e63b2e31d82f37cf04a8dac7d2f346023af 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 6267e6d2ee6ca9089a6e2655e44499078f72be16..c1b58b705625900b1c08bb6a7e8e3480fad001e9 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 296f756ec74104177e2228887c4e4f55639d5119..3ce0039e093b9151c38b48c7db1306db1cad13f3 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 61d8f1e5d65c6b1aa8d9a63b2b522f84da9d79a1..a06d3b981c2f28de2e8923d2055ae3d55117e44f 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 b43c90621949af6a17d0b5ca5c1cfd5c236208f4..d53f84b746547818c6392801fa2eb1e6b470eeaf 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 f9348a2fa3a23500b484d26b0372a86262fa3ccc..9589bb562992ebf4e407d57f5d6d82d3f1341f1f 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 95cbc183b4e40dc0f12c243efb7cc0603443663d..7a875cac7ef7fa36fc43f7552820d0f9c883a73d 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 59fc5a5862a26ea509de680fc4e116957e22238c..2bb105815cef03448a18c9eb6ddf7c04d061d7f1 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 ae08c341bc5772c22a3592c187d7b07a18ecc459..79d0940d16442095d41ebb54c2e52d2895abad10 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 1ff9526d1ccdcb04d60855351b2f2e285fcbd8dd..400493ea1c93b302f3297f5f8a6c626c4047d894 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 58373aba475d8538f2add22a51552bc0796a1fdf..31d3ff6c64b9052320e532b09ec52481194f70db 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 ddb2ee83f4e0bc2a8a94acbe2a5bb374b923fe04..5753614b30ff5f0df15c7bc8f6162b1fc0d7eca0 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 af824ad2ae13a24497ad274072382f35b0c7b2da..0000000000000000000000000000000000000000
--- 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 5c2d855c6de2278aa2754403e12b0ca973f5b3df..a0eb6b039279a2daadac95748f72ac1b8a3559d8 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 303f5e09cb5dc880cf4ba1d0cf4279875e1203e2..ed853b9e2b727de55c88bace5d584d76c1bf969f 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];
+ }
}